/**********************************************************************
 * 2-bit XOR problem.  
 *
 *  gcc -o xor xor.c -lm
 *  xor 5 0.2 0.8
 *
 * The first argument is the number of hidden units, the second is
 * the learning constant and the third is the momentum constant.  So 
 * the example runs with 5 hidden units, a learning constant
 * of 0.2 and momentum of 0.8.
 *
 * Set the  debugging  #define to control the output you get.
 **********************************************************************/

#include <math.h>
#include <stdio.h>
#include <sys/time.h>		/* for setting random number seed */

#define MAXNH 100			/* maximum number of hidden units */
#define debugging 2	/* set to 1 to see debugging messages */
			/* set to 2 to print just errors (for plotting) */

#define RN ((float) random() / (float)((1 << 31) - 1))

/*** function prototypes ***/

void calc_output(float x[2], int nh, float wh[MAXNH][3], float wo[MAXNH+1], 
		 float yh[MAXNH], float *yo);
void update_weights(float x[2], int nh, float wh[MAXNH][3], float wo[MAXNH+1], 
		    float yh[MAXNH], float yo, float error, float rate, 
		    float mom, float dwh[MAXNH][3], float dwo[MAXNH+1]);
void print_nets(float x[2], int nh, float wh[MAXNH][3], float wo[MAXNH+1], 
		float yh[MAXNH], float yo, float error);
float logistic(float x);
float Dlogistic(float y);

/**********************************************************************/

main(int argc, char *argv[])
{
  float wh[MAXNH][3];		/* hidden units' weights */
  float wo[MAXNH+1];		/* output unit's weights */
  float x[2];			/* input vector */
  float yh[MAXNH];			/* hidden units' outputs */
  float yo;			/* output unit's output  */
  int nh = 5;			/* number of hidden units */
  float rate = 0.2;		/* learning rate with default value*/
  float mom = 0.8;		/* momentum rate with default value*/
  float init_weight_range = 1.0 /*0.1*/; /* range of initial weights (+- this value) */
  float dwh[MAXNH][3];		/* previous weight change, for momentum */
  float dwo[MAXNH+1];		/* previous weight change, for momentum */

  /* Define input and target patterns */

  int patterns[4][2] = {0,0, 0,1, 1,0, 1,1};
  int targets[4] = {0, 1, 1, 0};

  int pattern, target, ncorrect, epoch, u, i;
  float error, epoch_error;

  /* Initialize the pseudo-random number generator by current time. */

  srandom(time(NULL));

  /* If arguments present, process them */

  if (argc > 1)
    nh = atoi(argv[1]);

  if (argc > 2)
    rate = atof(argv[2]);

  if (argc > 3)
    mom = atof(argv[3]);
  
  /* Initialize weights to small random values. */

  for (u=0; u<nh; u++)
    for (i=0; i<3; i++)
      wh[u][i] = init_weight_range * (RN * 2 - 1);

  for (i=0; i<nh+1; i++)
    wo[i] = init_weight_range * (RN * 2 - 1);

  /* Initialize previous weight changes to zero. */

  for (u=0; u<nh; u++)
    for (i=0; i<3; i++)
      dwh[u][i] = 0.;

  for (i=0; i<nh+1; i++)
    dwo[i] = 0.;

  /* Main loop. Run for this many repetitions of training set. */

  printf("%d hidden units, learning rate is %.2f, momentum rate is %.2f\n",
	 nh,rate,mom);

  for (epoch=0; epoch<2000; epoch++) {

    ncorrect = 0;
    epoch_error = 0;

    /* Loop for each pattern in the training set. */

    for (pattern=0; pattern<4; pattern++) {

      /* Load current training pattern */
      for (i=0; i<2; i++)
	x[i] = patterns[pattern][i];
      target = targets[pattern];

      calc_output(x,nh,wh,wo,yh,&yo); /* returns yo, network output */

      error = target - yo;
      epoch_error += error*error;

      update_weights(x,nh,wh,wo,yh,yo,error,rate,mom,dwh,dwo);

      if (debugging == 1) {
	printf("Just updated all weights for epoch %d pattern %d:\n",
	       epoch,pattern);
	print_nets(x,nh,wh,wo,yh,yo,error);
      }
      if (debugging == 2)		/* for plotting per pattern error */
	printf("pattern %d e %f\n",pattern,fabs(error));

      /* if output is within 0.2, count this pattern as correct */

      if (fabs(error) < 0.2) ncorrect++; 

    }				/* end of pattern loop */

    if (debugging==2)		/* for plotting per epoch error */
	printf("epoch %d e %f\n",epoch,sqrt(epoch_error/4));

    /* Stop epoch loop once all four patterns are correct */

    if (ncorrect == 4) break;
  }				/* end of epoch loop */

  /* Now for a little user feedback. */

  if (ncorrect == 4)
    printf("Converged after %d epochs.\n",epoch);
  else
    printf("Didn't converge within %d epochs.\n",epoch);

  for (pattern=0; pattern<4; pattern++) {

    /* Load current training pattern */
    for (i=0; i<2; i++)
      x[i] = patterns[pattern][i];
    target = targets[pattern];

    calc_output(x,nh,wh,wo,yh,&yo);
    error = target - yo;
    print_nets(x,nh,wh,wo,yh,yo,error);
  }
}

/**********************************************************************/

void calc_output(float x[2], int nh, 
		 float wh[MAXNH][3], float wo[MAXNH+1], 
		 float yh[MAXNH], float *yo)
{
  int u, i;
  float sum;

  for (u=0; u<nh; u++) {
    sum = wh[u][0];		/* bias weight (as if input of 1) */
    for (i=0; i<2; i++) {
      sum += wh[u][i+1] * x[i];
    }
    yh[u] = logistic(sum);
  }

  sum = wo[0];			/* bias weight (ditto) */
  for (i=0; i<nh; i++)
    sum += wo[i+1] * yh[i];
  *yo = logistic(sum);
}

/**********************************************************************/

void update_weights(float x[2], int nh, 
		    float wh[MAXNH][3], float wo[MAXNH+1], 
		    float yh[MAXNH], float yo, float error, 
		    float rate, float mom, 
		    float dwh[MAXNH][3], float dwo[MAXNH+1])
{
  float deltao;			/* delta value for output unit */
  float deltah[MAXNH];		/* delta values for hidden units */
  int u, i;

  /* Deltas */

  deltao = error * Dlogistic(yo);

  for (u=0; u<nh; u++)
    deltah[u] = deltao * wo[u+1] * Dlogistic(yh[u]);

  for (u=0; u<nh; u++) {
    wh[u][0] += dwh[u][0] = rate * deltah[u] + mom * dwh[u][0];
    for (i=0; i<2; i++) {
      wh[u][i+1] += dwh[u][i+1] = rate * deltah[u] * x[i] + mom * dwh[u][i+1];
    }
  }

  wo[0] += dwo[0] = rate * deltao  + mom * dwo[0];
  for (i=0; i<nh; i++) 
    wo[i+1] += dwo[i+1] = rate * deltao * yh[i] + mom * dwo[i+1];
}

/**********************************************************************/

void print_nets(float x[2], int nh, 
		float wh[MAXNH][3], float wo[MAXNH+1], 
		float yh[MAXNH], float yo, float error)
{
  int u, i;

  printf("\n ");
  for (u=0; u<nh; u++)
    printf("%6.2f",wh[u][0]);
  for (i=0; i<2; i++) {
    printf("\n%d",(int)x[i]);
    for (u=0; u<nh; u++)
      printf("%6.2f",wh[u][i+1]);
  }

  printf("\n ");
  for (u=0; u<nh; u++)
    printf("%6.2f",yh[u]);
  printf("\n ");
  for (u=0; u<nh; u++)
    printf("%6.2f",wo[u+1]);
  printf("%6.2f  %6.2f  (e=%6.2f)\n",wo[0],yo,error);
}

/**********************************************************************/

float logistic(float x)
{
  return (2. / (1. + exp(-x)) - 1);

  /** return (1. / (1. + exp(-x))); **/    /* For 0-1 logistic */
}

/**********************************************************************/

float Dlogistic(float y)
{
  return (1 - y*y);

  /** return (y * (1. - y)); **/    /* For 0-1 logistic */

}