#include "Intel_Compatibility.h"
#include <fstream>
#include <stdlib.h>
#include "dSFMT/dSFMT.h"
#include <math.h>
#include <iostream>

using namespace std;

/*  The purpose of this file is to serve as an extremely simple but reliable
	stochastic simulator.
*/

struct Reaction
{
	int r1, r2, p1, p2;
	bool isDoublet;
	double rate;
	double propensity;
	double sumPropensity;
};

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

int sort_function(const void *a, const void *b)
{
	Reaction *rxnA = (Reaction*)a;
	Reaction *rxnB = (Reaction*)b;
	
	if(rxnA->sumPropensity < rxnB->sumPropensity)
		return 1;
	if(rxnA->sumPropensity > rxnB->sumPropensity)
		return -1;
	return 0;
}

/* PARAMETERS ARE:
	1) INPUT REACTION FILE NAME
	2) INPUT SPECIES FILE NAME
	3) NUM ITERATIONS
	4) OUTPUT FILE NAME FOR CHEMICALS AND SORTED REACTIONS
	5) OUTPUT FILE NAME FOR SUM OF PROPENSITIES
	6) OUTPUT FILE NAME FOR CHEMICAL CONCENTRATIONS
*/
int main(int argc, char* argv[])
{
	dsfmt_t dsfmt;
	dsfmt_init_gen_rand(&dsfmt, 10235);
	 
	// BEGIN INPUT
	ifstream inputRxnFile(argv[1]);
	ifstream inputSpeciesFile(argv[2]);
	
	cout << "inputRxnFile = " << string(argv[1]) << std::endl;
	cout << "inputSpeciesFile = " << string(argv[2]) << std::endl;
	
	ofstream SteadyStateData("SteadyStateData.dat");
	
	long numIterations = atoi(argv[3]);
	int numChemicals, numReactions;
	
	inputSpeciesFile >> numChemicals;
	
	cout << "numChemicals = " << numChemicals;

	cout << "numChemicals = " << numChemicals << std::endl;
	
	// chems actually holds the concentrations during runtime
	int *chems = new int[numChemicals];

	// species just holds the standard data
	Chemical* species = new Chemical[numChemicals];
	
	for(int i = 0; i < numChemicals; i++)
	{
		inputSpeciesFile >> species[i].ID >> species[i].Name >> species[i].initNumber;
		chems[i] = species[i].initNumber;
	}

	inputRxnFile >> numReactions;
	
	Reaction *rxns = new Reaction[numReactions];

	cout << "numReactions = " << numReactions;
	
	for(int i = 0; i < numReactions; i++)
	{
		inputRxnFile >> rxns[i].r1 >> rxns[i].r2 >> rxns[i].p1 >> rxns[i].p2 >> rxns[i].rate;

		if(rxns[i].r1 == rxns[i].r2)
		{
			rxns[i].isDoublet = true;
			rxns[i].propensity = chems[rxns[i].r1] * (chems[rxns[i].r1] - 1) * rxns[i].rate / 2.0;
		} else {
			rxns[i].isDoublet = false;
			rxns[i].propensity = chems[rxns[i].r1] * chems[rxns[i].r2] * rxns[i].rate;
		}
		
		rxns[i].sumPropensity = 0.0;
	}

	// DONE WITH INPUT, BEGIN SIMULATION

	double a0 = 0.0;

	for(int j = 0; j < numReactions; j++)
	{
		if(rxns[j].isDoublet)
			rxns[j].propensity = chems[rxns[j].r1] * (chems[rxns[j].r1] - 1) * rxns[j].rate / 2.0;
		else
			rxns[j].propensity = chems[rxns[j].r1] * chems[rxns[j].r2] * rxns[j].rate;
	
		a0 += rxns[j].propensity;
	}
	
	cout << std::endl << "INITIAL SUM PROPENSITY = " << a0 << std::endl << std::endl;
	
	double simTime = 0.0;
	
	for(long i = 0; i < numIterations; i++)
	{
		a0 = 0.0;
		
		for(int j = 0; j < numReactions; j++)
		{
			if(rxns[j].isDoublet)
				rxns[j].propensity = chems[rxns[j].r1] * (chems[rxns[j].r1] - 1) * rxns[j].rate / 2.0;
			else
				rxns[j].propensity = chems[rxns[j].r1] * chems[rxns[j].r2] * rxns[j].rate;
			
			a0 += rxns[j].propensity;
			rxns[j].sumPropensity += rxns[j].propensity;
		}
		
		double randNum = a0 * dsfmt_genrand_close_open(&dsfmt);
		double runningSum = 0.0;
		
		int rxnToExecute = -1;
		
		// Now find the reaction we want
		for(int j = 0; j < numReactions; j++)
		{
			if(runningSum + rxns[j].propensity >= randNum)
			{
				rxnToExecute = j;
				break;
			}
			
			runningSum += rxns[j].propensity;
		}

		// if we didn't pick a reaction whatever, just do it again
		if(rxnToExecute == -1)
			continue;
		else
			if(i % 10000 == 0)
			{
				if(i % 1000000)
					cout << "Iteration: " << i << "	Executed reaction #: " << rxnToExecute << "	" << rxns[rxnToExecute].sumPropensity << std::endl;
				SteadyStateData << simTime << "	" << chems[197] << std::endl;
			}
		
		// Now execute the correct reaction
		chems[rxns[rxnToExecute].r1] -= 1;
		chems[rxns[rxnToExecute].r2] -= 1;
		chems[rxns[rxnToExecute].p1] += 1;
		chems[rxns[rxnToExecute].p2] += 1;
		simTime += -1.0/a0 * log(dsfmt_genrand_close_open(&dsfmt));
		
	}

	// Normalize the sum propensities!
	double sumOfSumPropensities = 0.0;
	
	for(int i = 0; i < numReactions; i++)
		sumOfSumPropensities += rxns[i].sumPropensity;
	
	cout << sumOfSumPropensities << std::endl;
	
	for(int i = 0; i < numReactions; i++)
		rxns[i].sumPropensity /= sumOfSumPropensities;
	
	// Sort the reactions based on the sum of propensties (from highest to lowest)
	qsort((void*)rxns, numReactions, sizeof(rxns[0]), sort_function);
	
	// Now lets do our outputting
	
	// First lets output our standard file with chemicals and reactions
	ofstream basicOutput(argv[4]);
	
	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].rate << std::endl;
	
	// Output a file containing the sum of propensities data in a format that is graphable
	ofstream sumPropsOutput(argv[5]);
	
	for(int i = 0; i < numReactions; i++)
		sumPropsOutput << i + 0.5 << "	" << rxns[i].sumPropensity << std::endl;
	
	// Now output the current chemical concentrations
	ofstream chemConcs(argv[6]);
	
	chemConcs << numChemicals << std::endl << std::endl;
	
	for(int i = 0; i < numChemicals; i++)
		chemConcs << species[i].ID << "	" << species[i].Name << "	" << chems[i] << std::endl;

	return 0;
}
