// On The Frequency of Chaotic Dynamics on Large Economies
// Author: Dave Albers
// e-mail: djalbers@students.wisc.edu

// note to users: this code must be recompiled for each combination of n, d,
// and s.  it was created this way in order to be able to look at any c
// combination of n, d, and s.  this code was written on the borland 4.0 c+++
// compiler.  to achieve the data format for dr. sprott's display program the
// data can be run through the convert program written for the power basic c
// compiler.

#include<math.h>
#include<iostream.h>
#include<fstream.h>
#include<stdlib.h>
#include<time.h>
#include<float.h>

// declarations for the random number generator
const long IM1 = 2147483563;
const double IM2 = 2147483399;
const double AM = (1.0/IM1);
const double IMM1 = (IM1-1);
const long IA1 = 40014;
const long IA2 = 40692;
const long IQ1 = 53668;
const long IQ2 = 52774;
const long IR1 = 12211;
const long IR2 = 3791;
const long NTAB = 32;
const double NDIV = (1+IMM1/NTAB);
const double EPS =  1.2e-7;
const double RNMX = (1.0-EPS);

// delcarations for # of neurons and the dimention of the system
const int n=64, d=64;

// number of iterations, number disregarded, and number of cases
const int thrown = 400;
const int iter = 2400;
const int cases = 1;

// array and variable declarations
long double B[n];
long double w[n][d+1];
long double y[d], dy[d];
long double ljmax, ljmin;
long double d0 = .00000001, ljsum=0.0, ljexp;
long double lypsum=0.0;

// function declarations
long double return_value(int a);
long double dreturn_value(int a);
void ljy(int a, long double f, int v);

float ran2 (long *idum);

ofstream output("testdata.dat");


//ifstream input("benchmrk.dat");


int main ()
{
	int i, j, k, l, a, u, v, m, z, x, ii, jj;
	long double f=0.0, b=0.0, ljexp, lexp;
	long double g=0.0, t=0.0, mm=0.0, qq=0.0;
	long double p;
	int s=8;
	long *idum;
	long kj;


	cout = output;

	//cin = input;

	srand(time(0));

	kj = rand();
	kj = -kj;
	idum = &kj;



	for (v=1; v<=cases; v++)
	{
		b=0;
		// assign arrays
		for (i=0; i<n; i++)
		{
			// for normal distribution of B
			for (m=1; m<=12; m++)
			{
				z = ran2(idum);
				x = x + z;
			}
			x = x - 6;
			B[i] = s*x;

			// just random B with no distribution
			B[i]  = ran2(idum);

			B[i] = B[i];
			if (B[i] >= 0)
			{
				B[i] = B[i];
			}
			else if (B[i]<0)
			{
				B[i] = -B[i];
			}

			b = b + B[i];


			for (j=0; j<=d; j++)
			{
				// initialize wij to have a normal distribution
				for (l=1; l<=12; l++)
				{
					t = ran2(idum);
					g = g + t;
				}
				g = g - 6;
				w[i][j] = s*g;
				g=0.0;
			}
		}

		// assign arrays for test
		//ii = 0;
		//for (i=0; i<128; i++)
		//{
		//	if (i<64)
		//	{
		//		input >> mm;
		//		B[i] = mm;
		//		//output << i << "  " << "next b = " << B[i] << endl;
		//		if (B[i] >= 0)
		//		{
		//			B[i] = B[i];
		//		}
		//		else if (B[i]<0)
		//		{
		//			B[i] = -B[i];
		//		}
		//
		//		b = b + B[i];
		//	}
		//	jj = 0;
		//	if (i == 64)
		//	{
		//		ii = 0;
		//	}
		//	if (i>=64)
		//	{
		//		for (jj=0; jj<=64; jj++)
		//		{
		//			input >> mm;
		//			w[ii][jj] = s*mm;
		//			//output << ii << "  " << "  " << jj << "next w = " <<
		//			//	w[ii][jj] << endl;
		//		}
		//		ii= ii + 1;
		//	}
		//}
		// end of test sequence //

		// to scale too n^2 instead of 1
		//b = b/(n*n);

		for (i=0; i<n; i++)
		{
			B[i] = B[i]/b;
		}

		for (j=0; j<d; j++)
		{
			// initialize y[j]
			y[j] = 2 * ran2(idum)-1;
			// for test
			//	y[j] = .5;
		}

		// for henon map only
		//y[0] = 2 * ran2(idum) - 1;

		for (a=0; a<=iter; a++)
		{
			f = return_value(a);

			for (u=0; u<(d-1); u++)
			{
				y[u] = y[u+1];
			}
			y[d-1] = f;
			ljy(a, f, v);
		}

	}

	return 0;

}



long double return_value(int a)
{
	int i, j;
	long double p, xn, xn1, x;
	long double fn = 0;

	// henon map iterations
	//fn = 1 - 1.4*y[1]*y[1] + .3*y[0];

	for (i=0; i<n; i++)
	{
		p = w[i][0];
		for (j=1; j<=d; j++)
		{
			p = p + w[i][j]*y[j-1];
		}
		fn = fn + B[i]*tanh(p);
	}

	// for test output first 10 y's for benchmark case
	//if (a<10)
	//{
	//	output << "f value = " << fn << endl;
	//}

	return fn;
}

long double dreturn_value(int a)
{
	int i, j;
	long double dp;
	long double dfn = 0; // set to zero when not iterating for henon map

	// henon map iterations
	//dfn = 1 - 1.4*dy[1]*dy[1] + .3*dy[0];

	for (i=0; i<n; i++)
	{
		dp = w[i][0];
		for (j=1; j<=d; j++)
		{
			dp = dp + w[i][j]*dy[j-1];
		}
		dfn = dfn + B[i]*tanh(dp);
	}

	return dfn;
}

void ljy(int a, long double f, int v)
{
	long double dstt, dst, df, rs, something;
	int i,j,t;
	if (a == 0)
	{
		dy[0] = y[0] + d0;
		for (i=1; i<d; i++)
		{
			dy[i] = y[i];
		}
		ljsum = 0;
		j = 0;
	}

	if (a==thrown)
	{
		ljsum = 0;
	}



	something = dreturn_value(a);



	for(i=0; i<(d-1); i++)
	{
		dy[i] = dy[i+1];
	}
	dy[d-1] = something;

	dstt = 0;

	for (i=0; i<d; i++)
	{
		dst = dy[i] - y[i];
		dstt = dstt + dst*dst;
	}


	if (dstt > 0)
	{
		df = dstt/(d0*d0);
		rs = 1/sqrt(df);
		for (i=0; i<d; i++)
		{
			dy[i] = y[i] + rs * (dy[i] - y[i]);
		}

		ljsum = ljsum + log(df);
	 }

	 else
	 {
		ljsum = -2e30 * (a-thrown);
		//output << "i like to eat" << endl;
	 }

	 if ( a > thrown )
	 {
		ljexp = ljsum/(a-thrown);

		if ( a == iter )
		{
			output << ljexp << endl;
			lypsum = lypsum + ljexp;
		}
	 }
	 
	  
}

float ran2 (long *idum)
{	
	int j;
	long k;
	static long idum2=123456789;
	static long iy=0;
	static long iv[NTAB];
	float temp;
	
	if (*idum <= 0) 
	{
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		idum2=(*idum);
		for (j=NTAB+7; j>=0; j--)
		{
			k=(*idum)/IQ1;
			*idum=IA1*(*idum-k*IQ1)-k*IR1;
			if (*idum < 0) *idum += IM1;
			if (j < NTAB) iv[j] = *idum;	
		}
		iy=iv[0];
	}
	
	k=(*idum)/IQ1;
	*idum=IA1*(*idum-k*IQ1)-k*IR1;
	if (*idum < 0) *idum += IM1;
	k=idum2/IQ2;
	idum2=IA2*(idum2-k*IQ2)-k*IR2;
	if (idum2 < 0) idum2 += IM2;
	j=iy/NDIV;
	iy=iv[j]-idum2;
	iv[j] = *idum;
	if (iy < 1) iy += IMM1;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp; 
}