// TurboEncabulator.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 sdgIndex, numElementsInSDG;
	int ddgIndex, numElementsInDDG;
	double rate;
	double realRate;
	double sumPropensity;
	double timeSinceLastUpdate;
};

vector<int> chems;
vector<Reaction> rxns;
vector<int> SingletDependencyGraph;
vector<int> DoubletDependencyGraph;

int numReactions;
int numChemicals;
int sizeOfDependencyGraph;

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

vector<Chemical> FinalChemicalConc;

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

    ifstream rxnFile(rxnFileName.c_str());
	ifstream speciesFile(speciesFileName.c_str());
    ifstream sdgFile(singletDependencyGraphFileName.c_str());
	ifstream ddgFile(doubletDependencyGraphFileName.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...";
	
	for(int i = 0; i < numReactions; i++) {
		sdgFile >> rxns[i].numElementsInSDG;
		ddgFile >> rxns[i].numElementsInDDG;
		rxns[i].sdgIndex = SingletDependencyGraph.size();
		rxns[i].ddgIndex = DoubletDependencyGraph.size();

		for(int j = 0; j < rxns[i].numElementsInSDG; j++) {
			int t;
			sdgFile >> t;
			SingletDependencyGraph.push_back(t);
			if(t >= numReactions) {
				cout << "BAD SDG Detected!!!" << std::endl;
				exit(0);
			}
		}

		for(int j = 0; j < rxns[i].numElementsInDDG; j++) {
			int t;
			ddgFile >> t;
			DoubletDependencyGraph.push_back(t);
			if(t >= numReactions) {
				cout << "BAD DDG Detected!!!" << std::endl;
				exit(0);
			}
		}
	}
	
	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) DoubletDependencyGraphFileName
	7) OUTPUT FILE NAME FOR SORTED REACTIONS
	8) OUTPUT FILE NAME FOR SUM OF PROPENSITIES
	9) OUTPUT FILE NAME FOR CHEMICAL CONCENTRATIONS
	10) ChemicalOutputFileName
	11) maxSimTime
*/
int main(int argc, char* argv[])
{
	cout << "********** TURBO-ENCABULATOR **********" << std::endl;
	cout << "Sagar Indurkhya, 2008" << std::endl << std::endl;
	
	InitializeSystem(string(argv[3]), string(argv[4]), string(argv[5]), string(argv[6]));
	
	int numIterations = atoi(argv[1]);
	int ChemicalToOutput = atoi(argv[2]);
	
	double simTime = 0.0, a0 = 0.0, maxSimTime = (double)(atoi(argv[11]));
	
	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[10]);
	
	// 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;
		}*/

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

		float 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));
				
				// Update chemical species
				chems.at(rxns.at(i).r1) -= 1;
				chems.at(rxns.at(i).r2) -= 1;
				chems.at(rxns.at(i).p1) += 1;
				chems.at(rxns.at(i).p2) += 1;
				
				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;
					
					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;
					exit(0);
				}
				
				if(chems.at(rxns.at(i).p1) < 0) {
					cout << "EPIC FAIL, Chemical Value is < 0 for p1 chem: " << rxns.at(i).p1 << std::endl;
					cout << "Iteration = " << numIterationsCompleted << std::endl;
					exit(0);
				}

				if(chems.at(rxns.at(i).p2) < 0) {
					cout << "EPIC FAIL, Chemical Value is < 0 for p2 chem: " << rxns.at(i).p2 << std::endl;
					cout << "Iteration = " << numIterationsCompleted << std::endl;
					exit(0);
				}

				// Reset the Null Chem
				chems[0] = 1;

				// Now recalculate all the appropriate SINGLET propensities
				//cout << "\t" << "Calculating SINGLET Propensities..." << std::endl;
				int numElementsInSDG = rxns.at(i).numElementsInSDG;
				for(unsigned int j = 0, pos = rxns.at(i).sdgIndex; j < numElementsInSDG; ++pos, j++) {
					//cout << "\t\t" << "Calculating Position in SDG for reaction " << i << " ..." << std::endl;
					int position = SingletDependencyGraph.at(pos);
					//cout << "\t\t\t" << "Calculating Difference in SDG ..." << std::endl;
					float diff = chems.at(rxns.at(position).r1) * chems.at(rxns.at(position).r2) * rxns.at(position).rate - propensities.at(position);

					rxns[position].sumPropensity += (simTime - rxns.at(position).timeSinceLastUpdate) * propensities.at(position);
					rxns[position].timeSinceLastUpdate = simTime;
					propensities.at(position) += diff;
					a0 += diff;
				}

				// Now recalculate all the appropriate DOUBLET propensities
				//cout << "\t" << "Calculating DOUBLET Propensities..." << std::endl;
				int numElementsInDDG = rxns.at(i).numElementsInDDG;
				for(unsigned int j = 0, pos = rxns.at(i).ddgIndex; j < numElementsInDDG; ++pos, j++) {
					//cout << "\t\t" << "Calculating Position in DDG " << i << " ..." << std::endl;
					int position = DoubletDependencyGraph.at(pos);
					//cout << "\t\t\t" << "Calculating Difference in DDG ..." << std::endl;
					float diff = chems.at(rxns.at(position).r1) * (chems.at(rxns.at(position).r1) - 1) * rxns.at(position).rate - propensities.at(position); // divide by two was calculated in previously
					rxns.at(position).sumPropensity += (simTime - rxns.at(position).timeSinceLastUpdate) * propensities.at(position);
					rxns.at(position).timeSinceLastUpdate = simTime;
					propensities.at(position) += diff;
					a0 += diff;
				}

				if(numIterationsCompleted % 10000 == 0)
				{
					if(numIterationsCompleted % 10000)
						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;
	
	//std::sort(SortedRxnMapping.begin(), SortedRxnMapping.end(), RxnComparison);
	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[7]);
	
	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[8]);
	
	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[9]);
	
	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;
}

