//neural.cpp


//Oguz Yetkin, Physics 499


//2 23 1998


//neural network implementation


//


//for now, I just want to write something that works.  We shall


//be expanding this as we go along


//


//array2d class works properly in allocating heightxwidth array


//as of 3/1/1998, 12:16 am


//


//purpose: generate http://sprott.physics.wisc.edu/neural/bench.dat





#include<math.h>


#include <iostream.h>


#include <stdio.h>


#include <stdlib.h>	//for rand(int)





#include <assert.h>





#include <limits.h>	//for RAND_MAX
#include<time.h>	//for srand 4/21/1998

#include "neural.h"





long double  squash(long double  value){





	return tanh(value);


}








//array2d implementation





array2d::array2d(int h, int w){







	if(debug){
		cout<<"\nCONSTRUCTING ARRAY2d, memory being allocated!!!"<<endl;
	}








	width=w;


	height=h;


	error=0;


	//base=new long double *[width];	//do these need to be reversed ???


	base=new long double * [height];	


	for(int i=0;i<height;i++){


		base[i]=new long double [width];


	}


}








array2d::~array2d(){




/* commented out 4/24/1998 --deletes are messing something up
	for(int i=0;i<height;i++){


                delete[]base[i];


        }	

  */

	delete[] base;


}





long double  array2d::getvalue(int i, int j){








	if(i>=height||j>=width){


		error=1;


		cerr<<"\ngetvalue: illegal subscripts!"<<i<<" "<<j<<endl;


		cerr<<"\nthe array size is: "<<height<<" "<<width<<endl;


		return 1;


	}


	return base[i][j];





}





int array2d::putvalue(int i, int j, long double  value){





 	if(i>=height||j>=width){


                error=1;


		cerr<<"\nputvalue: illegal subscripts!"<<i<<" "<<j<<endl; 


		cerr<<"\nthe array is: "<<height<<" "<<width<<endl;


		return -1;


        }


	


	base[i][j]=value;


	return 1;		


}





void array2d::dump(){


	


	//print the array


	printf("\nthe array is %d by %d\n",height, width);	


	for(int i=0;i<height;i++){


		for(int j=0;j<width;j++){


			


			//printf("%f ",base[i][j]);


			printf("%.15f ",getvalue(i,j));





		}


		printf("%c",'\n');


	}


}





void array2d::dump(int beg, int end){	//restricts width





  	for(int i=0;i<height;i++){


                for(int j=beg;j<end;j++){





                        //printf("%f ",base[i][j]);


                        printf("%f\t",getvalue(i,j));





                }


                printf("%c",'\n');


        }


		








}


void array2d::randomize(){





	for(int i=0;i<height;i++){


                for(int j=0;j<width;j++){





                        //base[i][j]=rand(0);





				//putvalue(i,j,rand(0));


				putvalue(i,j,rand());	//for windows


                }


       		//cout<<endl; 


	}	








}








Network::Network(int newd, int newn, long double  svalue){


	


	d=newd;


	n=newn;


	s=svalue;


	generation=0;




	if(debug){
		cout<<"\nCONSTRUCTING network, a HECK of a lot of memory allocated"<<endl;
	}


	w=new array2d(n+1,d+1);


	b=new array2d(n+1,d+1);


	input=new  array2d(1,d+1);


	output=new array2d(1,d+1);





	w->randomize();


	b->randomize();





	input->randomize();


	output->randomize();





}//end Network::Network





Network::~Network(){


	


	//deallocate memory


	delete w;


	delete b;


	delete input;


	delete output;





}//end Network::~Network





int Network::propagate(){	//returns the current generation





	long double outerpart=0;	


	long double innerpart=0;





	long double u;





	//could this be our memory leak???





	long double* phiu=new long double[n+1];	//from Prof. Sprott's code


    int i,j,k;











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





		u=w->base[i][0];





		for(j=1;j<=d;j++){





			u=u+(w->base[i][j])*(input->base[0][j]);





		}





		phiu[i]=squash(u);





	}





	for(j=1;j<=d;j++){





		output->base[0][j]=0;





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





			output->base[0][j]=output->base[0][j]+





				(b->base[i][j])*phiu[i];





		}





	}











	











    /* commented out 3/17/1998 by Oguz 











	for(int k=1;k<=d;k++){	


	  for(int i=1;i<=n;i++){


		innerpart=0;	


		 for(int j=1;j<d;j++){


			innerpart+= 


			  (w->base[i][j])*input->base[0][k];


		 }


		 //innerpart*=s; //commented out 3/17/1998 by Oguz


		 innerpart=innerpart+s*(w->base[i][j]);





		 innerpart=squash(innerpart);


		 outerpart+=((b->base[i][k])*innerpart);


	  }	


		


	  output->base[0][k]=outerpart;


	}











  */





	//unnecessary, here for extensibility purposes


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


	    input->base[0][i]=output->base[0][i];


	}








	delete[]phiu;	//hope this still works, this might fix the memory leak


	return generation++;





}








//overloaded assignment operator for array2d





//assumes left hand side is properly allocated and





//IS OF THE SAME SIZE





array2d &array2d::operator=(array2d &in){





	





	//cout<<"\noperator= for array2d got called!!!"<<endl;





	long i,j,k;





	i=j=k=0;





	for(i=0;i<in.height;i++){





		for(j=0;j<in.width;j++){





			base[i][j]=in.base[i][j];





			/* works, as of 3/31/1998, 3:13 am





			printf("\ni: %d j: %d base[i][j]: %f in.base[i][j]: %f ",





				i, j, base[i][j], in.base[i][j]);





			*/





		}





	}





	return *this;





}





	











//overloaded assignment operator for Network





//assumes networks are of the same size!!!





Network &Network::operator=(Network &in){





	//cout<<"\noperator= for Network got called!!!"<<endl;





	int i,j,k;





	i=j=k=0;





	generation=in.generation;





	//must dereference, because operator= will work on the





	//real object, but not the pointer.  WE want to copy, 





	//we DO NOT want a straight pointer assignment!!!











	d= in.d;	





	n= in.n;





	*w= *(in.w);	//relying on array2d's operator=





	*b= *(in.b);





	*input=*(in.input);





	*output=*(in.output);





	s=in.s;











	//rest to be completed





	//writing operator= for array2d now...























	return *this;











}













void Neuron::randomize(){
	
	input_weights->randomize();

}




void Neuron::dump(){











	cout<<"\ndumping neuron: "<<selfindex<<endl;





	cout<<"input weights"<<endl;





	input_weights->dump();





	cout<<"\nincoming neighbors: "<<endl;





	inputs->dump();





	cout<<"\noutput values"<<endl;





	//output_values->dump();





	//outvalues->dump();





	dumpoutput();





}











void Neuron::dumpoutput(){





	//cout<<endl;


	//cout<<"\noutput value for neuron: "<<selfindex<<endl;

	//output_values->dump();

	//outvalues->dump();

	//cout<<output<<endl;

	cout<<output;





}











void Neuron::commit(){











	output=temp_output;	//commit the changes made





}


List::List(){











	//create an empty list





	num_elements=0;





	L=NULL;





	current_pos=0;





}











List::List(int numelements){











	num_elements=numelements;





	L=new long double [num_elements];





	assert(L);





	current_pos=0;





	for(long i=0;i<num_elements;i++){	//zero everything out





		L[i]=0;





	}





}











List::~List(){





	delete[]L;





}

















int List::push(long double  num){





	





	if(current_pos>=num_elements){





		cout<<"\npush failed!"<<endl;





		return 0;





	}





	L[current_pos]=num;





	current_pos++;





	return 1;





}











int List::pop(long double & num){











	if(current_pos==0){





		cout<<"\ncan't pop past beginning of list!!!"<<endl;





		return 0;





	}





	





	current_pos-=1;





	num=L[current_pos];





	





	return 1;

















}











List &List::operator=(List &in){	//we also do inialization if necessary





	if(L){





		delete[]L;	//if there was anything, delete it





	}else{	//the input list was NULL





		cout<<"\nINPUT LIST WAS NULL, SOMETHING IS WRONG"<<endl;





		return *this;





	}





	num_elements=in.num_elements;





	current_pos=in.current_pos;





	L=new long double [num_elements];





	for(long i=0;i<num_elements;i++){





		L[i]=in.L[i];





	}



	return *this;

}




void List::dump(){




	long double value;


	int ok;


	if(num_elements==0){

		cout<<"\nlist is empty, can't dump it!"<<endl;


	}else{


		//cout<<"displaying list "<<endl;

		assert(L);

		for(long i=0;i<num_elements;i++){



			printf("%f ",L[i]);





	





		}





	}











	





}//end List::dump()



void List::randomize(){


	for(long i=0;i<num_elements;i++){

		L[i]=myrand();

	}


}








//useless default constructor





Neuron::Neuron(){




	//cout<<"\nthe DEFAULT CONSTRUCTOR FOR NEURON SHOULD NEVER BE CALLED! 4/24/1998 "<<endl;
	selfindex=-1;	//means invalid





	input_weights=NULL;




	//there is NO output weights now...4/17/1998
	//...but we still need the variable to compile.
	//WARNING!!!  Neuron::operator=() does NOT copy this value.
	//for now (around 4/17/1998) I'm (Oguz Yetkin) assuming that
	//output_weights will always stay NULL

	output_weights=NULL;




	outputs=NULL;





	//the value is arbitrary


	//will make this random!!!

	//output=0.1;	//changed to single output 4/3/1998
	//made random 4/17/1998
	output=squash(myrand());




}











Neuron::~Neuron(){

	cout<<"\ncomment the delete back in! 4/24/1998"<<endl;
	//delete [] input_weights;
	//delete[] output_weights;
	//delete [] inputs;
	//delete [] outputs;
	s=0;
	selfindex=0;




	//cout<<"\nDESTRUCTOR NOT YET IMPLEMENTED FOR NEURON!!!"<<endl;



}





//not really getting used now





Neuron::Neuron(long self, List* inweights, List* outweigths, List* newoutputs){ 







	selfindex=self;





	*input_weights=*inweights;	//using list's operator=(), it handles allocation





	*output_weights=*outweigths;





	*outputs=*newoutputs;





}











//simpler constructor--we'll be using this one





//THIS IS THE ONE WE USE for now





Neuron::Neuron(double in_s,long self, List* inweights, List* newinputs){











	//this isn't really getting called right now--we are doing the allocation manually





	selfindex=self;





	//*output_values=*inweights;





	//*outvalues=*inweights;





	*input_weights=*inweights;





	*inputs=*newinputs;





	//output_values=new List(input_weights->num_elements);





	//outvalues=new List(input_weights->num_elements);	





	//for now, we'll set output values equal to input weights, just to get it to compile





	//& not seg fault











	//just as a sanity check





	//output_values=new List(4);





	//assert(output_values);





	//assert(outvalues);





	s=in_s;





}





























//written 4/1/1998



void Neuron::propagate(Neuron* world){  //we need not pass anything in since the neuron knows who its inputs are!!!





	//changed 4/8/1998 to divide weights by sqrt(m)


	static int firstrun=1;


	static long double m;


	static long double sqrt_m;	//we don't want to do this more than once, or take a square root inside a loop





	if(firstrun){


		firstrun=0;


		m=inputs->num_elements;


		sqrt_m=sqrt(m);


	}





	//changed 4/3/1998 to take into account that we only
	//get one output from the incoming neighbor
	//write propagate function--imoprtant!!!


	long i,j,k;

	i=j=k=0;

	long double value; //generic variable
	long double neighindex;
	long double weight;
	long double in_value;	//incoming value from neuron
	long in_neighbor;
	int ok=1;//is everything OK?


	inputs->reset();





	





	input_weights->reset();	//get ready to pop





	





	//don't care about outvalues anymore 4/3/1998











	temp_output=0;





	value=0;





	//outvalues->current_pos=0;





	





	//we should have as many input weights as we have





	//incoming neighbors





	assert(input_weights->num_elements==inputs->num_elements);





	//for every incoming neighbor











	for(i=0;i<inputs->num_elements;i++){











		ok=inputs->pop(neighindex);  //get the index of a neighbor





		assert(ok);





		in_neighbor=(long)neighindex;





		in_value=world[in_neighbor].output;





		assert(in_neighbor>=0);





		ok=input_weights->pop(weight);	//get current input weight





		assert(ok);





		//value is an intermediate value





		value+=s*in_value*weight/sqrt_m;





	}





	temp_output=squash(value); //real output is assigned by commit() which is called by rNetwork::propagate() 4/3/1998











	/* commented out 4/3/1998











	for(j=0;j<world[selfindex].outvalues->num_elements;j++){ //for every component of output vector





			





		





		inputs->reset();			//reset the neighbor list





		value=0;





		for(i=0;i<inputs->num_elements;i++){	//for every incoming neighbor





				ok=inputs->pop(neighindex);  //get the index of a neighbor





				assert(ok);





		





				in_neighbor=(long)neighindex;





				assert(in_neighbor>=0);





	





				world[in_neighbor].outvalues->reset();	//so we can pop





		





				//value=0;





				//for every output value in the incoming neighbor's output vector





				input_weights->reset();





				for(k=0;k<world[in_neighbor].outvalues->num_elements;k++){





					





					//calculate one output value





					//get current weight





					ok=input_weights->pop(weight);





					assert(ok);	//for now, we should not run into problems--if we do, fix it!!!





					//ok=world[in_neighbor].outvalues->pop(in_value);





					in_value=world[in_neighbor].outvalues->L[j]; //assuming all output vectors the same size!!!





					assert(ok);





					value+=s*weight*in_value;





				}





			





		





		}//end middle for that has the purpose of handling output from a single incoming neighbor





			value=squash(value);	//apply squashing function





		//outvalues->current_pos=0;





		





			//ok=outvalues->push(value);





			//assert(ok);	//one of these is going to fail...I know it...





			outvalues->L[j]=value;





	}//end outer for, all neighbors done at this point





	*/











			





}
















//this is where we should get rid of dynamic allocation.  Causes bugs 
//when we keep track of all the outputs
//
//last update: 4/17/1998

Neuron & Neuron::operator=(Neuron& in){


	selfindex=in.selfindex;



	*input_weights=*(in.input_weights);

	//we'd forgotten to do this, added 4/17/1998, 3:09 am
	*inputs=*(in.inputs);



	//we no longer use output_weights


	//4/3/1998
	//4/17/1998, assuming following values are NULL

	//*output_weights=*(in.output_weights);

	//*outputs=*(in.outputs);





	//modified 4/33/1998


	output=in.output;	//is this working properly
	temp_output=in.temp_output;

	//this should NEVER fail!
	if(output!=in.output){
		cout<<"\n!!!!A VERY STRAIGHTFORWARD ASSIGNMENT JUST FAILED in neural.cpp Neuron::operator=().  GO BACK TO YOUR DEBUGGER!!!\nhalting program..."<<endl;

	}
	assert(output==in.output);





	return * this;

}

















//maybe at some point we want to have a FLAG for inputs that





//specifies random number of connections between neurons.





//just an idea 4/3/1998











rNetwork::rNetwork(long size, long double scale, long inputs){


	long i,j,k;
	i=j=k=0;


	generation=0;
	net_size=size;

	net=new Neuron[size]; //this _is_ calling the default constructor, unlike previously planned, so we deal with it.  Noticed : 4/24/1998


	//temp_net=new Neuron[size];





	s=scale;





	int ok;





	input_size=inputs;


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

		net[i].input_weights=getinweights(input_size);	//these are just ptr assignments, so we save time and memory

		//*temp_net[i].input_weights=*net[i].input_weights;
		net[i].output_weights=getoutweights();
		//*temp_net[i].output_weights=*net[i].output_weights;

		net[i].inputs=get_neigh_list();  //can be random or ordered

		//*temp_net[i].inputs=*net[i].inputs;
		net[i].outputs=new List(net[i].inputs->num_elements); //not used now

		//just for consistency's sake





		//*temp_net[i].outputs=*net[i].outputs;





		//net[i].outvalues=new List(net[i].input_weights->num_elements);





		//net[i].outvalues=get_output_seed();





		//*temp_net[i].outvalues=*net[i].outvalues;





	











		net[i].s=getscale();





		//temp_net[i].s=net[i].s;





		net[i].selfindex=i; //so we know who we are

		net[i].output=squash(myrand());



		//temp_net[i].selfindex=net[i].selfindex;











	}





	/*





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





		temp_net[i]=net[i];





	}





	*/

















}





rNetwork::~rNetwork(){

	delete[] net;
	//delete[] temp_net;



}

//added 5/13/1998
void rNetwork::randomize(){

	for(long i=0;i<net_size;i++){

		net[i].randomize();

	}

}

	





List* rNetwork::getinweights(long input_size){





	//static int firstrun=1;





	//static List* inlist;		//put anything there for now





	List* inlist;





	int ok;





/* commented out 4/3/1998





	if(firstrun){





		inlist=new List(input_size);	//ALL this can be VERY easily randomized





		/*





		inlist->push(rand()/RAND_MAX);





		inlist->push(rand()/RAND_MAX);





		inlist->push(rand()/RAND_MAX);





		inlist->push(rand()/RAND_MAX);





		//





		inlist->push(0.2321);





		inlist->push(0.6);





		inlist->push(0.9);





		inlist->push(0.1);





		firstrun=0;





	}





*/





	inlist=new List(input_size);





	assert(inlist);





	for(long i=0;i<input_size;i++){





		ok=inlist->push(myrand());





		//comment out later for optimization





		if(!ok){





			cout<<"\ngetinweights is pushing past end of list at "<<i<<endl;





		}





	}





	





	assert(inlist);





	return inlist;

}





List* rNetwork::getoutweights(){











		static int firstrun=1;





	static List* outlist;		//put anything there for now





	if(firstrun){





		outlist=new List(4);	//ALL this can be VERY easily randomized





		outlist->push(0.2);





		outlist->push(0.8);





		outlist->push(0.1);





		outlist->push(0.3);





		firstrun=0;





	}





		





	





	





	return outlist;











}











long double rNetwork::getscale(){











	





	return s;	//very arbitrary...

















}





List* rNetwork::get_neigh_list(){  //this list should actually contain ints











	List* neighlist;





	int max_neighbors=input_size;





	//int num_neighbors=rand()%max_neighbors +1;





	long num_neighbors;





	num_neighbors=max_neighbors;	//for now, for benchmarking purposes 4/1/1998





	int ok=1;




	if(debug){
		cout<<"\ndon't forget to delete list in real implementation!!!"<<endl;
	}





	neighlist=new List(num_neighbors);	//this will have to be deleted somehow





	assert(neighlist);





	for(long i=0;i<num_neighbors;i++){





		ok=neighlist->push( abs(rand()%net_size) );





		assert(ok);





	}





	return neighlist;





}











List* rNetwork::get_output_seed(){





	





	List* seed;




	if(debug){
		cout<<"\ndon't forget to delete list later!"<<endl;
	}





	//ASSUMING all output vectors the same size





	seed=new List(4);		//this will later be adjustable





	seed->push(0.2112);





	seed->push(0.3);





	seed->push(0.8);





	seed->push(0.5);





	return seed;

















}











int rNetwork::propagate(){











	long i,j,k;





	i=j=k=0;





	List* inputs;





	List* outputs;











	//initially, net and temp_net are the same network





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



			//inputs=net[i-1].outputs;

			//outputs=net[i].propagate(inputs); //outgoing neighbors of current neuron get returned

			//inputs=outputs;	//this should just be a pointer assignment

			net[i].propagate(net); //each neuron needs a pointer to the world

							   //a more C++ thing to do would be to pass in a pointer


							   // to 'this', but I don't like "this"



	} //end for





		





	//now all the outputs of the neurons in "net" have been updated, 





	//and temp_net can reflect the updated weights





	





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





		net[i].commit();	//swap temp_output with output. Thisis 





							//to ensure that we are advancing the whole network at once.





							//ASSUMPTION: weights do not change





	}





	

















	generation++;





	return 1;





}











void rNetwork::dump(){	//dump every single neuron in the network


	for(long i=0;i<net_size;i++){


		net[i].dump();


	}


}











void rNetwork::dumpoutputs(){








	cout<<endl;


	for(long i=0;i<net_size;i++){





		//net[i].dumpoutputs();		//I LOVE object oriented programming





		net[i].dumpoutput();





		cout<<" ";





	}


	//cout<<endl;


	





}











rNetwork& rNetwork::operator=(rNetwork &in){











	s=in.s;





	generation=in.generation;





	net_size=in.net_size;





	input_size=in.input_size;





	//modified 4/8/1998, becase we want the pointers to stay the same!


	//delete[] net;





	if(!net){	//this better not seg fault on me


		net=new Neuron[net_size];


	}





	for(long i=0;i<net_size;i++){





		net[i]=in.net[i];





	}











	return *this;





}




long double myrand(){



	static int firstrun=1;

	/* Seed the random-number generator with current time so that
    * the numbers will be different every time we run.
    */
   //cout<<"\ntime: "<<(unsigned)time(NULL);
	if(firstrun){
		srand( (unsigned)time( NULL ) );	//we only do this once in a given program
		firstrun=0;
	}





	//create Gaussian distribution,





	//add 12 uniform #s, subtract 6





	long n=12;





	long double thisrand;





	long double sum;





	long double returnval;





	





	sum=0;





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





		thisrand=(long double)rand()/(long double)RAND_MAX;





		//thisrand-=0.5;	//allow it to be negative





		sum+=thisrand;





	}





	//returnval=(sum-(n/2))/(long double)(n/2);


	//no need to divide


	returnval=(sum-(n/2));	//Standard Deviation is 1 when n==12





	return returnval;





}