with DC_Keyb;
with DC_Math;
with Text_IO; use Text_IO;
with DC_Scrn;

procedure Kohonen is
---------------------------------------------------------------------------
---------------------------------------------------------------------------
  Neuron_Count : constant long_integer := 400;
  Time_Max     : constant long_integer := 1_000_000;
  Neighborhood : constant long_integer := Neuron_Count;
  B : long_integer;
  L : float;
  type Float_Array is array ( positive range <> ) of float;
  Activity:  array ( 1..Neuron_Count ) of float;
  Winner  : long_integer;
  type Float_2 is array ( 1..2 ) of float;
  type Weights_Array is array ( 1..Neuron_Count ) of Float_2;
  Weights:  Weights_Array := ( others => ( 0.0, 0.0 ) );
  Inputs :  Float_2;
  Winners : array ( 0..3 ) of long_integer := ( others => 0 );
  Counts  : array ( 0..3 ) of long_integer := ( others => 1 );
  A : natural;
  Glix1,
  Glix2,
  Glix3 : long_integer;
  Glr : array ( long_integer range 1..97 ) of float;
  Dum, Sum : float := 0.0;
  package Float_IO is new Text_IO.Float_IO ( float );
  package Int_IO is new Text_IO.Integer_IO ( long_integer );
  Expect_Sum : array ( 0..3, 1..2 ) of float := ( others => ( others => 0.0 ) );
  Expect_Count : array ( 0..3 ) of long_integer := ( others => 0 );
  Learn_Ratio: float;

function Max ( A, B : long_integer ) return long_integer is
---------------------------------------------------------------------------
begin
  if A > B then
    return A;
  else
    return B;
  end if;
end Max;

function Min ( A, B : long_integer ) return long_integer is
---------------------------------------------------------------------------
begin
  if A < B then
    return A;
  else
    return B;
  end if;
end Min;

function "*" ( A:  float; B: Float_2 ) return Float_2 is
---------------------------------------------------------------------------
  Temp : Float_2;
begin
  for index in Float_2'range loop
    Temp ( index ) := A * B ( index );
  end loop;
  return Temp;
end "*";

function "+" ( A, B: Float_2 ) return Float_2 is
---------------------------------------------------------------------------
  Temp : Float_2;
begin
  for index in Float_2'range loop
    Temp ( index ) := A ( index ) + B ( index );
  end loop;
  return Temp;
end "+";


function "-" ( A, B: Float_2 ) return Float_2 is
---------------------------------------------------------------------------
  Temp : Float_2;
begin
  for index in Float_2'range loop
    Temp ( index ) := A ( index ) - B ( index );
  end loop;
  return Temp;
end "-";

function "Abs" ( A : in Float_2 ) return float is
---------------------------------------------------------------------------
  Sum: float := 0.0;
begin
  for index in A'range loop
    Sum := Sum + A ( index ) ** 2;
  end loop;
  return DC_Math.SqRt ( Sum );
end "Abs";

function Ran1 ( IDum: in long_integer := 1 ) return float is
---------------------------------------------------------------------------
-- Set IDum to a negative number to re-initialize.
---------------------------------------------------------------------------
  M1  : constant long_integer := 259_200;
  IA1 : constant long_integer :=   7_141;
  IC1 : constant long_integer :=  54_773;
  RM1 : constant float        := 1.0 / float ( M1 );
  M2  : constant long_integer := 134_456;
  IA2 : constant long_integer :=   8_121;
  IC2 : constant long_integer :=  28_411;
  RM2 : constant float        := 1.0 / float ( M2 );
  M3  : constant long_integer := 243_000;
  IA3 : constant long_integer :=   4_561;
  IC3 : constant long_integer :=  51_349;
  J   :          long_integer;
  Temp:          float;
begin
  if IDum < 0 then
    Glix1 := ( IC1 - IDum ) mod M1;
    Glix1 := ( IA1 * Glix1 + IC1 ) mod M1;
    Glix2 := Glix1 mod M2;
    Glix1 := ( IA1 * Glix1 + IC1 ) mod M1;
    Glix3 := Glix1 mod M3;
    for index in Glr'range loop
      Glix1 := ( IA1 * Glix1 + IC1 ) mod M1;
      Glix2 := ( IA2 * Glix2 + IC2 ) mod M2;
      Glr ( index ) := ( float ( Glix1 ) + float ( Glix2 ) * RM2 ) * RM1;
    end loop;
  end if;
  Glix1 := ( IA1 * Glix1 + IC1 ) mod M1;
  Glix2 := ( IA2 * Glix2 + IC2 ) mod M2;
  Glix3 := ( IA3 * Glix3 + IC3 ) mod M3;
  J := 1 + ( 97 * Glix3 ) / M3;
--  if ( J > 97 ) or ( J < 1 ) then
--    DC_Keyb.Pause ( "Pause in routine RAN1" );
--  end if;
  Temp := Glr ( J );
  Glr ( J ) := ( float ( Glix1 ) + float ( Glix2 ) * RM2 ) * RM1;
  return Temp;
end Ran1;

procedure Inputs_Set ( Inputs:  out Float_2; Mode : out natural ) is
---------------------------------------------------------------------------
  A : long_integer;
  Z, C, S : float;
  Pi : constant float := 3.14159;
  Random : float;
begin
  Random := Ran1;
  A := long_integer ( Random * 3.0 );
  if A = 0 then
    Z := ( 7.0 * Pi / 6.0 ) + ( 2.0 * Pi * Ran1 / 3.0 );
  elsif A > 0 then
    Z := 2.0 * Pi * Ran1;
  end if;
  C := DC_Math.Cos ( Z );
  S := DC_Math.Sin ( Z );
  case A is
    when 0 =>
      Inputs ( 1 ) :=        0.7 * C;
      Inputs ( 2 ) :=        0.7 * S;
    when 1 =>
      Inputs ( 1 ) :=  0.5 + 0.1 * C;
      Inputs ( 2 ) :=  0.7 + 0.2 * S;
    when 2 =>
      Inputs ( 1 ) := -0.5 + 0.1 * C;
      Inputs ( 2 ) :=  0.7 + 0.2 * S;
    when 3 =>
      Inputs ( 1 ) :=              C;
      Inputs ( 2 ) :=              S;
    when others =>
      raise Constraint_Error;
  end case;
  Mode := natural ( A );
end Inputs_Set;

---------------------------------------------------------------------------
---------------------------------------------------------------------------
begin
  Dum := Ran1 ( -36 );
  for Time in 1..Time_Max loop
    if Time mod 1000 = 0 then
      for index in 0..3 loop
	Winners ( index ) := Winners ( index ) / 2;
	Counts  ( index ) := Counts  ( index ) / 2;
      end loop;
    end if;
    Inputs_Set ( Inputs, A );
    Expect_Sum ( A, 1 ) := Expect_Sum ( A, 1 ) + Inputs ( 1 );
    Expect_Sum ( A, 2 ) := Expect_Sum ( A, 2 ) + Inputs ( 2 );
    Expect_Count ( A ) := Expect_Count ( A ) + 1;
    for Neuron in Activity'range loop
      Activity ( Neuron )
	:=  -Abs ( Inputs - Weights ( Neuron ) );
    end loop;
    Winner := Activity'first;
    for Neuron in Activity'range loop
      if Activity ( Neuron ) > Activity ( Winner ) then
	Winner := Neuron;
      end if;
    end loop;
    Winners ( A ) := Winners ( A ) + Winner;
    Counts  ( A ) := Counts  ( A ) + 1;
    L := 1.0 / ( float ( Time / 100 ) + 1.0 );
    B := Neuron_Count * ( Time_Max - Time ) / Time_Max;
    for Neuron in 1..Neuron_Count loop
      Learn_Ratio := L *
	DC_Math.Exp ( - float ( abs ( Winner - Neuron ) ) / float ( B  + 1  ) );
      Weights ( Neuron ) := Weights ( Neuron )
	  + Learn_Ratio * ( Inputs - Weights ( Neuron ) );
    end loop;
  end loop;
--  Int_IO.Put ( long_integer ( Time ), 7 );
  for Neuron in 1..Neuron_Count loop
    Float_IO.Put ( Weights ( Neuron ) ( 1 ), 2, 10, 0 );
    Put ( "  " );
    Float_IO.Put ( Weights ( Neuron ) ( 2 ), 2, 10, 0 );
    Put_Line ( "" );
  end loop;
  Put ( "  B:" );
  Int_IO.Put ( B, 4 );
  Put ( "  L:" );
  Float_IO.Put ( L, 1, 2 );
  Put ( "  A:" );
  Int_IO.Put ( long_integer ( A ), 1 );
  Put ( "  W:" );
  Int_IO.Put ( Winner, 4 );
  Put ( "  " );
  for index in Winners'range loop
    Float_IO.Put (
      float ( Winners ( index ) ) / float ( Counts ( index ) ), 4, 1, 0 );
    Put ( "  " );
  end loop;
  Put_Line ( "" );
  for index in 0..3 loop
    Float_IO.Put ( Expect_Sum ( index, 1 ) / float ( Expect_Count ( index ) ), 2, 2, 0 );
    Put ( "  " );
  end loop;
  Put_Line ( "" );
  for index in 0..3 loop
    Float_IO.Put ( Expect_Sum ( index, 2 ) / float ( Expect_Count ( index ) ), 2, 2, 0 );
    Put ( "  " );
  end loop;
  Put_Line ( "" );
end Kohonen;