#include "Intel_Compatibility.h"
#include <iostream>
#include <vector>
#include <sbml/SBMLTypes.h>
#include <string>
#include <map>
#include <list>
#include <fstream>

using namespace std;

struct SagarReaction
{
	int r1, r2, p1, p2;
	double rxnRate;
};

struct SagarParameter
{
	string name;
	bool evaluated;
	const ASTNode* expression;
	double value;
};

map<string, int> speciesNamesToID;
map<string, SagarParameter> paramNameToStruct;


// Function Prototypes
double getSagarValue(string);

double SagarProcessNode(const ASTNode* nd)
{
	if(nd->isNumber())
		return nd->getReal();
	
	if(nd->isName())
		return getSagarValue(nd->getName());
	
	if(nd->isOperator())
	{
		int numChildren = nd->getNumChildren();
		double sum = 0.0;
		double product = 1.0;
		
		switch(nd->getCharacter())
		{
		case '+':
			for(int i = 0; i < numChildren; i++)
				sum += SagarProcessNode(nd->getChild(i));
			return sum;
			break;
		case '-':
			return SagarProcessNode(nd->getLeftChild()) - SagarProcessNode(nd->getRightChild());
			break;
		case '*':
			for(int i = 0; i < numChildren; i++)
				product *= SagarProcessNode(nd->getChild(i));
			return product;
			break;
		case '/':
			return SagarProcessNode(nd->getLeftChild()) / SagarProcessNode(nd->getRightChild());
			break;
		default:
			cout << "ERROR: SINGLE OPERATOR CHARACTER UNDEFINED = " << nd->getCharacter() << std::endl;
			exit(0);
		}
	}

	cout << "ERROR: NonTrivial parsing Required = " << std::endl;
	exit(0);
}

double getSagarValue(string strParamName)
{	
	//cout << "Processing parameter: " << strParamName << "	Evaluated = ";
	
	if(paramNameToStruct.find(strParamName) == paramNameToStruct.end())
	{
			//cout << "paramNameToStruct query invalid: " << strParamName << std::endl;
			return 0.0;
	}
	
	SagarParameter &sp = paramNameToStruct[strParamName];
	
	if(sp.evaluated == true)
	{
		//cout << "True" << std::endl;
		return sp.value;
	}
	else {
		sp.value = SagarProcessNode(sp.expression);
		sp.evaluated = true;
		sp.expression = NULL;
		/*if(!sp.expression)
		{
			cout << "LOLCATS!!!" << std::endl;
			sp.value = 0.0;
		}
		else*/
			
		
		//cout << "False" << std::endl;
		return sp.value;
	}
}

void ProcessParameters(Model* theModel)
{
	// First read through all the basic parameter definitions,
	// and make assignments to those that are constant
	int numParameters = theModel->getNumParameters();
	
	cout << "Making first pass on parameters" << std::endl;
	
	vector<string> parameterNamesList;
	
	for(int i = 0; i < numParameters; i++)
	{
		string paramName = theModel->getParameter(i)->getId();
		paramNameToStruct[paramName];
		
		paramNameToStruct[paramName].name = paramName;
		
		if(theModel->getParameter(i)->getConstant()) {
			paramNameToStruct[paramName].evaluated = true;
			paramNameToStruct[paramName].value = theModel->getParameter(i)->getValue();
		} else
			paramNameToStruct[paramName].evaluated = false;
		
		paramNameToStruct[paramName].expression = NULL;
		//cout << paramNameToStruct[paramName].name << std::endl;
		parameterNamesList.push_back(paramNameToStruct[paramName].name);
	}

	// Now go through all the AssignmentRules and assign the root nodes
	// to the parameters.
	int numRules = theModel->getNumRules();
	list<string> paramsToProcess;
	
	cout << "Making second pass on parameters" << std::endl;
	
	for(int i = 0; i < numRules; i++)
	{
		AssignmentRule* ar = (AssignmentRule*)(theModel->getRule(i));

		if(!ar)
		{
			cout << "EPIC FAIL!" << std::endl;
			exit(0);
		}
		
		if(paramNameToStruct[ar->getVariable()].evaluated == false)
		{
			if(!(ar->getMath()))
			{
				cout << "FAILURE: ar->getMath() was not valid" << std::endl;
				exit(0);
			}
			
			if(paramNameToStruct.find(ar->getVariable()) == paramNameToStruct.end())
			{
				cout << "FAILURE: paramNameToStruct doesn't have " << ar->getVariable() << std::endl;
				exit(0);
			}
			
			paramNameToStruct[ar->getVariable()].expression = ar->getMath();
			paramsToProcess.push_back(ar->getVariable());
			//cout << paramsToProcess.back() << std::endl;
		}
		else
		{
			cout << "ERROR 2435234: EXITING NOW" << std::endl;
			exit(0);
		}
	}

	cout << "Making final pass on parameters" << std::endl;
	
	// Now we make a pass through each param to process
	while(paramsToProcess.empty() == false)
	{
		getSagarValue(paramsToProcess.front());
		//cout << "About the process the following variable according to assignment rule: " << paramsToProcess.front() << std::endl;
		//cout << "Processed Sagar Value: " << paramsToProcess.front() << " = " << getSagarValue(paramsToProcess.front());
		paramsToProcess.pop_front();
	}

	//for(int i = 0; i < (signed)parameterNamesList.size(); i++)
	//	cout << paramNameToStruct[parameterNamesList[i]].name << " = " << paramNameToStruct[parameterNamesList[i]].value << std::endl;

}

void ProcessModel(Model* theModel, string strOutputReactionFileName, string strOutputSpeciesFileName)
{
	int numReactions = theModel->getNumReactions();
	int numSpecies = theModel->getNumSpecies();
	
	cout << "Number Of Species = " << numSpecies << std::endl;
	cout << "Number Of Reactions = " << numReactions << std::endl;
	
	vector<int> initConcentrations;
	
	ofstream outputRxnFile(strOutputReactionFileName.c_str());
	ofstream outputSpeciesFile(strOutputSpeciesFileName.c_str());
	
	outputSpeciesFile << numSpecies + 1 << std::endl << std::endl;
	
	initConcentrations.push_back(1); // for "NULL"
	speciesNamesToID["NULL"] = 0;
	outputSpeciesFile << speciesNamesToID["NULL"] << "	" << "NULL	" << initConcentrations[0] << std::endl;
	
	cout << "Reading in species data" << std::endl;
	
	// Read in species data
	for(int i = 0; i < numSpecies; i++)
	{
		Species* z = theModel->getSpecies(i);
		speciesNamesToID[z->getName()] = i + 1;
		initConcentrations.push_back((int)z->getInitialConcentration());
		//cout << i << std::endl;

		outputSpeciesFile << speciesNamesToID[z->getName()] << "	" << z->getName() << "	" << initConcentrations.back() << std::endl;
	}
	
	// Read in Parameter Data
	cout << "Reading in parameter data" << std::endl;
	
	ProcessParameters(theModel);
	
	cout << "Beginning to process reactions" << std::endl;
	
	outputRxnFile << numReactions << std::endl << std::endl;
	
	// Now process reactions
	for(int i = 0; i < numReactions; i++)
	{		
		SagarReaction camel;
		
		Reaction* r = theModel->getReaction(i);
		
		// Process the reactants
		switch(r->getNumReactants())
		{
		case 0:
			// Shouldn't happen but ok...
			//cout << "ZERO CASE HAPPENING" << std::endl;
			camel.r1 = camel.r2 = speciesNamesToID["NULL"];
			break;
		case 1:
			//cout << "ONE CASE HAPPENING" << std::endl;
			if(r->getReactant(0)->getStoichiometry() > 1.0) 
			{
				// Assume it's 2 then
				camel.r1 = camel.r2 = speciesNamesToID[theModel->getSpecies(r->getReactant(0)->getSpecies())->getName()];
			} 
			else 
			{	
				camel.r1 = speciesNamesToID[theModel->getSpecies(r->getReactant(0)->getSpecies())->getName()];
				camel.r2 = speciesNamesToID["NULL"];
			}
			
			//cout << theModel->getSpecies(r->getReactant(0)->getSpecies())->getName() << std::endl;
			break;
		case 2:
			//cout << "TWO CASE HAPPENING" << std::endl;
			camel.r1 = speciesNamesToID[theModel->getSpecies(r->getReactant(0)->getSpecies())->getName()];
			camel.r2 = speciesNamesToID[theModel->getSpecies(r->getReactant(1)->getSpecies())->getName()];
			//cout << theModel->getSpecies(r->getReactant(0)->getSpecies())->getName() << std::endl;
			//cout << theModel->getSpecies(r->getReactant(1)->getSpecies())->getName() << std::endl;
			break;
		default:
			cout << "ERROR !!! TOO MANY REACTANTS!!!" << std::endl;
			exit(0);
			break;
		}
		
		//cout << "----" << std::endl;
		
		// Process the products
		switch(r->getNumProducts())
		{
		case 0:
			// Shouldn't happen but ok...
			//cout << "ZERO CASE HAPPENING" << std::endl;
			camel.p1 = camel.p2 = speciesNamesToID["NULL"];
			break;
		case 1:
			//cout << "ONE CASE HAPPENING" << std::endl;
			if(r->getProduct(0)->getStoichiometry() > 1.0) 
			{
				// Assume it's 2 then
				camel.p1 = camel.p2 = speciesNamesToID[theModel->getSpecies(r->getProduct(0)->getSpecies())->getName()];
			} 
			else 
			{	
				camel.p1 = speciesNamesToID[theModel->getSpecies(r->getProduct(0)->getSpecies())->getName()];
				camel.p2 = speciesNamesToID["NULL"];
			}
			
			//cout << theModel->getSpecies(r->getProduct(0)->getSpecies())->getName() << std::endl;
			break;
		case 2:
			//cout << "TWO CASE HAPPENING" << std::endl;
			camel.p1 = speciesNamesToID[theModel->getSpecies(r->getProduct(0)->getSpecies())->getName()];
			camel.p2 = speciesNamesToID[theModel->getSpecies(r->getProduct(1)->getSpecies())->getName()];
			break;
		default:
			cout << "ERROR !!! TOO MANY PRODUCTS!!!" << std::endl;
			exit(0);
			break;
		}
		
		// Now we need to get the rate coefficient
		//cout << "computing rate	" << r->getKineticLaw()->getNumParameters() << std::endl;
		
		// lol is the child-ish math node, and pol is the parent of that childish math node.
		const ASTNode* lol = r->getKineticLaw()->getMath();
		const ASTNode* pol = lol;
		camel.rxnRate = 1.0;
		
		if(lol->isOperator() == false)
		{
			if(lol->getChild(i)->isReal())
				camel.rxnRate = lol->getChild(i)->getReal();
			else
				camel.rxnRate = paramNameToStruct[lol->getName()].value;
		}
		else
		{
			while(lol->isOperator())
			{
				pol = lol;
				lol = lol->getChild(0);
			}

			// Now go through pol from left to right multiplying everything until we get to a species
	
			for(int i = 0; i < (signed)pol->getNumChildren(); i++)
			{
				if(pol->getChild(i)->isReal())
				{
					camel.rxnRate *= pol->getChild(i)->getReal();
					continue;
				}
		
				if(theModel->getSpecies(pol->getChild(i)->getName()))
					break;
				else
					camel.rxnRate *= paramNameToStruct[pol->getChild(i)->getName()].value;
			}
		}
		
		//cout << theModel->get
		//cout << lol->getName() << " = " << camel.rxnRate << std::endl;//<< camel.rxnRate << std::endl;
		//cout << camel.r1 << " + " << camel.r2 << " --> " << camel.p1 << " + " << camel.p2 << ",  k = " << camel.rxnRate << std::endl;
		outputRxnFile << camel.r1 << "	" << camel.r2 << "	" << camel.p1 << "	" << camel.p2 << "	" << camel.rxnRate << std::endl;
	}
}

/* PARAMETERS ARE:
	1) XML INPUT FILE
	2) OUTPUT REACTION FILE NAME
	3) OUTPUT SPECIES FILE NAME
*/

int main (int argc, char* argv[])
{
	SBMLDocument* document = readSBML(argv[1]);

	unsigned int errors = document->getNumErrors();

	cout << endl;
	cout << "  filename: " << argv[1] << endl;
	cout << "  error(s): " << errors  << endl;
	cout << endl;

	string strOutputReactionFileName(argv[2]);
	string strOutputSpeciesFileName(argv[3]);
	
	cout << "OutputReactionFile name = " << strOutputReactionFileName << std::endl;
	cout << "OutputSpeciesFile name = " << strOutputSpeciesFileName << std::endl;
	
	if (errors > 0) 
		document->printErrors(cerr);

	if(!document)
		cout << "TOTALLY SCREWED DATA." << std::endl;
	else
	{
		cout << "STILL A CHANCE" << std::endl;
		ProcessModel(document->getModel(), strOutputReactionFileName, strOutputSpeciesFileName);
	}
	
 	return errors;
}
