// lb-sim
// (c) A.C. Kaporis, L.K. Kirousis, and E.C. Stavropoulos, 2006 
// 
// Purpose: 
//
// Computes the expected size of the cut of a random graph G 
// by applying the differential equation technique. 
//
// Syntax: 
//
// lb-deq d step lbucket outfile
//		d is the average degree of G
//		step step measures the time increment at the end of which 
//		each traced parameter is updated according to its 1st derivative
//		lbucket measures the availability of any traced parameter during the process
//		The output is directed to output file name (also to the standand output)  
//
// Output format: 
//
// average degree and the expected size of the cut of graph G  
//
// Reference paper: 
//
// A.C. Kaporis, L.K. Kirousis, and E.C. Stavropoulos,
// Approximating almost all instances of Max-Cut within a ratio above the H{\aa}stad threshold. 
// In Proc of ESA 2006, LNCS 4168, pp. 432--443, 2006

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

#define MAXAVDEG 50

double ifact[MAXAVDEG+1];
double n[MAXAVDEG+1][MAXAVDEG+1][MAXAVDEG+1][2];
double ee[2];
double dcut=0.0;
double CUTout=0.0;
int flag = 0; 
int paint, selectdegree, BLUE, RED;


void ComputeIfactorial (int k) 
{
	int i;
	ifact[0] = 1.0;
	ifact[1] = 1.0;
	for (i=2; i<=k; i++)
		ifact[i] = ifact[i-1] / i;
	return ;
} 


double e (int dmax, int t)
{
	double Se=0.0;
	int s, b, r;
	for (s=0; s<=dmax; s++)
		for (b=0; b<=dmax; b++)
			for (r=0; r<=dmax; r++)
				Se = Se + (double)s*n[s][b][r][t%2];
	return Se; 
} //e


double ed (int dmax, int t)
{
	double Sed=0.0;
	int s;
	for (s=0; s<=dmax; s++)
		Sed = Sed + (double)s*n[s][0][0][t%2];
	return Sed; 
} //ed


int delta ( double x, double y, double z, double x0, double y0, double z0 )
{
	if (x==x0 && y==y0 && z==z0)
		return 1;
	else
		return 0;
} //delta


double Pr (int dmax, int s0, int s, int b, int r, int t)
{
	if ( 0 <= s && s <= dmax && 0 <= b && b <= dmax && 0 <= r && r <= dmax )
		return (double) s0 * (double) s * n[s][b][r][t%2] / ee[t%2];
	else
		return 0.0;
} //Pr



void Decimation (double d, int dmax, double step, double lbucket)
{
	dcut = 0.0; 
	int td=0;
	int s; 
	while (n[1][0][0][td%2] > lbucket) 
	{	
		for (s=0; s<=dmax; s++)
		{
			n[s][0][0][(td+1)%2] = n[s][0][0][td%2] + step*( Pr(dmax,1,s+1,0,0,td%2) 
									- Pr(dmax,1,s,0,0,td%2) - (double) delta(s,0,0,1,0,0) ) ;  
		}
		dcut = dcut + step; 
		ee[(td+1)%2] = ed(dmax,(td+1)%2);
		td = td + 1; 
	}

	for (s=0; s<=dmax; s++)
		n[s][0][0][0] = n[s][0][0][(td-1)%2];
	
	ee[0] = ee[(td-1)%2];
} //Decimation


void cBucket(int pnt, int deg, int bi, int re)
{
	paint = pnt;
	selectdegree = deg;
	BLUE = bi;
	RED = re;
} //cBucket



void OneHit (int dmax, int color, int s0, int b0, int r0, double step, int t)
{
	int s, b, r;
	if (color == 1)
		CUTout = CUTout + step * (double) b0;
	else
		CUTout = CUTout + step * (double) r0;
	
	for (s=0; s<=dmax; s++)
	{
		for (b=0; b<=dmax; b++)
		{
			for (r=0; r<=dmax; r++) 
			{
				if (color == 1)
				{
					n[s][b][r][(t+1)%2] = n[s][b][r][t%2] + step*( Pr(dmax,s0,s+1,b-1,r,t%2) 
										- Pr(dmax,s0,s,b,r,t%2) - (double)delta(s,b,r,s0,b0,r0) );
				}
				else 
				{
					n[s][b][r][(t+1)%2] = n[s][b][r][t%2] + step*( Pr(dmax,s0,s+1,b,r-1,t%2) 
										- Pr(dmax,s0,s,b,r,t%2) - (double)delta(s,b,r,s0,b0,r0) );

				}
				if ( n[s][b][r][(t+1)%2] < -0.000000001 ) 
				{ 
					printf("Abnormal termination! Check parameters \"step\" and/or \"lbucket\" \n");
					flag = 1; 
					exit (1); 
				}
			} //for r
		} // for b
	} //for s
	ee[(t+1)%2] = e(dmax, (t+1)%2);
}//OneHit


void InitDeg (double d, int dmax)
{
	int s, b, r;
	
	for (s=0; s<=dmax; s++)
		for (b=0; b<=dmax; b++)
			for (r=0; r<=dmax; r++)
				if (b==0 && r==0) 
				{
					n[s][b][r][0] = (double) (exp(-d)*pow(d,s)) * (double) ifact[s];
				}
				else
					n[s][b][r][0] = 0.0;

	ee[0] = e(dmax, 0); 
} //InitDeg


int Compute_dmax(double d)
{
	int s=1;
	double degree = 1.0; 
	while (degree > 0.00000001)
	{
		s++;
		degree = exp(-d)*pow(d,s)*ifact[s]; 
	}
	if (s < MAXAVDEG)
		return s; 
	else
	{
		return MAXAVDEG;
	}
}

 
//////////////////////////////////////////////////////////////////////////////
////////////////////////////////// MAIN //////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
int main (int argc, char *argv[])
{


	if (argc != 7) {
       	printf("Syntax: lb-deq-batch [dl] [du] [dstep] [step] [lbucket] [output file name] \n"); 
		exit(1);
	}

	srand( (unsigned)time( NULL ) );
	
	ComputeIfactorial(MAXAVDEG); 

	double dl = atof(argv[1]);	
	double du = atof(argv[2]);	
	double dstep = atof(argv[3]);	
	double step = atof(argv[4]);
	double lbucket =  atof(argv[5]);	
	double d = dl;
	int dmax, globaltime; 


	int s,b,r,discr;
	int Hit = 1;
	int HitTimes; 
	double SUM=0.0;
	FILE *fp;

	while (d < du + dstep)
	{
		printf("d: %f \n", d);
		dmax = Compute_dmax(d);
		InitDeg(d, dmax); 	
		dcut = 0.0; 
		if (d < 10.0)
		{
			Decimation(d, dmax, step, lbucket);
			printf("Decimation: dcut= %f \n",dcut);
		}
		CUTout = 0.0;
		globaltime = 0;
		while (globaltime <=30) 
		{
			cBucket(rand()%2,2,0,0); 
			OneHit(dmax, paint, selectdegree, BLUE, RED, step, globaltime);
			globaltime++; 	
		}

		do
		{
			for (discr=dmax; discr>=0; discr--)
				for (s=0; s <=dmax; s++)
					for (b=0; b<=dmax; b++)
						for (r=b; r<=dmax; r++) 
						{
							if ( (int) abs(b-r) == discr && (n[s][b][r][globaltime%2] >= lbucket || n[s][r][b][globaltime%2] >= lbucket)) 
							{
								Hit = 1;
								if (n[s][b][r][globaltime%2] > n[s][r][b][globaltime%2]) //lbucket) 
									cBucket(1,s,b,r);
								else 
									if (n[s][b][r][globaltime%2] < n[s][r][b][globaltime%2]) // >= lbucket)
										cBucket(0,s,r,b);
									else
										cBucket(rand()%2,s,b,r);
							
								for (HitTimes=0; HitTimes<60; HitTimes++) //
								{
									if ( n[selectdegree][BLUE][RED][globaltime%2] >= lbucket )
									{
										OneHit(dmax, paint, selectdegree, BLUE, RED, step, globaltime);
										globaltime ++;
									}
									else
										HitTimes=60; 
								} 
								r=dmax;
								b=dmax;
								s=dmax;
								discr = 0;
							}				
							else
							{
								Hit=0;
							}
						} 
		}while(Hit); 

		for (s=0; s <=dmax; s++)
			for (b=0; b<=dmax; b++)
				for (r=0; r<=dmax; r++)
					 SUM= SUM + (double) (s+(double)(b+r)/2.0)*n[s][b][r][(globaltime-1)%2];

		fp = fopen(argv[4], "a");
	
		printf ("d=%4.1f \t cut=%10.8f \n", d, 1.0 - CUTout * 2.0 / d);
		fprintf (fp, "%4.1f \t %10.8f \n", d, 1.0 - CUTout * 2.0 / d);
	
		fclose(fp);

		d += dstep;
	}

   	return 0; 
} //main