/****************************************************************************/
/*
/* File:		"rtrl_pc.c"
/*				Real-Time Recurrent Learning for an RMLP (PC version)
/*
/* Authors:		Scott A. Morrison, scott@cnel.ufl.edu
/*				Deniz Erdogmus, deniz@cnel.ufl.edu
/*				Justin Sanchez, justin@cnel.ufl.edu
/*				Dongho Han, dhan@cnel.ufl.edu
/*
/*				Computational NeuroEngineering Laboratory
/*				Department of Electrical and Computer Engineering
/*				The University of Florida
/*				Gainesville, FL, 32611, USA
/*				
/* Created:		26 November 2002
/*
/*****************************************************************************/

#define FILE_NAME_MAX		40
#define MAX_LABEL_LENGTH	40

/* the following will generate a random number between -0.8 and +0.8 */
#define	RANDOM_NUMBER	(rand()%2 ? 1 : -1)/(0.1f*(rand()%500) + 1)

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

/*********************** Parameters from "params.txt" ***********************/
char	PARAM_FILE[FILE_NAME_MAX] = "";
char	INPUT_FILE[FILE_NAME_MAX] = "";
char	DESIRED_FILE[FILE_NAME_MAX] = "";
int		TRAINING_START = 0;
int		TRAINING_END = 0;
int		TESTING_START = 0;
int		TESTING_END = 0;
int		TRAJ_LENGTH = 0;
int		TRAJ_PER_UPDATE = 0;
int		EPOCHS = 0;
int		N0 = 0;
int		N1 = 0;
int		N2 = 0;
float	SLOPE = 0.0f;
float	DERIVATIVE_OFFSET = 0.0f;
float	ETA1 = 0.0f;
float	ETA2 = 0.0f;
float	ETAF = 0.0f;
int		ADAPTIVE_STEP_SIZE = 0;
float	ETA_MIN = 0.0f;
float	ETA_MAX = 0.0f;
float	ETA_INCREASE = 0.0f;
float	MOMENTUM = 0.0f;
int		USE_INITS = 0;
char	W1_INIT_FILE[FILE_NAME_MAX*40] = "";
char	B1_INIT_FILE[FILE_NAME_MAX*40] = "";
char	W2_INIT_FILE[FILE_NAME_MAX] = "";
char	B2_INIT_FILE[FILE_NAME_MAX] = "";
char	WF_INIT_FILE[FILE_NAME_MAX] = "";
char	WEIGHT_OUT_FILE[FILE_NAME_MAX] = "";

int		TRAINING_LENGTH = 0;
int		UPDATES_PER_EPOCH = 0;

/*************************** Temporary Variables ***************************/
int		epoch	= 0;

float	d_power = 0.0f;	/* power of the desired signal, for normalized MSE */
float	e_power = 0.0f;	/* power of the error signal, for normalized MSE */
	
float	*W1,*b1,*W2,*b2,*Wf;			/* pointers to network weights */
float	*G_W1,*G_b1,*G_W2,*G_b2,*G_Wf;	/* pointers to gradients */
float	*U_W1,*U_b1,*U_W2,*U_b2,*U_Wf;	/* pointers to weight updates */
float	*c,*g,*h;						/* pointers to update parameters */
float	*y1p,*z1,*Y1,*y2,*e,*a,*temp;	/* pointers to network feedforward */
float	*input,*desired;				/* pointers to data */

float	*x_ptr = NULL;
float	*d_ptr = NULL;

float	MSE = 0.0f;
float	MSEp = 0.0f;

/*************************** Function Prototypes ***************************/
void displayMSE();
void updateGradient();
void doTraining();
void getParameters();
void showParameters();
void getWeights();
void putWeights();
float randomNumber(void);


/***************************** Main program **********************************/
int main(int argc, char **argv)
{

    int		i = 0;
	int		j = 0;
	int		lines = 0;

	FILE	* infile;
	FILE	* desfile;

	/* check for the proper command line arguments */
	if (argc != 2)
	{
		printf("\n*ERROR* Please specify only the parameter file name.\n");
		printf("Usage: bptt params.txt\n");
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);		
	}
	else
	{
		strcpy(PARAM_FILE,argv[1]);
	}

	/* get the parameters from the parameter file */
	getParameters();

	TRAINING_LENGTH = TRAINING_END - TRAINING_START;
	UPDATES_PER_EPOCH = floor(TRAINING_LENGTH / (TRAJ_LENGTH * TRAJ_PER_UPDATE));
	
	/* display title and parameters to screen */
	printf("\n**********************************************************************\n");
	printf("Training a Recurrent Neural Network using Backpropagation Through Time\n");
	showParameters(); /* display the parameters */
	printf("**********************************************************************\n");

	/* allocate the memory */
	W1		= calloc(N1*N0,sizeof(float));	/* first layer weights */
	b1		= calloc(N1*1,sizeof(float));	/* first layer bias */
	W2		= calloc(N2*N1,sizeof(float));	/* second layer weights */
	b2		= calloc(N2*1,sizeof(float));	/* second layer bias */
	Wf		= calloc(N1*N1,sizeof(float));	/* feedback weights */

	G_W1	= calloc(N1*N0,sizeof(float));	/* gradient updates for above */
	G_b1	= calloc(N1*1,sizeof(float));
	G_W2	= calloc(N2*N1,sizeof(float));
	G_b2	= calloc(N2*1,sizeof(float));
	G_Wf	= calloc(N1*N1,sizeof(float));

	U_W1	= calloc(N1*N0,sizeof(float));	/* weight updates for above */
	U_b1	= calloc(N1*1,sizeof(float));
	U_W2	= calloc(N2*N1,sizeof(float));
	U_b2	= calloc(N2*1,sizeof(float));
	U_Wf	= calloc(N1*N1,sizeof(float));

	c		= calloc(N1*N1,sizeof(float));	/* intermediate update parameters */
	g		= calloc(N1*N1*N0,sizeof(float));
	h		= calloc(N1*N1*N1,sizeof(float));

	y1p		= calloc(N1*1,sizeof(float));	/* previous Y1 */
	z1		= calloc(N1*1,sizeof(float));	/* input to hidden layer */
	Y1		= calloc(N1*1,sizeof(float));	/* output of hidden layer */
	y2		= calloc(N2*1,sizeof(float));	/* network output */
	e		= calloc(N2*1,sizeof(float));	/* network error */
	a		= calloc(N1*1,sizeof(float));	/* save sech(z1)^2 */
	temp	= calloc(N1*1,sizeof(float));

	input	= calloc(TRAINING_LENGTH*N0,sizeof(float));	/* input data */
	desired = calloc(TRAINING_LENGTH*N2,sizeof(float)); /* desired data */

	x_ptr = input;		/* initialize the pointer to the input data */
	d_ptr = desired;	/* initialize the pointer to the desired data */

	ETA1 = ETA1 / TRAJ_LENGTH;	/* adjust the step size according to the trajectory length */
	ETA2 = ETA2 / TRAJ_LENGTH;
	ETAF = ETAF / TRAJ_LENGTH;

	printf("Opening the input file...\n");

	/* open the input file */
	infile = fopen(INPUT_FILE,"r");
	if (infile == NULL)
	{
		printf("\n*ERROR* The input file \"%s\" could not be opened for reading.\n",INPUT_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	/* read from the input file */
	lines = 0;
	while( fscanf( infile,"%f",&input[lines]) != EOF)
	{	
		lines++;		
		if (lines >= TRAINING_LENGTH*N0) break;
	}

	printf("There were %d values read from the input file \"%s\".\n",lines,INPUT_FILE);

	/* close the input file */
    if (fclose(infile) != 0)
	{
		printf("\n*ERROR* The input file \"%s\" could not be closed properly.\n",INPUT_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	printf("Opening the desired file...\n");

	/* open the desired file */
	desfile = fopen(DESIRED_FILE,"r");
	if (infile == NULL)
	{
		printf("\n*ERROR* The desired file \"%s\" could not be opened for reading.\n",DESIRED_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	/* read from the desired file */
	lines = 0;
	while(fscanf( desfile,"%f",&desired[lines]) != EOF)
	{
		lines++;		
		if (lines >= TRAINING_LENGTH*N2) break;
	}

	printf("There were %d values read from the desired file \"%s\".\n",lines,DESIRED_FILE);

	/* close the desired file */
    if (fclose(desfile) != 0)
	{
		printf("\n*ERROR* The desired file \"%s\" could not be closed properly.\n",DESIRED_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	/* initialize weights */
	if (USE_INITS == 0)
	{
		printf("\nInitializing weights to small random numbers...\n");
		/* initialize to random numbers */	
		for (i=0;i<N0*N1;i++)
		{
			W1[i] = randomNumber()/2; 
		}

		for (i=0;i<N1;i++)
		{
			b1[i] = randomNumber()/2; 
		}

		for (i=0;i<N1*N1;i++)
		{
			Wf[i] = randomNumber()/2; 
		}

		for (i=0;i<N1*N2;i++)
		{
			W2[i] = randomNumber()/2; 
		}

		for (i=0;i<N2;i++)
		{
			b2[i] = randomNumber()/2; 
		}
	}
	else
	{
		/* initialize from files */
		getWeights();
	}

	/* train the network */
	printf("\nTraining started...\n\n");

	doTraining();
	
	printf("\nTraining finished.\n");
	printf("Saving the weights to files...\n");

	/* save weights */
	putWeights();

	/* free the memory */
	free(W1);
	free(b1);
	free(W2);
	free(b2);
	free(Wf);

	free(G_W1);
	free(G_b1);
	free(G_W2);
	free(G_b2);
	free(G_Wf);

	free(U_W1);
	free(U_b1);
	free(U_W2);
	free(U_b2);
	free(U_Wf);

	free(c);
	free(g);
	free(h);

	free(y1p);
	free(z1);
	free(Y1);
	free(y2);
	free(e);
	free(a);
	free(temp);

	free(input);
	free(desired);

	printf("\nExiting.\n");	/* we're done, so quit */
	exit(EXIT_SUCCESS);

} /* end of main() */


/*****************************************************************************/
void doTraining()
{

	int update	= 0;
	int traj	= 0;
	int i		= 0;
	int T		= 0;

	for(epoch=0;epoch<EPOCHS;epoch++)
	{
		/* reset the hidden state before staring epoch */
		for(i=0;i<N1;i++)
		{
			y1p[i] = 0.0f;
		}

		/* reset x_ptr */
		x_ptr = input;			/* double check this!!!! */
		d_ptr = desired;

		for(update=0;update<UPDATES_PER_EPOCH;update++)
		{

			/* reset gradients and momentum terms */
			for (i=0;i<N0*N1;i++)
			{
				U_W1[i] = 0.0f;
				G_W1[i] = 0.0f; 
			}

			for (i=0;i<N1;i++)
			{
				U_b1[i] = 0.0f;
				G_b1[i] = 0.0f; 
			}

			for (i=0;i<N1*N1;i++)
			{
				U_Wf[i] = 0.0f;
				G_Wf[i] = 0.0f; 
			}

			for (i=0;i<N1*N2;i++)
			{
				U_W2[i] = 0.0f;
				G_W2[i] = 0.0f; 
			}

			for (i=0;i<N2;i++)
			{
				U_b2[i] = 0.0f;
				G_b2[i] = 0.0f; 
			}


			/* for(B=0;B<TRAJ_LENGTH;B++) ORIGINALLY, I think it is a mistake */
			for(traj=0;traj<TRAJ_PER_UPDATE;traj++)
			{
				/* reset c, g, and h */
				for(i=0;i<N1*N1;i++)
				{
					c[i] = 0.0f;
				}
				for(i=0;i<N1*N1*N0;i++)
				{
					g[i] = 0.0f;
				}				
				for(i=0;i<N1*N1*N1;i++)
				{
					h[i] = 0.0f;
				}
				for(T=0;T<TRAJ_LENGTH;T++)
				{	
					updateGradient();

					/* save the state */
					for(i=0;i<N1;i++)
					{
						y1p[i] = Y1[i];
					}

					/* move to next input and desired */
					x_ptr = x_ptr + N0;
					d_ptr = d_ptr + N2;
				}

			} /* end of loop over all trajectories within this block */


			/* update weights */
			for (i=0;i<N0*N1;i++)
			{
				U_W1[i] = -1.0f * ETA1 * G_W1[i] + MOMENTUM * U_W1[i];
				W1[i] = W1[i] + U_W1[i]; 
			}

			for (i=0;i<N1;i++)
			{
				U_b1[i] = -1.0f * ETA1 * G_b1[i] + MOMENTUM * U_b1[i];
				b1[i] = b1[i] + U_b1[i]; 
			}

			for (i=0;i<N1*N1;i++)
			{
				U_Wf[i] = -1.0f * ETAF * G_Wf[i] + MOMENTUM * U_Wf[i];
				Wf[i] = Wf[i] + U_Wf[i]; 
			}

			for (i=0;i<N1*N2;i++)
			{
				U_W2[i] = -1.0f * ETA2 * G_W2[i] + MOMENTUM * U_W2[i];
				W2[i] = W2[i] + U_W2[i]; 
			}

			for (i=0;i<N2;i++)
			{
				U_b2[i] = -1.0f * ETA2 * G_b2[i] + MOMENTUM * U_b2[i];
				b2[i] = b2[i] + U_b2[i]; 
			}


		} /* end of loop over all update chunks */

		MSEp = MSE;

		/* display to screen the MSE after this epoch */
		displayMSE();

		if (ADAPTIVE_STEP_SIZE)
		{
			if (MSE < MSEp)
			{
				/* increase the step size */
				ETA1 = ETA1 * ETA_INCREASE;
				ETA2 = ETA2 * ETA_INCREASE;
				ETAF = ETAF * ETA_INCREASE;

				if (ETA1 > ETA_MAX)
				{
					ETA1 = ETA_MAX;
					printf("\n*WARNING* ETA1 reached ETA_MAX %f.\n",ETA_MAX);
				}
				if (ETA2 > ETA_MAX)
				{
					ETA2 = ETA_MAX;
					printf("\n*WARNING* ETA2 reached ETA_MAX %f.\n",ETA_MAX);
				}
				if (ETAF > ETA_MAX)
				{
					ETAF = ETA_MAX;
					printf("\n*WARNING* ETAF reached ETA_MAX %f.\n",ETA_MAX);
				}
			}
			else
			{
				/* decrease the step size */
				ETA1 = ETA1 / (ETA_INCREASE * 4);
				ETA2 = ETA2 / (ETA_INCREASE * 4);
				ETAF = ETAF / (ETA_INCREASE * 4);

				/* make sure the step size isn't too small */
				if (ETA1 < ETA_MIN)
				{
					ETA1 = ETA_MIN;
					printf("\n*WARNING* ETA1 reached ETA_MIN %f.\n",ETA_MIN);
				}
				if (ETA2 < ETA_MIN)
				{
					ETA2 = ETA_MIN;
					printf("\n*WARNING* ETA2 reached ETA_MIN %f.\n",ETA_MIN);
				}
				if (ETAF < ETA_MIN)
				{
					ETAF = ETA_MIN;
					printf("\n*WARNING* ETAF reached ETA_MIN %f.\n",ETA_MIN);
				}

			}
		}	/* end of ADAPTIVE_STEP_SIZE */


	} /* end of loop over epochs */

} /* end of doTraining() */



void displayMSE()
{
	int i = 0;
	int t = 0;
	int k = 0;
	int j = 0;
	int m = 0;
	float	*X_ptr = NULL;
	float	*D_ptr = NULL;

	float MSD = 0.0f;
	MSE = 0.0f;

	/* reset the state variables */
	for(i=0;i<N1;i++)
	{
		y1p[i] = 0.0f;
	}

	/* reset the pointers to the input and desired data */
	X_ptr = input;
	D_ptr = desired;

	/* go through the entire input/desired data */
	for(t=0;t<TRAINING_LENGTH;t++)
	{

		/* forward pass through the network */
		for(j=0;j<N1;j++)
		{
			z1[j]= 0;

			for (i=0;i<N0;i++)
			{
				z1[j] = z1[j] + X_ptr[t*N0+i]*W1[j*N0+i];
			}

			for (m=0;m<N1;m++)
			{
				z1[j] = z1[j] + y1p[m] * Wf[j*N1+m];
			}

			z1[j] = z1[j] + b1[j];	/* input to hidden layer */
			Y1[j] = tanh(SLOPE*z1[j]);	/* output from hidden layer */
			y1p[j] = Y1[j]; /* save the state for next time */
		}

		d_power = 0.0f;
		e_power = 0.0f;
		for(k=0;k<N2;k++)
		{
			y2[k]= 0;

			for (j=0;j<N1;j++)
			{
				y2[k] = y2[k] + Y1[j]*W2[k*N1+j];
			}

			y2[k] = y2[k] + b2[k];	/* linear output of the network */
			e_power = e_power + pow(D_ptr[t*N2+k] - y2[k],2);	/* error between desired and output */
			d_power = d_power + pow(D_ptr[t*N2+k],2);
		}
		MSE = MSE + (e_power);
		MSD = MSD + d_power;
	}

	MSE = MSE / (N2*TRAINING_LENGTH);
	MSD = MSD / (N2*TRAINING_LENGTH);

	printf("Epoch %d   \t MSE = %2.9f \t MSEN = %2.9f\n",epoch,MSE,MSE/MSD);

}


/*****************************************************************************/
void updateGradient()
{
	/* temporary variables used in this function */
	int		i = 0;
	int		j = 0;
	int		m = 0;
	int		k = 0;
	int		gamma = 0;
	int		alpha = 0;
	int		beta = 0;

	/* forward pass through the network */
	for(j=0;j<N1;j++)
	{
		z1[j]= 0;

		for (i=0;i<N0;i++)
		{
			z1[j] = z1[j] + x_ptr[i]*W1[j*N0+i];
		}

		for (m=0;m<N1;m++)
		{
			z1[j] = z1[j] + y1p[m] * Wf[j*N1+m];
		}

		z1[j] = z1[j] + b1[j];	/* input to hidden layer */
		Y1[j] = tanh(SLOPE*z1[j]);	/* output from hidden layer */
		a[j] = SLOPE*(1-pow(Y1[j],2)+DERIVATIVE_OFFSET);	/* function derivative at operating point */
	}

	for(k=0;k<N2;k++)
	{
		y2[k]= 0;

		for (j=0;j<N1;j++)
		{
			y2[k] = y2[k] + Y1[j]*W2[k*N1+j];
		}

		y2[k] = y2[k] + b2[k];	/* linear output of the network */
		e[k] = d_ptr[k] - y2[k];	/* error between desired and output */
	}


	/* The forward pass is done. Now backpropagate the updates. */

	/* update c */
	for(beta=0;beta<N1;beta++)
	{
		for(j=0;j<N1;j++)
		{
			temp[j] = 0.0f;
			for(i=0;i<N1;i++)
			{
				temp[j] = temp[j] + Wf[j*N1+i] * c[beta*N1+i];
			}

			if (j == beta)
			{
				c[beta*N1+j] = temp[j] + 1;
			}
			else
			{
				c[beta*N1+j] = temp[j];
			}
			
			c[beta*N1+j] = c[beta*N1+j] * a[j];
			
		}
	}

	/* update g */
	for(gamma=0;gamma<N0;gamma++)
	{
		for(beta=0;beta<N1;beta++)
		{
			for(alpha=0;alpha<N1;alpha++)
			{
				temp[alpha] = 0.0f;
				for (i=0;i<N1;i++)
				{
					//temp[alpha] = temp[alpha] + Wf[alpha*N1+i] * g[gamma*N0+beta*N1+i];			
					temp[alpha] = temp[alpha] + Wf[alpha*N1+i] * g[gamma*N1*N1+beta*N1+i];				

				}
				if (alpha == beta)
				{					
					//g[gamma*N0+beta*N1+alpha] = temp[alpha] + x_ptr[gamma];
					g[gamma*N1*N1+beta*N1+alpha] = temp[alpha] + x_ptr[gamma];
				}
				else
				{
					//g[gamma*N0+beta*N1+alpha] = temp[alpha];
					g[gamma*N1*N1+beta*N1+alpha] = temp[alpha];
				}
			
				//g[gamma*N0+beta*N1+alpha] = g[gamma*N0+beta*N1+alpha] * a[alpha];
				g[gamma*N1*N1+beta*N1+alpha] = g[gamma*N1*N1+beta*N1+alpha] * a[alpha];

			}
		}
	}

	/* update h */
	for(gamma=0;gamma<N1;gamma++)
	{
		for(beta=0;beta<N1;beta++)
		{
			for(alpha=0;alpha<N1;alpha++)
			{
				temp[alpha] = 0.0f;
				for (i=0;i<N1;i++)
				{
					temp[alpha] = temp[alpha] + Wf[alpha*N1+i] * h[gamma*N1+beta*N1+i];					
				}
				if (alpha == beta)
				{
					h[gamma*N1+beta*N1+alpha] = temp[alpha] + y1p[gamma];
				}
				else
				{
					h[gamma*N1+beta*N1+alpha] = temp[alpha];
				}
				
				h[gamma*N1+beta*N1+alpha] = h[gamma*N1+beta*N1+alpha] * a[alpha];

			}
		}
	}

	/* update the gradients */
	for (k=0;k<N2;k++)
	{
		G_b2[k] = G_b2[k] - e[k]; 

		for(j=0;j<N1;j++)
		{
			G_W2[k*N1+j] = G_W2[k*N1+j] - e[k]*Y1[j];
		}	

	}

	for (j=0;j<N1;j++)
	{
		temp[j] = 0.0f;
		for (i=0;i<N2;i++)
		{
			temp[j] = temp[j] + e[i]*W2[i*N1+j];
		}

		G_b1[j] = G_b1[j] - c[j*N1+j] * temp[j];

		for (i=0;i<N0;i++)
		{
			//G_W1[j*N0+i] = G_W1[j*N0+i] - g[i*N0+j*N1+j] * temp[j];
			G_W1[j*N0+i] = G_W1[j*N0+i] - g[i*N1*N1+j*N1+j] * temp[j];
		}
		
		for (i=0;i<N1;i++)
		{
			G_Wf[j*N1+i] = G_Wf[j*N1+i] - h[i*N1+j*N1+j] * temp[j];
		}

	}

} /* end of updateGradient() */


/*****************************************************************************/
void getParameters()
{

	char	label[MAX_LABEL_LENGTH];
	char	value[MAX_LABEL_LENGTH];
	FILE	*paramfile;

	/* open the parameter file */
	paramfile = fopen(PARAM_FILE,"r");
	if (paramfile == NULL)
	{
		printf("\n*ERROR* The parameter file \"%s\" could not be opened for reading.\n",PARAM_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	/* read from the parameter file */
	while( fscanf( paramfile,"%s %s",&label,&value) != EOF )
	{
		if (strcmp(label, "Input_Filename") == 0)
        {
			strcpy(INPUT_FILE,value);
		}
		else if (strcmp(label, "Desired_Filename") == 0)
		{
			strcpy(DESIRED_FILE,value);
		}
		else if (strcmp(label, "Training_Start") == 0)
		{
			TRAINING_START = atoi(value);
		}
		else if (strcmp(label, "Training_End") == 0)
		{
			TRAINING_END = atoi(value);
		}
		else if (strcmp(label, "Testing_Start") == 0)
		{
			TESTING_START = atoi(value);
		}
		else if (strcmp(label, "Testing_End") == 0)
		{
			TESTING_END = atoi(value);
		}
		else if (strcmp(label, "Trajectory_Length") == 0)
		{
			TRAJ_LENGTH = atoi(value);
		}
		else if (strcmp(label, "Trajectories_per_update") == 0)
		{
			TRAJ_PER_UPDATE = atoi(value);
		}
		else if (strcmp(label, "Epochs") == 0)
		{
			EPOCHS = atoi(value);
		}
		else if (strcmp(label, "N0") == 0)
		{
			N0 = atoi(value);
		}
		else if (strcmp(label, "N1") == 0)
		{
			N1 = atoi(value);
		}
		else if (strcmp(label, "N2") == 0)
		{
			N2 = atoi(value);
		}
		else if (strcmp(label, "Slope_of_Nonlinearity") == 0)
		{
			SLOPE = atof(value);
		}
		else if (strcmp(label, "Derivative_Offset") == 0)
		{
			DERIVATIVE_OFFSET = atof(value);
		}
		else if (strcmp(label, "Eta1") == 0)
		{
			ETA1 = atof(value);
		}
		else if (strcmp(label, "Eta2") == 0)
		{
			ETA2 = atof(value);
		}
		else if (strcmp(label, "Etaf") == 0)
		{
			ETAF = atof(value);
		}
		else if (strcmp(label, "ADAPTIVE_STEP_SIZE") == 0)
		{
			ADAPTIVE_STEP_SIZE = atoi(value);
		}
		else if (strcmp(label, "ETA_INCREASE") == 0)
		{
			ETA_INCREASE = atof(value);
		}
		else if (strcmp(label, "ETA_MAX") == 0)
		{
			ETA_MAX = atof(value);
		}
		else if (strcmp(label, "ETA_MIN") == 0)
		{
			ETA_MIN = atof(value);
		}
		else if (strcmp(label, "Momentum") == 0)
		{
			MOMENTUM = atof(value);
		}
		else if (strcmp(label, "Use_Inits_Below") == 0)
		{
			USE_INITS = atoi(value);
		}
		else if (strcmp(label, "W1init") == 0)
		{
			strcpy(W1_INIT_FILE,value);
		}
		else if (strcmp(label, "b1init") == 0)
		{
			strcpy(B1_INIT_FILE,value);
		}
		else if (strcmp(label, "W2init") == 0)
		{
			strcpy(W2_INIT_FILE,value);
		}
		else if (strcmp(label, "b2init") == 0)
		{
			strcpy(B2_INIT_FILE,value);
		}
		else if (strcmp(label, "Wfinit") == 0)
		{
			strcpy(WF_INIT_FILE,value);
		}
		else if (strcmp(label, "Weight_Output_Filename") == 0)
		{
			strcpy(WEIGHT_OUT_FILE,value);
		}		
		else if (strcmp(label, "end") == 0)
		{
			/* end of the parameter file */
			break;
		}
	}

	/* close the parameter file */
    if (fclose(paramfile) != 0)
	{
		printf("\n*ERROR* The parameter file \"%s\" could not be closed properly.\n",PARAM_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}
	
} /* end of getParameters() */


/*****************************************************************************/
void showParameters()
{
	printf("Parameter File: \"%s\"\n",PARAM_FILE);

	printf("   INPUT_FILE        = %s\n",INPUT_FILE);
	printf("   DESIRED_FILE      = %s\n",DESIRED_FILE);
	printf("   TRAINING_START    = %d\n",TRAINING_START);
	printf("   TRAINING_END      = %d\n",TRAINING_END);
	printf("   TESTING_START     = %d\n",TESTING_START);
	printf("   TESTING_END       = %d\n",TESTING_END);;
	printf("   TRAJ_LENGTH       = %d\n",TRAJ_LENGTH);
	printf("   TRAJ_PER_UPDATE   = %d\n",TRAJ_PER_UPDATE);
	printf("   EPOCHS            = %d\n",EPOCHS);
	printf("   N0                = %d\n",N0);
	printf("   N1                = %d\n",N1);
	printf("   N2                = %d\n",N2);
	printf("   SLOPE             = %f\n",SLOPE);
	printf("   DERIVATIVE_OFFSET = %f\n",DERIVATIVE_OFFSET);
	printf("   ETA1              = %2.8f\n",ETA1);
	printf("   ETA2              = %2.8f\n",ETA2);
	printf("   ETAF              = %2.8f\n",ETAF);
	printf("   ADAPTIVE_STEP_SIZE = %d\n",ADAPTIVE_STEP_SIZE);
	printf("   ETA_INCREASE       = %2.8f\n",ETA_INCREASE);
	printf("   ETA_MAX            = %2.8f\n",ETA_MAX);
	printf("   ETA_MIN            = %2.12f\n",ETA_MIN);
	printf("   MOMENTUM          = %f\n",MOMENTUM);
	printf("   USE_INITS         = %d\n",USE_INITS);
	printf("   W1_INIT_FILE      = %s\n",W1_INIT_FILE);
	printf("   B1_INIT_FILE      = %s\n",B1_INIT_FILE);
	printf("   W2_INIT_FILE      = %s\n",W2_INIT_FILE);
	printf("   B2_INIT_FILE      = %s\n",B2_INIT_FILE);
	printf("   WF_INIT_FILE      = %s\n",WF_INIT_FILE);
	printf("   WEIGHT_OUT_FILE   = %s\n",WEIGHT_OUT_FILE);
	printf("   TRAINING_LENGTH   = %d\n",TRAINING_LENGTH);
	printf("   UPDATES_PER_EPOCH = %d\n",UPDATES_PER_EPOCH);

} /* end of show Parameters */


/*****************************************************************************/
/* read W1 b1 W2 b2 Wf from separate files */
void getWeights()
{
	FILE	*initfile;
	int		i = 0;
	int		j = 0;
	int		k = 0;

	initfile = fopen(W1_INIT_FILE,"r");
	if (initfile == NULL)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not opened for reading.\n",W1_INIT_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (j=0;j<N1;j++)
	{
		for(i=0;i<N0;i++)
		{	
			fscanf(initfile,"%f",W1+(k++));
		}
	}

	/* close the initialization file */
    if (fclose(initfile) != 0)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be closed properly.\n",W1_INIT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}

	initfile = fopen(B1_INIT_FILE,"r");
	if (initfile == NULL)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be opened for reading.\n",B1_INIT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for(j=0;j<N1;j++)
	{
			fscanf(initfile,"%f",b1+(k++));
	}

	/* close the initialization file */
    if (fclose(initfile) != 0)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be closed properly.\n",B1_INIT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}

	initfile = fopen(W2_INIT_FILE,"r");
	if (initfile == NULL)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not opened for reading.\n",W2_INIT_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (i=0;i<N2;i++)
	{
		for(j=0;j<N1;j++)
		{
			fscanf(initfile,"%f",W2+(k++));
		}
	}

	/* close the initialization file */
    if (fclose(initfile) != 0)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be closed properly.\n",W2_INIT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}

	initfile = fopen(B2_INIT_FILE,"r");
	if (initfile == NULL)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be opened for reading.\n",B2_INIT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for(j=0;j<N2;j++)
	{
			fscanf(initfile,"%f",b2+(k++));
	}

	/* close the output file */
    if (fclose(initfile) != 0)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be closed properly.\n",B2_INIT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}

	initfile = fopen(WF_INIT_FILE,"r");
	if (initfile == NULL)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not opened for reading.\n",WF_INIT_FILE);
		printf("Exiting with nothing done.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (i=0;i<N1;i++)
	{
		for(j=0;j<N1;j++)
		{
			fscanf(initfile,"%f",Wf+(k++));
		}
	}

	/* close the initialization file */
    if (fclose(initfile) != 0)
	{
		printf("\n*ERROR* The initialization file \"%s\" could not be closed properly.\n",WF_INIT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}


} /* end of getWeights() */

/*****************************************************************************/
/* write W1 b1 W2 b2 Wf to separate files */
void putWeights()
{
	FILE	*outfile;
	char	OUT_FILE[FILE_NAME_MAX] = "";
	int		i = 0;
	int		j = 0;
	int		k = 0;

	strcpy(OUT_FILE,"w1"); strcat(OUT_FILE,WEIGHT_OUT_FILE);

	outfile = fopen(OUT_FILE,"w");
	if (outfile == NULL)
	{
		printf("\n*ERROR* The output file \"%s\" could not be created for writing.\n",OUT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (j=0;j<N1;j++)
	{
		for(i=0;i<N0;i++)
		{
			fprintf(outfile,"%e\t",W1[k++]);
		}
		fprintf(outfile,"\n");
	}

	/* close the output file */
    if (fclose(outfile) != 0)
	{
		printf("\n*ERROR* The output file \"%s\" could not be closed properly.\n",OUT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}

	printf("\nCreated the file \"%s\"\n",OUT_FILE);

	strcpy(OUT_FILE,"b1"); strcat(OUT_FILE,WEIGHT_OUT_FILE);

	outfile = fopen(OUT_FILE,"w");
	if (outfile == NULL)
	{
		printf("\n*ERROR* The output file \"%s\" could not be created for writing.\n",OUT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	for(j=0;j<N1;j++)
	{
		fprintf(outfile,"%e\n",b1[j]);
	}

	/* close the output file */
    if (fclose(outfile) != 0)
	{
		printf("\n*ERROR* The output file \"%s\" could not be closed properly.\n",OUT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}
	printf("Created the file \"%s\"\n",OUT_FILE);

	strcpy(OUT_FILE,"w2"); strcat(OUT_FILE,WEIGHT_OUT_FILE);

	outfile = fopen(OUT_FILE,"w");
	if (outfile == NULL)
	{
		printf("\n*ERROR* The output file \"%s\" could not be created for writing.\n",OUT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (i=0;i<N2;i++)
	{
		for(j=0;j<N1;j++)
		{
			fprintf(outfile,"%e\t",W2[k++]);
		}
		fprintf(outfile,"\n");
	}

	/* close the output file */
    if (fclose(outfile) != 0)
	{
		printf("\n*ERROR* The output file \"%s\" could not be closed properly.\n",OUT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}
	printf("Created the file \"%s\"\n",OUT_FILE);

	strcpy(OUT_FILE,"b2"); strcat(OUT_FILE,WEIGHT_OUT_FILE);

	outfile = fopen(OUT_FILE,"w");
	if (outfile == NULL)
	{
		printf("\n*ERROR* The output file \"%s\" could not be created for writing.\n",OUT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	for(j=0;j<N2;j++)
	{
		fprintf(outfile,"%e\n",b2[j]);
	}

	/* close the output file */
    if (fclose(outfile) != 0)
	{
		printf("\n*ERROR* The output file \"%s\" could not be closed properly.\n",OUT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}
	printf("Created the file \"%s\"\n",OUT_FILE);

	strcpy(OUT_FILE,"wf"); strcat(OUT_FILE,WEIGHT_OUT_FILE);

	outfile = fopen(OUT_FILE,"w");
	if (outfile == NULL)
	{
		printf("\n*ERROR* The output file \"%s\" could not be created for writing.\n",OUT_FILE);
		printf("Exiting without writing the file.\n");
		exit(EXIT_FAILURE);
	}

	k = 0;
	for (i=0;i<N1;i++)
	{
		for(j=0;j<N1;j++)
		{
			fprintf(outfile,"%e\t",Wf[k++]);
		}
		fprintf(outfile,"\n");
	}

	/* close the output file */
    if (fclose(outfile) != 0)
	{
		printf("\n*ERROR* The output file \"%s\" could not be closed properly.\n",OUT_FILE);
		printf("Exiting.\n");
		exit(EXIT_FAILURE);
	}
	printf("Created the file \"%s\"\n",OUT_FILE);

} /* end of putWeights() */


float randomNumber(void)
{
	static unsigned long jflone = 0x3F800000;
	static unsigned long jflmsk = 0x007FFFFF;

	static unsigned long idum = 0L;
	unsigned long itemp;
	float rand;

	if (idum==0L)
	{
		idum = time(0);
	}

	idum = 1664525L*idum + 1013904223L;
	itemp = jflone | (jflmsk & idum);
	rand = (*(float *)&itemp)-1.0;

	return(rand);
}