/************************************************************ 
;   Backpropagation with momentum                           * 
;   by Andres Perez-Uribe                                   *
;   Universidad del Valle, Cali, Colombia                   *
;   sep/93                                                  *
;                                                           *
;   Email : andres DOT perez-uribe AT heig-vd.ch            * 
;           REDS Institute @ HEIG-VD                        *
;           University of Applied Sciences of Western       *
            Switzerland (HES-SO)                            *
;************************************************************

   References :
   -  G. Hinton, "How neural networks learn from experience",
      Scientific American, sep 1992.
   -  P. Werbos,  "The Roots of Backpropagation: From ordered derivatives 
      to Neural Neworks and Political Forecasting", John Wiley and Sons, 
      New York, 1994

   Compile : gcc -o Bkprop Bkprop.c -lm
   Run     : see example at the end of the C code. 

    There is no guarantee that the code will do what you
    expect or that it is error free. It is simply meant
    to provide a useful way to experiment with the
    Backpropagation learning algorithm.

    Update Oct 7/99...thanks to Stephane Pouyet <pouyet@nist.gov>
    
    Last update March/2007 ... change of variable names for the sake of clarity
    The original names came from Hinton's paper.  
*/

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

#define sigm(x)    1/(1 + exp(-(double)x)) /* sigmoid function */
/* #define dsigm(x)   (float)(x)*(1.0-(float)x)) derivative of the sigmoid function */
#define IN         35    /* number of inputs */
#define HIDDEN     5     /* number of hidden units */
#define OUT        10    /* number of outputs */
#define EPSILON    0.05  /* maximum Mean Square Error to stop training */
#define NUMTRAIN   18    /* number of training Patterns */

float inhiddw[IN][HIDDEN];
float hidoutw[HIDDEN][OUT];
float deltaihw[IN][HIDDEN];
float deltahow[HIDDEN][OUT];
float x[IN];
float y[HIDDEN];
float z[OUT];

/* training patterns */
int   patterns[NUMTRAIN][IN] = { { 0,1,1,1,1,1,0,             
				  1,0,0,0,0,0,1,
				  1,0,0,0,0,0,1,
				  1,0,0,0,0,0,1,
				  0,1,1,1,1,1,0 },

				{ 0,0,0,0,0,0,0,
				  0,1,0,0,0,0,1,
				  1,1,1,1,1,1,1,
				  0,0,0,0,0,0,1,
				  0,0,0,0,0,0,0 },

				{ 0,1,0,0,0,0,1,
				  1,0,0,0,0,1,1,
				  1,0,0,0,1,0,1,
				  1,0,0,1,0,0,1,
				  0,1,1,0,0,0,1 },

				{ 1,0,0,0,0,1,0,
				  1,0,0,0,0,0,1,
				  1,0,0,1,0,0,1,
				  1,1,1,0,1,0,1,
				  1,0,0,0,1,1,0 },

				{ 0,0,0,1,1,0,0,
				  0,0,1,0,1,0,0,
				  0,1,0,0,1,0,0,
				  1,1,1,1,1,1,1,
				  0,0,0,0,1,0,0 },

				{ 1,1,1,0,0,1,0,
				  1,0,1,0,0,0,1,
				  1,0,1,0,0,0,1,
				  1,0,1,0,0,0,1,
				  1,0,0,1,1,1,0 },
				       
				{ 0,0,1,1,1,1,0,
				  0,1,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  0,0,0,0,1,1,0 },

				{ 1,0,0,0,0,0,0,             
				  1,0,0,0,0,0,0,
				  1,0,0,1,1,1,1,
				  1,0,1,0,0,0,0,
				  1,1,0,0,0,0,0 },

				{ 0,1,1,0,1,1,0,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  0,1,1,0,1,1,0 },

				{ 0,1,1,0,0,0,0,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,1,0,
				  0,1,1,1,1,0,0 },

				{ 1,1,1,1,0,0,0,   /* 4 */
				  0,0,0,1,0,0,0,
				  0,0,0,1,0,0,0,
				  0,0,0,1,0,0,0,
				  1,1,1,1,1,1,1 }, 

				{ 1,1,1,1,0,1,0,   /* 5 */
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,0,1,1,0 },

				{ 1,0,0,0,0,0,0,    /* 7 */  
				  1,0,0,0,0,0,0,
				  1,0,0,1,0,0,0,
				  1,1,1,1,1,1,1,
				  0,0,0,1,0,0,0 },

				{ 0,1,0,0,0,1,0,    /* 3 */
				  1,0,0,0,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,1,0,1,0,1,
				  0,1,1,0,1,1,0 },

				{ 1,0,0,0,0,1,1,   /* 2 */
				  1,0,0,0,1,0,1,
				  1,0,0,1,0,0,1,
				  1,0,1,0,0,0,1,
				  1,1,0,0,0,0,1 },

				{ 1,1,1,1,0,0,0,   /* 4 abierto */
				  0,0,0,1,0,0,0,
				  0,0,0,1,0,0,0,
				  1,1,1,1,1,1,1,
				  0,0,0,1,0,0,0 }, 

				{ 0,0,0,1,1,1,0,  /* 0 */
				  0,1,1,0,0,0,1,
				  1,0,0,0,0,0,1,
				  1,0,0,0,0,1,0,
				  1,1,1,1,1,0,0 },

				{ 0,1,1,0,0,0,1,     /* 9 */
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  1,0,0,1,0,0,1,
				  0,1,1,1,1,1,1 } };


/* desired outputs */
int   desout[NUMTRAIN][OUT] = { { 1,0,0,0,0,0,0,0,0,0 },
				{ 0,1,0,0,0,0,0,0,0,0 },
				{ 0,0,1,0,0,0,0,0,0,0 },
				{ 0,0,0,1,0,0,0,0,0,0 },
				{ 0,0,0,0,1,0,0,0,0,0 },
				{ 0,0,0,0,0,1,0,0,0,0 },
				{ 0,0,0,0,0,0,1,0,0,0 },
				{ 0,0,0,0,0,0,0,1,0,0 },
				{ 0,0,0,0,0,0,0,0,1,0 },
				{ 0,0,0,0,0,0,0,0,0,1 },
				{ 0,0,0,0,1,0,0,0,0,0 },
				{ 0,0,0,0,0,1,0,0,0,0 },
				{ 0,0,0,0,0,0,0,1,0,0 },
				{ 0,0,0,1,0,0,0,0,0,0 },
				{ 0,0,1,0,0,0,0,0,0,0 },
				{ 0,0,0,0,1,0,0,0,0,0 },
				{ 1,0,0,0,0,0,0,0,0,0 },
				{ 0,0,0,0,0,0,0,0,0,1 } };


float ehid[HIDDEN];
float eout[OUT];
int lrndpatr[NUMTRAIN]; /* learned pattern at a given time */
float sme[NUMTRAIN]; /* squared mean error */
float delta=0.5;     /* learning rate */
float alfa=0.1;      /* momentum */
long int itr;        /* iteration number */
int matrizin[35];    /* test input pattern */

void init();
void training(); 
void netanswer(int pattern[]); /* neural network output */ 
float se(int x[],float y[],int SIZE); /* squared error */
void backprop(int patternnum);
void error();

void init()
{
 int i,j;

  srand(time(NULL));
  for(i=0;i<IN;i++)
   for(j=0;j<HIDDEN;j++) {
    inhiddw[i][j] = -0.5 + ((float)rand() / RAND_MAX);
    deltaihw[i][j] = 0;
   }

  for(i=0;i<HIDDEN;i++)
   for(j=0;j<OUT;j++) {
    hidoutw[i][j] = -0.5 + ((float)rand() / RAND_MAX);
    deltahow[i][j] = 0;
   }

  for(i=0;i<NUMTRAIN;i++)
   lrndpatr[i] = 0;
   
  for(i=0;i<35;i++)
   matrizin[i]=0; 
}
  
void training()
{
 int i,l;           
 long int j;
 int t;
 int count;

 j=0;
 do {
  do {

/* select a random training pattern: i = (int)(NUMTRAIN*rnd), where 0<rnd<1 */
   i = (int)(NUMTRAIN * (float)rand() / RAND_MAX);
  } while(lrndpatr[i]);
   j++;
   netanswer(patterns[i]);   
   backprop(i);

  if(!(j%100)) /* show error every 100 training steps*/
   printf("\nAfter %ld steps ",j);
  error();
  l = 1;
  count=0;
  for(t=0;t<NUMTRAIN;t++) {
    if( sme[t] < EPSILON ) {
     lrndpatr[t] = 1;
     count++; 
    }
    else lrndpatr[t] = 0;
   l = l && (lrndpatr[t]);
  }
  if(!(j%100))
   printf("%d%% patterns learned",(int)(100 * (float)count/NUMTRAIN));
 } while(!l); /* iterate until sme of every pattern is less than the epsilon threshold */
   
 printf("\nAfter %ld steps ",j);
 printf("%d%% patterns learned",(int)(100 * (float)count/NUMTRAIN));
 printf("\n\n End of training\n");

}

void netanswer(int pattern[])
{
 int i,j;
 float totin;
   
 for(i=0;i<IN;i++)
  x[i] = (float)pattern[i];

 for(j=0;j<HIDDEN;j++) {
  totin = 0;
  for(i=0;i<IN;i++) 
   totin = totin + x[i]*inhiddw[i][j];
  y[j] = sigm(totin);
 }
  
 for(j=0;j<OUT;j++) {
  totin = 0;
  for(i=0;i<HIDDEN;i++)
   totin = totin + y[i]*hidoutw[i][j];
  z[j] = sigm(totin);
 }
}

float se(int a[],float b[],int SIZE)   /* Error measure */
{
 int i;
 float e=0;

 for(i=0;i<SIZE;i++)
  e = e + ((float)a[i] - b[i])*(a[i] - b[i]);
 e = 0.5 * e;
 return e;
}
  
void betaout(int patternnum)   /* error out */
{
 int k;
 for(k=0;k<OUT;k++)
  eout[k] = 0;

 for(k=0;k<OUT;k++)
  eout[k] = (z[k] - (float)desout[patternnum][k])*(z[k]*(1.0-z[k]));
}

void betahid()   /* error hidden */
{
 int j,k;
 for(j=0;j<HIDDEN;j++)
  ehid[j] = 0;

 for(j=0;j<HIDDEN;j++) 
  for(k=0;k<OUT;k++)
   ehid[j] = ehid[j] + hidoutw[j][k]*(y[j]*(1.0-y[j]))*eout[k];
}

void backprop(int patternnum)
{
 int i,j,k;
 float temp;

 betaout(patternnum); 
 betahid();
 
 for(j=0;j<HIDDEN;j++)
  for(k=0;k<OUT;k++) {
   temp = -delta*y[j]*eout[k];
   hidoutw[j][k] = hidoutw[j][k] + temp + alfa*deltahow[j][k];
   deltahow[j][k] = temp;
  }

 for(i=0;i<IN;i++)
  for(j=0;j<HIDDEN;j++) {
   temp = -delta*x[i]*ehid[j];
   inhiddw[i][j] = inhiddw[i][j] + temp + alfa*deltaihw[i][j];
   deltaihw[i][j] = temp;
  }
}   
void error()
{
 int i;

 for(i=0;i<NUMTRAIN;i++) {
  netanswer(patterns[i]);
  sme[i]=se(desout[i],z,OUT);
 }
}

void test()
{
 int i,j;

 for(;;) {
  printf("- Test -\n\n[");
  for(i=0;i<7;i++) { 
   for(j=0;j<5;j++) 
    scanf("%d",&matrizin[j*7+i]);
   printf("\n");
  }
  printf("]\n\n Output activations :\n");
  netanswer(matrizin);
  for(i=0;i<OUT;i++)
   printf("\nz[%d] = %f",i,z[i]);
  }
}

void main(int argc,char *argv[])
{
 init();
 training();
 test(); 
}

/*
Example :

gcc -o Bkprop Bkprop.c -lm

% ./Bkprop

102
204
306
408
510
612
714
816
918
1020
1122
1224
1326
1428
1530
1632
1734
1836
1938
2040
2142
2244
2346
2448
2550
2652
2754
2856
2958
3060
3162
3264
3366
3468
3570
3672
3774
3876
3978
4080

 End of training
- Test -

[0 0 1 0 0 

 0 0 1 0 0 

 0 0 1 0 0 

 0 0 0 0 0                            <-------- a '1' with some noise

 0 0 1 0 0 

 0 0 1 0 0 

 0 1 1 1 0

]

 Output activations :

z[0] = 0.073368
z[1] = 0.606160                        <------- the highest activation
z[2] = 0.101022
z[3] = 0.017971
z[4] = 0.101509
z[5] = 0.000393
z[6] = 0.014482
z[7] = 0.212412
z[8] = 0.003177
z[9] = 0.006917- Test -

*/