#include <stdio.h>
#include <math.h>

#define MXUNIT  16
#define MXTRSMPL 1000

float 
  eta,
  alpha,
  beta,
  ts,
  sys_err,
  maxe,
  wt[MXUNIT][MXUNIT],                   /* weight matrix            */
  dwt[MXUNIT][MXUNIT],                  /* delta weight             */
  yy[MXUNIT],                           /* node outputs             */
  oldppp[MXUNIT][MXUNIT][MXUNIT],       /* dynamics at time t       */
  newppp[MXUNIT][MXUNIT][MXUNIT],       /* dynamics at time t+1     */
  zz[MXUNIT],                           /* node output + ext. input */
  fs[MXUNIT],                           /* activation               */
  inp[MXTRSMPL][MXUNIT],                /* training inputs          */
  outpt[MXTRSMPL][MXUNIT],              /* computed outputs         */
  target[MXTRSMPL][MXUNIT],             /* target outputs           */
  ee[MXUNIT],                           /* error                    */
  err[MXTRSMPL];

int
  nsamples,                             /* number of training samples */
  nunits,                               /* number of units          */
  ninput=0,                             /* number of inputs         */
  noutput=0,                            /* number of outputs        */
  nhidden,                              /* number of hidden nodes   */
  mxniter,                              /* counter for no. itertn   */
  tot_iter,
  outindx[MXUNIT],                      /* unit numbers for output  */
  nrow,ncol,                            /* number of rows and columns */
  ntsdata, 
  ntsunits,
  nunf;

char 
  net_file[20];                         /* name of network data file */

FILE 
  *fp1,                                 /* training sample file   */
  *fp2,                                 /* net data file to write */
  *fp3,                                 /* net data file to read  */
  *fp4,                                 /* test data file         */
  *fp5,                                 /* test output file       */
  *fp6;                                 /* training output file   */

void set_up_train()
{
  int i,j;
  char tr_data_file[20], tofile[20];

  printf("\nLearning rate eta?: ");
  scanf("%f",&eta);
   
  printf("\nNumber of iterations in each free running?: ");
  scanf("%d",&nunf);
   
  printf("\nMax total error?: ");
  scanf("%f",&maxe);
   
  printf("\nMax number of iterations?: ");
  scanf("%d",&mxniter);

  printf("\nNumber of Units?: ");
  scanf("%d",&nunits);
  
  printf("\nNumber of Iutput Units?: ");
  scanf("%d",&ninput);
  
  printf("\nNumber of Output Units?: ");
  scanf("%d",&noutput);

  printf("\nEnter Unit Numbers for Output Units:");
  for (i=0; i<noutput; i++) {
    printf("\nUnit Number for %d th Output :", i);
    scanf("%d",&outindx[i]);
  }

  printf("\nTotal number of training samples?: ");
  scanf("%d",&nsamples);

  printf("\nEnter the Name of Training Data File: ");
  scanf("%s",tr_data_file);

  if ((fp1 = fopen(tr_data_file,"r")) == NULL) {
      perror("\n Cannot open data file ");
      exit(0);
    }

  for (i=0; i<nsamples; i++) {
    for (j=0; j<ninput; j++)
      fscanf(fp1,"%f",&inp[i][j]);
    for (j=0; j<noutput; j++)
      fscanf(fp1,"%f",&target[i][j]);
  }

  printf("\nEnter the Name of Net Data File: ");
  scanf("%s",net_file);

  if ((fp2 = fopen(net_file,"w")) == NULL) {
      perror("\n Cannot open data file ");
      exit(0);
    }

  printf("\nEnter the name of training output file: ");
  scanf("%s",tofile);

  if ((fp6 = fopen(tofile,"w")) == NULL) {
      perror("\n Cannot open data file ");
      exit(0);
    }

} /* end of SET_UP */

void initialize()
{
  int i,j,k;
  
  i=j=k=0;
  nrow=nunits;
  ncol=nunits+ninput+1;
  nhidden=nunits-noutput;

  srandom(time(0));
  for (i=0; i<nrow; i++) {
    for (j=0; j<ncol; j++) {
      wt[i][j]=random()/pow(2.0,31.0);
      dwt[i][j]=0.0;
    }
    yy[i]=0.1;
    ee[i]=0.0;
  }

  for (i=0; i<nrow; i++) 
    for (j=0; j<nrow; j++) 
      for (k=0; k<ncol; k++) 
	oldppp[i][j][k]=0.0;

} /* end of INITIALIZE */

void forward(ii)
     int ii;
{
  int   i,j,k;
  float net;

  zz[0]=1.0;
  for (i=1; i<=ninput; i++) {
    zz[i]=inp[ii][i-1];
    k=i;
  }
  k++;
  for (j=0; j<nunits; j++)
    zz[k+j] = yy[j];

  for (i=0; i<nrow; i++) {
    net=0.0;
    if (i != noutput-1) {
      for (j=0; j<ncol; j++)
	net += wt[i][j] * zz[j];
      yy[i] = 1/(1+exp(-net));
    }
    else {
      for (j=0; j<ncol; j++)
	net += wt[i][j] * yy[j];
      yy[i] = 1/(1+exp(-net));
    }
  }
  
  for (i=0; i<nunits; i++)
    outpt[ii][i]=yy[i];
} /* end of FORWARD */


void williams()
{
  int   h,i=0,j,k,m,p,q,r;
  float netp;

  do {

    sys_err=0.0;
    for (j=0; j<nsamples; j++) {            /* for each sample */
      for (m=0; m<nunf; m++) {

	forward(j);                         /* compute outputs */

	err[j]=0.0;
	for (k=0; k<noutput; k++) {            /* compute error   */
	  ee[outindx[k]]=target[j][k]-outpt[j][outindx[k]];
	  err[j] += ee[outindx[k]]*ee[outindx[k]];
	}
/*	printf("%f ",err[j]); */

/*	err[j] /= 2.0; */                      /* WHY? <---------------- */
	sys_err += err[j];

	for (k=0; k<nunits; k++) {           /* compute p(k,(i,j)) */
	  fs[k]=yy[k]*(1-yy[k]);
	  for (p=0; p<nrow; p++)
	    for (q=0; q<ncol; q++) {
	      netp=0.0;
	      for (r=ninput+1; r<ncol; r++)
		  netp += oldppp[r][p][q] * wt[k][r];  
	      newppp[k+ninput+1][p][q] = netp;
	      if (k==p)
		newppp[k+ninput+1][p][q] += zz[q];

	      newppp[k+ninput+1][p][q] *= fs[k];
	      oldppp[k+ninput+1][p][q] = newppp[k+ninput+1][p][q];
	    }
	}
	for (p=0; p<nrow; p++)             /* adjust w[i,j] */
	  for (q=0; q<ncol; q++) {
	    dwt[p][q]=0.0;
	    for (r=ninput+1; r<ncol; r++)
	      dwt[p][q] += newppp[r][p][q] * ee[r-ninput-1];
	    wt[p][q] += eta * dwt[p][q];
	  }
      }
    }
    
    for (m=0; m<nunits; m++)
      yy[m] = 0.01;

    
    sys_err /= nsamples;
    i++;
    printf("%d============ Mean Square Error : %f\n",i,sys_err);
    
  } while ((i < mxniter) && (sys_err > maxe));

  tot_iter=i;

} /* end of WILLIAMS */


void wrap_up_train()
{
  int   i,j;

  char  key;

    /* save weights in a file for later use */

  if (( fp2 = fopen(net_file,"w")) == NULL) {
      perror("\n Cannot open network data file ");
      exit(0);
    }

  fprintf(fp2,"%d %d %d %d\n",nunits,ninput,noutput,nhidden);

  for (i=0; i<noutput; i++) 
    fprintf(fp2,"%d ",outindx[i]);
  fprintf(fp2,"\n");
  
  for (i=0; i<nrow; i++) {
    for (j=0; j<ncol; j++)
      fprintf(fp2, "%f ",wt[i][j]);
    fprintf(fp2, "\n");
  }

     /* store in tofile the output generated in training phase */

  for (i=0; i<nsamples; i++) {
    for (j=0; j<noutput; j++)
      fprintf(fp6,"%f ", outpt[i][j]);
    fprintf(fp6,"\n");
  }

     /* print out the output from the net and target output */

  printf("\nTraining has been finished with the system error : ");
  printf("%f\n",sys_err);
  fprintf(fp2,"%f\n",sys_err);
  fprintf(fp2,"%d\n",tot_iter);

  printf("\n   Target Output           Computed Output      \n");
  for (i=0; i<nsamples; i++) {
    for (j=0; j<noutput; j++)
      printf("   %f                %f\n",target[i][j],outpt[i][j]);
    printf("\n");
    if ((i % 10) == 9) {
      printf("\nContiue to next outputs? ");
      do {
	printf("Waiting for you to type 'y'... ");
	scanf("%c",&key);
      }	while (key != 'y');
    }
  }
  close(fp1);
  close(fp2);
  close(fp6);
}

void set_up_test()
{
  int i,j;

  char
    testfile[20],                         /* test data file name  */
    gdfile[20];

  printf("\nNumber of test data?: ");
  scanf("%d",&ntsdata);

  printf("\nEnter the name of the file ");
  printf("that contains the network informations: ");
  scanf("%s",net_file);
  
  if (( fp3 = fopen(net_file,"r")) == NULL) {
      perror("\n Cannot open network data file ");
      exit(0);
    }

  fscanf(fp3,"%d%d%d%d",&nunits,&ninput,&noutput,&nhidden);

  for (i=0; i<noutput; i++) 
    fscanf(fp3,"%d",&outindx[i]);

  nrow=nunits;
  ncol=nunits+ninput+1;

  for (i=0; i<nrow; i++) 
    for (j=0; j<ncol; j++) 
      fscanf(fp3,"%f",&wt[i][j]);

  printf("\nEnter the name of the test data file: ");
  scanf("%s",testfile);
  if (( fp4 = fopen(testfile,"r")) == NULL) {
      perror("\n Cannot open test data file ");
      exit(0);
    }
  
  for(i=0; i<ntsdata; i++) {
    for (j=0; j<ninput; j++)
      fscanf(fp4,"%f",&inp[i][j]);
    for (j=0; j<noutput; j++)
      fscanf(fp4,"%f",&target[i][j]);
  }

  printf("\nEnter name of file where generated output will be stored: ");
  scanf("%s",gdfile);
  if (( fp5 = fopen(gdfile,"w")) == NULL) {
      perror("\n Cannot open output data file ");
      exit(0);
    }
  
  
} /* end of SET_UP_TEST */

void testforward(ii)
     int ii;
{
  int   i,j,k;
  float net;

  zz[0]=1.0;
  for (i=1; i<=ninput; i++) {
    zz[i]=inp[ii][i-1];
    k=i;
  }
  k++;
  for (j=0; j<nunits; j++)
    zz[k+j]=yy[j];

  for (i=0; i<nrow; i++) {
    net=0.0;
    for (j=0; j<ncol; j++)
      net += wt[i][j] * zz[j];
    yy[i] = 1/(1+exp(-net));
  }
  
  for (i=0; i<nunits; i++)
    outpt[ii][i]=yy[i];
} /* end of TESTFORWARD */



void testnet()
{
  int i,j,k;

  for (i=0; i<nrow; i++) 
    yy[i]=0.1;

  for (i=0; i<ntsdata; i++) {
    testforward(i);

  /* what should be the initial state of output units? */
  /* whether same forward function can be used depend on this. */
  /* no iterrations in test also? */
  
    for (j=0; j<noutput; j++) {
      printf("%f ",target[i][j]);
      fprintf(fp5,"%f ",target[i][j]);
    }
    printf("     ");
    fprintf(fp5, "     ");
    for (j=0; j<noutput; j++) {
      printf("%f ",outpt[i][j]);
      fprintf(fp5,"%f ",outpt[i][j]);
    }
    printf("\n");
    fprintf(fp5,"\n");
  }
  printf("\nTest finished and gererated output stored...\n");
}
  
main()
{
  char  ans[5];

  printf("\n\n... Starting Recurrent Net Simulator ...\n\n");
  printf("Begin with Training Session? (y/n): ");
  scanf("%s",ans);
  if (ans[0] == 'y') {
    set_up_train();
    initialize();
    williams();
    wrap_up_train();
  }
  printf("\n\nReady To Begin The Testing Session...\n");
  printf("Start It Now or Quit? (s/q): ");
  scanf("%s",ans);
  if (ans[0] == 's') {
    set_up_test();
    testnet();
    close(fp3);
    close(fp4);
    close(fp5);
  }
  else
    printf("\n\nAll finished...\n");
      
}