#include "Simulation.h"
#include <iostream>
#include <fstream>

using namespace std;
// Initialize Static Variables
bool Simulation::isInstantiated = false;
Simulation * Simulation::sim;


//////////////////////////////////////////////////////////////////////
// Get the instance of the Cell Simulation
// ->	Only allow one instance of the cell to exist at any one time
// 		If not yet exist, create a new Cell Simulation
// 		otherwise return the pointer to the current existing Cell Simulation
//
//Simulation * Simulation::getSimulation(){}


//////////////////////////////////////////////////////////////////////
// Constructor for the Simulation
// ->	Create and initialize properties of the Simulation space
Simulation:: Simulation(){

	// Cell Simulation Environment Manager
	lEnv = LatticeEnvironment::getLatticeEnv();

	// Build the Interaction table
	buildInteractionNetworkTable();

	// Initialize status before simulation
	simulationCompleted = 0;
	numCycles = 0;
	cutSitesActive = 0;
}



//////////////////////////////////////////////////////////////////////
// Destructor for the Simulation
//
Simulation:: ~Simulation(){}



//////////////////////////////////////////////////////////////////////
// Add Molecule to the Simulation
// ->	Add the molecule to the Cell Simulation Space
//
void Simulation:: add(Molecule * m){
	// add one to the number of molecules in the Simulation
	numMolecules ++;
	
	// add the molecule to the molecule list within the Simulation
	// Find the appriate molecule type and add it to its list

	molTable[m->getType()].push_back(m);
}



//////////////////////////////////////////////////////////////////////
//  Remove the Molecule from the Simulation
//  ->	Remove the molecule from the Simulation Space
//
void Simulation :: remove(Molecule *m){

	// Get the molecule type
	moltype type = m->getType();
	
	// Number of elements of the type 'type'
	int numElements = molTable[m->getType()].size();
	
	// temporary vector to store non-removing molecules
	vector <Molecule *> temp;
	
	for(int i = 0 ; i < numElements ; i ++){
	
		// Keep non-removing molecules
		if(molTable[type][i] != m){
			temp.push_back(molTable[type][i]);
		}
	}
	
	// decrement number of molecules in the Simulation Space
	numMolecules --;
	// Update the list of existing molecules
	molTable[type] = temp;
}



//////////////////////////////////////////////////////////////////////
//  Display all Molecules in the Simulation
//  ->	Display position and status of all molecules within the Simulation Space
//
void Simulation:: display(){

	printf("\n# types: %d / # molecules: %d\n" , molTable.size(), numMolecules);
	
	map<moltype, vector <Molecule *> >::const_iterator itr = molTable.begin();

	for (itr = molTable.begin(); itr != molTable.end(); itr ++) 
	{
		for(unsigned int i = 0 ; i < (itr->second).size() ; i ++){
			Molecule * m = (itr->second)[i];
			
			printf("TYPE NAME: %s, TYPE: %d, Position: %s, Activated: %d \n",
				(m->getTypeName()).c_str(), 
				itr->first, 
				(m->position()).c_str(),
				m->isActivated());
		}

	}
}

//////////////////////////////////////////////////////////////////////
//  Display a summary of the Molecules in the Simulation
//
void Simulation:: summary(){

        map<moltype, vector <Molecule *> >::const_iterator itr = molTable.begin();
	
	int numA=0; int Aactive=0;
	int numB=0; int Bactive=0;

        for (itr = molTable.begin(); itr != molTable.end(); itr ++)
        {
                for(unsigned int i = 0 ; i < (itr->second).size() ; i ++){
                        Molecule * m = (itr->second)[i];
			if ( m->getType() == MOL_ENZYME_A ) {
				numA++;
				if( m->isActivated() ) {
					Aactive++;
				}
				//printf("A at %s\n",m->position().c_str());
				//printf("A has mobility: %f\n",m->getMobility());
			}
			if ( m->getType() == MOL_ENZYME_B ) {
                                numB++;
                                if( m->isActivated() ) {
                                        Bactive++;
                                }
				//printf("B at %s\n",m->position().c_str());
				//printf("B has mobility: %f\n",m->getMobility());
                        }
                }

        }
	printf("#A: %d (%d active), #B: %d (%d active), cuts: %d\n",numA,Aactive,numB,Bactive,cutSitesActive);
}

//////////////////////////////////////////////////////////////////////
//  Create Molecule where following are specified
//  num: number of molecules to create
//  type: the type of molecule to create
//  env: the environment to create in
//


void Simulation::createMolecule(int num, int type, int env){
	// Get list of available locations to place molecule
	vector <Point *> * list = lEnv->getAvailableSites(type);
	int numAvailableSites = list->size();

	// For each molecule to be created
	for(int molNum = 0 ; molNum < num ; molNum ++){
	
		// Randomly choose an allowable location
		int randomIndex = RandomNG::randInt(0,numAvailableSites-1);
		Point * location = ((*list)[randomIndex]);
		Molecule * m = createNewMolecule(type);
		// Add the created molecule to the molTable
		add(m);
		m->setPosition(*location);

		int env = lEnv->getEnvironment(*location);
		if(!((env!=-1) && lEnv->isCompatibleEnvironment(type,env)))
		{
		  cout << "Badly placed molecule: type " << type << " at ( " << location->x << " " << location->y << " " << location->z << " ) env " << env << endl;
		  exit (2);
		}
	}
}



//////////////////////////////////////////////////////////////////////
//  Build the molecule activation Table
//
void Simulation::buildInteractionNetworkTable(){

	// Build Interaction Network of Molecules
	string line;
	string curLine;
	string delimiter = "\t";
	
	// Open input data file
	ifstream inFile (INPUT_FILE_ACTIVATION_NETWORK);

	while (!inFile.eof())
	{	
		getline(inFile,line,'\n');

		if(line.length() > 0){
			moltype activator, activatee;
		
			int keyPos = 0;
			
			// Get Activator Molecule type
			keyPos = line.find(delimiter);	
			activator = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Get Activatee Molecule type
			keyPos = line.find(delimiter);	
			activatee = atoi(line.substr(0,keyPos).c_str());
			line = line.substr(keyPos+delimiter.length());
			
			// Build Activator <-> Activatee relationship and add to 
			// the Interaction Network Table
			interactionNetworkTable[activator].push_back(activatee);

		}
	}
	
	// Close input data file
	inFile.close();
}



//////////////////////////////////////////////////////////////////////
//  Molecules randomly move within the Simulation
//
void Simulation:: move(){

	// Step through the Molecule table
	map<moltype, vector <Molecule *> >::const_iterator itr;
	
	// Iterate through molecules of specified TYPE
	for(itr = molTable.begin() ; itr != molTable.end() ; itr ++){
		int numMol = (itr->second).size();
		// Allow molecules of the specified TYPE to move randomly
		for(int molNum = 0 ; molNum < numMol ; molNum ++){
			// Get the next molecule in the list
			Molecule * mol = (itr->second)[molNum];
			
			// Move molecule
			if(mol->isDiffusible())	mol->move();
		}
	}		
}

void Simulation :: transportAcrossMembranes() {
	#define THRESHOLD	0.0025
	#define	NUM_CYCLES_OPEN	150

	// Iterate the table of molecules that are in the lattice
	map<int, vector <MembranousMolecule *> >::const_iterator mItr = memMolTable.begin();

	for (mItr = memMolTable.begin(); mItr != memMolTable.end(); mItr++)
	{
		for (int i = 0; i < (signed)(mItr->second).size() ; i++) {
			MembranousMolecule * m = (mItr->second)[i];

			// this is largely HARDCODED, beware!

			Point p_well(m->pos.x, m->pos.y, 5.0);
			Point p_cyto(m->pos.x, m->pos.y, 3.0);

			int id_well = lEnv->getLatticeSite(p_well)->getID();
			int id_cyto = lEnv->getLatticeSite(p_cyto)->getID();

			Tuple * well = lEnv->getStateMatrix()->getRow(id_well);
			Tuple * cytosol = lEnv->getStateMatrix()->getRow(id_cyto);

			if (m->getType() == MOL_MEMBRANOUS_FIXED_GATE) {
				if (m->openedAt < 0 && cytosol->tuple[0] > THRESHOLD) {
					m->state = GATE_IS_OPEN;
					m->openedAt = numCycles;
				}

				if (m->state == GATE_IS_OPEN) {
					float add_to_well = (1.0 / 26) * cytosol->tuple[0];
					float add_to_cyto = (1.0 / 26) * well->tuple[0];

					well->tuple[0] += add_to_well - add_to_cyto;
					cytosol->tuple[0] += add_to_cyto - add_to_well;

					if (m->openedAt < numCycles - NUM_CYCLES_OPEN)
						m->state = GATE_IS_CLOSED;
				}
			}
			else {  // MOL_MEMBRANOUS_FIXED_PUMP
			  //				float ratio = cytosol->tuple[0] / well->tuple[0];
				
				float add_to_well = (cytosol->tuple[0] - 0.001);
//				if (add_to_well > 0.01) add_to_well = 0.01;

				if (add_to_well > 0.0) {
					well->tuple[0] += add_to_well;
					cytosol->tuple[0] -= add_to_well;
				}
			}
		}
	}
}


//////////////////////////////////////////////////////////////////////
//  Check for Interaction between molecules
//
void Simulation :: checkMoleculeInteractions(){

	// Create Iterator for the molecule table
	map<moltype, vector <Molecule *> >::const_iterator itr = molTable.begin();

	// Iterate through the list of molecules of differnet types
	for (itr = molTable.begin(); itr != molTable.end(); itr ++) {
		int numMol = (itr->second).size();
		for(int molIndex = 0 ; molIndex < numMol ; molIndex ++){
			Molecule * mol = (itr->second)[molIndex];
			
			// Check interactions with neigboring molecules
			interactNeighborMolecules(mol);
		}
	}
}



//////////////////////////////////////////////////////////////////////
//  Check for interating molecules
//
void Simulation::interactNeighborMolecules(Molecule * m){

	// Get the list of possible interacting molecules that are near
	vector <Molecule *> interactionNeighborList = buildNeighborList(m);

	int numInteractingMolecules = interactionNeighborList.size();

	// For each of the neigboring molecules
	for(int molIndex = 0 ; molIndex < numInteractingMolecules ; molIndex ++){

		// Interact with the molecule
		Molecule * neighbor= interactionNeighborList[molIndex];
		// Newly activated is false
		bool newAct = false;
		
		if(!neighbor->isActivated()){
			newAct = true;			
			m->activates(neighbor);

			if (neighbor->getType() == MOL_ENZYME_B)
			{
			  cutSitesActive++;
			}
		}		

		if(newAct){
			// Display critical events during the simulation
			//temp
			//printInteractionDetails(m,neighbor,1);
			
			actEventTable[neighbor->getTypeName()] ++;
		}
				
		return;
	}
}



//////////////////////////////////////////////////////////////////////
//  Build a list of possible interacting molecules
//
vector <Molecule *> Simulation::buildNeighborList(Molecule * m){

	// Get the list of molecule types that molecule m interacts with
	vector <moltype> typeList = getInteractionTypeList(m->getType());
	
	// Initialize a list that holds molecules near the molecule m
	vector <Molecule *> neighborList;
	
	// Go through the list of molecules that molecule m interacts with
	for(unsigned int molType = 0 ; molType < typeList.size() ; molType ++){
	
		// Get the list of molecules with the specified type
		vector <Molecule *> molList = molTable[typeList[molType]];
		
		// Add the molecule to the neighbor list if it is near to molecule m
		for(unsigned int molNum = 0 ; molNum < molList.size() ; molNum ++){
			
			Molecule * neighbor = molList[molNum];
			
			// Conditions to add to the neigbor list
			if(m->isActivated() && /*!(n->isActivated(m)) &&*/ m->isNear(neighbor))
			{
			  neighborList.push_back(neighbor);
			}
		}	
	}
	
	// Return a list of possible interacting molecules that are near
	return neighborList;
}



//////////////////////////////////////////////////////////////////////
//  Build a list of interacting molecules for the specified type
//
vector <int> Simulation::getInteractionTypeList(moltype type){
	// return a list of interacting molecules for the specified type
	return interactionNetworkTable[type];
}


//added by Adam
//////////////////////////////////////////////////////////////////////
//  Inject molecules in a fixed location at a given rate
//

void Simulation :: inject(){
	int gen_rate = getI("GEN_ENZYME_RATE");
	int FINE=100;
	int GEN_RATE_PER=10000000;
	//make sure that FINE*GEN_RATE_PER fits into an int
	int randomInt = RandomNG::randInt(0,GEN_RATE_PER*FINE); 
	if ( getCycles() == 1 || randomInt < gen_rate*FINE ) {
        	Molecule * m = createNewMolecule(getI("GEN_TYPE"));
		add(m);
                Point loc(getI("GEN_LOC_X"),getI("GEN_LOC_Y"),getI("GEN_LOC_Z"));
		m->setPosition(loc);
		m->activate();

          }
}


//////////////////////////////////////////////////////////////////////
//  Print details of activation event to output stream
//
void Simulation::printInteractionDetails(Molecule * activator, Molecule * activatee, int act){

	// '>' = activates , '|' = inhibits/ubitiquates
	string arrow = (act == 1) ? ">" : "|";
	
	// Output details
	printf(">\t%d\t%s\t%s\t--%s\t%s\t%s\n",
			numCycles,
			(activator->getTypeName()).c_str(),
			(activator->position()).c_str(),
			arrow.c_str(),
			(activatee->getTypeName()).c_str(),
			(activatee->position()).c_str());
}



//////////////////////////////////////////////////////////////////////
//  Return the current number of cycles executed in the simulation
//  ->	A cycle is the basic time unit in the simulation
//
int Simulation :: getCycles(){
	// Return the current number of cycles
	return numCycles;
}



//////////////////////////////////////////////////////////////////////
//  Get number of molecules in the Simulation
//
int Simulation :: getNumMolecules(){
	// Return the number of active/existing molecules within the Simulation Space
	return numMolecules;
}



//////////////////////////////////////////////////////////////////////
//  Get number of molecules with the specified type in the Simulation
//
int Simulation :: getNumMolecules(moltype type){
	// Return the number of molecules of the specified type in the Simulation Space
	return molTable[type].size();
}



//////////////////////////////////////////////////////////////////////
//  Is simulation completed
//  ->	Is the termination condition reached?
//
bool Simulation :: isSimulationCompleted(){
	// Return whether the simulation is completed
	return simulationCompleted;
}



//////////////////////////////////////////////////////////////////////
//  Get access to the Molecule table
//
map <moltype, vector <Molecule *> > * Simulation :: getMoleculeTable(){
	// Return a pointer to the Molecule table that manages all 
	// existing molecules in the Simulation Space
	return &(molTable);
}



//////////////////////////////////////////////////////////////////////
//  Get access to the Membranous Molecule table
//
map <moltype, vector <MembranousMolecule *> > * Simulation :: getMembranousMoleculeTable(){
	// Return a pointer to the Molecule table that manages all 
	// existing molecules in the Simulation Space
	return &(memMolTable);
}


//////////////////////////////////////////////////////////////////////
//  Perform periodic modifications to the Simulation Space
//  ->	Temporal events updating behaviour/Properties of the Simulation
//
void Simulation :: updateSimulation(){

}



//////////////////////////////////////////////////////////////////////
//  Start Simulation
//
void Simulation:: simulate(){

  // While the simulating pathway has not yet complete	
  while(!simulationCompleted){
    // Continue Simulation
    continueSimulation();
  }
}

//////////////////////////////////////////////////////////////////////
//  Continue simulation
//  ->	Events involved in each cycle
//
void Simulation:: continueSimulation(){
  const static int EQUILIBRIUM_TIME = 0;
  //  const static int METABOLITE_UPTAKE_TIME = 100;
  
  // Increment number of cycles (Time Units)
  numCycles ++;
  
  checkMoleculeInteractions();     //This line only used for signalling simulation
/* This block not needed for signalling simulation
  if (numCycles == 1)
  {
    lEnv->pumpInMetabolite(1.0);
  }
  
  lEnv->diffuse();
  */

  // inject enzyme A - added by Adam
  inject(); 

  transportAcrossMembranes();
  
  // Move the Molecules
  move();
  
  if (numCycles > EQUILIBRIUM_TIME)
  {
    // Check interaction between molecules
    checkEnzymaticReactions();
  }
  
  const static int MAX_CYCLES = 10000000;
  const static int NUM_OUTPUTS = 10000;

  // Periodically updating the Simulation
  //updateSimulation();
  if(numCycles % (MAX_CYCLES / NUM_OUTPUTS) == 0){
    Tuple * curStat = lEnv->getStateMatrix()->collapse();
    lEnv->getStatisticsMatrix()->appendRow(curStat);
  }

  if (numCycles % (MAX_CYCLES / NUM_OUTPUTS ) == 0)
  { //summary(); //debug
    //cout << numCycles / (MAX_CYCLES / NUM_OUTPUTS) << endl;
  }

  if(numCycles==MAX_CYCLES){
    printf("MAX_CYCLES reached");
    simulationCompleted = true;
  }else if(cutSitesActive == getI("CUT_SITES_TERM")){ 
    printf("%d\n",numCycles);
    simulationCompleted = true;
  }

  if (simulationCompleted)
  {
    // These stats are for the diffusible small molecules.
    //Matrix * stats = lEnv->getStatisticsMatrix();       line not needed for signalling simulation
    //stats->display();       line not needed for signalling simulation
    exit(1);
  }
}




//////////////////////////////////////////////////////////////////////
//  Create and return Molecule of the type TYPE
//
Molecule * Simulation :: createNewMolecule(moltype type){
	return NULL;
}



//////////////////////////////////////////////////////////////////////
//  Input of molecules to the Simulation Space
//
void Simulation::initialDataInput(){}



//////////////////////////////////////////////////////////////////////
//  Display data collected from the Simulation
//
void Simulation :: displayResults(){}



void Simulation :: checkEnzymaticReactions(){

        map<moltype, vector <Molecule *> >::const_iterator itr = molTable.begin();

        // Iterate through the list of molecules of differnet types
        for (itr = molTable.begin(); itr != molTable.end(); itr ++) {
                int numMol = (itr->second).size();
                for(int molIndex = 0 ; molIndex < numMol ; molIndex ++){
                        Enzyme * enz = (Enzyme *)(itr->second)[molIndex];

			int reactant = enz->reactant;
			int product = enz->product;
			int enzType = enz->getType();
			float forward_rate = enz->forwardReactionRate;
			float reverse_rate = enz->reverseReactionRate;
			float equilibrium_constant = enz->equilibriumConstant;
			float maximum_rate = enz->maximumRate;

			lEnv->reactWithMetabolite(enz->pos, reactant, product, enzType, forward_rate, reverse_rate, equilibrium_constant, maximum_rate);
                }
        }

}

