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

Title  : 	Linear Network (Classification)
Inputs : 	1. Input File Containing Training Data         (Infile)
				2. No. of Inputs                               (N)
				3. Possible No. of Classes                     (Nc)

Output :    classification error %

Notations Used:
	Nv   --> No. of patterns
	N    --> No. of inputs
	Nout --> No. of outputs
	X    --> Input Vector
	Xa   --> Augmented input vector
	Wo	  --> Weight vector connecting input to the output
	to	  --> Output threshold
	ty	  --> Desired output
	y	  --> Actual output
	R	  --> Auto correlation Matrix
	C	  --> Cross correlation matrix
	dstd --> desired mean
	dstd --> desired standard deviation


********************************************************************************
! IMPORTANT: Significant changes have been made to the program. Have a
				 thorough look at it.

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



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

#define PI 3.1415926
#define dmean 0.0
#define dstd 0.5


int IX=3;
int IY=4009;
int IZ=234;

char Infile[14];

void main(void)
{

	int i, j;
	int N, Nc, Nv;
	float *X, *Xa, *ty;
	float **Wo,*to, **W, **R, **C;
	int ClassId;
	FILE *ifs;

	void CGrad(int ,float * ,float **,float *);
	void error(float **, int, int);
	void classification_error(float **, int, int);
	float slete(float, float);


	printf("Enter Input data file : ");
	gets(Infile);

	printf("Enter No. of Inputs, N : ");
	scanf("%d",&N);

	printf("Enter Possible no. of classes, Nc : ");
	scanf("%d",&Nc);


	ifs=fopen(Infile,"r");

	if (ifs == NULL ) {
		perror(Infile);
		exit(1);
	}


	/*	Allocating memory dynamically for the variables used */

	X= (float *)malloc(sizeof(float)*(N) );
	Xa= (float *)malloc(sizeof(float)*(N+1) );
	ty= (float *)malloc(sizeof(float)*Nc );
	W= (float **)malloc(sizeof(float*)*Nc);
	Wo= (float **)malloc(sizeof(float*)*Nc);
	to= (float *)malloc(sizeof(float)*Nc );
	R= (float **)malloc(sizeof(float*)*(N+1));
	C= (float **)malloc(sizeof(float*)*Nc);

	for (i=0;i<N+1;i++) {
		R[i]= (float*)malloc(sizeof(float)*(N+1));
	}

	for (i=0;i<Nc;i++) {
		C[i]=(float* )malloc(sizeof(float)*(N+1));
	}

	for (i=0;i<Nc;i++) {
		Wo[i]= (float*)malloc(sizeof(float)*N);
	}

	for (i=0;i<Nc;i++) {
		W[i]= (float*)malloc(sizeof(float)*(N+1));
	}


	/* Initializations */


	/* Initialzing weights and thresholds to gaussian random nos. */

	for (i=0;i<Nc;i++){

		to[i]=slete(dstd, dmean);

		for (j=0;j<N;j++){
			Wo[i][j]=slete(dstd, dmean);

		}
	}

	/* W --> Weight vector containing Wo and to */

	for (i=0;i<Nc;i++){

		for (j=0;j<N;j++){
			W[i][j]=Wo[i][j];
		}
		W[i][N]= to[i];

	}



	/* Initializing R and C with zeroes */

	for (i=0;i<N+1;i++) {
		for (j=0;j<N+1;j++) {
			R[i][j]=0.0;
		}
	}

	for (i=0;i<Nc;i++) {
		for (j=0;j<N+1;j++) {
			C[i][j]=0.0;
		}
	}

	Nv=0;



	/*	Here we find the auto and cross correlation matrices */

	while (!feof(ifs))
	{
		for (i=0;i<N && !feof(ifs);i++) {
			fscanf(ifs,"%f",&X[i]);
		}

		if (!feof(ifs) ) {
			fscanf(ifs,"%d", &ClassId);
      }

		for (i=0;i<N;i++) {
			Xa[i]= X[i];
      }

		Xa[N]=1.0; 	

		for (i=0;i<Nc;i++) {
			ty[i]=0.0;
		}

		ty[ClassId-1]=1.0;

		/* The o/p node corr. to the class id is made equal to 1 */
		/* all other nodes have desired o/p = 0 */

		/* Since C allocates memory starting from zero and the classes are usually
		named from 1 onwards, the node ty[Class_Id -1] is made equal to 1 instead
		of just ty[Class_Id] */


		/* Finding the auto-correlation */

		for (i=0;i<N+1;i++) {
			for (j=0;j<N+1;j++) {
				R[i][j]+=Xa[i]*Xa[j];
			}
		}


		/* Finding the cross-correlation */

		for (i=0;i<Nc;i++) {
			for (j=0;j<N+1;j++) {
				C[i][j]+=ty[i]*Xa[j];
			}
		}

		Nv++;

	}  									/* end of while loop */

	fclose(ifs);
	Nv--;

	/*	 Normalizing auto and cross correlations */

	for (i=0;i<N+1;i++) {
		for (j=0;j<N+1;j++) {
			R[i][j]/=Nv;
		}
	}

	for (i=0;i<Nc;i++) {
		for (j=0;j<N+1;j++) {
			C[i][j]/=Nv;
		}
	}


	/*	Calling Conjugate gradient Subroutine */

	for (i=0;i<Nc;i++) {
		CGrad(N+1, W[i], R, C[i]);
	}


	/* calling "error" subroutine that calculates mapping error */
	error(W, N, Nc);

	/* Calling subroutine that calculates the classification error */

	classification_error(W, N, Nc);

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


}  // End of main






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

Title    : Subroutine- Conjugate gradient
Comments : Follow the derivation done in the class for better understanding of
			  this routine

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



void CGrad(int Nu,float *w,float **r,float *c)
{

	float XD, XN, Num, Den, tempg, B1, B2;
	float *p, *g;
	int l, m , iter, pass;


	XD = 1;


	p= (float *)malloc(sizeof(float)*Nu);
	g= (float *)malloc(sizeof(float)*Nu);


	/* passing through the CG routine two times for improved results */
	for (pass =0;pass<2;pass++)
	{

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

		for(iter= 0; iter<Nu; iter++)	  					/* start of Iteration loop */
		{
			for(l=0;l<Nu;l++) {  		  					/* start of l -loop */
				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;
			}                                         /* end of l -loop*/

			XN = 0;

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

			B1 = XN/XD;
			XD = XN;

			/* direction vector */
			for(l=0;l<Nu;l++) {
				p[l] = -g[l]+B1*p[l];
			}

			Den = Num = 0.0;

			for(l=0;l<Nu;l++)										/* start of l - loop */
			{
				/* 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];
				}

			}                                   			/* end of   l - loop */

			B2 = Num/Den;


			for(l=0;l<Nu;l++) {
				w[l]+= B2*p[l];								/* Updating the weights */
			}

		}													/* end of Iteration loop */

	}												  /* end of pass - loop */

}                              /* end of CjGrad */





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

This subroutine determines the difference between the desired and obtained
output. This requires passing through the data once.

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

void error(float **W, int N2, int Nout2)
{

	float *Error, y ;
	float *x, *xa, *ty;
	int i, j;
	int classid;
	FILE *ifs1;
	float MSE =0.0;
	int Nv=0;;


	/* memory allocation */

	x= (float *)malloc(sizeof(float)*(N2));
	ty= (float *)malloc(sizeof(float)* Nout2);
	xa= (float *)malloc(sizeof(float)*(N2+1));
	Error= (float *)malloc(sizeof(float)* Nout2);


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

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

	/* Initializing error to zeroes */
	for (i=0;i<Nout2;i++) {
		Error[i]=0.0;
	}

	while (!feof(ifs1) )
	{

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

		if ( ! feof(ifs1) ) {
			fscanf(ifs1,"%d",&classid);
		}

		for (i=0;i<N2;i++) {
			xa[i] =x[i];
		}

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

		ty[classid-1]=1.0;

		xa[N2]=1.0; 							/* Threshold */


		/* for each set of input features, calculating and accumulating error */
		for (i=0;i<Nout2;i++)
		{

			y = 0.0;

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

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

	}										/* end of while loop */
	fclose(ifs1);
	Nv--;

	/* printing errors at each node */

	printf("\nMapping error:\n");
	printf("---------------------------\n");

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

	printf("---------------------------\n");
	printf("Total MSE =      %f\n",MSE);
	printf("---------------------------\n\n");

}







/******************************************************************************
Title : subroutine- Classification error
Comments : This subroutine calculates the classification error percentage 

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


void classification_error(float **W1, int N1, int Nc1)
{
	float *y;
	float Max_y;
	int Class_obtained;
	int Class_err=0;
	float *X;
	int ClassId;
	int i, j;
	FILE *ifsc;
	int Nv;
	float percentage_error;
	Nv=0;

	/* Memory allocation */
	X=( float *)malloc(sizeof(float)*N1);
	y= (float *) malloc ( sizeof(float) *Nc1);

	for (i=0;i<Nc1;i++)
		y[i]=0.0;

	ifsc = fopen(Infile,"r");

	while ( !feof(ifsc) )
	{

		Nv++;

		for (i=0;i<N1 && !feof(ifsc);i++) {
			fscanf(ifsc,"%f", &X[i]);
      }

		X[N1]=1.0;

		if ( !feof(ifsc) ) {
			fscanf(ifsc,"%d",&ClassId);
      }

		for (i=0;i<Nc1;i++) {
			y[i]=0.0;
      }

		for (i=0;i<Nc1;i++) {
			for (j=0;j<N1+1;j++) {
				y[i]+= X[j]*W1[i][j];
			}
      }

/*		Here we find the maximum of the outputs:
		The node with the maximm o/p represents the desired class */

		Max_y= y[0];
		Class_obtained=1;

		for (i=1;i<Nc1;i++)
		{
			if(y[i]>Max_y) {
				Max_y= y[i];
				Class_obtained= i+1;
			}
		}

/*	   Checking whether the obtained class is equal to the desired class or not.
		if they are not equal it means that an error has occurred. so we increment
		the classification error count by 1 */

		if ( ClassId != Class_obtained ) {
			Class_err += 1;
		}

	}


	/* calculating classification error% */
	percentage_error=((float)Class_err/(float)Nv)*100.0;

	printf("\nclassification error = %d\n",Class_err);
	printf("\nclassification error percentage = %f \n",percentage_error);

	fclose(ifsc);

} // end of  classification error subroutine





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

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

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


float rand1(int *ix, int *iy, int *iz)
{
	int ixx, iyy, izz;
	float itemp;
	float 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=(float)(*ix)/30629.0+(float)(*iy)/30307.0+(float)(*iz)/30323.0;
	itemp=floor(temp);
	return (temp-itemp);

}





/* Gaussian Random No. generator */

float slete(float std, float xmean)
{
	float rand1(int *, int *, int *);
	return (xmean+std*cos(2*PI*rand1(&IX, &IY, &IZ))*sqrt(-2.0*log(rand1(&IX, &IY, &IZ))));

}