// NextReactionMethod.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"
#include <utility>
#include <limits>

using namespace std;

struct Reaction {
	int r1, r2, p1, p2;
	int posInHeap;
	double rate;
	double propensity;
	double previousPropensity;
	double tau;
	double previousTau;
	double previousNonZeroTime;
};

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

vector< pair<int, vector<int> > > dependencyGraph;

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...";
	
	dependencyGraph.reserve(2 * numReactions);
	
	for(int i = 0; i < numReactions; i++) {
		int numRxnsInSDG = -1, numRxnsInDDG = -1;
		
		sdgFile >> numRxnsInSDG;
		ddgFile >> numRxnsInDDG;
		
		vector<int> a, b;
		pair<int, vector<int> > p1(numRxnsInSDG, a);
		pair<int, vector<int> > p2(numRxnsInDDG, b);

		dependencyGraph.push_back(p1);
		
		for(int j = 0; j < numRxnsInSDG; j++) {
			int t;
			sdgFile >> t;
			
			if(t != i)
				dependencyGraph.back().second.push_back(t);
			if(t >= numReactions) {
				cout << "BAD SDG Detected!!!" << std::endl;
				exit(0);
			}
		}
		
		dependencyGraph.back().first = dependencyGraph.back().second.size();

		dependencyGraph.push_back(p2);
		
		for(int j = 0; j < numRxnsInDDG; j++) {
			int t;
			ddgFile >> t;
			if(t != i)
				dependencyGraph.back().second.push_back(t);
			if(t >= numReactions) {
				cout << "BAD DDG Detected!!!" << std::endl;
				exit(0);
			}
		}
		
		dependencyGraph.back().first = dependencyGraph.back().second.size();
	}
	
	cout << " FINISHED" << std::endl << std::endl;
}

inline void SwapElementsInHeap(int heap[], int n1, int n2) {
	int temp = heap[n1];
	
	heap[n1] = heap[n2];
	heap[n2] = temp;
	
	rxns[heap[n1]].posInHeap = n1;
	rxns[heap[n2]].posInHeap = n2;
}

// Functions that manipulate the heap
void CreateMinHeap(int heap[], int node) {
	if(node < 1) {
		cout << "CreateMinHeap Failed... node < 1" << std::endl;
		exit(0);
	}
	
	if(2 * node <= numReactions) {
		CreateMinHeap(heap, 2 * node);
	
		if(rxns[heap[node]].tau > rxns[heap[2 * node]].tau)
			SwapElementsInHeap(heap, node, 2 * node);
		
		if(2 * node + 1 <= numReactions) {
			CreateMinHeap(heap, 2 * node + 1);

			if(rxns[heap[node]].tau > rxns[heap[2 * node + 1]].tau)
				SwapElementsInHeap(heap, node, 2 * node + 1);
		}
	}
}

// Written iteratively for slight speed boost. Sorry no tail-recursion in C!
void MoveDownElementInHeap(int heap[], int node) {
	// Assume the node being passed in is valid
	while(true) {
		int maxChild = -1;
	
		if(node * 2 <= numReactions)
			maxChild = node * 2;
		
		if(2 * node + 1 <= numReactions && rxns[heap[2 * node + 1]].tau > rxns[heap[2 * node]].tau)
			maxChild = 2 * node + 1;
	
		if(maxChild != -1 && rxns[heap[node]].tau > rxns[heap[maxChild]].tau) {
			SwapElementsInHeap(heap, node, maxChild);
			node = maxChild;
		} else
			return;
	}
}

void MoveUpElementInHeap(int heap[], int node) {
	while(true) {
		if(node > 1 && rxns[heap[node]].tau < rxns[heap[node / 2]].tau) {
			SwapElementsInHeap(heap, node, node / 2);
			node /= 2;
		}
		else
			return;
	}
}

void SiftElementInHeap(int heap[], int node) {
	// See if we want to sift up in the heap first
	MoveUpElementInHeap(heap, node);
	
	// See if we want to instead sift down in the heap first
	if(node <= numReactions)
		MoveDownElementInHeap(heap, node);
}

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

	// 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)
			rxns[i].propensity = chems[rxns[i].r1] * chems[rxns[i].r2] * rxns[i].rate;
		else
			rxns[i].propensity = chems[rxns[i].r1] * (chems[rxns[i].r1] - 1) * rxns[i].rate;

		rxns[i].tau = (1.0 / rxns[i].propensity) * log(1.0 / dsfmt_genrand_close_open(&dsfmt));
	}
	
	int *heap = new int[numReactions + 1];

	heap[0] = -22;

	for(int i = 1; i <= numReactions; ++i) {
		heap[i] = i - 1;
		rxns[i - 1].posInHeap = i;
	}
	
	CreateMinHeap(heap, 1);
	
	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;
		}*/

		if(rxns[heap[1]].tau <= 0.0)
			break;
		
		double oldSimTime = simTime;
		simTime += rxns[heap[1]].tau;
		
		// Update chemical species
		chems[rxns[heap[1]].r1] -= 1;
		chems[rxns[heap[1]].r2] -= 1;
		chems[rxns[heap[1]].p1] += 1;
		chems[rxns[heap[1]].p2] += 1;
		
		// Reset the Null Chem
		chems[0] = 1;
		
		// Compute the propensity of the next reaction
		if(rxns[heap[1]].r1 == chems[rxns[heap[1]].r2])
			rxns[heap[1]].propensity = chems[rxns[heap[1]].r1] * (chems[rxns[heap[1]].r1] - 1) * rxns[heap[1]].rate;
		else
			rxns[heap[1]].propensity = chems[rxns[heap[1]].r1] * chems[rxns[heap[1]].r2] * rxns[heap[1]].rate;
		
		// Update the tau of the next reaction
		if(rxns[heap[1]].propensity <= 0.0)
			rxns[heap[1]].tau = numeric_limits<float>::max();
		else
			rxns[heap[1]].tau = (1.0 / rxns[heap[1]].propensity) * log(1.0 / dsfmt_genrand_close_open(&dsfmt));
		
		// See if the next reaction needs to be sifted down the heap
		MoveDownElementInHeap(heap, 1);
		
		// Now recalculate all the appropriate SINGLET propensities
		int numElementsInSDG = dependencyGraph[2 * heap[1]].first;
		
		for(vector<int>::iterator iter = dependencyGraph[2 * heap[1]].second.begin(); iter != dependencyGraph[2 * heap[1]].second.end(); iter++) {
			int position = *iter;
			double newPropensity = chems[rxns[position].r1] * chems[rxns[position].r2] * rxns[position].rate;
			if(rxns[position].propensity <= 0.0) {
				if(newPropensity > 0.0) {
					rxns[position].tau = rxns[position].propensity / newPropensity * (rxns[position].previousTau - simTime) + simTime;
					rxns[position].propensity = newPropensity;
				}
			} else {
				if(newPropensity <= 0.0) {
					rxns[position].previousTau = rxns[position].tau;
					rxns[position].tau = numeric_limits<float>::max();
					rxns[position].previousNonZeroTime = oldSimTime;
					rxns[position].previousPropensity = rxns[position].propensity;
					rxns[position].propensity = 0.0;
				} else {
					rxns[position].tau = rxns[position].propensity / newPropensity * (rxns[position].tau - simTime) + simTime;
					rxns[position].propensity = newPropensity;
				}
			}

			SiftElementInHeap(heap, rxns[position].posInHeap);
		}

		// Now recalculate all the appropriate DOUBLET propensities
		int numElementsInDDG = dependencyGraph[2 * heap[1] + 1].first;

		for(vector<int>::iterator iter = dependencyGraph[2 * heap[1] + 1].second.begin(); iter != dependencyGraph[2 * heap[1] + 1].second.end(); iter++) {
			int position = *iter;
			//cout << "\t\t\t" << "Calculating Difference in DDG ..." << std::endl;
			double newPropensity = chems[rxns[position].r1] * (chems[rxns[position].r1] - 1) * rxns[position].rate; // divide by two was calculated in previously
			
			if(rxns[position].propensity <= 0.0) {
				if(newPropensity > 0.0) {
					rxns[position].tau = rxns[position].propensity / newPropensity * (rxns[position].previousTau - simTime) + simTime;
					rxns[position].propensity = newPropensity;
				}
			} else {
				if(newPropensity <= 0.0) {
					rxns[position].previousTau = rxns[position].tau;
					rxns[position].tau = numeric_limits<float>::max();
					rxns[position].previousNonZeroTime = oldSimTime;
					rxns[position].previousPropensity = rxns[position].propensity;
					rxns[position].propensity = 0.0;
				} else {
					rxns[position].tau = rxns[*iter].propensity / newPropensity * (rxns[position].tau - simTime) + simTime;
					rxns[position].propensity = newPropensity;
				}
			}
	
			SiftElementInHeap(heap, rxns[position].posInHeap);
		}
	}
	
	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;

	ofstream outputTime(argv[9]);
	outputTime << (EndTime - StartTime) << " seconds." << std::endl;
	
	cout << std::endl;
	cout << "Press any key to exit simulation.";

	//cin >> numIterations;
	
	return 0;
}

