program Recurrent_Auto_Associate_Memory;

const p_max       = 10;
      q_max       = 10;
      lambda      = 0.05;   {learning rate}
      mu          = 0.7;   {momentum rate}
      epoch_max   = 30000; {max. number of adaptation epochs}
      N_train_max = 50;    {max. number of training set patterns}
      length_max  = 5;    {max. length of training set patterns}
      N_token_max = 20;    {max. number of tokens in paterns}
      print_freq  = 100;
      epsilon     = 0.1;   {precision of reconstruction}

type vector_Left     = array [1..p_max] of real;
     vector_Right    = array [1..q_max] of real;
     vector_Left_Long= array [1..p_max] of longint;
     vector_Right_Long=array [1..q_max] of longint;
     vector_Hidden   = array [1..p_max] of real;
     matrix_Left_OH  = array [1..p_max,1..p_max] of real;
     matrix_Right_OH = array [1..q_max,1..p_max] of real;
     matrix_Left_HI  = array [1..p_max,1..p_max] of real;
     matrix_Right_HI = array [1..p_max,1..q_max] of real;
     vector_train    = array [1..N_train_max] of integer;
     matrix_train    = array [1..N_train_max,1..length_max] of byte;
     matrix_token    = array [1..N_token_max,1..q_max] of byte;
     vector_string   = array [1..length_max] of integer;
     matrix_train_y  = array [1..N_train_max,1..p_max] of real;


var data_in       : text;
    i,j,k         : integer;
    p,q           : integer;
    time          : integer;
    index         : integer;
    index_prime   : integer;
    lengthx       : integer;
    N_train       : integer;
    N_token       : integer;
    sum           : real;
    theta_Left_O  : vector_Left;
    theta_Right_O : vector_Right;
    theta_H       : vector_Hidden;
    w_Left_OH     : matrix_Left_OH;
    w_Right_OH    : matrix_Right_OH;
    w_Left_HI     : matrix_Left_HI;
    w_Right_HI    : matrix_Right_HI;

    delta_theta_Left_O  : vector_Left;
    delta_theta_Right_O : vector_Right;
    delta_theta_H       : vector_Hidden;
    delta_w_Left_OH     : matrix_Left_OH;
    delta_w_Right_OH    : matrix_Right_OH;
    delta_w_Left_HI     : matrix_Left_HI;
    delta_w_Right_HI    : matrix_Right_HI;

    train_set     : matrix_train;
    train_length  : vector_train;
    token         : matrix_token;
    sequence      : vector_string;
    x_L,x_prime_L : vector_Left;
    diff_L        : vector_Left;
    diff_R        : vector_Right;
    diff_L_tot    : vector_Left;
    diff_R_tot    : vector_Right;
    x_R,x_prime_R : vector_Right;
    y,y_code      : vector_Hidden;
    y_req         : vector_hidden;
    grad_theta_Left_O : vector_Left;
    grad_theta_Right_O: vector_Right;
    grad_theta_H      : vector_Hidden;

    grad_theta_Left_O_tot : vector_Left;
    grad_theta_Right_O_tot: vector_Right;
    grad_theta_H_tot      : vector_Hidden;
    grad_w_Left_OH_tot    : matrix_Left_OH;
    grad_w_Right_OH_tot   : matrix_Right_OH;
    grad_w_Left_HI_tot    : matrix_Left_HI;
    grad_w_Right_HI_tot   : matrix_Right_HI;

    counter_L : vector_Left_Long;
    counter_R : vector_Right_Long;

    pop_y     : matrix_train_y;

procedure Data_input;
var i,j : integer;
begin write('RandSeed='); readln(RandSeed);

      readln(data_in,N_train);
      if N_train>N_train_max then
      begin write('N_train is too large');
           halt;
      end;
      readln(data_in,N_token,p,q);
      if N_token>N_token_max then
      begin write('N_token is too large');
           halt;
      end;
      if p>p_max then
      begin write('p is too large');
           halt;
      end;
      if q>q_max then
      begin write('q is too large');
           halt;
      end;
      for i:=1 to N_token do
      begin for j:=1 to q do read(data_in,token[i,j]);
            readln(data_in);
      end;
      for i:=1 to N_train do
      begin read(data_in, train_length[i]);
            for j:=1 to train_length[i] do
            read(data_in,train_set[i,j]);
            readln(data_in);
      end;
end;


procedure Ini_coefficient;
var i,j : integer;
begin for i:=1 to p do theta_Left_O[i]:=1-2*random;
      for i:=1 to q do theta_Right_O[i]:=1-2*random;
      for i:=1 to p do theta_H[i]:=1-2*random;
      for i:=1 to p do
      for j:=1 to p do w_Left_OH[i,j]:=1-2*random;
      for i:=1 to q do
      for j:=1 to p do w_Right_OH[i,j]:=1-2*random;
      for i:=1 to p do
      for j:=1 to p do w_Left_HI[i,j]:=1-2*random;
      for i:=1 to p do
      for j:=1 to q do w_Right_HI[i,j]:=1-2*random;

      for i:=1 to p do delta_theta_Left_O[i]:=0.0;
      for i:=1 to q do delta_theta_Right_O[i]:=0.0;
      for i:=1 to p do delta_theta_H[i]:=0.0;
      for i:=1 to p do
      for j:=1 to p do delta_w_Left_OH[i,j]:=0.0;
      for i:=1 to q do
      for j:=1 to p do delta_w_Right_OH[i,j]:=0.0;
      for i:=1 to p do
      for j:=1 to p do delta_w_Left_HI[i,j]:=0.0;
      for i:=1 to p do
      for j:=1 to q do delta_w_Right_HI[i,j]:=0.0;
end;


procedure Calcul_grad_tot(    index                 : integer;
                          var grad_theta_Left_O_tot : vector_Left;
                          var grad_theta_Right_O_tot: vector_Right;
                          var grad_theta_H_tot      : vector_Hidden;
                          var grad_w_Left_OH_tot    : matrix_Left_OH;
                          var grad_w_Right_OH_tot   : matrix_Right_OH;
                          var grad_w_Left_HI_tot    : matrix_Left_HI;
                          var grad_w_Right_HI_tot   : matrix_Right_HI;
                              theta_Left_O  : vector_Left;
                              theta_Right_O : vector_Right;
                              theta_H       : vector_Hidden;
                              w_Left_OH     : matrix_Left_OH;
                              w_Right_OH    : matrix_Right_OH;
                              w_Left_HI     : matrix_Left_HI;
                              w_Right_HI    : matrix_Right_HI;
                          var diff_L        : vector_Left;
                          var diff_R        : vector_Right;
                          var y_code        : vector_hidden);

var i,j,k,lengthx     : integer;
    sum               : real;
    counter_L         : vector_Left_Long;
    Counter_R         : vector_Right_Long;
    grad_theta_Left_O : vector_Left;
    grad_theta_Right_O: vector_Right;
    grad_theta_H      : vector_Hidden;
    x_L,x_prime_L     : vector_Left;
    x_R,x_prime_R     : vector_Right;
    y                 : vector_Hidden;

begin lengthx:=train_length[index];
      for i:=1 to lengthx do sequence[i]:=train_set[index,i];

      for i:=1 to p do
      begin diff_L[i]:=0.0; counter_L[i]:=0 end;
      for i:=1 to q do
      begin diff_R[i]:=0.0; counter_R[i]:=0 end;

      for i:=1 to p do grad_theta_Left_O_tot[i]:=0.0;
      for i:=1 to q do grad_theta_Right_O_tot[i]:=0.0;
      for i:=1 to p do grad_theta_H_tot[i]:=0.0;
      for i:=1 to p do
      for j:=1 to p do grad_w_Left_OH_tot[i,j]:=0.0;
      for i:=1 to q do
      for j:=1 to p do grad_w_Right_OH_tot[i,j]:=0.0;
      for i:=1 to p do
      for j:=1 to p do grad_w_Left_HI_tot[i,j]:=0.0;
      for i:=1 to p do
      for j:=1 to q do grad_w_Right_HI_tot[i,j]:=0.0;

      for index_prime:=1 to lengthx do
      begin if index_prime=1 then
            for j:=1 to p do x_L[j]:=0.5 else
            for j:=1 to p do x_L[j]:=y[j];
            for j:=1 to q do
            x_R[j]:=token[sequence[index_prime],j];

            for j:=1 to p do
            begin sum:=theta_H[j];
                  for k:=1 to p do
                  sum:=sum+w_Left_HI[j,k]*x_L[k];
                  for k:=1 to q do
                  sum:=sum+w_Right_HI[j,k]*x_R[k];
                  y[j]:=1/(1+exp(-sum));
            end;
            for j:=1 to p do
            begin sum:=theta_Left_O[j];
                  for k:=1 to p do
                  sum:=sum+w_Left_OH[j,k]*y[k];
                  x_prime_L[j]:=1/(1+exp(-sum));
            end;
            for j:=1 to q do
            begin sum:=theta_Right_O[j];
                  for k:=1 to p do
                  sum:=sum+w_Right_OH[j,k]*y[k];
                  x_prime_R[j]:=1/(1+exp(-sum));
            end;

            for j:=1 to p do
            begin diff_L[j]:=diff_L[j]+abs(x_L[j]-x_prime_L[j]);
                  counter_L[j]:=counter_L[j]+1;
            end;
            for j:=1 to q do
            begin diff_R[j]:=diff_R[j]+abs(x_R[j]-x_prime_R[j]);
                  counter_R[j]:=counter_R[j]+1;
            end;

            for j:=1 to p do
            grad_theta_Left_O[j]:= (x_prime_L[j]-x_L[j])*
                                   x_prime_L[j]*(1-x_prime_L[j]);
            for j:=1 to q do
            grad_theta_Right_O[j]:=(x_prime_R[j]-x_R[j])*
                                   x_prime_R[j]*(1-x_prime_R[j]);
            for j:=1 to p do
            begin sum:=0.0;
                  for k:=1 to p do
                  sum:=sum+grad_theta_Left_O[k]*w_Left_OH[k,j];
                  for k:=1 to q do
                  sum:=sum+grad_theta_Right_O[k]*w_Right_OH[k,j];
                  grad_theta_H[j]:=sum*y[j]*(1-y[j]);
            end;

            for i:=1 to p do
            grad_theta_Left_O_tot[i]:=grad_theta_Left_O_tot[i]+
                                      grad_theta_Left_O[i];
            for i:=1 to q do
            grad_theta_Right_O_tot[i]:=grad_theta_Right_O_tot[i]+
                                       grad_theta_Right_O[i];
            for i:=1 to p do
            grad_theta_H_tot[i]:=grad_theta_H_tot[i]+
                                 grad_theta_H[i];
            for i:=1 to p do
            for j:=1 to p do
            grad_w_Left_OH_tot[i,j]:=grad_w_Left_OH_tot[i,j]+
                                     grad_theta_Left_O[i]*y[j];
            for i:=1 to q do
            for j:=1 to p do
            grad_w_Right_OH_tot[i,j]:=grad_w_Right_OH_tot[i,j]+
                                      grad_theta_Right_O[i]*y[j];
            for i:=1 to p do
            for j:=1 to p do
            grad_w_Left_HI_tot[i,j]:=grad_w_Left_HI_tot[i,j]+
                                     grad_theta_H[i]*x_L[j];
            for i:=1 to p do
            for j:=1 to q do
            grad_w_Right_HI_tot[i,j]:=grad_w_Right_HI_tot[i,j]+
                                      grad_theta_H[i]*x_R[j];
      end;
      for i:=1 to p do y_code[i]:=y[i];
      {for j:=1 to p do
      diff_L[j]:=diff_L[j]/counter_L[j];
      for j:=1 to q do
      diff_R[j]:=diff_R[j]/counter_R[j];}
end;



procedure Coefficient_update(var delta_theta_Left_O    : vector_Left;
                             var delta_theta_Right_O   : vector_Right;
                             var delta_theta_H         : vector_Hidden;
                             var delta_w_Left_OH       : matrix_Left_OH;
                             var delta_w_Right_OH      : matrix_Right_OH;
                             var delta_w_Left_HI       : matrix_Left_HI;
                             var delta_w_Right_HI      : matrix_Right_HI;
                                 grad_theta_Left_O_tot : vector_Left;
                                 grad_theta_Right_O_tot: vector_Right;
                                 grad_theta_H_tot      : vector_Hidden;
                                 grad_w_Left_OH_tot    : matrix_Left_OH;
                                 grad_w_Right_OH_tot   : matrix_Right_OH;
                                 grad_w_Left_HI_tot    : matrix_Left_HI;
                                 grad_w_Right_HI_tot   : matrix_Right_HI;
                             var theta_Left_O          : vector_Left;
                             var theta_Right_O         : vector_Right;
                             var theta_H               : vector_Hidden;
                             var w_Left_OH             : matrix_Left_OH;
                             var w_Right_OH            : matrix_Right_OH;
                             var w_Left_HI             : matrix_Left_HI;
                             var w_Right_HI            : matrix_Right_HI);
var j,k : integer;
begin for j:=1 to p do
      delta_theta_Left_O[j]:=-lambda*grad_theta_Left_O_tot[j]+
                              mu*delta_theta_Left_O[j];
      for j:=1 to q do
      delta_theta_Right_O[j]:=-lambda*grad_theta_Right_O_tot[j]+
                               mu*delta_theta_Right_O[j];
      for j:=1 to p do
      delta_theta_H[j]:=-lambda*grad_theta_H_tot[j]+
                         mu*delta_theta_H[j];
      for j:=1 to p do
      for k:=1 to p do
      delta_w_Left_OH[j,k]:=-lambda*grad_w_Left_OH_tot[j,k]+
                             mu*delta_w_Left_OH[j,k];
      for j:=1 to q do
      for k:=1 to p do
      delta_w_Right_OH[j,k]:=-lambda*grad_w_Right_OH_tot[j,k]+
                              mu*delta_w_Right_OH[j,k];
      for j:=1 to p do
      for k:=1 to p do
      delta_w_Left_HI[j,k]:=-lambda*grad_w_Left_HI_tot[j,k]+
                             mu*delta_w_Left_HI[j,k];
      for j:=1 to p do
      for k:=1 to q do
      delta_w_Right_HI[j,k]:=-lambda*grad_w_Right_HI_tot[j,k]+
                              mu*delta_w_Right_HI[j,k];

      for j:=1 to p do
      theta_Left_O[j]:=theta_Left_O[j]+delta_theta_Left_O[j];
      for j:=1 to q do
      theta_Right_O[j]:=theta_Right_O[j]+delta_theta_Right_O[j];
      for j:=1 to p do
      theta_H[j]:=theta_H[j]+delta_theta_H[j];
      for j:=1 to p do
      for k:=1 to p do
      w_Left_OH[j,k]:=w_Left_OH[j,k]+delta_w_Left_OH[j,k];
      for j:=1 to q do
      for k:=1 to p do
      w_Right_OH[j,k]:=w_Right_OH[j,k]+delta_w_Right_OH[j,k];
      for j:=1 to p do
      for k:=1 to p do
      w_Left_HI[j,k]:=w_Left_HI[j,k]+delta_w_Left_HI[j,k];
      for j:=1 to p do
      for k:=1 to q do
      w_Right_HI[j,k]:=w_Right_HI[j,k]+delta_w_Right_HI[j,k];
end;



procedure create_string(     y_ini         : vector_hidden;
                             theta_Left_O  : vector_Left;
                             theta_Right_O : vector_Right;
                             w_Left_OH     : matrix_Left_OH;
                             w_Right_OH    : matrix_Right_OH;
                         var length        : integer;
                         var sequence      : vector_string);
var index,i,j   : integer;
    long        : integer;
    sum         : real;
    dist_min    : real;
    iterat_stop : Boolean;
    alpha       : vector_Right;
    y,y_aux     : vector_Hidden;

label 1;

function index_token(alpha : vector_Right) : integer;
var i,j       : integer;
    tok       : vector_string;
    finishing : Boolean;
begin for i:=1 to q do
      if alpha[i]<0.2 then tok[i]:=0 else
      if alpha[i]>0.8 then tok[i]:=1 else tok[i]:=2;
      i:=0;
      finishing:=true;
      while (i<N_token) and (finishing) do
      begin i:=i+1;
            sum:=0;
            for j:=1 to q do
            sum:=sum+abs(tok[j]-token[i,j]);
            if sum=0 then finishing:=false;
            index_token:=i;
      end;
      if finishing then index_token:=3;
end;

begin for i:=1 to p do y[i]:=y_ini[i];
      length:=0;
      for index:=1 to length_max do
      begin for i:=1 to p do
            begin sum:=theta_Left_O[i];
                  for j:=1 to p do
                  sum:=sum+w_Left_OH[i,j]*y[j];
                  y_aux[i]:=1/(1+exp(-sum));
            end;
            for i:=1 to q do
            begin sum:=theta_Right_O[i];
                  for j:=1 to p do
                  sum:=sum+w_Right_OH[i,j]*y[j];
                  alpha[i]:=1/(1+exp(-sum));
            end;
            j:=index_token(alpha);
            if j=3 then goto 1;
            sequence[index]:=j;
            length:=index;
            for i:=1 to p do y[i]:=y_aux[i];
      end;
1:end;



procedure Print_check;
var i,j,pq,length_max : integer;
    length            : integer;
    index,index_prime : integer;
    y_req,y_ini       : vector_Hidden;
    sequence          : vector_string;
begin writeln('epoch=',time);
      if p<q then pq:=q else pq:=p;
      for i:=1 to pq do
      if (i<=p) and (i<=q) then
      writeln('# ',i:2,'  diff_L=',diff_L_tot[i]:5:3,'  diff_R=',
              diff_R_tot[i]:5:3) else
      if (i<=p) and (i>q) then
      writeln('# ',i:2,'  diff_L=',diff_L_tot[i]:5:3) else
      if (i>p) and (i<=q) then
      writeln('# ',i:2,'                diff_R=',diff_R_tot[i]:5:3);

      length_max:=0;
      for index:=1 to N_train do
      if Train_length[index]>length_max then
      length_max:=Train_length[index];


      for index:=1 to N_train do
      begin lengthx:=train_length[index];
            for i:=1 to lengthx do
            sequence[i]:=train_set[index,i];

            write('pattern # ',index:2,'   (');
            for i:=1 to lengthx do write(sequence[i]);
            write(')');
            for i:=1 to length_max-lengthx do write(' ');
            write('  y=(');
            for i:=1 to p do
            begin write(pop_y[index,i]:3:2);
                  y_ini[i]:=pop_y[index,i];
                  if i<p then write(',');
            end;
            write(')    alpha=(');
            create_string(y_ini,theta_Left_O,theta_Right_O,
                          w_Left_OH,w_Right_OH,length,sequence);
            for i:=1 to length do write(sequence[length-i+1]:1);

            writeln(')');

      end;
      writeln;
end;


begin assign(data_in,'data_in'); reset(data_in);
      Data_input;
      Ini_coefficient;

      time:=0;
      while time<epoch_max do
      begin time:=time+1;
            for i:=1 to p do diff_L_tot[i]:=0.0;
            for i:=1 to q do diff_R_tot[i]:=0.0;

            for index:=1 to N_train do
            begin Calcul_grad_tot
                  (index,
                   grad_theta_Left_O_tot,grad_theta_Right_O_tot,
                   grad_theta_H_tot,
                   grad_w_Left_OH_tot,grad_w_Right_OH_tot,
                   grad_w_Left_HI_tot,grad_w_Right_HI_tot,
                   theta_Left_O,theta_Right_O,theta_H,
                   w_Left_OH,w_Right_OH,
                   w_Left_HI,w_Right_HI,
                   diff_L,diff_R,
                   y_code);

                  for i:=1 to p do pop_y[index,i]:=y_code[i];

                  for i:=1 to p do diff_L_tot[i]:=diff_L_tot[i]+diff_L[i];
                  for i:=1 to q do diff_R_tot[i]:=diff_R_tot[i]+diff_R[i];

                  Coefficient_update
                  (delta_theta_Left_O,delta_theta_Right_O,
                   delta_theta_H,
                   delta_w_Left_OH,delta_w_Right_OH,
                   delta_w_Left_HI,delta_w_Right_HI,
                   grad_theta_Left_O_tot,grad_theta_Right_O_tot,
                   grad_theta_H_tot,
                   grad_w_Left_OH_tot,grad_w_Right_OH_tot,
                   grad_w_Left_HI_tot,grad_w_Right_HI_tot,
                   theta_Left_O,theta_Right_O,
                   theta_H,
                   w_Left_OH,w_Right_OH,
                   w_Left_HI,w_Right_HI);
            end;

            for i:=1 to p do diff_L_tot[i]:=diff_L_tot[i]/N_train;
            for i:=1 to q do diff_R_tot[i]:=diff_R_tot[i]/N_train;

            if (time=1) or (time mod print_freq = 0) then Print_check;
      end;
end.