/*===========================================================================*/
/*                      Feedforward neural network                           */
/*                           D.J.Cowan, 2001.                                */
/*===========================================================================*/
/* A set of platform independant routines to implement a two tier feed-
	forward neural net with biased threshold and momentum terms. This file to
   be comiled as object code and linked with calling funnction.              */

/*Standard C header files*/
#include <stdio.h>                                         /*printf() scanf()*/
#include <stdlib.h>                                        /*rand() exit()   */
#include <math.h>                                          /*pow() exp()     */
#include <time.h>                                          /*clock()         */
#define BIAS 1
/*===========================================================================*/

enum {FALSE, TRUE};

/*Define data structure of network.*/
typedef struct fnet
{
	int    i                                                    /*architecture*/
	int    hid;
	int    out;
	float  error;
	float  rate;
	float  mom;
	float* p_in;                                      /*pointers to the arrays*/
	float* p_hid;
	float* p_out;
	float* p_targ;
	float* p_hid_w;
	float* p_out_w;
	float* p_out_mom;
	float* p_hid_mom;
} FFNET;

/*===========================================================================*/
/*Function declarations                                                      */

/*Public functions, methods*/
FFNET* init_ffnet(int, int, int, float, float);
FFNET* open_ffnet(char*);
FFNET* morph_ffnet(int*, int, int, int, int, int);
int    save_ffnet(FFNET*, char*);
int    del_ffnet(FFNET*);
float  train_ffnet(FFNET*, float*, float*);
float* test_ffnet(FFNET*, float*);
int    print_ffnet(FFNET*);

/*Private neural engine and utility funcions*/
static void   forwards(float*, float*, float*, int, int);
static void   forward_pass(FFNET*);
static void   backward_pass(FFNET*);
static float  sum_squared_error(FFNET*);
static float  sum_dEdxi(FFNET*, int);
static float  delta_w(FFNET*, int, int);
static float  delta_v(FFNET*, float, int, int);
static float* init(int, int);
static float  ff_bin_dec(int*, int);
static void   print_values(float*, int);
static void   fferr(int);

/*===========================================================================*/
/*Public interface functions                                                 */

/*Initialise net and all global variables, return true*/
FFNET* init_ffnet(int inputs, int hidden, int outputs, float r, float m)
{
   	FFNET* p_this = (FFNET*)malloc(sizeof(FFNET));
   	if (p_this == NULL)
   	{                                         /*error checking*/
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->in   = inputs;
   	p_this->hid  = hidden;                              /*transcribe architecture*/
   	p_this->out  = outputs;
   	p_this->rate = r;
   	p_this->mom  = m;
   	p_this->p_in = init(inputs, 0);                     /*initialise arrays with*/

   	if (p_this->p_in == NULL)
   	{                                 /*...error checks*/
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_hid = init(hidden, 0);
   	if (p_this->p_hid == NULL)
   	{
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_out = init(outputs, 0);
   	if (p_this->p_out == NULL)
   	{
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_targ = init(outputs, 0);
   	if (p_this->p_targ == NULL)
   	{
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_hid_w = init((inputs+BIAS) * hidden, 1);
   	if (p_this->p_hid_w == NULL){
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_out_w = init((hidden+BIAS) * outputs, 1);
   	if (p_this->p_out_w == NULL)
   	{
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_hid_mom = init((inputs+BIAS) * hidden, 0);
   	if (p_this->p_hid_mom == NULL)
   	{
   		fferr(1);
      	return(FALSE);
   	}

   	p_this->p_out_mom = init((hidden+BIAS) * outputs, 0);
   	if (p_this->p_out_mom == NULL)
   	{
   		fferr(1);
    	return(FALSE);
   	}

   	return(p_this);
}

/*Take a pointer to a filename, open and write net contents*/
int save_ffnet(FFNET* p_this, char* filename)
{
   	FILE* out_file;                                              /*open stream*/

   	if(p_this == NULL)
   	{
   		fferr(7);
      	return(FALSE);
   	}

   	out_file=fopen(filename, "wb");
   	if (out_file == NULL)
   	{
   		fferr(3);
      	return(FALSE);
   	}

   	fputc(p_this->in, out_file);
   	fputc(p_this->hid, out_file);
   	fputc(p_this->out, out_file);
   	p_this->rate*=1000;                              /*rescue to 4dp befor....*/
   	fputc(p_this->rate, out_file);                   /*..implicit cast to char*/
   	p_this->mom*=1000;
   	fputc(p_this->mom, out_file);                        /*write binary chunks*/
   	fwrite(p_this->p_in, sizeof(float), p_this->in, out_file);
   	fwrite(p_this->p_hid, sizeof(float), p_this->hid, out_file);
   	fwrite(p_this->p_out, sizeof(float), p_this->out, out_file);
   	fwrite(p_this->p_hid_w, sizeof(float), (p_this->in+1)*p_this->hid, out_file);
   	fwrite(p_this->p_out_w, sizeof(float), (p_this->hid+1)*p_this->out, out_file);
   	fwrite(p_this->p_hid_mom, sizeof(float), (p_this->in+1)*p_this->hid, out_file);
   	fwrite(p_this->p_out_mom, sizeof(float), (p_this->hid+1)*p_this->out, out_file);
   	fclose(out_file);

   	return(TRUE);
}

/*Take a pointer to a filename, open and initialise new from file
  contents returning pointer to FFNET*/
FFNET* open_ffnet(char* filename)
{
	int in, hid, out, rate, mom;
   	FILE * in_file;
	FFNET* p_this;

   	in_file = fopen(filename, "rb");                            /*open as binary*/
   	if (in_file == NULL)
   	{
   		fferr(4);
      	return(FALSE);
   	}
   	                                               /*read in architecture*/
   	in    = (int)fgetc(in_file);
   	hid   = (int)fgetc(in_file);
   	out   = (int)fgetc(in_file);
   	rate  = fgetc(in_file);
   	mom   = fgetc(in_file);
   	rate /=1000;                                         /*recover float to 4dp*/
   	mom  /=1000;

   	p_this = init_ffnet(in, hid, out, rate, mom);       /*initialise new network*/
   	if (p_this == NULL)
   	{
   		fferr(6);
      	return(FALSE);
   	}
   	                                                 /*read in data from file*/
   	fread(p_this->p_in, sizeof(float), p_this->in, in_file);
    fread(p_this->p_hid, sizeof(float), p_this->hid, in_file);
    fread(p_this->p_out, sizeof(float), p_this->out, in_file);
    fread(p_this->p_hid_w, sizeof(float), (p_this->in+1)*p_this->hid, in_file);
   	fread(p_this->p_out_w, sizeof(float), (p_this->hid+1)*p_this->out, in_file);
   	fread(p_this->p_hid_mom, sizeof(float), (p_this->in+1)*p_this->hid, in_file);
   	fread(p_this->p_out_mom, sizeof(float), (p_this->hid+1)*p_this->out, in_file);

   	fclose(in_file);
   	return(p_this);
}

/*Takes a pointer to a net and returns memory to calloc.*/
int del_ffnet(FFNET* p_this)
{
   	if(p_this == NULL)
   	{                                         /*net to delete?*/
   		fferr(6);
    	return(FALSE);
   	}
	free(p_this->p_in);                               /*return all allocations*/
	free(p_this->p_hid);
   	free(p_this->p_out);
   	free(p_this->p_targ);
   	free(p_this->p_hid_w);
   	free(p_this->p_out_w);
   	free(p_this->p_hid_mom);
   	free(p_this->p_out_mom);
   	return(TRUE);
}

/*Take a pointer to a genotype and initialis net as phenotype.*/
FFNET* morph_ffnet(int* gene, int length, int depth, int in, int hid, int out)
{
	int i, rate = 0, mom = 0;
   	FFNET* p_this;
   	int units = ((in+BIAS)*hid)+((hid+BIAS)*out);             /*enough genotype?*/

   	if (length != units*depth)
   	{
		fferr(5);
      	return(FALSE);
   	}

   	p_this = init_ffnet(in, hid, out, rate, mom);          /*initialise and test*/
   	if(p_this == NULL)
   	{
   		fferr(6);
      	return(FALSE);
   	}                                                /*translate weight arrays*/

   	for (i=0; i<(p_this->in+BIAS)*p_this->hid; i++)
   	{
   		*(p_this->p_hid_w + i) = ff_bin_dec(gene, depth);
      	gene += depth;
   	}

   	for (i = 0; i < (p_this->hid + BIAS)*p_this->out; i++)
   	{
   		*(p_this->p_out_w + i) = ff_bin_dec(gene, depth);
      	gene += depth;
   	}

   	return(p_this);
}

/*Take pointers to input and target vectors and train the net*/
float train_ffnet(FFNET* p_this, float* in_vect, float* targ_vect)
{
	int i;
   	if(p_this == NULL)
   	{                                         /*error check*/
   		fferr(7);
    	return(FALSE);
   	}

   	for (i = 0; i < p_this->in; i++)
   	{                          /*read data to p_in*/
      	*(p_this->p_in + i) = *(in_vect + i);
   	}

   	for (i = 0; i < p_this->in; i++)
   	{                                 /*and p_targ*/
      	*(p_this->p_targ + i) = *(targ_vect + i);
   	}

   	forward_pass(p_this);                                            /*run net*/
   	backward_pass(p_this);

   	return(sum_squared_error(p_this));
}

/*Take a pointer to an input vector and return a pointer to
  an output vector*/
float* test_ffnet(FFNET* p_this, float* in_vect)
{
	int i;
	if (p_this == NULL)
	{
   		fferr(7);                                                 /*error check*/
      	return(FALSE);
   	}

	for (i = 0; i < p_this->in; i++)
	{
      *(p_this->p_in + i) = *(in_vect + i);                             /*read data*/
   	}

   	forward_pass(p_this);                                                /*run*/

   	return(p_this->p_out);
}

/*Print node and weight data to screen.*/
int print_ffnet(FFNET* p_this)
{
	if (p_this == NULL)
	{                                            /*error check*/
   		fferr(7);
      	return(FALSE);
    }

   	printf("Error:%6.4f\nOutput:", sum_squared_error(p_this));
   	print_values(p_this->p_out, p_this->out);
   	printf("   Target:");
   	print_values(p_this->p_targ, p_this->out);
   	printf("\n");
   	printf("O Weights:");
   	print_values(p_this->p_out_w, (p_this->hid + BIAS)*p_this->out);
   	printf("\n");
   	printf("Hidden:");
   	print_values(p_this->p_hid, p_this->hid);
   	printf("\n");
   	printf("H Weights:");
   	print_values(p_this->p_hid_w, (p_this->in + BIAS)*p_this->hid);
   	printf("\n");
   	printf("Input:");
   	print_values(p_this->p_in, p_this->in);
   	printf("\n=====================================\n");

   	return(TRUE);
}

/*==========================================================================*/
/*Neural engine functions                                                   */

/*Take pointers to two layers and thier weights and sizes*/
/*and calculate the forward activation of the net*/
static void forwards(float* p_in, float* p_out, float* p_weight,
							int inputs, int outputs){
   int i, j;
   float sum;
   for (j=0; j<outputs; j++){                            /*for each output, j*/
      sum=0;
      for (i=0; i<inputs; i++){                           /*for each input, i*/
      	sum+=(*(p_weight+(j*inputs+i+j))*(*(p_in+i)));    /*sum input.weight*/
      }
      sum+=(*(p_weight+(j*inputs+inputs+j)));           /*and the bias weight*/
      *(p_out+j)=1/(1+exp(-sum));                         /*logistic function*/
   }
}

/*Run the net forwards*/
static void forward_pass(FFNET* p_this){
   forwards(p_this->p_in,
   			p_this->p_hid,
            p_this->p_hid_w,
            p_this->in,
            p_this->hid);
   forwards(p_this->p_hid,
   			p_this->p_out,
            p_this->p_out_w,
            p_this->hid,
            p_this->out);
}

/*Calculate and return the sum squared error of the net*/
static float sum_squared_error(FFNET* p_this){
   int j;
   p_this->error=0;
   for (j=0; j<p_this->out; j++){
       p_this->error+=pow((*(p_this->p_out+j)-*(p_this->p_targ+j)),2);
   }
   return(p_this->error);
}

/*Calculate the effect of hidden unit i on sum squared error*/
static float sum_dEdxi(FFNET* p_this, int i){
   int j;
   float sum=0;
   for (j=0; j<p_this->out; j++){                           /*for each output*/
       sum+=2*(*(p_this->p_out+j)-*p_this->p_targ+j)*(*p_this->p_out+j)*
       	(1-(*p_this->p_out+j))*(*(p_this->p_out_w+j*p_this->hid+i+j));
   }
   return(sum);
}

/*Calculate the scaled change in output weight j.i.*/
static float delta_w(FFNET* p_this, int i, int j){
   float delta;
   delta=(2*(*(p_this->p_out+j)-*(p_this->p_targ+j))*
   	(*(p_this->p_out+j))*(1-*(p_this->p_out+j))*(*(p_this->p_hid+i)));
   return(delta);
}

/*Calculate the scaled change in hidden weight i.h.*/
static float delta_v(FFNET* p_this, float deltaX, int h, int i){
   float delta;
   delta=(*(p_this->p_hid+i)*(1-*(p_this->p_hid+i))*
    (*(p_this->p_in+h))*deltaX);
   return(delta);
}

/*Propagate the error backwards*/
static void backward_pass(FFNET* p_this){
   int h, i, j, offset;
   float deltaX;
   float* p_temp_w=(float*)calloc((p_this->hid+BIAS)*p_this->out,
   	sizeof(float));
   float* p_temp_v=(float*)calloc((p_this->in+BIAS)*p_this->hid,
   	sizeof(float));
   if(p_temp_w==NULL || p_temp_v==NULL){
   	fferr(1);
   }
   for (j = 0; j < p_this->out; j++){                      /*for each output node*/
     offset = (j*p_this->hid) + j;
       for (i = 0; i < p_this->hid; i++)
       {                  /*for each hidden node*/
      	   *(p_temp_w + offset + i) = *(p_this->p_out_w + offset + i) - (p_this->rate*delta_w(p_this, i, j)) + (*(p_this->p_out_mom+offset+i)*p_this->mom);          /*adjust*/
       }
       offset = (j*p_this->hid + p_this->hid + j);
       *(p_temp_w + offset) = *(p_this->p_out_w + offset) - (p_this->rate*(2*(*(p_this->p_out + j)-*(p_this->p_targ + j))*(*(p_this->p_out + j))*(1 - *(p_this->p_out + j)))) + (*(p_this->p_out_mom + offset)*p_this->mom);            /*for the bias*/
   }
   for (i=0; i<p_this->hid; i++){                      /*for each hidden node*/
       deltaX=sum_dEdxi(p_this, i);                    /*get its contribution*/
       offset=(i*p_this->in+i);
       for (h=0; h<p_this->in; h++){                         /*for each input*/
           *(p_temp_v+offset+h)=*(p_this->p_hid_w+offset+h)-
           	(p_this->rate*delta_v(p_this, deltaX, h, i))
           	+(*(p_this->p_hid_mom+offset+h)*p_this->mom);            /*adjust*/
       }
       offset=(i*p_this->in+p_this->in+i);
       *(p_temp_v+offset)=*(p_this->p_hid_w+offset)-
        (p_this->rate*(*(p_this->p_hid+i)*(1-*(p_this->p_hid+i))*deltaX))+
        (*(p_this->p_hid_mom+offset)*p_this->mom);              /*adjust bias*/
   }
   for (j=0; j<(p_this->hid+BIAS)*p_this->out; j++){      /*weights->momentum*/
       *(p_this->p_out_mom+j)=*(p_this->p_out_w+j);
   }
   for (j=0; j<(p_this->in+BIAS)*p_this->hid; j++){
     *(p_this->p_hid_mom+j)=*(p_this->p_hid_w+j);
   }
   for (j=0; j<(p_this->hid+BIAS)*p_this->out; j++){     /*insert new weights*/
       *(p_this->p_out_w+j)=*(p_temp_w+j);
   }
   for (j=0; j<(p_this->in+BIAS)*p_this->hid; j++){
       *(p_this->p_hid_w+j)=*(p_temp_v+j);
   }
   free(p_temp_w);                                  /*return memory to calloc*/
   free(p_temp_v);
}

/*===========================================================================*/
/*Internal utility functions.                                                */

/*Take two ints and return pointer to an initialised array of floats*/
static float* init(int length, int dim){
   int i, temp;
   float ftemp;
   float* array=(float*)calloc(length, sizeof(float));
   if(array==NULL) return(FALSE);                               /*error check*/
   if (dim){
       for (i=0; i<length; i++){
   	   temp=rand()%10000;
      	   ftemp=1-(float)temp/5000;                         /*range -1<n<1 */
   	   array[i]=ftemp;
       }
   }
   return (array);
}

/*Take a pointer to a binary string and return a float.*/
float ff_bin_dec(int* gene, int depth){
	int i;
   float phen=0;
   for (i=0; i<depth; i++){                          /*translate 'depth' bits*/
   	if(*(gene+i)){
	   	phen+=pow(2,depth-1-i);
      }
   }
   phen=1-(phen/pow(2,depth-2));                        /*condition to -1<n<1*/
   return(phen);
}

/*Take a pointer to an array and its size, print the output values*/
static void print_values(float* p_array, int size){
   int i;
   for (i=0; i<size; i++){
       printf(" %6.4f .", *(p_array+i));
   }
}

/*Take an error signal and print message to stderr*/
static void fferr(signal){
	if(signal==1)fprintf(stderr, "Unable to allocate enough memory.");
   if(signal==2)fprintf(stderr, "Weights have exceeded platform limits.");
   if(signal==3)fprintf(stderr, "Unable to open output stream.");
   if(signal==4)fprintf(stderr, "Unable to open input stream.");
   if(signal==5)fprintf(stderr, "Unit numbers and bit length do not agree.");
   if(signal==6)fprintf(stderr, "Unabble to initialise net parameters.");
   if(signal==7)fprintf(stderr, "Net parameers must be initialised.");
}

/*===========================================================================*/
/*Last modified 31.12.2001*/