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

// FLNTESTING------Using the Conjugate Gradient Menehod 
/*
	Nv   --> No. of patterns
	N    --> No. of inputs
	Nout --> No. of outputs
	x    --> Input Vector
	xa   --> Augmented input vector including all the degree of the input vector
	W	  --> I/P to O/P Weights		
	to	  --> Output threshold          
	ty	  --> Desired output
	y	  --> Actual output
	R	  --> Auto-correlation matrix
	C    --> Cross-correlation matrix
	d   --> the degree of the network
	L   --> the demention of the augmented input vector
*/

/* The bigest d = 3*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <direct.h>
#include <string.h>
                 
#define dmean 0.0
#define dstd 0.5
#define WITHOUTPUT  1
#define NOOUTPUT    0
int IX=3;
int IY=4009;
int IZ=234;



void main(void)
{

	char Infile[100],Infile1[100],str[100],str1[100];
	FILE *ifs;
	int i, j,jj, d, L;
	int Nout, N, Nv;
	double *XD;
	double *x, *xa, *ty,*y, *MSE,*Error;
	double *to, **W;
	double *p, *g;
	double TotalError;

	void CalculateGradient(int, double *, double *, double **, double *);
	void CalculateDirection(int, double *, double *, double *);
	void CalculateB2(int ,double *, double *, double **, double *);
	void error(int, int, int, int, char *, double **, double *);
	void Getaugmentedvector(double *, double *,int,int, int);
	double slete(double, double);
	double rand1(int *, int *, int *);
	int factorial(int, int);

	printf("Enter the Input File Name: ");
	scanf("%s", Infile);

	printf("Enter the Weights File Name : ");
	scanf("%s", Infile1);

	
	/* Opening the input file in "read" mode */
	ifs= fopen(Infile1,"r");


	/* Checking for the existance of the file */
	if (ifs == NULL )
	{
		perror(Infile1);
		exit(1);
	}
	fscanf(ifs,"%s",&str);
	fscanf(ifs,"%d",&N);
	fscanf(ifs,"%d",&Nout);
	fscanf(ifs,"%d",&d);
	
	// calculate the dimension of the input vector
	L = factorial(N,d);
 	/* Allocating Memory dynamically*/
	

	x= (double *)malloc(sizeof(double)*(N));
	y = (double *)malloc(sizeof(double)*Nout);
	xa= (double *)malloc(sizeof(double)*L);
	MSE = (double *)malloc(sizeof(double)*1);
	Error = (double *)malloc(sizeof(double)*Nout);
	p= (double *)malloc(sizeof(double)*L);
	g= (double *)malloc(sizeof(double)*L);
	ty= (double *)malloc(sizeof(double)*Nout );
	to= (double *)malloc(sizeof(double)*Nout );
	XD = (double *)malloc(sizeof(double)*1);

	W= (double **)malloc(sizeof(double *)*Nout);

	for (i=0;i<Nout;i++) {
		W[i]= (double *)malloc(sizeof(double )*L);
	}

	/* Initialize the weights */
	for(i = 0 ; i < Nout; i ++)
	{
		for( j = 0 ; j < L; j ++)
		{
			fscanf(ifs,"%s",str1);
			W[i][j] = atof(str1);
			
		}
	
	}


	Nv= 0;
	ifs= fopen(Infile,"r");


	/* Checking for the existance of the file */
	if (ifs == NULL )
	{
		perror(Infile);
		exit(1);
	}

	while(!feof(ifs))
	{
		Nv ++;
		for(i = 0; i < N ; i++)
		{
			fscanf(ifs,"%s",str);
			x[i] = atof(str);
			xa[i] = x[i];
		}
		for(i = 0; i < Nout; i ++)
		{
			fscanf(ifs,"%s",str);
			ty[i] = atof(str);
		}
		
	}
	fclose(ifs);
	Nv--;


	printf("\n\n\n\t\tFLN  Process  Program");
	printf("\n\n\tTo be processed File Name:%s\t\t\t",Infile);
	printf("\n\n\tWeights File Name:\t%s",Infile1);
	printf("\n\n\tNo of Inputs:\t%d",N);
	printf("\n\n\tNo of Outputs:\t%d",Nout);
	
	for(i = 0; i < Nout; i ++)
	{
		Error[i] = 0.0;
	}
	TotalError = 0.0;
	ifs= fopen(Infile,"r");


	/* Checking for the existance of the file */
	if (ifs == NULL )
	{
		perror(Infile);
		exit(1);
	}
	for(jj=0;jj<Nv;jj++)
	{
		for(i = 0; i < N ; i++)
		{
			fscanf(ifs,"%s",str);
			x[i] = atof(str);
			xa[i] = x[i];
		}
		for(i = 0; i < Nout; i ++)
			{
				fscanf(ifs,"%s",str);
				ty[i] = atof(str);
			}
		Getaugmentedvector(x, xa, d, N, L);
		for(i = 0; i < Nout; i ++)
		{
			y[i] = 0.0;
			for(j = 0;j < L; j ++)
			{
				y[i] += W[i][j]*xa[j];
			}
			Error[i] += (ty[i]-y[i])*(ty[i]-y[i]);
		}
	}
	fclose(ifs);

	error(N, Nout, L, d, Infile, W, MSE);
	printf("\n\n\tThe total Error is: \t%f",MSE[0]);

	printf("\n\nTraining complete !\n");

}																	/* end of main */










/******Sub-routine: Conjugate Gradient****************************************/



void CalculateB2(int Nu,double *p, double *g,double **r,double *w)
{

	int l, m;
	double Num, Den, B2;

	Num = Den = 0.0;
	B2 = 0;

	for(l=0;l<Nu;l++)								
	{
		/* Num --> numerator of B2 */
		Num += p[l]*g[l]/-2.0;
		/* Den --> denominator of B2 */
		for(m=0;m<Nu;m++) 
		{
			Den += p[m]*p[l]*r[m][l];
		}
	}
		B2 = Num/Den;
		//update the weights
	
		for(l=0;l<Nu;l++) 
		{
			w[l]+= B2*p[l];								/* Updating the weights */
		}                                   			
	             							
}                             












/******** Subroutine :error ***************************************************

This subroutine calculates the mean square error at each node and prints them.

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

void error(int N2, int Nout2, int L, int d, char *Infile, double **W, double * MSE)
{

	  FILE  *ifs;
	  double *Error, y;
	  double *x, *xa, *ty;
	  int i, j;
      int Nv=0;
	   MSE[0]=0.0;

	  /* Memory allocation */

	  x = (double *) malloc(sizeof(double)*(N2));
	  xa = (double *) malloc(sizeof(double)*(L));
	  ty = (double *) malloc(sizeof(double)*(Nout2));
	  Error = (double *) malloc(sizeof(double)*(Nout2));


	  /* opening the input file in read mode */	
	  ifs= fopen(Infile,"r");


	  for (i=0;i<Nout2;i++){
		  Error[i]=0.0;
	  }

	  while(!feof(ifs))
	  {
     		Nv++;
			for(i=0;i<N2 && !feof(ifs);i++) {
				fscanf(ifs,"%lf",&x[i]);
			}

			for(i=0;i<Nout2 && !feof(ifs);i++) {
				fscanf(ifs,"%lf",&ty[i]);
			}

			Getaugmentedvector(x, xa, d, N2, L);


		/*here we calculate the error for all the output nodes at the same time*/

			for (i=0;i<Nout2;i++)
			{
				y = 0.0;

				for(j=0;j<L;j++) {
					y += xa[j]*W[i][j];
				}

				Error[i]+= (ty[i]-y)*(ty[i]-y);
			}

			/* y= actual output  */
			/* t= desired output */


	}
	fclose(ifs);
   Nv--;

	for (i=0;i<Nout2;i++) {
		MSE[0]+=Error[i]/Nv;
		//printf("\nError at node %d : %f",(i+1),Error[i]/Nv);
	}
 
//	printf("\n\nTotal MSE : %f",MSE);

} 	/*End of error() */


/******Sub-routine: Conjugate Gradient****************************************/



void CalculateGradient(int Nu,double * g,double *w,double **r,double *c)
{
	
	int l,m;
	double tempg;

	for(l=0;l<Nu;l++) 
	{
		g[l] = 0.0;
	}

	
	for(l=0;l<Nu;l++)
	{  		  					
		tempg = 0.0;
		for(m=0;m<Nu;m++)
		{
			tempg += w[m]*r[m][l];
			
		}
		/* gradient vector */
		g[l] = -2.0*c[l]+2.0*tempg;
	}
	

}


/********Get the Direction Vector***********************************************/

void CalculateDirection(int Nu, double *XD, double *p, double *g)
{
	
	int i;
	double XN = 0.0, B1;
/*	for(i = 0;i<Nu;i++)
	{
		p[i] = 0.0;
	}
*/	
	for(i = 0; i < Nu; i ++)
	{
		XN += (g[i]*g[i]);
	}

	B1 = XN/XD[0];
	XD[0] = XN;
	
	for(i = 0; i < Nu; i ++)
	{
		p[i] = -g[i] + B1*p[i];
	}

}


/************************************************************************************/
 factorial(int N, int d)
{
	int result = 1;
	int i =0;

	if (d == 0)
		return result;
	else{
			for (i = N+d; i > N; i --)
			{
				result *= i;
			}
			for (i = 1; i <= d; i ++)
			{
				result /= i;
			}
			return result;
		}
}

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

	Subroutine : Random no. Generator
		-the random nos. generated are uniformly distributed
		 between 0 and 1.
		-IX, IY, IZ are the SEEDS

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


double rand1(int *ix, int *iy, int *iz)
{
	int ixx, iyy, izz;
	double itemp;
	double temp;


	ixx=(*ix)/177;
	*ix=171*(*ix%177)-2*ixx;

	if(*ix < 0)
		*ix+=30269;

	iyy=(*iy)/176;
	*iy=176*(*iy%176)-2*iyy;

	if(*iy < 0)
		*iy+=30307;

	izz=(*iz)/178;
	*iz=170*(*iz%178)-2*izz;

	if(*iz < 0)
		*iz+=30323;

	temp=(double)(*ix)/30629.0+(double)(*iy)/30307.0+(double)(*iz)/30323.0;
	itemp=floor(temp);
	return (temp-itemp);

}




/********* Gaussian Random No. generator **************************************/


double slete(double std, double xmean)
{
	double cc=3.1415926;
	return xmean+std*cos(2*cc*rand1(&IX, &IY, &IZ))*sqrt(-2.0*log(rand1(&IX, &IY, &IZ)));

}

/*******************************************************************************/
void Getaugmentedvector(double *x, double *xa, int d, int N, int L)
{
	int i,j,k,m;
	for(i = 0; i < N ; i ++)
		xa[i] = x[i];//first degree
	m = N;
	if(d==1){xa[N] = 1.0; return;}
	for (i = 0; i < N; i++){
		for(j = 0; j <= i; j ++){
			xa[m] = x[i]*x[j];
			m++;}}//second degree
	if(d==2){xa[m]=1.0;return;}
	for(i = 0; i < N; i++){
		for(j = 0; j <= i; j ++){
			for(k = 0; k <= j; k++){
				xa[m] = x[i]*x[j]*x[k];
				m++;}}}//third degree
	xa[m]=1.0;
	return;
	//Right now I only calculate the vector when the d =2
/*	int i, j, k, l, index;
	double temp,ljl;
	int temp1;
	int flag = 1;
	for(i = 0; i < N; i ++){
		xa[i] = x[i];
		}
	if(d==1){
		xa[L-1] =1.0;
		return;
	}
	else{
		for(i = 1; i <d; i ++)
		{
			index = factorial(N, i)-1;
			temp1 = factorial(N,i)-factorial(N,i-1);//the No. of the most newer vector
			for(k = 0; k < N; k ++)
				{
				for(j = 0; j < temp1; j ++)
				{
				
					temp = xa[j+factorial(N,i-1)-1]*x[k];
					for(l = factorial(N, i)-2; l < index; l ++)
					{
						ljl=(temp - xa[l]);
						if(ljl==0.0)
						{
							flag = 0;
							break;
						}
						if(flag == 0)
							break;
					}
					if(flag == 1)
					{
						xa[index] = temp;
						index ++;
					}
					flag = 1;

				}
			}
			
		}
		xa[L-1] = 1.0;
	}
*/
 }

