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

double nonlinear(double x);
double inverse_nonlinear(double x);

double trainingfunction(double x, double y, double z);

int main()
{
#define NUM_INPUTS 3
#define NUM_NODES NUM_INPUTS
#define NUM_LAYERS 3
#define NUM_OUTPUTS 1

double inputs[NUM_INPUTS];
double weights[NUM_LAYERS][NUM_NODES][NUM_INPUTS];
double internal[NUM_LAYERS][NUM_NODES];
double output[NUM_OUTPUTS];
double outputweights[NUM_NODES][NUM_OUTPUTS];
double mu;
double traininput[NUM_INPUTS];
double trainoutput[NUM_OUTPUTS];
double outputerror[NUM_OUTPUTS];
double internalerror[NUM_LAYERS][NUM_NODES];
double collector;

int i,j,k,l;
mu=0.1;
srand(10);

// training part

for(l=0;l<1000;l++)
{
// select the first training vector pair of input and output
  for(i=0;i<NUM_INPUTS;i++)
   traininput[i]=0.1+rand()-0.1;
  trainoutput[0]=trainingfunction(traininput[0],traininput[1],traininput[2]); // HACK!!

// use the input vector as the output from the input layer
  for(i=0;i<NUM_INPUTS;i++)
   inputs[i]=traininput[i];
  for(i=0;i<NUM_INPUTS;i++)
   internal[0][i]=inputs[i];

// compute the activation of each unit on the next layer
  for(i=1;i<NUM_LAYERS;i++)
   for(j=0;j<NUM_NODES;j++)
   {
     internal[i][j]=0;
    for(k=0;k<NUM_INPUTS;k++)
     internal[i][j]+=weights[i][j][k]*internal[i-1][k];
    internal[i][j]=nonlinear(internal[i][j]);
   }

// mark down the outputs
  i=NUM_LAYERS;
  for(j=0;j<NUM_OUTPUTS;j++)
  {
    output[j]=0;
    for(k=0;k<NUM_INPUTS;k++)
     output[j]+=outputweights[j][k]*internal[i-1][k];
    output[j]=nonlinear(output[j]);
  } 

// calculate the error for the pattern across the output layer:
// error = (yk - ok)f'(netok)
  for(j=0;j<NUM_OUTPUTS;j++)
  {
    outputerror[j]=(trainoutput[j]-output[j])*inverse_nonlinear(output[j]);
  }

// computer the error for all hidden layer units:
// error h pj = f'(net h j) sum (output error)*weight(kj)
  for(i=NUM_LAYERS;i>-1;i--)
   for(j=0;j<NUM_NODES;j++)
   {
     collector=0;
     for(k=0;k<NUM_OUTPUTS;k++)
      collector += outputerror[k]*outputweights[j][k];
     internalerror[i-1][j]=inverse_nonlinear(internal[i-1][j])*collector;
   }

// update the connection weight values as follows:
// new_weight = old_weight + mu*internalerror*input
  for(i=0;i<NUM_LAYERS;i++)
   for(j=0;j<NUM_NODES;j++)
    for(k=0;k<NUM_NODES;k++)
     weights[i][j][k]+=mu*internalerror[i][j]*internal[i][j];
}

// calculation output part
for(l=0;l<10;l++)
{
// make up a random puzzle to play
  for(i=0;i<NUM_INPUTS;i++)
   traininput[i]=0.1+rand()-0.1;
  trainoutput[0]=trainingfunction(traininput[0],traininput[1],traininput[2]); // HACK!!

// use the input vector as the output from the input layer
  for(i=0;i<NUM_INPUTS;i++)
   inputs[i]=traininput[i];
  for(i=0;i<NUM_INPUTS;i++)
   internal[0][i]=inputs[i];

// compute the activation of each unit on the next layer
  for(i=1;i<NUM_LAYERS;i++)
   for(j=0;j<NUM_NODES;j++)
   {
     internal[i][j]=0;
    for(k=0;k<NUM_INPUTS;k++)
     internal[i][j]+=weights[i][j][k]*internal[i-1][k];
    internal[i][j]=nonlinear(internal[i][j]);
   }

// mark down the outputs
  i=NUM_LAYERS;
  for(j=0;j<NUM_OUTPUTS;j++)
  {
    output[j]=0;
    for(k=0;k<NUM_INPUTS;k++)
     output[j]+=outputweights[j][k]*internal[i-1][k];
    output[j]=nonlinear(output[j]);
  } 
 printf("input %d %d %d\n",traininput[0],traininput[1],traininput[2]);
 printf("output %d\n",output[0]);  // HAAAAACK
 printf("correct answer %d\n\n",trainoutput[0]); // HACK
}

 return(0);
}

double nonlinear(double x)
{
  if (x<-M_PI_2) return 0;
  if (x>M_PI_2) return 1;
  return (sin(x)+1)/2;
}

double inverse_nonlinear(double x)
{ 
  if (x<-M_PI_2) return 0;
  if (x>M_PI_2) return 0;
  return (cos(x))/2;
}
double trainingfunction(double x, double y, double z)
{
  return x+y+z;
}