/***************************************************************************
                            -  description
                             -------------------
    begin                : Sun Mar 12 2000
    copyright            : (C) 2000 by 
    email                : 
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
/***************************************************************************
                          Backpropagation.cpp  -  description
                             -------------------
    begin                : Sun Mar 12 2000
    copyright            : (C) 2000 by 
    email                : 
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
/***************************************************************************

                          Backpropagation.cpp  -  description
                             -------------------
    begin                : Fri Mar 10 2000
    copyright            : (C) 2000 by 
    email                : 
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
/* <c>1998 Thomas Hasenknopf */
/* Listing zu "Künstliche Intelligenz, Teil 3 */
/* Elektronik, Heft 7/2000 */

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

#define MAX_NEURON_INPUTS 64
#define MAX_NEURON_PER_LAYER 512
#define MAX_LAYER_PER_NET 16
#define MAX_EPOCH 10000

typedef struct NF{          /* Definition eines Neurons */
  int NumberInputs;         /* Anzahl der Neuroneneingänge */
  double w[MAX_NEURON_INPUTS];     /* Gewichtsfaktoren (max. 16) */
  double dw[MAX_NEURON_INPUTS];    /* letzte Gewichtsaenderungen */
  double b;                 /* Schwellwert */
  double db;                /* letzte Schwellwertaenderung */
  double n;                 /* Neuronenaktivierung */  
  double a;                 /* Ausgang des Neurons */
  double ebp;               /* Rückvermittlungs-Fehler */
  double (*F)(double);      /* Zeiger auf Aktivierungsfunktion */
  double (*Fd)(struct NF*); /* Zeiger auf Ableitungsfunktion */
} Neuron;

typedef struct {            /* Definition einer Neuronenschicht */
  int NumberNeurons;        /* Anzahl der Neuronen in der Schicht */
  Neuron N[MAX_NEURON_PER_LAYER];             /* Neuronen der Schicht (max. 16) */
} NeuronLayer;

typedef struct {            /* Definition eines Neuronennetzwerks */
  int NumberLayers;         /* Anzahl der Schichten im Netzwerk */
  NeuronLayer L[MAX_LAYER_PER_NET];         /* Schichten im Netzwerk (max. 8) */
} NeuronNet;

typedef struct {            /* Trainingsergebnis */
  double sse;               /* Summierter quadratischer Fehler */
  unsigned long epochs;     /* Trainingsepochen */
  int overflow;             /* Ziel des Trainings verfehlt? */
} TrainingResult;

double LinearFunction(double x) {
  return x;
}

double LinearDerivative(Neuron *N) {
  return 1;
}

double TanhDerivative(Neuron *N) {
 return (1.0-pow(N->a,2));
}

/* Initialisieren einer Neuronenschicht mit Zufallszahlen */
void InitLayer(int NumberInputs, int NumberNeurons, NeuronLayer *L,
	       double(*ActFct)(double), double(*DerivFct)(Neuron*)) {
  int i,n;
  L->NumberNeurons=NumberNeurons;
  srand(1);                              /* Init. des Generators */
  for(n=0;n<NumberNeurons;n++) {            /* Init. der Neuronen */
    L->N[n].NumberInputs=NumberInputs;
    L->N[n].F=ActFct;                       /* Zuweisung von F */  
    L->N[n].Fd=DerivFct;                    /* Zuweisung von F' */
    for(i=0;i<NumberInputs;i++) {           /* Init. der Gewichte */
      L->N[n].w[i]=rand()*2.0/RAND_MAX-1.0;
    }
    L->N[n].b=rand()*2.0/RAND_MAX-1.0;      /* Init. des Schwellwertes */
  }
}
  Neuron N;

/* Simulation einer Neuronenschicht */
void SimuLayer(double *Input, NeuronLayer *L) {
  int m,i;
  double n;
  for(m=0;m<L->NumberNeurons;m++) {
    N=L->N[m];
    n=0.0;
    for(i=0;i<N.NumberInputs;i++) { /* Summieren d. gewichteten Eingänge */
      n+=Input[i]*N.w[i];
    }
    n+=N.b;                         /* Summieren der Schwellwerte */
    N.n=n;                          /* Neuronenaktivierung */ 
    N.a=N.F(n);                     /* Aufruf der Aktivierungsfunktion */
    L->N[m]=N;
  }
}

void DisplayNeuron(Neuron *N) {
  int i;
  for(i=0;i<N->NumberInputs;i++) {
    printf("w[%d] = %lf\n",i,N->w[i]); /* Ausgabe der Eingangsgewichte */
  }
  printf("b = %lf\n",N->b);           /* Ausgabe des Schwellwertes */
}

void DisplayNet(NeuronNet* Net) {
  int i,m;
  for(i=0;i<Net->NumberLayers;i++) {
    printf("\n>>>> Layer %d <<<<\n",i);
    for(m=0;m<Net->L[i].NumberNeurons;m++) {
      printf("Neuron %d \n",m);
      DisplayNeuron(&Net->L[i].N[m]);
    }
  }
}

void DisplayTrainingResult(TrainingResult Result) {
  printf("\nTraining %ld epochs, SSE=%e:\n",Result.epochs,Result.sse);
  if(Result.overflow) {
    printf("Training not succeeded. Change training parameters!\n");
  }
}
  Neuron *pthis;

TrainingResult TrainNet(int NumberTargets, double *Input, double *Target,
                 NeuronNet *Net, double Lrate, double Mf, double SseMax) {
  long layer,epoch,inputs,outputs,set,numSets,numLayers,numNeurons,i,m,n;
  double sse,inp[MAX_NEURON_INPUTS],err;
  TrainingResult result;
  numLayers=Net->NumberLayers;
  outputs=Net->L[numLayers-1].NumberNeurons;/* Anz. d. Netzwerkausgaenge*/
  numSets=NumberTargets/outputs;            /* Trainings-Wertepaare*/
  inputs=Net->L[0].N[0].NumberInputs;       /* Anz. d. Netzwerkeingaenge*/

  /* Backpropagation Algorithmus */
  for(epoch=0;epoch<MAX_EPOCH;epoch++)
  {
    sse=0.0;
    printf("Epoch %d\n",epoch);
	  for(set=0;set<numSets;set++)
    {
      /* Vorwaertsvermittlung */
      for(layer=0;layer<numLayers;layer++)
		  {
				if(layer==0)
        {
	  			SimuLayer(&Input[set*inputs],&Net->L[0]);
				}
				else
        {
	  			for(i=0;i<Net->L[layer-1].NumberNeurons;i++)
          {
	    			inp[i]= Net->L[layer-1].N[i].a;
          }
	  		  SimuLayer(&inp[0],&Net->L[layer]);
	      }
      }
      /* Fehlerrückvermittlung */ 
      for(layer=numLayers-1;layer>=0;layer--)
      {
		 		numNeurons=Net->L[layer].NumberNeurons;
				for(m=0;m<numNeurons;m++)
        {
			    pthis=&Net->L[layer].N[m];
	 			  if(layer== (numLayers-1))
          {       /* Ausgangsschicht */
	   			  err=Target[set*outputs+m] - pthis->a;
	   		    sse+=pow(err,2);
	 			  }
				  else
					{ /* verborgene Schicht */
	    			err=0;
	    			for(n=0;n<Net->L[layer+1].NumberNeurons;n++)
				  	{
	      			err+=Net->L[layer+1].N[n].w[m]*Net->L[layer+1].N[n].ebp;
	   				}
	  			}
          pthis->ebp=pthis->Fd(pthis)*err;
				}
      }
      /* Gewichte berechnen */
      for(layer=0;layer<numLayers;layer++)
		  {
				numNeurons=Net->L[layer].NumberNeurons;
				for(m=0;m<numNeurons;m++)
			  {
	  			pthis=&Net->L[layer].N[m];
	  			for(i=0;i<pthis->NumberInputs;i++)
				  {
            if(layer==0)
						{ /* Eingangsschicht */
	      			pthis->dw[i]=(1.0-Mf)*pthis->ebp*Input[set*inputs+i]*Lrate + Mf*pthis->dw[i];
	    			}
	   			  else
            {
	      			pthis->dw[i]=(1.0-Mf)*pthis->ebp*Net->L[layer-1].N[i].a*Lrate + Mf*pthis->dw[i];
	    			}
	   				pthis->w[i]+=pthis->dw[i];
	 			  }
	 			  pthis->db=(1.0-Mf)*pthis->ebp*Lrate + Mf*pthis->db;
	 			  pthis->b+=pthis->db;
			  }
      }
    }
    if(sse<=SseMax) break;
  }
  result.sse=sse;
  result.epochs=epoch;
  result.overflow=(epoch==MAX_EPOCH);
  return result;
}

double P[]={-1.0,
	    -0.5,
	     0.0,
	     0.2,
	     0.5,
	     1.0,
			 -0.6};
double T[]={ 1.0,
	     0.5,
	     0.0,
	     0.0,
	     0.3,
	     0.2,
			 0.1};

NeuronNet Net;

void main(void) {
  TrainingResult result;
  unsigned long epochs;
  /* Definition eines 3-lagigen Netzwerks mit 1 Eingang und 1 Ausgang,
     2 verborgenen nichtlin. Schichten und lin. Ausgangsschicht
  */
  Net.NumberLayers=3;
  /* 1 Input 4 Neuronen */
  InitLayer(1,4,&Net.L[0],tanh,TanhDerivative);
  /* 4 Input 2 Neuronen */
  InitLayer(4,2,&Net.L[1],tanh,TanhDerivative);
  /* 2 Inputs 1 Neuron */   
  InitLayer(2,1,&Net.L[2],LinearFunction,LinearDerivative);

  /* Lrate=0.7,  Momentfaktor=0.8,  SSE=0.005 */
  result=TrainNet(sizeof(T)/sizeof(double),P,T,&Net,0.7,0.8,0.005);
  DisplayNet(&Net);
  DisplayTrainingResult(result);
}