// RetroEncabulator.cpp - Sagar Indurkhya 2008
#include "Intel_Compatibility.h"
#include <iostream>
#include <fstream>
#include <math.h>
#include <string>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <list>
#include "dSFMT/dSFMT.h"

using namespace std;

struct Reaction
{
	int r1, r2, p1, p2;
	int chemIndex, numElementsInChemDG;
	double rate;
	double realRate;
	double sumPropensity;
	double timeSinceLastUpdate;
};

vector<int> chems;
vector<Reaction> rxns;
vector<int> RxnToAffectedChem;
vector<int> ChemQuad;
vector<int> ChemToAffectedRxns;

int numReactions;
int numChemicals;
int sizeOfDependencyGraph;

struct Chemical
{
	int ID;
	string Name;
	int initNumber;
};

vector<Chemical> FinalChemicalConc;

void InitializeSystem(string rxnFileName, string speciesFileName, string dgFileName)
{
    cout << "Processing input files...";

    ifstream rxnFile(rxnFileName.c_str());
	ifstream speciesFile(speciesFileName.c_str());
    ifstream dgFile(dgFileName.c_str());

    speciesFile >> numChemicals;
	cout << "Num Chemicals = " << numChemicals << std::endl;

    for(int i = 0; i < numChemicals; i++)
	{
		Chemical c;
	
		int blah, chemVal;
		string blahblah;
        speciesFile >> c.ID >> c.Name >> c.initNumber;
		chems.push_back(c.initNumber);
		FinalChemicalConc.push_back(c);
	}

	cout << std::endl << "Processed Species...";

	rxnFile >> numReactions;
	cout << "Num Reactions = " << numReactions << std::endl;
	
	for(int i = 0; i < numReactions; i++) {
		Reaction r;
		rxns.push_back(r);
		rxnFile >> rxns[i].r1 >> rxns[i].r2 >> rxns[i].p1 >> rxns[i].p2 >> rxns[i].realRate;
		rxns[i].rate = rxns[i].realRate;
		rxns[i].sumPropensity = 0.0;
		
		if(rxns[i].r1 == rxns[i].r2)
			rxns[i].rate /= 2.0;

		rxns[i].sumPropensity = rxns[i].timeSinceLastUpdate = 0.0;
		
		bool badReaction = false;

		if(rxns[i].r1 >= numChemicals)
			badReaction = true;
		if(rxns[i].r1 >= numChemicals)
			badReaction = true;
		if(rxns[i].p1 >= numChemicals)
			badReaction = true;
		if(rxns[i].p1 >= numChemicals)
			badReaction = true;

		if(badReaction)	{
			cout << "BAD REACTION DETECTED" << std::endl;
			exit(0);
		}
	}
	
	cout << std::endl << "Processed Reactions..." << std::endl;
	
	cout << "Input File Processing ... FINISHED." << std::endl;

	cout << "Processing Dependency Graph...";
	cout << "\t" << "Beginning Stage 1" << std::endl;
	
	// First for each reaction read which chemicals are affected
	for(int i = 0; i < numReactions; i++) {
		rxns[i].chemIndex = RxnToAffectedChem.size();
		dgFile >> rxns[i].numElementsInChemDG;
		
		for(int j = 0; j < rxns[i].numElementsInChemDG; j++) {
			int a, b;
			dgFile >> a >> b;
			RxnToAffectedChem.push_back(a);
			RxnToAffectedChem.push_back(b);
		}
	}
	
	cout << "\t" << "Beginning Stage 2" << std::endl;
	// Second for each chemical give the quadruple with CSM and DSM Info
	for(int i = 0; i < numChemicals; i++) {
		int a, b;
		dgFile >> a >> b; // CSM and DSM sizes
		ChemQuad.push_back(-1);
		ChemQuad.push_back(a);
		ChemQuad.push_back(-1);
		ChemQuad.push_back(b);
	}
	
	cout << "\t" << "Beginning Stage 3" << std::endl;
	// Third for each chemical we need to know the affected reactions
	// Now for each chemical we need to know 
	for(int i = 0; i < numChemicals; i++) {
		ChemQuad[4*i] = ChemToAffectedRxns.size();

		for(int j = 0; j < ChemQuad[i*4 + 1]; j++) {
			int Temp;
			dgFile >> Temp;
			ChemToAffectedRxns.push_back(Temp);
		}

		ChemQuad[4*i + 2] = ChemToAffectedRxns.size();

		for(int j = 0; j < ChemQuad[i*4 + 3]; j++) {
			int Temp;
			dgFile >> Temp;
			ChemToAffectedRxns.push_back(Temp);
		}
	}
	
	for(int i = 0; i < 4*numChemicals; i++) {
		if(ChemQuad[i] < 0)	{
			cout << "EXTREME FAILURE!!!" << std::endl;
			exit(0);
		}
	}
	
	cout << "ChemQuad.size() = " << ChemQuad.size() << std::endl;
	cout << "ChemToAffectedRxns.size() = " << ChemToAffectedRxns.size() << std::endl;
	cout << "FINISHED" << std::endl << std::endl;
}

bool RxnComparison(const Reaction &a, const Reaction &b) { 
	return (a.sumPropensity >= b.sumPropensity); 
}

/*
	Arguments = 
	1) maxNumIterations
	2) OutputChemNumber
	3) ReactionFileName
	4) SpeciesFileName
	5) SingletDependencyGraphFileName
	6) OUTPUT FILE NAME FOR SORTED REACTIONS
	7) OUTPUT FILE NAME FOR SUM OF PROPENSITIES
	8) OUTPUT FILE NAME FOR CHEMICAL CONCENTRATIONS
	9) ChemicalOutputFileName
	10) maxSimTime
*/
int main(int argc, char* argv[])
{
	cout << "********** RETRO-ENCABULATOR **********" << std::endl;
	cout << "Sagar Indurkhya, 2008" << std::endl << std::endl;
	
	InitializeSystem(string(argv[3]), string(argv[4]), string(argv[5]));
	
	int numIterations = atoi(argv[1]);
	int ChemicalToOutput = atoi(argv[2]);
	
	double simTime = 0.0, a0 = 0.0, maxSimTime = (double)(atoi(argv[10]));
	
	vector<double> propensities;

	for(int i = 0; i < numReactions; i++)
		propensities.push_back(0.0);

	// Start up our random number generator
	dsfmt_t dsfmt;
	dsfmt_init_gen_rand(&dsfmt, 6426);

	for(int i = 0; i < numReactions; i++) {
		if(rxns[i].r1 != rxns[i].r2)
			propensities[i] = chems[rxns[i].r1] * chems[rxns[i].r2] * rxns[i].rate;
		else
			propensities[i] = chems[rxns[i].r1] * (chems[rxns[i].r1] - 1) * rxns[i].rate;
		a0 += propensities[i];
	}
	
	cout << std::endl << std::endl;
	cout << "a0 = " << a0 << std::endl;
	cout << "maxSimTime = " << maxSimTime << std::endl;
	cout << "numIterations = " << numIterations << std::endl << std::endl;
	cout << "chemicalToOutput = " << ChemicalToOutput << std::endl;
	cout << "The program is ready to run. All code from here down has been ";
	cout << "carefully written to run extremely efficiently. Enter any key ";
	cout << "to begin the simulation, and the program will stay silet until ";
	cout << "the simulation is done." << std::endl;

	//int lol;

	//cin >> lol;
	
	// OPTIMIZED CODE BELOW!!!
	
	unsigned startTime = 0; // Set start time here.

	ofstream chemicalOutputSS(argv[9]);
	
	// Do simulation here
	for(int numIterationsCompleted = 0; numIterationsCompleted < numIterations && simTime < maxSimTime; ) {
		// We should now record which chemcicals we need to output
		// to the GNUPlot file
		/*if(numOutputChems > 0)
		{
			outputFile << simTime << "	";
			
			for(int i = 0; i < numOutputChems; i++)
				outputFile << chems[outputChems[i]] << "	";
			
			outputFile << std::endl;
		}*/

		double randNum1 = a0*(double)dsfmt_genrand_close_open(&dsfmt);

		double runningSum = 0.0;
		
		for(int i = 0; i < numReactions; i++) {
			runningSum += propensities.at(i);

			if(randNum1 > runningSum)
				continue;
			else {
				//cout << "Reaction chosen = " << i << std::endl;
				// time to execute the reaction!
				simTime += -1.0/a0 * log(dsfmt_genrand_close_open(&dsfmt));

				//cout << "Reaction selected = " << i << std::endl;
				
				for(int j = 0, zilch = rxns.at(i).chemIndex; j < rxns.at(i).numElementsInChemDG; j++, zilch += 2) {
					// First handle singlets
					int curChem = RxnToAffectedChem[zilch];
					int delta = RxnToAffectedChem[zilch + 1];

					chems[curChem] += delta;
					
					//cout << "\t" << "curChem = " << curChem << " delta = " << delta << std::endl;
					
					//cout << "\t\t" << "ChemQuad.at(4*curChem) = " << ChemQuad.at(4*curChem) << std::endl;
					//cout << "\t\t" << "ChemQuad.at(4*curChem + 1) = " << ChemQuad.at(4*curChem + 1) << std::endl;
					
					for(int k = 0, pos = ChemQuad.at(4 * curChem); k < ChemQuad.at(4 * curChem + 1); k++, pos++) {
						int rxnPos = ChemToAffectedRxns.at(pos);
						
						double diff = chems.at(rxns.at(rxnPos).r1) * chems.at(rxns.at(rxnPos).r2) * rxns.at(rxnPos).rate - propensities.at(rxnPos);
						
						rxns.at(rxnPos).sumPropensity += (simTime - rxns.at(rxnPos).timeSinceLastUpdate) * propensities.at(rxnPos);
						rxns.at(rxnPos).timeSinceLastUpdate = simTime;
						propensities.at(rxnPos) += diff;
						a0 += diff;
					}
					
					//cout << "\t\t" << "Singlets completed" << std::endl;
					
					for(int k = 0, pos = ChemQuad.at(4 * curChem + 2); k < ChemQuad.at(4 * curChem + 3); k++, pos++) {
						int rxnPos = ChemToAffectedRxns.at(pos);
	
						double diff = chems.at(rxns.at(rxnPos).r1) * (chems.at(rxns.at(rxnPos).r1) - 1) * rxns.at(rxnPos).rate - propensities.at(rxnPos);
						
						rxns.at(rxnPos).sumPropensity += (simTime - rxns.at(rxnPos).timeSinceLastUpdate) * propensities.at(rxnPos);
						rxns.at(rxnPos).timeSinceLastUpdate = simTime;
						propensities.at(rxnPos) += diff;
						a0 += diff;
					}
		
					//cout << "\t\t" << "Doublets completed" << std::endl;
				}
				
				//cout << "\t" << "UPDATED REACTIONS" << std::endl;
				
				if(chems.at(rxns.at(i).r1) < 0) {
					cout << "EPIC FAIL, Chemical Value is < 0 for r1 chem: " << rxns.at(i).r1 << std::endl;
					cout << "Iteration = " << numIterationsCompleted << std::endl;
					cout << "Reaction = " << i << std::endl;
					cout << "runningSum = " << runningSum << std::endl;
					cout << "randNum1 = " << randNum1 << std::endl;
					cout << "a0 = " << a0 << std::endl;
					exit(0);
				}
				
				if(chems.at(rxns.at(i).r2) < 0) {
					cout << "EPIC FAIL, Chemical Value is < 0 for r2 chem: " << rxns.at(i).r2 << std::endl;
					cout << "Iteration = " << numIterationsCompleted << std::endl;
					cout << "Reaction = " << i << std::endl;
					cout << "runningSum = " << runningSum << std::endl;
					cout << "randNum1 = " << randNum1 << std::endl;
					cout << "a0 = " << a0 << std::endl;
					exit(0);
				}
				
				if(numIterationsCompleted % 100 == 0) {
					if(numIterationsCompleted % 10000 == 0)
						cout << "Iteration: " << numIterationsCompleted << "	Executed reaction #: " << i << "	" << rxns.at(i).sumPropensity << std::endl;
					if(ChemicalToOutput >= 0)
						chemicalOutputSS << simTime << "	" << chems.at(ChemicalToOutput) << std::endl;
					else {
						chemicalOutputSS << simTime;
						for(int z = 0; z < numChemicals; z++)
							chemicalOutputSS << "	" << chems.at(z);
						chemicalOutputSS << std::endl;
					}
				}
		
				numIterationsCompleted++;
				
				break;
			}
		}

		// Now see if a0 is > 0
		if(a0 <= 0.0000000001) {
			cout << "Simulation Terminated because a0 < 0.00000001, a0 = " << a0 << std::endl;
			break;
		}
	}
	
	unsigned int stopTime = 0; // Put finish time here

	for(int i = 0; i < numReactions; i++)
		rxns.at(i).sumPropensity += (simTime - rxns.at(i).timeSinceLastUpdate) * propensities.at(i);
	
	cout << std::endl << "SIMULATION FINISHED!" << std::endl;
	
	cout << "Normalizing Propensities..." << std::endl;
	// Normalize the sum propensities!
	double sumOfSumPropensities = 0.0;
	
	for(int i = 0; i < numReactions; i++)
		sumOfSumPropensities += rxns.at(i).sumPropensity;
	
	cout << "sumOfSumPropensities = " << sumOfSumPropensities << std::endl;
	
	for(int i = 0; i < numReactions; i++)
		rxns.at(i).sumPropensity /= sumOfSumPropensities;
	
	// Sort the reactions based on the sum of propensties (from highest to lowest)
	cout << "Sorting Propensities..." << std::endl;

	
	list<Reaction> tempRxns;
	
	for(int i = 0; i < numReactions; i++)
		tempRxns.push_back(rxns[i]);
	
	cout << "Num Reactions = " << numReactions << " and rxns.size() = " << rxns.size() << std::endl;
	
	tempRxns.sort(RxnComparison);
	
	rxns.clear();

	for(list<Reaction>::iterator iter = tempRxns.begin(); iter != tempRxns.end(); iter++)
		rxns.push_back(*iter);
	
	// Now lets do our outputting
	
	// First lets output our standard file with chemicals and reactions
	cout << "Outputting Reactions..." << std::endl;
	ofstream basicOutput(argv[6]);
	
	basicOutput << numReactions << std::endl << std::endl;
	
	for(int i = 0; i < numReactions; i++)
		basicOutput << rxns[i].r1 << "	" << rxns[i].r2 << "	" << rxns[i].p1 << "	" << rxns[i].p2 << "	" << rxns[i].realRate << std::endl;
	
	// Output a file containing the sum of propensities data in a format that is graphable
	cout << "Outputting Sum of Propensities..." << std::endl;
	
	ofstream sumPropsOutput(argv[7]);
	
	for(int i = 0; i < numReactions; i++)
		sumPropsOutput << i + 0.5 << "	" << rxns[i].sumPropensity << std::endl;
	
	// Now output the current chemical concentrations
	cout << "Outputting Chemical Concentrations..." << std::endl;
	
	ofstream chemConcs(argv[8]);
	
	chemConcs << numChemicals << std::endl << std::endl;
	
	for(int i = 0; i < numChemicals; i++)
		chemConcs << FinalChemicalConc[i].ID << "	" << FinalChemicalConc[i].Name << "	" << chems[i] << std::endl;

	cout << std::endl;
	cout << "Press any key to exit simulation.";

	cin >> numIterations;
	
	return 0;
}

