// OptimizedDirectMethod.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 <sys/time.h>
#include "dSFMT/dSFMT.h"

using namespace std;

struct Reaction
{
        int r1, r2, p1, p2;
        int sdgIndex, numElementsInSDG;
		int ddgIndex, numElementsInDDG;
        double rate;
};

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

int numReactions;
int numChemicals;
int sizeOfDependencyGraph;

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++)
	{
		int blah, chemVal;
		string blahblah;
        speciesFile >> blah >> blahblah >> chemVal;
		chems.push_back(chemVal);
	}

	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].rate;
		if(rxns[i].r1 == rxns[i].r2)
			rxns[i].rate /= 2.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;
}

/*
	Arguments = 
	1) numIterations
	2) JobName
	3) ReactionFileName
	4) SpeciesFileName
	5) SingletDependencyGraphFileName
	6) DoubletDependencyGraphFileName
	7) randNumSeed
	8) maxSimTime
*/
int main(int argc, char* argv[])
{
	cout << "STOCHASTIC SIMULATION ALGORITHM - OPTIMIZED DIRECT METHOD" << 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]);
	string JobName(argv[2]);

	double simTime = 0.0, a0 = 0.0;
	double maxSimTime = atoi(argv[8]);
	
	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, atoi(argv[7]));

	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 << "numIterations = " << numIterations << std::endl << 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!!!
	
	timeval timVal;
	gettimeofday(&timVal, NULL);
	double StartTime = timVal.tv_sec + (timVal.tv_usec/1000000.0);

	// Do simulation here
	for(; --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[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[rxns[i].r1] -= 1;
				chems[rxns[i].r2] -= 1;
				chems[rxns[i].p1] += 1;
				chems[rxns[i].p2] += 1;
				
				// Reset the Null Chem
				chems[0] = 1;

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

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

		// Now see if a0 is > 0
		if(a0 <= 0.0)
			break;
	}
	
	gettimeofday(&timVal, NULL);
	double EndTime = timVal.tv_sec + (timVal.tv_usec/1000000.0);
	
	cout << "Time Taken = " << EndTime - StartTime << " seconds." << std::endl;

	cout << std::endl << "SIMULATION FINISHED!" << std::endl;
	cout << "\tNumIterations remaining = " << numIterations << std::endl;

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

	//cin >> numIterations;
	
	return 0;
}

