//main.cpp
//neural network implementation, 2/23/1998
//Oguz Yetkin
//
//working so far: benchmark for feedforward network
//				  Lyapunov exponent for feedforward network
//				  recurrent network implementation (last correction 4/8/1998, to correct weights with sqrt(num_inputs))
//				  TO DO: is Lyapunov working correctly?  Converges @ incorrect value.  Why?  Question outstanding as of 4/8/1998
//
//last modified: 4/8/1998 to add recurrent Lyapunov calculatoin
//

//BUG DISCOVERED: 4/17/1998 1:06 AM.

//					setting d0=0 on recurrent network calculation still yields

//					a nonzero d1.  Suspected: rNet::operator=() is messing up

//





#include<iostream.h>
#include<fstream.h>
#include<stdio.h>


#include<math.h>	//for square root
#include<stdlib.h>	//for rand()
#include<assert.h>


#include<limits.h>	//let's know what our limits are!!!
#include <string.h>
#include "neural.h"
#include "Henon.h"

#ifdef WIN32
#include "plotter.h"
#endif

const int SUICIDE_VALUE=2;	//exit program after this many iterations in "gather data" mode.  4/29/1998.
								//this is so that NT can do our garbage cleaning, and the program can be respawned

extern int debug;




union bytetool{	//to avoid byte order confusion
		unsigned char bytes[4];
		float number;
} bt;







const long double DELTA=1E-15;	//very small value




//handles all lyapunov calculations for given network





long double Lyapunov(Network* net, long iter, long displayfreq);





//written 4/8/1998
//slightly modified version of Lyapunov routine, handles
//recurrent networks





long double Lyapunov_r(rNetwork* net, long iter=500, long displayfreq=20);

void Lyapunov_r_driver(); //does all its own prompting
void Lyapunov_r_benchmarkdriver();
int SaveSettings(char* filename, double min_s, double max_s,double  s_increment,double  min_size, double max_size, double size_increment,
				 long min_inputs_per_neuron, long max_inputs_per_neuron, long inputs_per_neuron_increment, long iter,
				 int screenoutput, int ascii);


int LoadSettings(char* filename, double& min_s, double& max_s,double&  s_increment,double&  min_size, double& max_size, double& size_increment,
				 long& min_inputs_per_neuron, long& max_inputs_per_neuron, long& inputs_per_neuron_increment, long& iter,
				 int& screenoutput, int& ascii);


				 
				 
		

//based very closely on Lyapunov_r, but uses 2 decoupled Henon maps, written
//to test Lyapunov calculator.  Written 4/8/1998, just like everything else that matters


long double Lyapunov_henon(long iter, long displayfreq, long map_size);




//given 2 vectors and their dimensions, calculates distance in N-space





long double distanceN(long double* vector1,long double* vector2, long n);


//same thing, but takes a pointer to an array of ptrs to long double
//written 4/8/1998


long double distanceNptr(long double** vector1, long double** vector2, long n);



long double sqr(long double num);	//square a num


long double log2(long double num);


void recurrent_benchmark();


enum choices{BENCHMARK=1,LYAPUNOV,RECURRENT_BENCHMARK, RECURRENT_LYAPUNOV, HENON_LYAPUNOV, GATHER_DATA, PLOT, TRAIN};





int main(int argc, char** argv){


	int choice;
	long i,j;
	long iter;
	long map_size;



	long displayfreq;
	long double lyapunov;








	long double s;	//this block of variables required for recurrent Lyapunov calculation 4/8/1998
	long size;
	long inputs_per_neuron;



	//variables for Lyapunov Exponent calculation


	Network * net;
	rNetwork * rnet;
	Network * net2;


	long double d0=0.00000001;

	
	if(argc>1){
		if(strcmp(argv[1],"debug")==0){
				cout<<"\nentering debug mode";
				debug=1;
		}else if(strcmp(argv[1],"benchmark")==0){
				//cout<<"\nbenchmarking!"<<endl;
				debug=0;
				Lyapunov_r_driver();
				exit(1);
		
		}else if(strcmp(argv[1], "data")==0){
			cout<<"\ngathering data in batch mode"<<endl;
			Lyapunov_r_benchmarkdriver();
			return 1;	//return an error code from main, so that they can respawn us
		} else if(strcmp(argv[1], "lyapunov")==0){
			if(argc!=5){
				cout<<"\nuse: main lyapunov net_size inputs_per_neuron scale_factor"<<endl;
				exit(1);
			}

			long n;	//# of neurons
			long m;	//inputs /neuron 
			double s;	//scale factor
			double lyapunov;
			n=strtol(argv[2],NULL,10);
			m=strtol(argv[3],NULL, 10);
			s=strtod(argv[4],NULL);
			cerr<<"\nn m s: "<<n<<" "<<m<<" "<<s<<endl;
			rnet=new rNetwork(n,s,m);
			lyapunov=Lyapunov_r(rnet);
			cout<<n<<"  "<<m<<"  "<<s<<"  "<<lyapunov<<endl;
			exit(1);
		}else if(strcmp(argv[1], "tweak")==0){
			if(argc!=6){
				cout<<"\nuse: main tweak net_size inputs_per_neuron scale_factor iter"<<endl;
				exit(1);
			}

			long n;	//# of neurons
			long m;	//inputs /neuron 
			double s;	//scale factor
			double lyapunov;
			long tweakiter;
			n=strtol(argv[2],NULL,10);
			m=strtol(argv[3],NULL, 10);
			s=strtod(argv[4],NULL);
			tweakiter=strtol(argv[5], NULL, 10);
			cerr<<"\nn m s: "<<n<<" "<<m<<" "<<s<<endl;
			rnet=new rNetwork(n,s,m);
			lyapunov=Lyapunov_r(rnet);
			cout<<n<<"  "<<m<<"  "<<s<<"  "<<lyapunov<<endl;
			for(int i=0;i<tweakiter;i++){
				rnet->randomize();
				lyapunov=Lyapunov_r(rnet);
				cout<<lyapunov<<endl;
			}
			exit(1);
		}


	
		else{
			printf("\nuse: %s [debug] [benchmark] [data] [lyapunov n m s] [tweak n m s iter]",argv[0]);
		
		
		}
		
		
	}else{
		debug=0;
			printf("\ntype %s [debug] [benchmark] [data] [lyapunov n m s] [tweak n m s iter]\n",argv[0]);
	
	}


	


		

	

	/*
	Henon* testhenon;
	testhenon=new Henon();
	testhenon->randomize();
	testhenon->discard();
	cout<<"\ntesting henon!"<<endl;


    for(i=0;i<100;i++){
		testhenon->advance();
		printf("\nx: %f y: %f ",testhenon->x, testhenon->y);
	}
	*/




	if(debug){
		cout<<"\nmaking new network!"<<endl;
	}
	net=new Network(4, 4, 4); //n=4 d=4 s=4
	net2=new Network(4,4,4);	//parallel universe?


	//long iter;


	printf("\nWARNING: memory deallocation temporarily disabled, running this program for an extended period of time may crash your machine");
	printf("\n4/24/1998\n");


	printf("\n%d)Benchmark \n%d)Lyapunov \n%d)Recurrent Benchmark \n%d)Recurrent Lyapunov\n%d)Henon Lyapunov\n%d)Gather Data\n",


		BENCHMARK, LYAPUNOV, RECURRENT_BENCHMARK, RECURRENT_LYAPUNOV,HENON_LYAPUNOV,GATHER_DATA);



	cin>>choice;



	//I'm somewhat confused about benchmark.dat
	//if there's no 0 for b, then are all values 
	//initialized to 0, or does the array start from 1???
	//I'm putting everything @ 0 to make things work
	//OY 3/2/1998 5:47 pm





	net->w->base[0][0]=0;
        net->b->base[0][0]=0;
        net->w->base[0][1]=0;
        net->b->base[0][1]=0;
        net->w->base[0][2]=0;
        net->b->base[0][2]=0;
        net->w->base[0][3]=0;
        net->b->base[0][3]=0;





	net->w->base[1][0]=-1.161204;
	net->b->base[1][0]=0;
	net->w->base[1][1]=-.561551;
	net->b->base[1][1]=.1636758;
	net->w->base[1][2]=5.066455;
	net->b->base[1][2]=.2757667;
        net->w->base[1][3]=7.139109; 
	net->b->base[1][3]=9.563043E-2;
 	net->w->base[1][4]=-5.131825; 
        net->b->base[1][4]=8.334725E-2;
        net->w->base[2][0]=9.370738E-2; 
        net->b->base[2][0]=0;
        net->w->base[2][1]=-5.116127;
        net->b->base[2][1]=2.233952E-2;
        net->w->base[2][2]=4.733435; 
        net->b->base[2][2]=.3650303;
	net->w->base[2][3]=-8.787042 ; 
        net->b->base[2][3]=.355971;
        net->w->base[2][4]= -.8605409; 
        net->b->base[2][4]=.4066904; 
        net->w->base[3][0]=4.888542; 
        net->b->base[3][0]=0; 
        net->w->base[3][1]=-4.576682; 
        net->b->base[3][1]=.5034996; 
 	net->w->base[3][2]=7.925086  ;
        net->b->base[3][2]=.1963863 ;
        net->w->base[3][3]=-5.740646 ;
        net->b->base[3][3]=.2648118;
        net->w->base[3][4]=3.809366 ;
        net->b->base[3][4]=.4218021 ;
        net->w->base[4][0]=2.468891;
        net->b->base[4][0]=0;
       
        net->w->base[4][1]=2.116171 ;
        net->b->base[4][1]= .3104851;
        net->w->base[4][2]=-3.889002E-2;
        net->b->base[4][2]=.1628166;
        net->w->base[4][3]=-2.105793E-2;
        net->b->base[4][3]=.2835867 ;
        net->w->base[4][4]= .5917256 ;
        net->b->base[4][4]= 8.816015E-2;


       





	//now attempting to verify benchmark.dat

		if(debug){
			printf("\nattempting to initialize to weights at http://sprott.physics.wisc.edu/neural/bench.dat\n");
		}


	/* the basic code is:
	  for i = 1 to n
       		 for j = 0 to d
            		print w(i, j)
            		print b(i, j)
        	next j
    	  next i
	*/


  switch(choice){


   case BENCHMARK:


		cout<<"\nenter times to iterate: ";


	cin>>iter;



	  for(i=1;i<=(net->n);i++){


		for(j=0;j<=(net->d);j++){


			printf("\nnet->w->base[%d][%d]=%f",i,j,


				net->w->base[i][j]);


			//printf("\n%f",net->w->base[i][j]);


			printf("\nnet->b->base[%d][%d]=%f",i,j,


				net->b->base[i][j]);


			//printf("\n%f",net->b->base[i][j]);


		}


	}


	


	





	for(i=0;i<=net->d;i++){
		net->input->base[0][i]=0.5;
		net->output->base[0][i]=0.5;


	}//initialize values





	cout<<"\ninput: "<<endl;




	net->input->dump(1,5);
	cout<<"\noutput: "<<endl;
	net->output->dump(1,5);


	


	cout<<"\nbrace for protection fault, first use of operator= !!!"<<endl;




	*net2= *net;


	net2->input->base[0][1]+=d0;	//the first input will always be 0





	cout<<"\nnow iterating "<<iter<<" times"<<endl;
	printf("\nx(1)\t       x(2)\t       x(3)\t       x(4)\t       \n");


	for(i=0;i<iter;i++){
		net->input->dump(1,5); //restrict width of array dumped



		//net2->input->dump(1,5);


		net->propagate();
		//net2->propagate();


	}






	/*


	cout<<"\nnow doing the same thing for net2"<<endl;
	cout<<"\ninput: "<<endl;
	net2->input->dump(1,5);
	cout<<"\noutput: "<<endl;
	net2->output->dump(1,5);



	cout<<"\nnow iterating "<<iter<<" times"<<endl;


	printf("\nx(1)\t       x(2)\t       x(3)\t       x(4)\t       \n");


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





		net2->input->dump(1,5); //restrict width of array dumped





		net2->propagate();



	

	  }





	*/





   break;





  case LYAPUNOV:



	cout<<"\niterations for lyapunov calculation? ";
	cin>>iter;
	cout<<"\nenter number of iters to wait before displaying a result: "<<endl;
	cin>>displayfreq;
	lyapunov=Lyapunov(net, iter,displayfreq);
	cout<<"\nThe Exponent after "<<iter<<" iterations is: "<<lyapunov<<endl;



  break;











  case RECURRENT_BENCHMARK:





	  cout<<"\ntesting recurrent randomly connected network!!!"<<endl;
	  recurrent_benchmark();


  break;











  case RECURRENT_LYAPUNOV:





	  //cout<<"\nlyap calculation for recurrent randomly connected network not yet implemented"<<endl;





	cout<<"\nLyapunov exponent for recurrent network, first implemented 4/8/1998"<<endl;
	cout<<"\nscale factor? "<<endl;
	cin>>s;
	cout<<"\nnumber of neurons? "<<endl;
	cin>>size;
	cout<<"\ninputs per neuron? "<<endl;
	cin>>inputs_per_neuron;
	//cout<<"\niterations? "<<endl;
	//cin>>iter;




	cout<<"\ncreating network!"<<endl;
	rnet=new rNetwork(size, s,inputs_per_neuron);
  	cout<<"\niterations for recurrent network lyapunov calculation? ";
	cin>>iter;
	cout<<"\nenter number of iters to wait before displaying a result: "<<endl;
	cin>>displayfreq;
	lyapunov=Lyapunov_r(rnet, iter,displayfreq);


	cout<<"\nThe Exponent after "<<iter<<" iterations is: "<<lyapunov<<endl;


	  break;





  case HENON_LYAPUNOV:


	  cout<<"\nbenchmarking Lyapunov calculator with decoupled henon maps"<<endl;
	cout<<"\niterations for henon map  lyapunov calculation? ";
	cin>>iter;
	cout<<"\nenter number of iters to wait before displaying a result: "<<endl;
	cin>>displayfreq;
	cout<<"\nhow many decoupled maps?"<<endl;

	cin>>map_size;




    lyapunov=Lyapunov_henon(iter,displayfreq,map_size);





	cout<<"\nThe Exponent after "<<iter<<" iterations is: "<<lyapunov<<endl;


    //Lyapunov_henon(iter, displayfreq);
	break;

  case GATHER_DATA:
	Lyapunov_r_benchmarkdriver();

	break;

  case PLOT:{

#ifdef WIN32
	rNetwork* rn;
	rn=new rNetwork(64, 64, 2);
	Plotter* pl;
	pl=new Plotter(100, 100);
	for(int i=0;i<64;i++){




	}



#endif


			}
	  break;

   
  case TRAIN:
		cout<<"\ntraining not yet implemented";

	  break;




  }//end switch








	cout<<"\ndestroying net"<<endl;
	delete net;


	cout<<"\ndone"<<endl;
	return 0;


}














long double Lyapunov(Network* net, long iter, long displayfreq){




	long interval=displayfreq;	//reporting interval for lyapunov


	//input and output vectors of networks
	//are array2d's of size 1xd+1




	long n,d;
	long double s; 
	long discard=500;	//put on attractor, should actually be an adjustable parameter
	long double d0, d1;





	d0=0.00000001;


	long i,j,k;	//counters
	long double sum, lyapunov;




	sum=0;
	lyapunov=0;


	i=j=k=0;





	Network* net1;


	Network* net2;


	Network* net1new;


	Network* net2new;





	





	n=net->n;
	d=net->d;
	s=net->s;





	//now allocate new temporary networks





	net1=new Network(n, d, s);
	net1new=new Network(n, d, s);
	net2=new Network(n, d, s);
	net2new=new Network(n, d, s);




	//dereferencing done in order to use copying 
	//(initialization) operator=




	*net1=*net;	//now we are done with the net passed into us



	//iterate net1 until on attractor, then we'll copy into 
	//other temporary networks





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





		net1->propagate();





	}





	//now assign to other networks





	sum=0;
	lyapunov=0;





	*net1new=*net1;
	*net2=*net1;


   	net2->input->base[0][1]+=d0;	//the first input will always be 0





	*net2new=*net2;




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





		//*net1new=*net1;
		//*net2=*net1;
		//net2->input[0]+=d0;


		//*net2new=*net2;



		//now propagate the "new" versions of both networks.
		//this way, we have one "old" version and one "new" version





		//of each






		net1new->propagate();
		net2new->propagate();




		/*




		cout<<"iter: "<<i<<endl;
		cout<<"net1new: "<<endl;
		net1new->input->dump(1,5);
		cout<<"net2new: "<<endl;
	




		net2new->input->dump(1,5);





	*/





		//calculate d1 (the inputs are updated in both networks)





		//we're passing the ADDRESS of net1new->output->base[0][1]


		//which is treated as a 1-D array beginning at base[0][1]
		//messy pointer arithmetic, that's all...





		d1=distanceN(&(net1new->output->base[0][1]), &(net2new->output->base[0][1]), d);





		//cout<<"\nd1: "<<d1<<endl;
		//LEFT OFF HERE 3/31/1998, 3:45 am





		//readjust orbit of net2 in the same dir and distance as net1
		//this SHOULD work in the N dimensional case





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



			net2->input->base[0][j]=


				net1new->input->base[0][j]+


				(d0/d1)*((net2new->input->base[0][j]-net1new->input->base[0][j]));


		}





		//*net1=*net1new; superfluous 3/31/1998 9:02 p




		/* commented out 3/31/1998


		cout<<"net1: "<<endl;
		net1->input->dump(1,5);
		cout<<"net2: "<<endl;
		net2->input->dump(1,5);


		*/





	





		//the following assertis inefficient, but we want it there
		//assertion fails, of course




		/*


		printf("\nreadjusted distance should be: %f, but it is: %f",





					d0, distanceN(&(net1->input->base[0][1]), &(net2->input->base[0][1]), d));





		*/


		/*	I think the numbers are close enough that I don't have to worry about the assertion failing...


		assert(
			distanceN(&(net1->input->base[0][1]), &(net2->input->base[0][1]), d)==d0);


		*/




		sum+= log (d1/d0);	//base e


		lyapunov=sum/i;


		if(i%interval==0){
			//printf("\niter: %d lyapunov: %f",i, lyapunov);

			printf("\niter: %d lyapunov: %f d0: %f d1: %f",i, lyapunov,d0,d1);




		}



		//*net1new=*net1;
		//*net2new=*net2;
		//*net1new=*net1;	//superfluous, 3/31/1998 9:03p


		*net2new=*net2;


	}//end for



	return lyapunov;





}











//written 4/8/1998
//slightly modified version of Lyapunov routine, handles
//recurrent networks





long double Lyapunov_r(rNetwork* net, long iter, long displayfreq){





	long interval=displayfreq;	//reporting interval for lyapunov
	//input and output vectors of networks


	//are array2d's ofsize 1xd+1


	long net_size, num_inputs;
	long double s;	//the scale factor  4/8/1998, but it was there before :)
	//long discard=500;	//put on attractor, should actually be an adjustable parameter
	//modified: 4/21/1998
	long discard=250;
	long double d0, d1;

	long double dreal;	//the actual readjusted distance, for comparison purposes





	long double** net1vector;	//original values
	long double** net2vector;	//values of parallel universe
	long double** net1vectornew;
	long double** net2vectornew;




	d0=0.00000001;

	//sanity check:

	//d0=0;

	if(debug){
		printf("\nd0 should be zero, it is : %E",d0);
	}


	long i,j,k;	//counters
	long double sum, lyapunov;
	sum=0;
	lyapunov=0;








	i=j=k=0;
	rNetwork* net1;
	rNetwork* net2;
	rNetwork* net1new;
	rNetwork* net2new;




	net_size=net->net_size;
	//d=net->d;
	s=net->s;
	num_inputs=net->input_size;





	net1vector=new long double*[net_size]; //pointers to outputs of all the neurons in the network
	net2vector=new long double*[net_size];
	net1vectornew=new long double*[net_size];
	net2vectornew=new long double*[net_size];


	//cout<<"\nsanity check!  do we crash now? "<<endl;
	//delete[] net1vector;
	//cout<<"\nSO DID WE CRASH?"<<endl;




	//now allocate new temporary networks





	net1=new rNetwork(net_size, s, num_inputs);
	net1new=new rNetwork(net_size, s, num_inputs);
	net2=new rNetwork(net_size, s, num_inputs);
	net2new=new rNetwork(net_size, s, num_inputs);


	//dereferencing done in order to use copying 
	//(initialization) operator=



	*net1=*net;	//now we are done with the net passed into us





	//iterate net1 until on attractor, then we'll copy into 
	//other temporary networks





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


		net1->propagate();


	}





	//now assign to other networks



	sum=0;
	lyapunov=0;


	*net1new=*net1;
	*net2=*net1;
	//net2->net[0].output+=d0;	//we're only modifying one of the neurons.  This may have to change in the future 4/8/1998

	//modified 4/15/1998

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

		net2->net[i].output+=d0;  //perturb all values

	}
	*net2new=*net2;





	//now get our input and output vectors


	//added 4/8/1998
	//now these values should stay w/ us forever.
	//otherwise we run into super-messy pointer problems


	


	for(i=0;i<net_size;i++){	//!!! Make sure we're not changing locations of outputs in memory when we reasisgn!


		net1vector[i]=&(net1->net[i].output);
		net2vector[i]=&(net2->net[i].output);
		net1vectornew[i]=&(net1new->net[i].output);
		net2vectornew[i]=&(net2new->net[i].output);
	}


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





		//now propagate the "new" versions of both networks.
		//this way, we have one "old" version and one "new" version
		//of each


		net1new->propagate();
		net2new->propagate();




		//calculate d1 (the inputs are updated in both networks)
		//we're passing in the 1-d array of long double* which points
		//to every single output in a given network


		d1=distanceNptr(net1vectornew, net2vectornew, net_size);	//net_size is also vector size

		
		


		//readjust distance?? 4/8/1998
		for(j=0;j<net_size;j++){


			//net2->input->base[0][j]=


			*(net2vector[j])=				//once again, messy ptr arithmetic, 4/8/1998





				//net1new->input->base[0][j]+


				*(net1vectornew[j])+





				//(d0/d1)*((net2new->input->base[0][j]-net1new->input->base[0][j]));


				(d0/d1)*


					(*(net2vectornew[j])-*(net1vectornew[j]));


		}



		//calculate readjusted distance 4/17/1998

		dreal=distanceNptr(net1vectornew, net2vector, net_size);	//net_size is also vector size

		


		sum+= log (d1/d0);	//base e
		lyapunov=sum/i;



		if(i%interval+1==0){


			printf("\niter: %d lyapunov: %E d0: %E d1: %E dreal: %E",i, lyapunov,d0,d1,dreal);


		}

		//cout<<"\ndreal - d0="<<dreal-d0<<endl;

		//printf("\n||dreal-d0||= %E",(float)dreal-(float)d0);
		if(debug){
			printf("\n||dreal-d0||= %E",sqrt(sqr((float)dreal)-sqr((float)d0)));
		}
		//assert(sqrt(sqr(dreal)-sqr(d0))<DELTA);  //the distance is usually adjusted to within 1e-17
		//assert((float)dreal==(float)d0); //float precision is good enough...

		*net2new=*net2;


	}//end for




	//delete stuff that we allocated
	//we crash here, but we do NOT crash if we delete before we do a 
	//whole bunch of stuff

		//delete[] net1vector;
		//delete[] net2vector;
		//delete[] net1vectornew;
		//delete[] net2vectornew;


		//delete[] net1;
		//delete[] net1new;
		//delete[] net2;
		//delete[] net2new;


	return lyapunov;


}











//based very closely on Lyapunov_r, but uses 2 decoupled Henon maps, written
//to test Lyapunov calculator.  Written 4/8/1998, just like everything else that matters


long double Lyapunov_henon(long iter, long displayfreq, long map_size){




	//int map_size=1;  //how many henon maps?
	long interval=displayfreq;	//reporting interval for lyapunov


	//input and output vectors of networks
	//are array2d's ofsize 1xd+1




	long net_size, num_inputs;
	long double s;	//the scale factor  4/8/1998, but it was there before :)
	long discard=500;	//put on attractor, should actually be an adjustable parameter
	long double d0, d1;
	long double dreal;	//actual calculated distance




	DecoupledMap* dm1;
	DecoupledMap* dm2;
	DecoupledMap* dm1new;
	DecoupledMap* dm2new;


	
	dm1=new DecoupledMap(map_size);
	dm2=new DecoupledMap(map_size);
	dm1new=new DecoupledMap(map_size);
	dm2new=new DecoupledMap(map_size);





	d0=0.00000001;

	//changed 4/15/1998 to look at changes when it's 10^-3 smaller

	//d0=0.00000000001;

	//d0=0.00001;	//let's see what happens when d0 is large

	//d0=2;
	//COMMENT:  changing d0 as above makes no difference with the henon maps







	long i,j,k;	//counters
	long double sum, lyapunov;
	sum=0;
	lyapunov=0;


	i=j=k=0;


	
	//rNetwork* net1;
	//rNetwork* net2;
	//rNetwork* net1new;
	//rNetwork* net2new;


	//net_size=net->net_size;
	//d=net->d;
	//s=net->s;
	//num_inputs=net->input_size;


	long double** map1vector; 
	long double** map2vector;
	long double** map1vectornew;
	long double** map2vectornew;



	/*


	map1vector=new long double*[map_size]; //pointers to outputs of all the neurons in the network
	map2vector=new long double*[map_size];
	map1vectornew=new long double*[map_size];
	map2vectornew=new long double*[map_size];


	*/


	


	map1vector=dm1->map;
	map2vector=dm2->map;
	map1vectornew=dm1new->map;
	map2vectornew=dm2new->map;


	


	//now allocate new temporary networks





	//net1=new rNetwork(net_size, s, num_inputs);
	//net1new=new rNetwork(net_size, s, num_inputs);
	//net2=new rNetwork(net_size, s, num_inputs);
	//net2new=new rNetwork(net_size, s, num_inputs);




	//dereferencing done in order to use copying 
	//(initialization) operator=











	//*net1=*net;	//now we are done with the net passed into us



	//iterate net1 until on attractor, then we'll copy into 
	//other temporary networks


	/* not needed for Henon map, 4/8/1998


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


		net1->propagate();


	}





  */


	//DECOUPLE the maps!
	//already decoupled!!!


	/*


	map2->randomize();
	for(i=0;i<map_size;i++){
		map1[i].discard();
		map2[i].discard();
	}


	*/


	//now assign to other networks


	sum=0;
	lyapunov=0;


	//for(i=0;i<map_size;i++){


		*dm1new=*dm1;
		*dm2=*dm1;
	//perturb orbit

		//(*(map2vector[0]))+=d0;	//we're only modifying one of the neurons.  This may have to change in the future 4/8/1998

		//time for the change!!!

		//we're modifying EVERYTHING by d0



		cout<<"\nd0 is : "<<d0<<endl;

		for(i=0;i<map_size*2;i++){



			(*(map2vector[i]))+=d0;

		}








		*dm2new=*dm2;


	//}





	//now get our input and output vectors




	//added 4/8/1998
	//now these values should stay w/ us forever.
	//otherwise we run into super-messy pointer problems


	


	/* commented out for henon map test case
	for(i=0;i<map_size;i++){	//!!! Make sure we're not changing locations of outputs in memory when we reasisgn!
		map1vector[i]=&(net1->net[i].output);
		net2vector[i]=&(net2->net[i].output);
		net1vectornew[i]=&(net1new->net[i].output);
		net2vectornew[i]=&(net2new->net[i].output);


	}


	*/


		/*


	map1vector[0]=&(map1[0]->x);	//take address, same idea as network
	map1vector[1]=&(map1[0]->y);
	map1vector[2]=&(map1[1]->x);
	map1vector[3]=&(map1[1]->y);


	map2vector[0]=&(map2[0]->x);	//take address, same idea as network
	map2vector[1]=&(map2[0]->y);
	map2vector[2]=&(map2[1]->x);
	map2vector[3]=&(map2[1]->y);


	map1vectornew[0]=&(map1new[0]->x);	//take address, same idea as network
	map1vectornew[1]=&(map1new[0]->y);
	map1vectornew[2]=&(map1new[1]->x);
	map1vectornew[3]=&(map1new[1]->y);


	map2vectornew[0]=&(map2new[0]->x);	//take address, same idea as network
	map2vectornew[1]=&(map2new[0]->y);
	map2vectornew[2]=&(map2new[1]->x);
	map2vectornew[3]=&(map2new[1]->y);


  */


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





		//now propagate the "new" versions of both networks.
		//this way, we have one "old" version and one "new" version
		//of each




		//for(k=0;k<map_size;k++){
			dm1new->advance();
			dm2new->advance();


		//}





	





	


		//calculate d1 (the inputs are updated in both networks)
		//we're passing in the 1-d array of long double* which points
		//to every single output in a given network


		//we do the same thing with the map vectors as we did with the net vectors
		d1=distanceNptr(map1vectornew, map2vectornew, map_size*2);	//net_size is also vector size
		//assert(d1!=0);



		if(d1==0){

			cout<<"\nWARNING! d1 is 0, calculation will be in error!"<<endl;

		}else{



		//readjust distance?? 4/8/1998
		for(j=0;j<map_size*2;j++){


			//net2->input->base[0][j]=
			*(map2vector[j])=				//once again, messy ptr arithmetic, 4/8/1998
				//net1new->input->base[0][j]+
				*(map1vectornew[j])+
				//(d0/d1)*((net2new->input->base[0][j]-net1new->input->base[0][j]));
				(d0/d1)*


					(*(map2vectornew[j])-*(map1vectornew[j]));



		}


		//actual readjusted distance, should be equal to d0 (4/17/1998)
		dreal=distanceNptr(map1vectornew, map2vector, map_size*2);	//net_size is also vector size
		


		
		sum+= log (d1/d0);	//base e


		lyapunov=sum/i;


		}//end else !!!



		if(i%interval==0){


			printf("\niter: %d lyapunov: %E d0: %E d1: %E dreal: %E",i, lyapunov,d0,d1,dreal);

				if(debug){
					printf("\n||dreal-d0||= %E",sqrt(sqr((float)dreal)-sqr((float)d0)));
				}


		}




		//*map2new=*map2;
		*dm2new=*dm2;


		//*(map2new[0])=*(map2[0]);
		//*(map2new[1])=*(map2[1]);




	}//end for
	return lyapunov;


}






//given 2 vectors and their dimensions, calculates distance in N-space
//vector1 and vector2 MUST be arrays of the right size!!!





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




	long i,j;
	i=j=0;
	long double innersum=0;
	long double distance=0;




	/*


	cout<<"\ndistance between:"<<endl;
	for(i=0;i<n;i++){


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


	}




	cout<<endl;


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


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





	}


	*/




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


		innersum+=sqr(vector1[i]-vector2[i]);


		//cout<<"\ninnersum: "<<innersum<<endl;




	}



	distance=sqrt(innersum);


	//cout<<distance<<endl;


	return distance;





}








//added 4/8/1998


long double distanceNptr(long double** vector1, long double** vector2, long n){



	long i,j;


	i=j=0;



	long double innersum=0;
	long double distance=0;




	/*
	cout<<"\ndistance between:"<<endl;


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


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


	}




	cout<<endl;
	for(i=0;i<n;i++){
		printf("%f ",vector2[i]);


	}




	*/


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


		innersum+=sqr(*(vector1[i])-*(vector2[i]));


		//cout<<"\ninnersum: "<<innersum<<endl;


	}



	distance=sqrt(innersum);


	//cout<<distance<<endl;


	return distance;




}







long double sqr(long double num){

	return num*num;



}





long double log2(long double num){




	static long double ln2=log(2);	//so we only have to do it once!
	//return 1.4427*log(num);	//return log base 2
	return log(num)/ln2;


}





void recurrent_benchmark(){


	char ch; //dummy value




	rNetwork* rnet;
	rNetwork* rnet2;


	long double s;
	long inputs_per_neuron=2;
	long iter;


	long size;


	long i,j,k;


	i=j=k=0;


	long testiter=50;
	long double value;


	cout<<"\nfirst testing list!"<<endl;



	List* testlist;


	testlist=new List(testiter);


	while(testlist->push(myrand())); //convoluted expression...



	while(testlist->pop(value)){


			cout<<value<<endl;


	}





	cout<<"\nrand_max is: "<<RAND_MAX;
	cout<<"\nscale factor? "<<endl;
	cin>>s;
	cout<<"\nnumber of neurons? "<<endl;
	cin>>size;
	cout<<"\ninputs per neuron? "<<endl;
	cin>>inputs_per_neuron;
	cout<<"\niterations? "<<endl;
	cin>>iter;





	cout<<"\ncreating network!"<<endl;
	rnet=new rNetwork(size, s,inputs_per_neuron);
	cout<<"\nnetwork looks like: "<<endl;
	rnet->dump();
	cout<<"\npress something and then press enter to continue commence iterations!"<<endl;
	cin>>ch;





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


		rnet->propagate();	//should use 2 networks
		rnet->dumpoutputs();	



	}


}


//to dump stuff to file
void Lyapunov_r_driver(){
	
	
	char filename[1000];
	long discards=250;
	long num_cases=900;
	long iter=1000;
	long displayfreq=2000; //never display
	long size;
	double s;
	double inputs_per_neuron;
	long double lyapunov;
	rNetwork* rn;

	/*

	ofstream outfile;
	cout<<"\nfilename to dump lyapunov values to: "<<endl;
	cin>>filename;
	outfile.open(filename);


	
	cout<<"\nLyapunov exponent for recurrent network, first implemented 4/8/1998"<<endl;
	cout<<"\nscale factor? "<<endl;
	cin>>s;
	cout<<"\nnumber of neurons? "<<endl;
	cin>>size;
	cout<<"\ninputs per neuron? "<<endl;
	cin>>inputs_per_neuron;

	*/
	//cout<<"\niterations? "<<endl;
	//cin>>iter;
	//commented out (temporarily) 4/24/1998
	s=2;
	size=100;
	inputs_per_neuron=10;
	

	
	//cout<<"\ncreating network!"<<endl;
	//rn=new rNetwork(size, s,inputs_per_neuron);
  	//cout<<"\niterations for recurrent network lyapunov calculation? ";
	//cin>>iter;
	//cout<<"\nenter number of iters to wait before displaying a result: "<<endl;
	//cin>>displayfreq;

	rn=new rNetwork(size, s,inputs_per_neuron);
	lyapunov=Lyapunov_r(rn, iter,displayfreq);
	cout<<lyapunov<<endl;
	
	//commented back in 4/24/1998, 5:15 am
	/* commented back out 4/24/1998 6:19 am
	for(long i=0;i<num_cases;i++){
		
		rn=new rNetwork(size, s,inputs_per_neuron);
		lyapunov=Lyapunov_r(rn, iter,displayfreq);
		//cout<<"\nThe Exponent after "<<iter<<" iterations is: "<<lyapunov<<endl;
		outfile<<lyapunov<<endl;
		
		
		//delete rn;
		//rn=NULL;
	}
	outfile.close();
	*/
	

}


void Lyapunov_r_benchmarkdriver(){ //to create datafile in Sprott & Albers format

	int suicide_count=1;	//values to bail out after
	char datafilename[]="neural.ini";
	ifstream indatafile;
	ofstream outdatafile;
	int fileok=0;
	
	//commented out 4/29/1998, because we now have LoadSettings() and SaveSettings()
	indatafile.open(datafilename);
	if(indatafile.peek()==EOF){

		fileok=0;
		indatafile.close();
		//outdatafile.open(datafilename);	//assuming this will work

	}else{
		fileok=1;
		indatafile.close();	//will be handled by LoadSettings()
	}
	long i,j,k;
	i=j=k=0;

	int screenoutput=1;
		char filename[1000];
	long discards=250;
	long num_cases=900;
	long iter=1000;
	long displayfreq=2000; //never display
	
	double size;
	double min_size;
	double max_size;
	double size_increment=10;

	double s;  //leave s alone for now
	double min_s;
	double max_s;

	double inputs_per_neuron;
	long min_inputs_per_neuron;
	long  max_inputs_per_neuron;
	long inputs_per_neuron_increment=10;
	double s_increment;

	int ascii=0;
	long iter_lyap=1000;
	 
	long double lyapunov;
	rNetwork* rn;


	char ch;
	

	ofstream outfile;	//for ascii
	


	if(!fileok){
		cout<<"\nfilename to dump data to: "<<endl;
		cin>>filename;
	

	
		cout<<"\nLyapunov generator, 4/24/1998"<<endl;
		cout<<"\nmin scale factor? "<<endl;
		cin>>min_s;
		cout<<"\nmax scale factor? "<<endl;
		cin>>max_s;
		cout<<"\nscale factor increment multiplier? "<<endl;
		cin>>s_increment;
		cout<<"\nmin number of neurons? "<<endl;
		cin>>min_size;
		cout<<"\nmax number ofneurons? "<<endl;
		cin>>max_size;
		cout<<"\nsize increment factor? "<<endl;	//modified to increment factors 4/29/1998 12:48 am
		cin>>size_increment;
		cout<<"\nmin inputs per neuron? "<<endl;
		cin>>min_inputs_per_neuron;
		cout<<"\nmax inputs per neuron? "<<endl;
		cin>>max_inputs_per_neuron;
		cout<<"\ninputs per neuron increment? "<<endl;
		cin>>inputs_per_neuron_increment;
		cout<<"\nuse ascii? (y/n)"<<endl;
		cin>>ch;
		if(ch=='n'||ch=='n'){
			ascii=0;
		}else{
			ascii=1;
		}


		cout<<"\niterations per case? "<<endl;
		cin>>iter;
		cout<<"\ndisplay results on screen? (y/n)"<<endl;;
		cin>>ch;
		if(ch=='n'||ch=='n'){
			screenoutput=0;
		}else{
		screenoutput=1;
		}
		//now write the values
		
		SaveSettings(filename, min_s, max_s, s_increment,min_size, max_size, size_increment,
				 min_inputs_per_neuron, max_inputs_per_neuron, inputs_per_neuron_increment, iter,
				 screenoutput, ascii);
		/* commented out and replaced by SaveSettings() function 4/29/1998
 		outdatafile<<filename<<endl;
		outdatafile<<min_s<<endl;
		outdatafile<<max_s<<endl;
		outdatafile<<s_increment<<endl;
		outdatafile<<min_size<<endl;
		outdatafile<<max_size<<endl;
		outdatafile<<size_increment<<endl;
		outdatafile<<min_inputs_per_neuron<<endl;
		outdatafile<<max_inputs_per_neuron<<endl;
		outdatafile<<inputs_per_neuron_increment<<endl;
		outdatafile<<iter<<endl;
		outdatafile<<screenoutput<<endl;
		outdatafile<<ascii<<endl;
		*/

	}else{//end of user input


			LoadSettings(filename, min_s, max_s, s_increment,min_size, max_size, size_increment,
				 min_inputs_per_neuron, max_inputs_per_neuron, inputs_per_neuron_increment, iter,
				 screenoutput, ascii);

		/* commented out 4/29/1998 1:07 am
		indatafile>>filename;
		indatafile>>min_s;
		indatafile>>max_s;
		indatafile>>s_increment;
		indatafile>>min_size;
		indatafile>>max_size;
		indatafile>>size_increment;
		indatafile>>min_inputs_per_neuron;
		indatafile>>max_inputs_per_neuron;
		indatafile>>inputs_per_neuron_increment;
		indatafile>>iter;
		indatafile>>screenoutput;
		indatafile>>ascii;
		*/
	}

		//we can use file input!!!




	outfile.open(filename, ios::app);	//open file for append, added 4/29/1998
	//outfile.open(filename, ios::binary);

	//commented out (temporarily) 4/24/1998
	//s=2;
	//size=100;
	//inputs_per_neuron=10;
	
	//no neuron will get more inputs than there are nodes in the network!!!
	//we'll make sure of that here 4/24/1998
	if(screenoutput){
		printf("\nneurons        inputs/neuron        Lyapunov      s\n");
	}
	
	if(ascii){
		if(debug){
		  outfile<<"neurons "<<"       inputs/neuron"<<"        s"<<"    lyapunov"<<endl;
		}

	}
	//for(size=min_size;size<max_size;size+=size_increment){ 
	for(size=min_size;size<max_size;size*=size_increment){	//changed 4/29/1998 to increment by MULTIPLYING
															//(APPLIES TO ALL FOLLOWING FOR LOOPS)

			  

		
		for(j=min_inputs_per_neuron;j<max_inputs_per_neuron;j*=inputs_per_neuron_increment){
		  for(s=min_s;s<max_s;s*=s_increment){
			
			  
				
			  if(j>size){	//constrain size of inputs (does this mess up the stats at all?
				inputs_per_neuron=size; 
			}else{
				inputs_per_neuron=j;
		
				
			}
				
			



			//now iterate
			
			for(i=0;i<iter;i++){
				
			
				rn=new rNetwork(size, s,inputs_per_neuron);
				lyapunov=Lyapunov_r(rn, iter_lyap,displayfreq);
				if(screenoutput){
					//printf("\n1%d   %d   %f ",size,inputs_per_neuron,lyapunov);
					cout<<(int)size<<"        "<<(float)inputs_per_neuron<<"        "<<(float)s<<"        "<<(float)lyapunov<<endl;


				}
				if(ascii){
					outfile<<(int)size<<"        "<<(float)inputs_per_neuron<<"        "<<(float)s<<"        "<<(float)lyapunov<<endl;
				}else{
					outfile.put((char)((int)log2(size)));
					outfile.put((char)((int)log2(inputs_per_neuron)));
					outfile.put((char)((int)log2(s)));
					bt.number=float(lyapunov);	//bt is a union of type bytetool
					//outfile<<(float)lyapunov;
					//outfile<<bt.bytes;
					//hopefully this is the correct byte ordering
					outfile.put((char)bt.bytes[0]);
					outfile.put((char)bt.bytes[1]);
					outfile.put((char)bt.bytes[2]);
					outfile.put((char)bt.bytes[3]);

				}


				
				//decide if we should commit suicide yet
				
			//enhancement for batch mode, program "commits suicide" after SUICIDE_VALUE iterations to be respawned.
				//this is a temporary hack for garbage collction, implemented 4/29/1998, and meant to be run under Windows NT
						//Phoenix-like wishes: 4/29/1998
					
			  if(suicide_count>SUICIDE_VALUE){
					//save current minima as starting values and bail out
						
				 
					SaveSettings(filename, s, max_s, s_increment,size, max_size, size_increment,
						inputs_per_neuron, max_inputs_per_neuron, inputs_per_neuron_increment, iter,
						screenoutput, ascii);
					cout<<"\ncommiting suicide after "<<SUICIDE_VALUE<<" iterations, respawn with [data] argument"<<endl;
					outfile<<endl<<"EOF"<<endl;
					outfile.close();
					exit(1);	//very inelegant, but NT will collect our garbage and hopefully we'll be reborn from the ashes
				}	
				
			}//end innermost for
		


		}
	}

	
	//cout<<"\ncreating network!"<<endl;
	//rn=new rNetwork(size, s,inputs_per_neuron);
  	//cout<<"\niterations for recurrent network lyapunov calculation? ";
	//cin>>iter;
	//cout<<"\nenter number of iters to wait before displaying a result: "<<endl;
	//cin>>displayfreq;

	//rn=new rNetwork(size, s,inputs_per_neuron);
	//lyapunov=Lyapunov_r(rn, iter,displayfreq);
	//cout<<lyapunov<<endl;
	
	//commented back in 4/24/1998, 5:15 am
	/* commented back out 4/24/1998 6:19 am
	for(long i=0;i<num_cases;i++){
		
		rn=new rNetwork(size, s,inputs_per_neuron);
		lyapunov=Lyapunov_r(rn, iter,displayfreq);
		//cout<<"\nThe Exponent after "<<iter<<" iterations is: "<<lyapunov<<endl;
		outfile<<lyapunov<<endl;
		
		
		//delete rn;
		//rn=NULL;
	}
	outfile.close();
	*/


	
		
			  //exit code used to be here, right above increment--4/29/1998, 3:16 am
			  suicide_count++;
	}
		
	
	outfile.close();

	
	SaveSettings(filename, s, max_s, s_increment,size, max_size, size_increment,
						inputs_per_neuron, max_inputs_per_neuron, inputs_per_neuron_increment, iter,
						screenoutput, ascii);
	cout<<"\ndone, press a key to continue!!! "<<endl;
	cin>>ch;
	exit (1);


}



int SaveSettings(char* filename, double min_s, double max_s,double  s_increment,double  min_size, double max_size, double size_increment,
				 long min_inputs_per_neuron, long max_inputs_per_neuron, long inputs_per_neuron_increment, long iter,
				 int screenoutput, int ascii){



	ofstream outdatafile;


	outdatafile.open("neural.ini");
	




		//now write the values
		outdatafile<<filename<<endl;
		outdatafile<<min_s<<endl;
		outdatafile<<max_s<<endl;
		outdatafile<<s_increment<<endl;
		outdatafile<<min_size<<endl;
		outdatafile<<max_size<<endl;
		outdatafile<<size_increment<<endl;
		outdatafile<<min_inputs_per_neuron<<endl;
		outdatafile<<max_inputs_per_neuron<<endl;
		outdatafile<<inputs_per_neuron_increment<<endl;
		outdatafile<<iter<<endl;
		outdatafile<<screenoutput<<endl;
		outdatafile<<ascii<<endl;

		outdatafile.close();



		return 1;

}


int LoadSettings(char* filename, double& min_s, double& max_s,double&  s_increment,double&  min_size, double& max_size, double& size_increment,
				 long& min_inputs_per_neuron, long& max_inputs_per_neuron, long& inputs_per_neuron_increment, long& iter,
				 int& screenoutput, int& ascii){




	ifstream indatafile;
	indatafile.open("neural.ini");
	


		indatafile>>filename;
		indatafile>>min_s;
		indatafile>>max_s;
		indatafile>>s_increment;
		indatafile>>min_size;
		indatafile>>max_size;
		indatafile>>size_increment;
		indatafile>>min_inputs_per_neuron;
		indatafile>>max_inputs_per_neuron;
		indatafile>>inputs_per_neuron_increment;
		indatafile>>iter;
		indatafile>>screenoutput;
		indatafile>>ascii;



		return 1;
}