#include "Intel_Compatibility.h"
#include <fstream>
#include <iostream>
#include <vector>
#include <list>
#include <map>
#include <set>

using namespace std;

struct Reaction
{
	int r1, r2, p1, p2;
	double rate;
	double normalizedSumProp;
	int numAffectedRxns;
};

void ProcessReaction(Reaction &r, int z, list<int> &affectedChemicals) {
	if(z == 0)
		return;
	
	int one = 0;
	if(r.r1 == z) one--;
	if(r.r2 == z) one--;
	if(r.p1 == z) one++;
	if(r.p2 == z) one++;

	if(one != 0) {
		for(list<int>::iterator iter = affectedChemicals.begin(); iter != affectedChemicals.end(); iter++)
			if(*iter == z)
				return;

		affectedChemicals.push_back(z);
	}
}

// Compute Affected Chemicals for a Reaction
void ComputeAffectedChemicalsForReaction(Reaction &r, list<int> &affectedChemicals) {
	affectedChemicals.clear();
	ProcessReaction(r, r.r1, affectedChemicals);
	ProcessReaction(r, r.r2, affectedChemicals);
	ProcessReaction(r, r.p1, affectedChemicals);
	ProcessReaction(r, r.p2, affectedChemicals);
}

/* ARGUMENTS:
	1) Sorted Reaction file
	2) Sum Prop file
	3) Output File Name
*/
int main(int argc, char* argv[])
{
	cout << "@@@ CAEK PROGRAM - Sagar Indurkhya (2008) &&&" << std::endl;
	cout << "Sorted Reaction File = " << argv[1] << std::endl;
	cout << "Sum Prop File = " << argv[2] << std::endl;
	cout << "Output File Name = " << argv[3] << std::endl;
	
	int numReactions;
	vector<Reaction> rxns;
	
	ifstream rxnInput(argv[1]);
	ifstream sumPropData(argv[2]);
	
	rxnInput >> numReactions;
	
	double sumOfPropensities = 0.0;
	
	for(int i = 0; i < numReactions; i++) {
		Reaction r;
		int temp;
		rxnInput >> r.r1 >> r.r2 >> r.p1 >> r.p2 >> r.rate;
		r.numAffectedRxns = 0;
		sumPropData >> temp >> r.normalizedSumProp;
	
		if(r.normalizedSumProp < 0.0)
		{
			cout << "ERROR: r.normalizedSumProp < 0.0. r.normalizedSumProp = " << r.normalizedSumProp << std::endl;
			return -1;
		}
		
		sumOfPropensities += r.normalizedSumProp;
		
		rxns.push_back(r);
	}

	for(int i = 0; i < numReactions; i++)
		rxns[i].normalizedSumProp /= sumOfPropensities;
	
	map<int, list<int> > chemToRxn;
	
	for(int i = 0; i < numReactions; i++) {
		if(rxns[i].r1 != 0)
			chemToRxn[rxns[i].r1].push_back(i);
		if(rxns[i].r2 != 0 && rxns[i].r2 != rxns[i].r1)
			chemToRxn[rxns[i].r2].push_back(i);
	}
	
	int highestDegree = 0;
	
	// Now for each reaction
	for(int i = 0; i < numReactions; i++) {
		list<int> affectedChemicals;
		ComputeAffectedChemicalsForReaction(rxns[i], affectedChemicals);
		
		set<int> affectedRxns;
		
		for(list<int>::iterator iter = affectedChemicals.begin(); iter != affectedChemicals.end(); iter++) {
			for(list<int>::iterator jiter = chemToRxn[*iter].begin(); jiter != chemToRxn[*iter].end(); jiter++)
				affectedRxns.insert(*jiter);
		}
		
		rxns[i].numAffectedRxns = affectedRxns.size();
		
		if(rxns[i].numAffectedRxns > highestDegree)
			highestDegree = rxns[i].numAffectedRxns;

		//cout << "Finished Reaction " << i << std::endl;
	}

	cout << "highest Degree = " << highestDegree << std::endl;
	
	double *graphToOutput = new double[highestDegree + 1];
	
	for(int i = 0; i < (highestDegree + 1); i++)
		graphToOutput[i] = 0.0;
	
	// Now add the sumProp of each reaction
	for(int i = 0; i < numReactions; i++)
		graphToOutput[rxns[i].numAffectedRxns] += rxns[i].normalizedSumProp;
	
	ofstream output(argv[3]);
	
	double runningSum = 0.0;
	
	// Finally output this to a file
	for(int i = 0; i < (highestDegree + 1); i++) {
		if(graphToOutput[i] < 0.0) {
			cout << "ERROR: graphToOutput[i] < 0.0, i = " << i << ", graphToOutput[i] = " << graphToOutput[i] << std::endl;
			return -1;
		}
		
		runningSum += graphToOutput[i];
	
		output << i << "	" << graphToOutput[i] << "	" << runningSum << std::endl;
	}
	
	return 0;
}
