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



//		SCIENTIFIC RESULTS:
// Percent of plasmids duplicated in S-Phase (1 or 2 digit INTEGER)
#define SPHASEDUP 84
// Percent of duplicated plasmids distributed evenly to daughter cells (1 or 2 digit INTEGER)
#define SPHASEDIST 78
// Number of plasmids equal to or below which a cell dies, -1 is no cell death
#define PLASMIDDELETE 5
// If cell death should be based on the number of plasmids (i.e. 100% of cells with
// 0 plasmids die while only 50% of cells with plasmid 1 die), set this equal to 1.
// If PLASMIDDELETE should be a cutoff, set the following number to 0.
#define PERCENTDEATH 1
// Zerodouble can be set to 1 or 2. If the 
#define ZERODOUBLE 2


//		SCIENTIFIC PROCESS:
// Set the number of generations (maximum of 1000)
#define GENERATIONS 25
// Set the name of the input file
#define INPUTFILE "1.txt"


// 		ADVANCED SETUP:
// Set the number of runs to average for a general result: note 40 is a good number
// (maximum of 1000)
#define AVRUNNUM 40



int gens[100], maxsize, totalmax=1, daughter1, daughter2, seednum, tens;
long ncount=0, ncountres, betweentonext[2], zerocount=0;
long daughterps[15][26000000], tempcells[26000000], avdist[AVRUNNUM][1000];
char *outputname;

void assignation();
void betweenruns();
void calcrepeats();
void input();
void integrate();
void nextinputs();
void output();
int rnd();
int testparams();

// *****************************************************************************
// function main calls the integration function								MAIN
// The program will only run if the parameters are valid.

main(int argc, char *argv[])
{
	int a, b, h, i, writepercent, runprog;
	float percentdone;

	if (argc > 1) {
		outputname = argv[1];
	}
	else {
		outputname = "averages.txt";
	}

	runprog = testparams();

	if (runprog == 1) {
	
	//	seed C's random number generator with the time
		srand((unsigned)time(NULL));
		
		calcrepeats();
	//	--------------
		
	//	run the main body of the program- the first loop, h, goes through the number
	//	of runs necessary to make an average while the second loop, i, steps through
	//	the number of sets of approximately ten generations (from 5-14)
		for (h = 0; h < AVRUNNUM; h++) {
			for (i = 0; i <= tens; i++) {
				if (i == 0) input();
				else nextinputs(h);
				
				integrate(gens[i]);
	//			-------------------
				
				percentdone = (((float) h*(tens + 1) + i)/((float) (tens + 1)*AVRUNNUM))*100. + .5;
				writepercent = ((int) percentdone);
				printf("%3i %%\n", writepercent);
				
				betweenruns(h, i);
	//			--------------
			}
		}
	
	//	print out the results to the averages.txt file	
		output();
	//	---------
		printf("100 %%\n");
	}
	
	if (runprog == 0) {
		printf("     ERROR: program not run\n");
	}
}

// *****************************************************************************
// function assignation assigns the plasmids to the daughter cells	 ASSIGNATION
// a random integer (0-1000) is picked and based on the SPHASEDUP
// and SPHASEDIST from the pars.h file, the duplication and distribution
// of each cell are randomly calculated

void assignation()
{
	int random, rand2;
	
	random = rnd();
	rand2 = rnd();
	
	if (random <= (SPHASEDUP*10)) {
		if (rand2 <= (SPHASEDIST*10)) {
			daughter1 = 1;
			daughter2 = 1;
		}
		else if ((SPHASEDIST*10) < rand2 < (SPHASEDIST*10 + 5*(100-SPHASEDIST))) {
			daughter1 = 2;
			daughter2 = 0;
		}
		else {
			daughter1 = 0;
			daughter2 = 2;
		}
	}
	else {
		if (rand2 <= 500) {
			daughter1 = 1;
			daughter2 = 0;
		}
		else {
			daughter1 = 0;
			daughter2 = 1;
		}
	}
}

// *****************************************************************************
// function betweenruns calculates the average for each run/prints   BETWEENRUNS

void betweenruns(int i, int l)
{
	int a, b, j, k, numberofcells, maxnumofplasmids=1;

	for (a = 0; a < AVRUNNUM; a++) {
		for (b = 0; b < 1000; b++) {
			avdist[a][b] = 0;
		}
	}

//	numberofcells is the entire number of cells in a generation. the number of
//	cells with zero plasmids is calculated separately (zerocount) because any
//	cell with zero plasmids can only have daughter cells with 0 plasmids.
	numberofcells = zerocount + ncount;
	avdist[i][0] = zerocount;
	
//	convert each cell with x plasmids into a list (avdist) where the number
//	of cells with each number of plasmids is counted. daughterps is the
//	computational approximation of real life while avdist is the count thereof
	for (j = 0; j < ncount; j++) {
		avdist[i][daughterps[gens[l]][j]]++;
		if (daughterps[gens[l]][j] > maxnumofplasmids) maxnumofplasmids = daughterps[gens[l]][j];
	}
	
	if (totalmax < maxnumofplasmids) totalmax = maxnumofplasmids;

	betweentonext[0] = numberofcells;
	betweentonext[1] = maxnumofplasmids;

//	ncount must be reset as well as zerocount and ncountres is used to step
//	through the integration function for the next time that it is used
	ncountres = ncount;
	ncount = 0;
	zerocount = 0;
}

// *****************************************************************************
// function calcrepeats calculates the number of times the code must CALCREPEATS
// be repeated so as not to exceed the maximum array size. This means
// that an actual simulation of division will only be run up to fourteen
// generations until it is converted into percents (to the tenths place)
// and used to seed the next generations. Generations are split into lengths
// of ten with the final number of consecutive generations being no smaller 
// than 5 and no greater than 14

void calcrepeats()
{
	int a, i;
	
	a = GENERATIONS%10;
	
	if (GENERATIONS < 15) {
		tens = 0;
		gens[0] = GENERATIONS;
	}
	else {
		if (a != 0 && a < 5) {
			tens = ((int) ((float) GENERATIONS)/10.) - 1;
			gens[tens] = a + 10;
		}
		else {
			tens = (int) ((float) GENERATIONS)/10.;
			gens[tens] = a;
		}
		
		for (i = 0; i < tens; i++) {
			gens[i] = 10;
		}
	}
}

// *****************************************************************************
// function input takes the numbers of plasmids per cell from a file	   INPUT
// the input file should be formatted as the number of plasmids followed
// by a space and then the number of cells containing the number of plasmids
// i.e. 25 cells with 4 plasmids each would be written as "4 25" (without
// quotation marks) in the input file specified in the pars.h file

void input()
{
	int i, j, counter=0, number[100], numplasmids[100];
	FILE *readin;
	
	readin = fopen(INPUTFILE,"r");

	while (!feof(readin)) {
		fscanf(readin,"%i %i", &number[counter], &numplasmids[counter]);
		counter++;
	}
	
	fclose(readin);

//	the counts of plasmids are then converted into a computational representation
//	of the cells. Each number up to ncount in an array is a cell with the value
//	of the position of the array (number[j]) being the number of plasmids the cell
//	has
	for (j = 0; j <= counter - 1; j++) {
		for (i = 0; i < numplasmids[j]; i++) {
			daughterps[0][ncount] = number[j];
			ncount++;
		}
	}
	
	ncountres = ncount;
	ncount = 0;
	maxsize = counter - 1;
}

// *****************************************************************************
// function integrate steps through the generations and doubles		   INTEGRATE
// each cell through the assignation function. integrate() steps
// through each plasmid in the previous generation (i-1) and
// calculates on a plasmid by plasmid basis what the distribution
// to the daughter cells is

void integrate(int l)
{
	int i, j, k, m, n, cell1, cell2, tempcount, random, acceptint;
	float acceptance, plasmiddelete;
	
	plasmiddelete = PLASMIDDELETE;
	
	for (i = 1; i <= l; i++) {
		ncount = 0;
		tempcount = 0;
		zerocount = zerocount*ZERODOUBLE;
		
		for (j = 0; j < ncountres; j++) {
			cell1 = 0;
			cell2 = 0;
			
			for (k = 0; k < daughterps[i-1][j]; k++) {
				assignation();
				cell1 = cell1 + daughter1;
				cell2 = cell2 + daughter2;
			}
			
			if (cell1 == 0) {
				zerocount++;
			}
			else {
				daughterps[i][ncount] = cell1;
				ncount++;

				if (daughterps[i][ncount] > maxsize) {
					maxsize = daughterps[i][ncount];
				}
			}
			
			if (cell2 == 0) {			
				zerocount++;
			}
			else {
				daughterps[i][ncount] = cell2;
				ncount++;			

				if (daughterps[i][ncount] > maxsize) {
					maxsize = daughterps[i][ncount];
				}
			}
		}
		
		ncountres = ncount;
		
//		Everything above is for the standard version of the code- there are no
//		cell deaths. Below is the option of deleting cells with 0 plasmids. By
//		definition, in all cases, all cells with 0 plasmids will be deleted if
//		cell death is turned on.
		if (PLASMIDDELETE > -1) {
			zerocount = 0;
		}
		
//		Below the cell deaths can either have a definite cutoff (PERCENTDEATH == 0)
//		or have a percent based cutoff depending on the PERCENTDEATH setting.
		if (PLASMIDDELETE > 0) {
			if (PERCENTDEATH == 0) {
				for (m = 0; m < ncount; m++) {
					if (daughterps[i][m] > PLASMIDDELETE) {
						tempcells[tempcount] = daughterps[i][m];
						tempcount++;
					}
				}
			}
			
			if (PERCENTDEATH == 1) {
				for (m = 0; m < ncount; m++) {
					if (daughterps[i][m] > PLASMIDDELETE) {
						tempcells[tempcount] = daughterps[i][m];
						tempcount++;
					}
					
					if (daughterps[i][m] <= PLASMIDDELETE) {
						random = rnd();
						acceptance = ((float) daughterps[i][m])/plasmiddelete;
						acceptint = (int) (acceptance*1000);
						
						if (random < acceptint) {
							tempcells[tempcount] = daughterps[i][m];
							tempcount++;
						}
					}
				}
			}
			
			for (n = 0; n < tempcount; n++) {
				daughterps[i][n] = tempcells[n];
			}
			
			ncountres = tempcount;
		}
	}
}

// *****************************************************************************
// function nextinputs creates the inputs for gens > 10. It converts  NEXTINPUTS
// the values of the cells to percents of cells with each number of 
// plasmids. The percents are then multiplied by 1000 and rounded off
// to the nearest integer and are then used to seed the next generations

void nextinputs(int h)
{
	int j, k, l, distribution[1000], numberofcells, maxnumofplasmids;
	float percentdistribution;
	
	numberofcells = betweentonext[0];
	maxnumofplasmids = betweentonext[1];
	
	for (j = 0; j <= maxnumofplasmids; j++) {
		percentdistribution = (((float) avdist[h][j])/((float) numberofcells))*1000 + .5;
		distribution[j] = ((int) percentdistribution);
	}
	
	for (k = 0; k <= maxnumofplasmids; k++) {
		for (l = 0; l < distribution[k]; l++) {
			daughterps[0][ncount] = k;
			ncount++;
		}
	}
	
	ncountres = ncount;
}

// *****************************************************************************
// function output converts the counts to percentages (much like the	  OUTPUT
// function nextinputs) and prints them to a file titled "averages.txt"

void output()
{
	int i, j, k, l, m, finaldist=0;
	float printtotals[1000], totalcells=0, printpercent=0;
	FILE *fj;
	
	fj = fopen(outputname, "w");
	
	for (i = 0; i <= totalmax; i++) {
		finaldist = 0;
		
		for (j = 0; j < AVRUNNUM; j++) {
			finaldist = finaldist + avdist[j][i];
		}
		
		printtotals[i] = (((float) finaldist)/((float) AVRUNNUM))*100.;
		totalcells = totalcells + printtotals[i];
	}

	k = totalmax;

//	This section makes it so that only those values greater than 0 are printed, even
//	if there are intermediate values that are 0.
	while (printpercent < .05) {
		printpercent = printtotals[k]/totalcells*100;
		k--;
	}
	
	for (m = 0; m <= k+1; m++) {
		printpercent = printtotals[m]/totalcells*100;
		fprintf(fj, "%3i     %3.1f\n", m, printpercent);
	}

	fclose(fj);
}

// *****************************************************************************
// function rnd creates a random number between 1-1000 using C's			 RND
// built-in random number generator. It's accuracy is adequate for
// our use.

int rnd()
{
	return (rand() % 1000) + 1;
}

// *****************************************************************************
// function testparams makes sure all parameters are valid. This	  TESTPARAMS
// prevents bus errors and segmentation faults due to user error.
// In fact, with luck it will plain old prevent user error.

int testparams()
{
	int errors=0;
	FILE *testfp;
	
	if (SPHASEDUP < 0 || SPHASEDUP > 100) {
		errors++;
		printf("   Duplication in S-phase (SPHASEDUP) must be between 0-100\n");
	}
	
	if (SPHASEDIST < 0 || SPHASEDIST > 100) {
		errors++;
		printf("   Distribution in S-phase (SPHASEDIST) must be between 0-100\n");
	}
	
	if (GENERATIONS > 999) {
		errors++;
		printf("   Generations must be less than 1000\n");
	}
	
	if (AVRUNNUM < 5) {
		printf("     WARNING: More runs should be averaged for improved accuracy\n");
	}
	
	testfp = fopen(INPUTFILE,"r");
	
	if (testfp == 0) {
		errors++;
		printf("   Input file not present-- check file/filename\n");
	}
	
	if (errors > 0) return 0;
	else return 1;
}