// 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)
{
	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 != 5) {
       	printf("Syntax: lb-deq [d] [step] [lbucket] [output file name] \n"); 
		exit(1);
	}

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

	double d = atof(argv[1]);	
	printf("d: %f \n", d); 

	int dmax = Compute_dmax(d);
	printf("dmax: %d \n", dmax); 
	
	InitDeg(d, dmax); 	

	double step = atof(argv[2]);
	printf("step: %e \n", step); 

	double lbucket =  atof(argv[3]);	
	printf("lbucket: %e \n", lbucket); 

	if (d<10.0)
	{
		dcut=0.0;
		Decimation(d, dmax, step, lbucket);
		printf("Decimation: dcut= %f \n",dcut);
	}
	else
		printf("No decimation\n");

	int globaltime = 0;
	while (globaltime <=30) 
	{
		cBucket(rand()%2,2,0,0); 
		OneHit(dmax, paint, selectdegree, BLUE, RED, step, globaltime);
		globaltime++; 	
	}

	int s,b,r,discr;
	int Hit = 1;
	int HitTimes; 

	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); 

	double SUM=0.0;
	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];

	FILE *fp;
	fp = fopen(argv[4], "a");

	printf ("d=%4.1f \t cut=%10.8f \n", d, 1.0 - CUTout * 2.0 / d);
	fprintf (fp, "d=%4.1f \t cut=%10.8f \n", d, 1.0 - CUTout * 2.0 / d);
	
	fclose(fp);
		
   	return 0; 
} //main