//========================================
//Info on loading of the metabolic network
//========================================
//network metabolites are accessed by the array spec[]
//network reactions are accessed by the array reac[]
//reac[] = reaction array
	//*.nreactants = number of reactants
	//*.nproducts = number of products
	//*.reversible = reversible reaction (=1) or not (=0)
	//*.k = normalized transcript level at current experimental time point
	//*.RFrp = reaction factor for the direction  R -> P
	//*.RFpr = reaction factor for the direction  P -> R
	//*.r[] = reactant array
		//*.index = position of metabolite in spec[]
		//*.smetry = stoichiometry
	//*.p[] = product array
		//*.index = position of metabolite in spec[]
		//*.smetry = stoichiometry
	//*.kc[] = array for experimental time points (n=ExpIter)
		//*.knorm = contains normalized transcript levels
//spec[] = network metabolite (species) array 
	//*.id = KEGG identifier of network metabolite
	//*.plus = variable that contains additions for every metabolite
	//*.minus = variable that contains subtraktions for every metabolite
	//*.SF = temporary variable for saving ratios of plus and minus
	//*.calc[] = array for all iterations (n=iter)
		//*.y = network metabolite level at each iteration
	//*.kc[] = array for experimental time points (n=ExpIter)
		//*.y = integral (=CMV)
//=========	
//Variables
//=========
unsigned int ExpIter; //number of time points of experimental dataset
unsigned int iter; //number of desired iterations
double K; //constant of the algorithm (S/(S+K)), e.g. K=5
double MaxFlow; //maximum desired metabolite transfer at each iteration
unsigned int nspec; //number of network metabolites
unsigned int nreac; //number of network reactions
double valueRP=valuePR=0;
double valuer=valuep=0; 
double minRP, minPR;
unsigned int minRPi,  minPRi;
double tempSF;
double S=0;
double Simpson;
double fa, fb, fab2;
double a, b, ab6;
for (unsigned int t=0; t<=(ExpIter); t++) //outer loop
{		
	for (unsigned int i=1; i<=nspec; i++) spec[i].calc[0].y=1;	//set level of every metabolite to 1
	for (unsigned int i=0; i<nreac; i++) reac[i].k=reac[i].kc[t].knorm;	//annotate normalized transcript-value of current experimental time point to reac[i].k
	a=b=0;
	unsigned int x=0; //counter for integral calculation (counts up to 2)
	for (unsigned int xi=0; xi<iter; xi++) //inner loop
	{
		for (unsigned int i=1; i<=nspec; i++) spec[i].minus=spec[i].plus=0; 
		for (unsigned int i=0; i<nreac; i++) 
		{		
		//===================================
		//Calculation of reaction factor RF (L)
			//direction R -> P
			valuer=valuep=1;
			minRP=minPR=nspec*2+100000; //this is just a high value
			for (unsigned int r=0; r<reac[i].nreactants; r++) 
			{
				S=spec[reac[i].r[r].index].calc[x].y / reac[i].r[r].smetry;
				if (S<=minRP)
				{
					minRP=S;
					minRPi=reac[i].r[r].index;
				}
				valuer=valuer * (S / (S + K)); 
			}
			reac[i].RFrp=valuer; 
			//direction P -> R
			if (reac[i].reversible)
			{
				valuer=valuep=1;
				for (unsigned int p=0; p<reac[i].nproducts;  p++) 
				{
					S=spec[reac[i].p[p].index].calc[x].y / (reac[i].p[p].smetry);
					if (S<=minPR)
					{
						minPR=S;
						minPRi=reac[i].p[p].index;
					}
					valuep=valuep * (S / (S + K)); 
				}			
				reac[i].RFpr=valuep;		
			}
			//===================================
			//Calculation of plus and minus for the metabolites in this reaction
			//direction R -> P
			for (unsigned int r=0; r<reac[i].nreactants; r++) 
			{
				valuer= reac[i].RFrp * reac[i].k * minRP * reac[i].r[r].smetry;
				spec[reac[i].r[r].index].minus=spec[reac[i].r[r].index].minus + valuer;
				valuep=valuep + valuer;		
			}
			valuer=0;
			for (unsigned int p=0; p<reac[i].nproducts;  p++) valuer=valuer + reac[i].p[p].smetry;
			for (unsigned int p=0; p<reac[i].nproducts;  p++) spec[reac[i].p[p].index].plus=spec[reac[i].p[p].index].plus + ((valuep/ valuer) * reac[i].p[p].smetry);
			//direction P -> R
			if (reac[i].reversible)
			{
				valuer=valuep=0;
				for (unsigned int p=0; p<reac[i].nproducts; p++) 
				{
					valuer= reac[i].RFpr * reac[i].k * minPR * reac[i].p[p].smetry;
					spec[reac[i].p[p].index].minus=spec[reac[i].p[p].index].minus + valuer;
					valuep=valuep + valuer;	
				}
				valuer=0;
				for (unsigned int r=0; r<reac[i].nreactants;  r++) valuer=valuer + reac[i].r[r].smetry;
				for (unsigned int r=0; r<reac[i].nreactants;  r++) spec[reac[i].r[r].index].plus=spec[reac[i].r[r].index].plus + ((valuep/ valuer) * reac[i].r[r].smetry);
			}
		}
		//== Determine the step size of iteration for the case minus>plus
		for (unsigned int i=1; i<=nspec; i++) {if ((spec[i].plus-spec[i].minus)<0) {spec[i].SF=(spec[i].minus-spec[i].plus)/spec[i].calc[x].y;} else {spec[i].SF=0;}}
		tempSF=spec[1].SF;
		for (unsigned int i=2; i<=nspec; i++)
		{
			if (spec[i].SF>tempSF) tempSF=spec[i].SF;
		}
		//== Consider maximum flow if desired
		IterFac=MaxFlow / tempSF;
		//== Determine the step size of iteration for the case minus<plus
		for (unsigned int i=1; i<=nspec; i++) {if ((spec[i].plus-spec[i].minus)>0) {spec[i].SF=(spec[i].plus-spec[i].minus)/spec[i].calc[x].y;} else {spec[i].SF=0;}}
		tempSF=spec[1].SF;
		for (unsigned int i=2; i<=nspec; i++)
		{
			if (spec[i].SF>tempSF) tempSF=spec[i].SF;	
		}
		//== Consider maximum flow if desired			
		if ((MaxFlow / tempSF)<IterFac) IterFac=(MaxFlow / tempSF);
		if (IterFac>1) IterFac=pow(IterFac,0.1); //non-linear step size calculation;  empirical trade-off between analysis speed and numerical stability
		b=b+Iterfac; //save step size for later integral calculation
		//== Calculation of new network metabolite levels after the iteration 
		for (unsigned int i=1; i<=nspec; i++) spec[i].calc[x+1].y=spec[i].calc[x].y + (IterFac * (spec[i].plus-spec[i].minus));
		x++;
		if (x==2)
		{
		//=================================================================
		// Calculation of integral of the curve using Simpson's rule => CMV
		//=================================================================				
			b = b+IterFac;
			ab6=(b-a)/6;
			for (unsigned int i=1; i<=nspec; i++)
			{
				fa = spec[i].calc[x-2].y;
				fab2 = spec[i].calc[x-1].y;
				fb = spec[i].calc[x].y;
				Simpson= ( ab6 * (fa + (4* fab2) + fb) );
				spec[i].kc[t].y=spec[i].kc[t].y + Simpson;
			}			
			//== Set value at timepoint x (2) back to zero to start next iterations
			for (unsigned int i=1; i<=nspec; i++) spec[i].calc[0].y=spec[i].calc[x].y;
			a=b;
			x=0;
		}
	} //end of loop for iterations (xi)
} //end of loop for experimental time points (t)
//==========
//Print CMV
//==========	
ofstream CMV("CMV.txt", ios::app);
CMV<<"ID";
for (unsigned int x=0; x<=ExpIter; x++) CMV<<";t"<<x;
CMV<<"\n"; 
for (unsigned int i=1; i<=nspec; i++)
{
	CMV<<spec[i].id;
	for (unsigned int x=0; x<=ExpIter; x++) CMV<<";"<<(spec[i].kc[x].y);
	CMV<<"\n";
}
CMV.close();