//neural.cpp, the "reincarnated" version.

//(c) 1999, Oguz Yetkin

//NO OBJECTS IN THIS ONE.

//the header was last modified 1/11/1999

//when the program was last modified is anyone's guess :)

//

//WARNING! CURRENTLY PROGRAM DOESN'T TEST EXISTING NETS FOR ERROR, AND IS NOT

//SUITEDFOR DISTRIBUTED TRAINING!



#include<iostream.h>

#include<fstream.h>

#include<stdio.h>

#include<stdlib.h>

#include<string.h>

#include<math.h>

#include<time.h>



#include "array2d.h"

#include "neural.h"





#define PI 3.141592654



//for export of data to outside applications

double** trainset;

double** trainnew;	//newset



const int MAX=64;

int num=5;



//we'll use globals for speed and simplicity



long n=32;

long m=32;

double s=1;

double d0=0.000001;

long gaussnum=32;





array2d* wts;

double* x;

double* y;

double* xe;

double* ye;



//globals



long Gimprovements=0;

long Gsteps=0;

double Gerror=0;

double Gle=0;

long   Gleiters=200;	//iterations for lyapunov calculation

double Gmisses=0;

double Gconstraint=0;

double Gcalculate_lyapunov=1;

double Gtemperature=1;

//added 2/8/1999
double Glower=0;	//lower LE constraint
double Gupper=0;	//upper LE constraint






//1 for passing constraint, 0 for failing it.

int constrain(double constraint, double le){

	if(Glower!=0||Gupper!=0){
		//ignore others, work on global constraints
		if(le>Gupper){
			return 0;
		}
		if(le<Glower){
			return 0;
		}
		return 1;

	}
	//we'll never get here if the above is executed

	if(constraint==0){

		return 1;

	}



	

	// the following should work, but the opposite does...go figure

	if(constraint==1 && le>0){

		return 1;

	}

	if(constraint==-1 && le<=0){

		return 1;

	}

	if(le==0){

		return 1; //just to get things going...

	}

	

	return 0;

}



void report(void){	//output to screen

	static int firstrun=1;

	if(firstrun==1){

		firstrun=0;

		cout<<"%improvements steps error le temp misses constraint:"<<Gconstraint<<endl;

	}



	cout<<endl<<Gimprovements<<" "<<Gsteps<<" "<<Gerror<<" "<<Gle<<" "<<Gtemperature<<" "<<Gmisses;



}



//assumes input and output are of size n and wts is n by m

void niterate(double* output, double* input, double** wts, long n, long m, double s){

	long i,j,k;

	i=j=k=0;

	double u=0;

	//cout<<"\n* ";

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

	   u=0;

	   for(j=0;j<m;j++){

	     u=u+wts[i][j]*input[j];

	   }

	   output[i]=tanh(s*u);

	  // cout<<output[i]<<" ";



	}

//	output[i]=tanh(s*u);

}



//d0 is global

double lyapunov(long n, long m, double s, double** wts, long iters){



	long i,j,k;

	i=j=k=0;   



	double* x=new double[n];

	double* xe=new double[n];

	double* y=new double[n];

	double* ye=new double[n];

	double dsum=0;

	double lsum=0;

	double le=0;

	double d1=0;

	//d0 is global



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

		x[i]=y[i]=1;

		xe[i]=ye[i]=x[i]+d0/(float)sqrt(m);

	}

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

		for(j=0;j<n;j++){

			x[j]=y[j];

			xe[j]=ye[j];

		}

		

		

		niterate(y,x, wts, n, m,s);

		niterate(ye, xe, wts, n, m, s);

		dsum=0;

		for(j=0;j<n;j++){

		  dsum=dsum+(y[j]-ye[j])*(y[j]-ye[j]);

		 }

		 d1=sqrt(dsum);

		 if(d1==0){

			 cerr<<"\n**LE too small, returning 0";

			 return 0;

		 }

		 for(j=0;j<n;j++){ //readjust ye for separation d0

			ye[j]=y[j]+(d0/d1)*(ye[j]-y[j]); 

		 }

		 

		lsum=lsum+log(fabs(d1/d0));

		le=lsum/(double)i;  //running average

	}

	//patched 8/7/1999 to fix memory leak
	delete[] x;
	delete[] xe;
	delete[] y;
	delete[] ye;

	return le;

}





inline double distanceN(double* vector1,double* vector2, long n){



	long i,j;

	i=j=0;

	long double innersum=0;

	long double distance=0;



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

		innersum+=pow((vector1[i]-vector2[i]),2);//square

	}



	distance=sqrt(innersum);

	//cout<<distance<<endl;

	return distance;

}



//error is used to calculate the error between "width" timeseries

//containing "height" points.  

//specifically, This is the RMS error of 2 matrices containing width column

//vectors of length height.  The distance of each columnvector in array1

//and array2 is computed, and the square root of this value is taken



inline double error(double** array1, double** array2, int width, int height){



	double sum=0;

	for(int i=0;i<width;i++){//for each column vector

		sum+=distanceN(array1[i], array2[i], height); //sum the dist between each vector

	}

	sum=sqrt(sum);

	return sum;

}



//train()

//trains until it finds a solution better than error

//returns: new error

//action: trains for one "improvement"

//inputs: tset, a matrix consisting of N timeseries (represented as column vectors)

//		        each containing tlength points

//		  wts,	an n by m network of weights (to be modified)

//		  n, m,s (number of neurons, number of connections, scaling factor)

//		  tlength (number of time points in training set, old error)

//		  temp, cooldown(temperature to start with, cooldown factor);



double train(double** tset, double** wts, int n, int m, double s, double old_error,

			 int tlength, double& temp, double cooldown){



			//we'll assume wts is disposable.  We'll let the caller

			//worry about saving a copy

			

			static firstrun=1;

			long i,j,k;

			i=j=k=0;

			static array2d* newset;	//new set of time series generated

			static array2d* scratch;	//temporary copy of network

			static double* output;

			static double* input;

			static oldn=0;

			static oldm=0;

			static oldtlength=0;



			if(firstrun==1){

					newset=new array2d(n, tlength);

					scratch=new array2d(n,m);

					output=new double[n];//the output vector (temporary variables)

					input=new double[n];//the vector to be used as input



					oldn=n;

					oldm=m;

					oldtlength=tlength;



		

				firstrun=0;

			}

			//don't allocate unless we need to 

			if(oldn!=n||oldm!=m||oldtlength!=tlength){//reallocate everything

					cerr<<"\n**reallocating everything!"<<endl;

					delete newset;

					//patch added 8/7/1999 to fix mem leak
					delete scratch;

					delete[] output;	//these could be static, saving some time

					delete[] input;

					newset=new array2d(n, tlength);

					scratch=new array2d(n,m);

					output=new double[n];//the output vector (temporary variables)

					input=new double[n];//the vector to be used as input



					oldn=n;

					oldm=m;

					oldtlength=tlength;



			}



			//double old_error=0;

			double new_error=old_error;



		

		

			

			//input will remain the same during this training run

			//no it won't! realized: 1/10/1999

			for(j=0;j<n;j++){

				input[j]=tset[j][0];	//get the first value in each column

			}

			



			//get the initial error

			//cout<<"\ngetting initial error"<<endl;



			/* commented out 1/4/1998

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

			  niterate(output,input, wts, n,m,s);

				  for(j=0;j<n;j++){

					newset->base[j][i]=output[j];

				  }

				

			}

			*/



			//old_error=error(tset, newset->base, n, tlength);

			//new_error=old_error;

			//cout<<"\ninit error is: "<<old_error<<endl;

			//now try to make it better





			//optimization added 1/11/1999. 

			//Determine old_error from net, and then do something

			//about it...

			/* start optimization */

			if(old_error==0){

				cerr<<"\ndetermining init error!"<<endl;

				//initialize the input to be 

				//the same as tset[][0]

				for(j=0;j<n;j++){

					input[j]=tset[j][0];	//get the first value in each column

				}



				//get newoutput to compare with training set

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

				  //output is what we get, input is always tset[][0]

				  niterate(output,input, wts, n,m,s);

				  					

				  for(j=0;j<n;j++){

						newset->base[i][j]=output[j];

						input[j]=output[j];	//added 1/6/1998, reassign!

				  }

				 

			

				}

				//calculate new error

				new_error=error(tset, newset->base, n, tlength);

				cerr<<"\noriginalerror determined to be: "<<new_error<<endl;

			/* end optimization */



			old_error=new_error;	

			}

			while(new_error>=old_error){

				Gsteps++;	//global variable for stats tracking

				

				//printf("\n**new_error:%f>old_error:%f temp: %f",new_error, old_error,temp); 		

				//create copy of network

				



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

					for(j=0;j<m;j++){

						//scratch->base[i][j]=wts[i][j];

						//changed 1/4/1998

						//changed 1/6/1998--let caller handle zeroing out!

						//scratch->base[i][j]=0;

						scratch->base[i][j]=wts[i][j];

					}

					

				}

				

			

				//perturb the weights

				//cout<<"\nperturbing weights!"<<endl;

			lyapunov_bailout:

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

					//cout<<".";		

					for(j=0;j<m;j++){

						

						//scratch->base[i][j]+=((temp*(rand()/(double)RAND_MAX))/(float)sqrt(m));

						//division of gauss() by sqrt(m) added 1/10/1999

						scratch->base[i][j]+=temp*gauss()/(float)sqrt(m); 



					}

					

				}

			

				//cout<<"\nweights perturbed!"<<endl;



				//initialize the input to be 

				//the same as tset[][0]

				for(j=0;j<n;j++){

					input[j]=tset[j][0];	//get the first value in each column

				}



					//get newoutput to compare with training set

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

				  //output is what we get, input is always tset[][0]

				  niterate(output,input, scratch->base, n,m,s);

				  					

				  for(j=0;j<n;j++){

						newset->base[i][j]=output[j];

						input[j]=output[j];	//added 1/6/1998, reassign!

				  }

				 

			

				}

				/*

				cout<<"\n***dumping input: "<<endl;

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

					cout<<input[i]<<" ";

				}

				cout<<"\n***dumping output: "<<endl;

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

					cout<<output[i]<<" ";

				}

				cout<<"\n***dumping newset[i][n-1]: "<<endl;

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

					cout<<newset->base[i][n-1]<<" ";

				}

				*/





			::trainset=tset;	//set global variables for comm w/ outside

			:: trainnew=newset->base;





				//calculate new error

				new_error=error(tset, newset->base, n, tlength);

				//printf("\n--t: %f e:%f !<goal e: %f",temp, new_error, old_error);

				

				

				//lower the temperature  [this could go outside the loop]

				if(new_error<old_error){ 

					if(Gcalculate_lyapunov){

						Gle=lyapunov(n,m,s,scratch->base,Gleiters);

					}

					if(constrain(Gconstraint, Gle)==0){

						

						//the way this works is backwards, 

						//don't ask me why Oguz 1/12/1999 2:39 am



						//cout<<"\nfailed!"<<endl;

						//cout<<Gconstraint<<endl;

						//cout<<Gle<<endl;

						Gmisses++;

						//cout<<"\n**bailing out, le: "<<Gle<<" err: "<<new_error<<endl;

						new_error=old_error+1;	//just so we can get out ASAP

						//cout<<"\n**bailing out, le: "<<Gle<<" err: "<<new_error<<endl;

						//goto lyapunov_bailout;

						goto lyapunov_bailout;

					}

					//copy scratch back to wts

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

						for(j=0;j<m;j++){

							wts[i][j]=scratch->base[i][j];

						}

					}

						//since we've found something better, double the temperature

					if(temp<1){

						temp*=2;

					}

				}else{

					/* label fail */

					//means scratch will not be copied to wts either because

					//of bad error or lyapunov constraint failure

fail:

					//temp*=1.1;

					temp*=cooldown; //cooldown must be less than one

				}



			







			}//end while

			





		

		

			return new_error;	

}



//I/O functions.  We assume user knows file size in network.

//this may change in the future



//written 1/10/1999 3:40 am

//neural net file has the form

int save(double** wts, int n, int m, double s, char* filename, char* header){

	int returnval=0,i,j;

	ofstream outfile;

	if(header==NULL){

		header=new char[MAX];

		sprintf(header,"%d %d %f",n,m,s); //write a minimal amt of info

	}

	outfile.open(filename);

	if(!outfile){

		returnval=0;

		goto end;	//bail out

	}



	outfile<<header<<endl;

	for( i=0;i<n;i++){	//for every neuron

		for( j=0;j<m;j++){//write a list of weights

			outfile<<wts[i][j]<<endl;

		}

	}

	outfile.close();

	returnval=1;

end:



	delete [] header;
	return returnval;

}





//caller may ignore header at will, assuming wts is allocated

int load(double** wts, int n, int m, char* filename, char* header){



	int returnval=0,i,j;

	ifstream infile;

	infile.open(filename);

	if(!infile){

		returnval=0;

		goto end;	//bail out

	}

	//header=new char[MAX];

	infile.getline(header,MAX);

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

		for(j=0;j<m;j++){

			infile>>wts[i][j];

			if(infile.eof()){

				goto end;	//bail out with error

			}



		}

	}





	returnval=1;

end:

	return returnval;

}



//takes training input from a file, and trains for "improvements" steps

//runs the trained net ever run_freq improvements for run_len iters if run_len

//is >0.  If run_len==0, then outfilename is ignored



//tset file looks like:

//n1 n2 n3...nN

//t0 t0 t0...t0

//t1 t1 t1...t1

//...

//tTlength...

//where every column is a desired time series entry for a neuron



int trainonfile(char* trainfilename, char*netfilename, char* outfilename,  

				 int n, int m, double s, long goal_improvements, 

				 long run_freq, long run_len, double& error){



	/* set the globals */

	Gimprovements=0;

	Gsteps=0;

	Gerror=0;

	Gle=0;

	Gmisses=0;

	if(Gconstraint!=0){

		Gcalculate_lyapunov=1;

	}

	



	long i,j,k;

	long improvements=0;

	int ok=0;

	ifstream trainfile;

	ofstream outfile;

	ofstream netfile;

	array2d* tset;

	array2d* wts;

	double* input;

	double* output;

	char header[MAX];	//neural file header, may come in handy

	long tlength=n;	//length of training set

	double temp=1;

	double cooldown=0.9995;

	cerr<<"\n**allocating "<<n<<" by "<<tlength<<" array for tset"<<endl;
	tset=new array2d(n,tlength);

	wts=new array2d(n,m);

	output=new double[n];//the output vector (temporary variables)

	input=new double[n];//the vector to be used as input



	//static double error=10;	//initially set to a high value

	char ch;

/* begin initial file I/O operations */

	ok=load(wts->base, n, m, netfilename, header);

	//delete[] header;


	if(!ok){//there's no file there, we have to create our own net

		cerr<<"**zeroing out weights"<<endl;

		wts->bzero();

		/* disabled because this messes up sine positive 1/13/1999 2:43 am */

		//error=0;//if error is 0, train() automatically determines it for you

		

		

		//wts->randomize();

	}

	

	ok=save(wts->base, n, m, s, netfilename,NULL);

	//ok=load(wts->base, n, m, netfilename, header);	//redundant, but will make it work



	if(!ok){

		cerr<<"\ncouldn't save: "<<netfilename<<endl;

		goto end;

	}

	//cerr<<"\nsaved the net we just loaded!"<<endl;

	//cin>>ch;

	//load training set

	trainfile.open(trainfilename);

	if(!trainfile){

		cerr<<"\ncouldn't open training file"<<endl;

		ok=0;

		goto end;

	}

	for(i=0;i<tlength;i++){	//length of trainset.  Should actually be 'n'

		if(trainfile.eof()){

			ok=0;

			cerr<<"\nerror reading from training file!"<<endl;

			goto end;

		}

		for(j=0;j<n;j++){

			trainfile>>tset->base[j][i];	//yes, stuff looks reversed, but this should be OK 1/10/1999

			

		}

	}

/* end of initial file I/O operations */

	//remember that tlength=n for now (1/10/1999)

	for(improvements=1;improvements<=goal_improvements;improvements++){

		//printf("\nimp: %d n: %d m: %d s: %f err: %f tlength: %d temp: %f cooldown: %f",	improvements, n, m, s, error, tlength, temp, cooldown);

		Gimprovements=improvements;

		//train handles lyapunov constraint checking and LE calculation through globals

		error=train(tset->base, wts->base, n, m, s, error, tlength, temp, cooldown);

		Gerror=error;

		Gtemperature=temp;

		report();	//output stuff to screen

		if(improvements%run_freq==0){

			ok=save(wts->base, n, m, s, netfilename,NULL);

			if(run_len>0){

				/* begin run code */

				outfile.open(outfilename);

				outfile<<"% improvements: "<<Gimprovements<<

					" steps: "<<Gsteps<<" error: "<<Gerror<<" LE: "<<Gle<<" constraint: "<<Gconstraint<<endl;

					for(j=0;j<n;j++){

					input[j]=tset->base[j][0];	//get the first value in each column

				}

				

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

					niterate(output, input, wts->base, n, m, s);

					for(j=0;j<n;j++){

						outfile<<output[j]<<" ";

						//output[j]=input[j];

						input[j]=output[j];

					}

					outfile<<endl;

				}

				outfile.close();

			}//end inner if

			/* end run code */

		}//end outer if

	}//end for



			



end:

		return ok;

	



}



//utility functions



//generates sine training set file

void gensineset(char* filename, int n, int tlength){

	long i,j,k;

	array2d* tset;

	ofstream outfile;

	tset=new array2d(n,tlength);

	outfile.open(filename);

	//outfile<<endl<<"test"<<endl;

	//outfile.close();

	if(!outfile){
		cerr<<"\ncan't open file for writing!"<<endl;
		exit(1);
		

	}
			

	

	cerr<<"\nmaking assignments!"<<endl;

	for(i=0;i<n;i++){	//for each neuron in the row

		for(j=0;j<m;j++){	//do the timecourse

				tset->base[i][j]=0.5*sin(((2*PI*(j+i))/(float)n));

		}

		

	}



	cerr<<"\nwriting to outfile!"<<endl;

	//outfile<<endl<<"test"<<endl;	

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

		for(j=0;j<n;j++){

			//cout<<"i j tset->base[i][j]"<<endl;

			//cout<<i<<" "<<j<<" "<<tset->base[i][j]<<endl;

			outfile<<(tset->base[i][j])<<" ";

		}

		outfile<<endl;

	}

	outfile.close();



	delete tset;

}