//============================================================================
// Name        : NBOS.cpp
// Author      : Sagar Indurkhya
// Version     : 0.1
// Copyright   : 2008
// Description : n-Birds-One-Stone algorithm
//============================================================================

#include "Intel_Compatibility.h"
#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include "dSFMT/dSFMT.h"

using namespace std;

typedef double NUMBER;

#define TROUGH_LENGTH 64
#define VERBOSE
//#define HALT_BEFORE_STARTING_ALGORITHM
//#define DEBUG_ON
#ifdef DEBUG_ON
	#define VERBOSE_DURING_ALGORITHM
	#define VERBOSE_DURING_ALGORITHM_LEVEL_2
	#define ULTRA_SPECIAL_DEBUG_MODE
	#define NORMAL_DEBUG_MODE
	#define EPSILON_CLOUD_UPDATE_TEST
#endif
#define EPSILON_VALUE 0.00000001
//#define PRINT_OUT_SPECIES_CONCENTRATIONS

vector<NUMBER> speciesConc;

enum VM_INSTRUCTION {
	UPDATE_SINGLE_RXN_IN_SUPERCACHE,	// 0
	UPDATE_DOUBLET_RXN_IN_SUPERCACHE,	// 1
	SET_CLOUD,							// 2
	UPDATE_RXN_IN_CLOUD,				// 3
	UPDATE_SUBTREE_IN_CLOUD,			// 4
	UPDATE_CLOUD_NO_PRIMARY,			// 5
	UPDATE_CLOUD_WITH_PRIMARY,			// 6
	SET_AND_UPDATE_CLOUD_ONE_PASS,		// 7
	SC_URIC_UCNP,						// 8
	SC_URIC_UCNP_UP_ONE,				// 9
	SC_URIC_UCNP_UP_TWO,				// 10
	SC_URIC_UCNP_UP_THREE				// 11
};

struct Reaction {
	int r1, r2, p1, p2;
	NUMBER rate;
};

class MySubTree {
public:
	unsigned int offset;
	unsigned int reBoost;
	unsigned int deMultiplyingChemical;
	unsigned int depth;
	NUMBER *heap;
	NUMBER *rxnRates;

	MySubTree() {}

	int ChooseReaction(NUMBER value) {
		value /= speciesConc[deMultiplyingChemical];

		int chosenRxn = 2;

		for(unsigned char i = 0; i < depth; ++i)
			for(int j = 0; j < TROUGH_LENGTH; ++j, ++chosenRxn)
				if(value <= heap[chosenRxn]) {
#ifdef ULTRA_SPECIAL_DEBUG_MODE
					if(heap[chosenRxn] <= 0.0)
						return -1;
#endif
					chosenRxn = TROUGH_LENGTH * (chosenRxn + j - 1) + 2;
					break;
				}
				else
					value -= heap[chosenRxn];

		return chosenRxn - offset + reBoost;
	}

	void UpdateRxnInTree(int rxn, NUMBER value) {
		rxn += offset;

		while(rxn > 1) {
			heap[rxn] += value;
			rxn = (rxn - 2)/TROUGH_LENGTH + 1;
		}

		heap[1] += value;
	}
};

class MyCloud {
public:
	NUMBER total, temp;

	unsigned int deMultiplyingChemical;

	// Main tree stuff
	unsigned int offset, reBoost, depth;
	NUMBER *heap;
	NUMBER *rxnRates;
#ifdef ULTRA_SPECIAL_DEBUG_MODE
	vector<Reaction> debugReactions;
#endif

	// Sub-Trees
	unsigned short numTrees;
	vector<NUMBER> subTreeWeights;
	vector<NUMBER> subTreeTotals;
	vector<MySubTree> subTrees;

	MyCloud() {}

	int ChooseReaction(NUMBER value) {
#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t\t" << "value before deMultiplying = " << value << std::endl;
		cout << "\t\t" << "speciesConc[deMultiplyingChemical] = " << speciesConc[deMultiplyingChemical] << std::endl;
		cout << "\t\t" << "deMultiplyingChemical = " << deMultiplyingChemical << std::endl;
		cout << "\t\t" << "heap[1] = " << heap[1] << std::endl;
#endif

		value /= speciesConc[deMultiplyingChemical];

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t\t" << "value = " << value << std::endl;
#endif

		if(value <= heap[1]) {

#ifdef VERBOSE_DURING_ALGORITHM
			cout << "\t\t" << "value is in main-tree." << std::endl;
			cout << "\t\t" << "depth = " << depth << std::endl;
#endif

			int chosenRxn = 2;

			for(unsigned char i = 0; i < depth; ++i)
				for(int j = 0; j < TROUGH_LENGTH; ++j, ++chosenRxn)
				{
#ifdef VERBOSE_DURING_ALGORITHM
					cout << "\t\t\t" << "chosenRxn = " << chosenRxn << "	heap[chosenRxn] = " << heap[chosenRxn] << std::endl;
#endif
					if(value <= heap[chosenRxn]) {
#ifdef ULTRA_SPECIAL_DEBUG_MODE
						if(heap[chosenRxn] <= 0.0)
							return -1;
#endif

						chosenRxn = (TROUGH_LENGTH)*(chosenRxn - 1) + 2;
						break;
					}
					else
						value -= heap[chosenRxn];
				}

			chosenRxn = (chosenRxn - 2)/TROUGH_LENGTH + 1;

#ifdef VERBOSE_DURING_ALGORITHM
			cout << "\t\t" << "chosenRxn = " << chosenRxn << std::endl;
			cout << "\t\t" << "offset = " << offset << std::endl;
			cout << "\t\t" << "reBoost = " << reBoost << std::endl;
#endif
			return chosenRxn - offset + reBoost;
		}

		value -= heap[1];

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t\t" << "value = " << value << std::endl;
#endif

		for(int i = 0; i < numTrees; i++) {
			if(value <= subTreeWeights[i])
				return subTrees[i].ChooseReaction(value);
			value -= subTreeWeights[i];
		}

		return -1;
	}

	void UpdateRxnInCloud(int rxn, NUMBER delta) {
		rxn += offset;

		while(rxn > 1) {
			heap[rxn] += delta;
			rxn = (rxn - 2)/TROUGH_LENGTH + 1;
		}

		heap[1] += delta;
		temp += delta;
	}

	// Optimized function for updating a reaction in a cloud with
	// no branch mispredictions! (Depth = 1)
	void UpdateRxnInCloud_UP_ONE(int rxn, NUMBER delta) {
		heap[rxn + offset] += delta;
		heap[1] += delta;
		temp += delta;
	}

	// Optimized function for updating a reaction in a cloud with
	// no branch mispredictions! (Depth = 2)
	void UpdateRxnInCloud_UP_TWO(int rxn, NUMBER delta) {
		rxn += offset;
		heap[rxn] += delta;
		rxn = (rxn - 2)/TROUGH_LENGTH + 1;
		heap[rxn] += delta;
		heap[1] += delta;
		temp += delta;
	}

	// Optimized function for updating a reaction in a cloud with
	// no branch mispredictions! (Depth = 3)
	void UpdateRxnInCloud_UP_THREE(int rxn, NUMBER delta) {
		rxn += offset;
		heap[rxn] += delta;
		rxn = (rxn - 2)/TROUGH_LENGTH + 1;
		heap[rxn] += delta;
		rxn = (rxn - 2)/TROUGH_LENGTH + 1;
		heap[rxn] += delta;
		heap[1] += delta;
		temp += delta;
	}
};

void EpsilonPreCheck(NUMBER theTree[], int locInTree, MyCloud *curCloud) {
	cout << "\t\t\t\t\t" << "Now performing an Epsilon test with eps = " << EPSILON_VALUE << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->heap[1] = " << curCloud->heap[1] << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->total = " << curCloud->total << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->temp = " << curCloud->temp << std::endl;
	cout << "\t\t\t\t\t" << "mainTree[mainTreeOffset + opCode] = " << theTree[locInTree] << std::endl;
	cout << "\t\t\t\t\t" << "speciesConc[curCloud->deMultiplyingChemical] = " << speciesConc[curCloud->deMultiplyingChemical] << std::endl;
	cout << "\t\t\t\t\t" << "mainTree[mainTreeOffset + opCode]/speciesConc[curCloud->deMultiplyingChemical] = " << theTree[locInTree]/speciesConc[curCloud->deMultiplyingChemical] << std::endl;
}

void EpsilonCheckTest(NUMBER theTree[], int locInTree, MyCloud *curCloud) {
	cout << "\t\t\t\t\t" << "Now performing an Epsilon test with eps = " << EPSILON_VALUE << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->numSubTrees = " << curCloud->subTrees.size() << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->heap[1] = " << curCloud->heap[1] << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->total = " << curCloud->total << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->temp = " << curCloud->temp << std::endl;
	cout << "\t\t\t\t\t" << "mainTree[mainTreeOffset + opCode] = " << theTree[locInTree] << std::endl;
	cout << "\t\t\t\t\t" << "curCloud->deMultiplyingChemical = " << curCloud->deMultiplyingChemical << std::endl;
	cout << "\t\t\t\t\t" << "speciesConc[curCloud->deMultiplyingChemical] = " << speciesConc[curCloud->deMultiplyingChemical] << std::endl;
	cout << "\t\t\t\t\t" << "mainTree[mainTreeOffset + opCode]/speciesConc[curCloud->deMultiplyingChemical] = " << theTree[locInTree]/speciesConc[curCloud->deMultiplyingChemical] << std::endl;

	NUMBER difference = theTree[locInTree]/speciesConc[curCloud->deMultiplyingChemical] - curCloud->total;
	cout << "\t\t\t\t\t" << "difference = " << difference << std::endl;
	if(speciesConc[curCloud->deMultiplyingChemical] > 0.0 && (difference > EPSILON_VALUE || difference < -1.0 * EPSILON_VALUE)) {
		cout << "ERROR: (difference > epsilon) or (difference < -1.0 * epsilon)" << std::endl;
		cout << "curCloud.offset = " << curCloud->offset << std::endl;
		exit(0);
	}
	else
		cout << "\t\t\t\t\t" << "Epsilon test passed!" << std::endl;

	if(curCloud->subTrees.size() == 0) {
		NUMBER diff2 = curCloud->total - curCloud->heap[1];
		if(diff2 > EPSILON_VALUE || diff2 < -1.0 * EPSILON_VALUE) {
			cout << "ERROR: (diff2 > epsilon) or (diff2 < -1.0 * epsilon)" << std::endl;
			cout << "diff2 = " << diff2 << std::endl;
			exit(0);
		}
	}
}

// Assumes that entries will have at least depth of 2. This is not a
// normal tree-depth measurement as most people would think.
int getSpecialDepth(int numEntries) {
	int depth = 2, k = 64;
	while(k < numEntries) {
		depth += 1;
		k *= TROUGH_LENGTH;
	}

	return depth - 1;
}

void ValidateMainTree(NUMBER theTree[], NUMBER epsilon, int TotalSizeOfTree)
{
	for(int parent = 1; parent < TotalSizeOfTree; ++parent) {
		// Now we want to show that the sum of the parents
		NUMBER valueOfParent = theTree[parent];

		NUMBER sumOfChildren = 0.0;

		// Now establish boundaries
		int firstChild = TROUGH_LENGTH*(parent - 1) + 2;

		bool done = false;

		for(int i = 0; i < TROUGH_LENGTH; i++) {
			if(firstChild + i < TotalSizeOfTree) {

				if(theTree[firstChild + i] < 0.0) {
					cout << "\t\t\t\t" << "theTree[firstChild + i] = " << theTree[firstChild + i] << std::endl;
					exit(0);
				}
				sumOfChildren += theTree[firstChild + i];
			} else {
				done = true;
				break;
			}
		}

		// Now check to make sure and report the violation if one exists
		NUMBER difference = valueOfParent - sumOfChildren;

		cout << "\t\t" << "For parent = " << parent << std::endl;

		if(difference < 0.0 || difference > epsilon) {
			cout << "\t\t\t" << "valueOfParent = " << valueOfParent << std::endl;
			cout << "\t\t\t" << "sumOfChildren = " << sumOfChildren << std::endl;
		}

		cout << "\t\t\t" << "difference = " << difference << std::endl;

		if(done)
			return;
	}
}

#ifdef UMT_DEPTH_PLAIN
void UpdateMainTree(NUMBER mainHeap[], int startNode, NUMBER increment) {
	while(startNode > 1) {
		mainHeap[startNode] += increment;
#ifdef VERBOSE_DURING_ALGORITHM
		if(mainHeap[startNode] < -0.000000001) {
			cout << "ERROR: mainHeap[startNode] < 0.0" << std::endl;
			cout << "\t\t" << "increment = " << increment << std::endl;
			cout << "\t\t" << "startNode = " << startNode << std::endl;
			cout << "\t\t" << "mainHeap[startNode] = " << mainHeap[startNode] << std::endl;
			exit(0);
		}
#endif
		startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	}

	mainHeap[1] += increment;
}
#endif

#ifdef UMT_DEPTH_ONE
void UpdateMainTree(NUMBER mainHeap[], int startNode, NUMBER increment) {
	mainHeap[startNode] += increment;
	mainHeap[1] += increment;
}
#endif

#ifdef UMT_DEPTH_TWO
void UpdateMainTree(NUMBER mainHeap[], int startNode, NUMBER increment) {
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	mainHeap[1] += increment;
}
#endif

#ifdef UMT_DEPTH_THREE
void UpdateMainTree(NUMBER mainHeap[], int startNode, NUMBER increment) {
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	mainHeap[1] += increment;
}
#endif

#ifdef UMT_DEPTH_FOUR
void UpdateMainTree(NUMBER mainHeap[], int startNode, NUMBER increment) {
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	startNode = (startNode - 2)/TROUGH_LENGTH + 1;
	mainHeap[startNode] += increment;
	mainHeap[1] += increment;
}
#endif

/* PARAMETERS ARE:
 * 	1) numIterations
 *  2) BasicInputFileName
 *  3) Species File
 *  4) Reaction File
 *	5) randNumGenSeed
 *	6) OutputTimeFileName
 */
int main(int argc, char* argv[]) {

	cout << setprecision(10);

	dsfmt_t dsfmt;
	dsfmt_init_gen_rand(&dsfmt, atoi(argv[5]));

#ifdef VERBOSE
	cout << "Beginning NBOS PROGRAM" << std::endl;
#endif

	// Load stuff in
	NUMBER simulationTime = 0.0;
	NUMBER* mainTree = NULL;
	NUMBER SuperCacheTotalWeight = 0.0;

	int numSpecies, numReactions, numSuperCacheRxns, numClouds;
	int mainTreeOffset = -1, mainTreeDepth = -1;

	vector<MyCloud> myClouds;
	vector<Reaction> superCacheRxns;
	vector<NUMBER> superCacheWeights;

	vector<Reaction> SystemReactions;

	vector<unsigned int> vmTapeIndex; // Has size of numChemicals
	vector<unsigned int> vmTape;

	vector<unsigned int> rxnTapeIndex; // Has size of numRxns
	vector<int> rxnTape;

	MyCloud dudCloud;

	int numIterations = atoi(argv[1]);

	// Load up our basic input file
	ifstream basicInput(argv[2]);
	ifstream speciesInput(argv[3]);
	ifstream rxnInput(argv[4]);

#ifdef VERBOSE
	cout << "Read in input files, preparing to load data." << std::endl;
	cout << "argv[2] = " << argv[2] << std::endl;
	cout << "argv[3] = " << argv[3] << std::endl;
	cout << "argv[4] = " << argv[4] << std::endl;
	cout << std::endl;
#endif

	// Load in the numSpecies and numReactions
	basicInput >> numSpecies;
	basicInput >> numReactions;
	basicInput >> numSuperCacheRxns;
	basicInput >> numClouds;
	basicInput >> mainTreeOffset;

#ifdef VERBOSE
	cout << "\t" << "numIterations = " << numIterations << std::endl;
	cout << "\t" << "numSpecies = " << numSpecies << std::endl;
	cout << "\t" << "numReactions = " << numReactions << std::endl;
	cout << "\t" << "numSuperCacheRxns = " << numSuperCacheRxns << std::endl;
	cout << "\t" << "numClouds = " << numClouds << std::endl;
	cout << "\t" << "mainTreeOffset = " << mainTreeOffset << std::endl;
#endif

	mainTreeDepth = getSpecialDepth(numClouds);

	// Now load in the species
	for(int i = 0; i < numSpecies; i++) {
		NUMBER c;
		speciesInput >> c;
		speciesConc.push_back(c);
	}

	// Now load the super-cached reactions
	for(int i = 0; i < numSuperCacheRxns; i++) {
		Reaction r;
		rxnInput >> r.r1 >> r.r2 >> r.p1 >> r.p2 >> r.rate;
		if(r.r1 == r.r2)
			r.rate /= 2.0;
		superCacheRxns.push_back(r);
		SystemReactions.push_back(r);

		if(r.r1 == r.r2)
			superCacheWeights.push_back(speciesConc[r.r1] * (speciesConc[r.r2] - 1.0) * r.rate);
		else
			superCacheWeights.push_back(speciesConc[r.r1] * speciesConc[r.r2] * r.rate);
	}

#ifdef VERBOSE
	cout << "Loaded Super-Cached Reactions" << std::endl;
#endif

	// Now load in the basic clouds
	for(int i = 0; i < numClouds; i++) {

#ifdef VERBOSE
		cout << "\t" << "Loading in a basic cloud # " << i << std::endl;
#endif
		MyCloud mc;
		myClouds.push_back(mc);

		int numRxnsInThisMainTree;
		int mainHeapSizeForThisTree;

		basicInput >> myClouds.back().deMultiplyingChemical;
		basicInput >> myClouds.back().offset;
		basicInput >> myClouds.back().reBoost;
		basicInput >> myClouds.back().numTrees;
		basicInput >> numRxnsInThisMainTree;

		mainHeapSizeForThisTree = myClouds.back().offset + numRxnsInThisMainTree;
		myClouds.back().depth = getSpecialDepth(numRxnsInThisMainTree);

#ifdef VERBOSE
		cout << "\t\t" << "deMultiplyingChemical = " << myClouds.back().deMultiplyingChemical << std::endl;
		cout << "\t\t" << "reBoost = " << myClouds.back().reBoost << std::endl;
		cout << "\t\t" << "numSubTrees = " << myClouds.back().numTrees << std::endl;
		cout << "\t\t" << "numRxnsInThisMainTree = " << numRxnsInThisMainTree << std::endl;
		cout << "\t\t" << "mainHeapSizeForThisTree = " << mainHeapSizeForThisTree << std::endl;
		cout << "\t\t" << "depth = " << myClouds.back().depth << std::endl;
		cout << "\t\t---" << std::endl;
		cout << "\t\t" << "Setting up main tree... " << std::endl;
		cout << "\t\t\t" << "myClouds.size() = " << myClouds.size() << std::endl;
#endif

		// Set up the main tree
		myClouds.back().heap = new NUMBER[mainHeapSizeForThisTree];

#ifdef VERBOSE
		cout << "\t\t\t" << "Was able to allocate an array." << std::endl;
#endif

		for(int j = 0; j < mainHeapSizeForThisTree; j++)
			myClouds.back().heap[j] = 0.0;

		myClouds.back().rxnRates = new NUMBER[numRxnsInThisMainTree];

		for(int j = 0; j < numRxnsInThisMainTree; j++)
			myClouds.back().rxnRates[j] = 0.0;

#ifdef VERBOSE
		cout << "\t\t\t" << "COMPLETE." << std::endl;
		cout << "\t\t" << "Loading cloud's main tree's reactions..." << std::endl;
#endif

		// Now load in the main tree's reactions
		for(int j = 0; j < numRxnsInThisMainTree; j++) {

#ifdef VERBOSE
			cout << "\t\t\t\t" << "Loading reaction #" << j << std::endl;
#endif

			Reaction r;
			rxnInput >> r.r1 >> r.r2 >> r.p1 >> r.p2 >> r.rate;
			if(r.r1 == r.r2)
				r.rate /= 2.0;

			SystemReactions.push_back(r);
			// Now see if we need to correct this reaction
			if(r.r1 != myClouds.back().deMultiplyingChemical) {
				cout << "HELLO WORLD!!! ERROR: r.r1 != myClouds.back().deMultiplyingChemical" << std::endl;
				cout << "myClouds.back().deMultiplyingChemical = " << myClouds.back().deMultiplyingChemical << std::endl;
				cout << "r.r1 = " << r.r1 << std::endl;
				cout << "r.r2 = " << r.r2 << std::endl;
				exit(0);
			}

			myClouds.back().rxnRates[j] = r.rate;

			if(r.r1 == r.r2)
				myClouds.back().UpdateRxnInCloud(j, (speciesConc[r.r1] - 1.0) * r.rate);
			else
				myClouds.back().UpdateRxnInCloud(j, speciesConc[r.r2] * r.rate);

#ifdef ULTRA_SPECIAL_DEBUG_MODE
			myClouds.back().debugReactions.push_back(r);
#endif
		}

		myClouds.back().total = myClouds.back().heap[1];

#ifdef VERBOSE
		cout << "\t\t\t" << "COMPLETE." << std::endl;
		cout << "\t\t" << "Setting up the subTrees..." << std::endl;
#endif

		// Now setup the subTrees
		for(int j = 0; j < myClouds.back().numTrees; j++) {

#ifdef VERBOSE
			cout << "\t\t\t" << "Loading subTree #" << j << std::endl;
#endif

			MySubTree mst;
			myClouds.back().subTrees.push_back(mst);

			int numRxnsInThisSubTree;

			basicInput >> myClouds.back().subTrees.back().offset >> myClouds.back().subTrees.back().reBoost;
			basicInput >> myClouds.back().subTrees.back().deMultiplyingChemical >> numRxnsInThisSubTree;

			myClouds.back().subTrees.back().heap = new NUMBER[myClouds.back().subTrees.back().offset + numRxnsInThisSubTree];
			myClouds.back().subTrees.back().rxnRates = new NUMBER[numRxnsInThisSubTree];
			myClouds.back().subTrees.back().depth = getSpecialDepth(numRxnsInThisSubTree);

			for(int k = 0; k < numRxnsInThisSubTree; k++) {
				Reaction r;
				rxnInput >> r.r1 >> r.r2 >> r.p1 >> r.p2 >> r.rate;
				if(r.r1 == r.r2)
					r.rate /= 2.0;
				SystemReactions.push_back(r);
				myClouds.back().subTrees.back().UpdateRxnInTree(k, r.rate);
			}
		}

		for(int j = 0; j < myClouds.back().numTrees; j++)
			myClouds.back().total += myClouds.back().subTrees[j].heap[1] * speciesConc[myClouds.back().subTrees[j].deMultiplyingChemical];

		myClouds.back().temp = 0.0;
	}

#ifdef VERBOSE
	cout << "Loaded Clouds and SubTrees." << std::endl;
	cout << "Num Species = " << numSpecies << std::endl;
#endif

	// Now load in the virtual machine tape.
	for(int i = 0; i < numSpecies; i++) {
		int numEntriesForTape;
		basicInput >> numEntriesForTape;

		vmTapeIndex.push_back(vmTape.size());
		vmTape.push_back(numEntriesForTape);

		for(int j = 0; j < numEntriesForTape; j++) {
			int s, t, u;
			basicInput >> s >> t;

			if(s == SC_URIC_UCNP)
			{
				basicInput >> u;
				// Handle the case of optimized SC_URIC_UCNP's.
				if(myClouds.at(t).depth == 1)
					s = SC_URIC_UCNP_UP_ONE;
				else if(myClouds.at(t).depth == 2)
					s = SC_URIC_UCNP_UP_TWO;
				else if(myClouds.at(t).depth == 3)
					s = SC_URIC_UCNP_UP_THREE;
				vmTape.push_back(s);
				vmTape.push_back(t);
				vmTape.push_back(u);
			}
			else
			{
				vmTape.push_back(s);
				vmTape.push_back(t);
			}

		}
	}

#ifdef VERBOSE
	cout << "Loaded Virtual Machine Tape." << std::endl;
	cout << "numReactions = " << numReactions << std::endl;
#endif

	// Now load in the reaction tape.
	for(int i = 0; i < numReactions; i++) {
		int numEntriesForTape;
		basicInput >> numEntriesForTape;

		rxnTapeIndex.push_back(rxnTape.size());
		rxnTape.push_back(numEntriesForTape);

		for(int j = 0; j < numEntriesForTape; j++) {
			int s, t;
			basicInput >> s >> t;
			rxnTape.push_back(s);
			rxnTape.push_back(t);
		}
	}

#ifdef VERBOSE
	cout << "Loaded Reaction Tape." << std::endl;
	cout << "mainTreeOffset = " << mainTreeOffset << "  and numClouds = " << numClouds << std::endl;
#endif

	// Now we have to create and initialize the main tree.
	mainTree = new NUMBER[mainTreeOffset + numClouds];

	for(int i = 0; i < mainTreeOffset + (signed)myClouds.size(); i++)
		mainTree[i] = 0.0;

	for(int i = 0; i < (signed)myClouds.size(); i++)
	{
		MyCloud *curCloud = &(myClouds[i]);
		UpdateMainTree(mainTree, mainTreeOffset + i, curCloud->total * speciesConc.at(curCloud->deMultiplyingChemical));
		cout << "EpsilonCheckTest, i = " << i << std::endl;
		EpsilonCheckTest(mainTree, mainTreeOffset + i, curCloud);
	}

#ifdef VERBOSE
	cout << "Validating Main Tree... " << std::endl;
	ValidateMainTree(mainTree, 0.00000001, mainTreeOffset + numClouds);

	cout << "Created Main Tree structure." << std::endl;
#endif

	for(int i = 0; i < superCacheWeights.size(); i++)
		SuperCacheTotalWeight += superCacheWeights[i];

#ifdef PRINT_OUT_SPECIES_CONCENTRATIONS
	cout << "Creating file for species concentrations" << std::endl;
	ofstream outputSpecies("NBOS_SPECIES_OUTPUT.dat");
#endif

#ifdef VERBOSE
	cout << std::endl;
	cout << "-------------------------------------------------" << std::endl;
#endif

#ifdef HALT_BEFORE_STARTING_ALGORITHM
	int tempHalt = 0;
	cin >> tempHalt;
#endif

	timeval timVal;
	gettimeofday(&timVal, NULL);
	double StartTime = timVal.tv_sec + (timVal.tv_usec/1000000.0);

	// Algorithm starts here!
	for (int i = 0; i < numIterations; ) {

		START_OF_ALGORITHM:

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "Starting Iteration #" << i << "." << std::endl;
		cout << "\t" << "mainTree[1] = " << mainTree[1] << std::endl;
		cout << "\t" << "SuperCacheTotalWeight = " << SuperCacheTotalWeight << std::endl;
#endif

		// See if simulation is done
		if(mainTree[1] + SuperCacheTotalWeight < 0.00000000001)
		{
			cout << "SIMULATION TERMINATED BECAUSE a0 < 0.000000000001" << std::endl;
			cout << "Simulation Time = " << simulationTime << std::endl;
			exit(0);
			break;
		}

		NUMBER randNum = (mainTree[1] + SuperCacheTotalWeight) * dsfmt_genrand_close_open(&dsfmt);

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t" << "randNum = " << randNum << std::endl;
#endif

		int rxnToExecute = -1;

		if(randNum <= SuperCacheTotalWeight)
		{
			// Check if it is in the super-cached weights
			for (int j = 0; j < numSuperCacheRxns; j++) {
				if (randNum <= superCacheWeights[j]) {
					rxnToExecute = j;
#ifdef VERBOSE_DURING_ALGORITHM
					cout << "\t" << "Chose a reaction from Super-Cache." << std::endl;
#endif
#ifdef ULTRA_SPECIAL_DEBUG_MODE
					if(superCacheWeights[j] <= 0.0)
						goto START_OF_ALGORITHM;
#endif
					goto RXN_TO_EXECUTE_HAS_BEEN_DECIDED;
				}

				randNum -= superCacheWeights[j];
			}
		} else
			randNum -= SuperCacheTotalWeight;

		// Search through the main cloud tree
		int selectedCloud = 2;

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t\t" << "Main Tree Depth = " << mainTreeDepth << std::endl;
		cout << "\t\t" << "mainTreeOffset = " << mainTreeOffset << std::endl;
#endif

		for(int j = 0; j < mainTreeDepth; ++j)
			for(int k = 0; k < TROUGH_LENGTH; ++k)
			{
#ifdef VERBOSE_DURING_ALGORITHM
				cout << "\t\t\t" << "*** randNum = " << randNum << std::endl;
				cout << "\t\t\t" << "selectedCloud + k = " << (selectedCloud + k) << "	mainTree[selectedCloud + k] = " << mainTree[selectedCloud + k] << std::endl;
#endif
				if(randNum <= mainTree[selectedCloud + k]) {
#ifdef ULTRA_SPECIAL_DEBUG_MODE
					if(mainTree[selectedCloud + k] <= 0.0)
						goto START_OF_ALGORITHM;
#endif
					selectedCloud = (TROUGH_LENGTH)*((selectedCloud + k)-1) + 2;
					break;
				}
				else
					randNum -= mainTree[selectedCloud + k];
			}

		selectedCloud = (selectedCloud - 2)/TROUGH_LENGTH + 1;
		selectedCloud -= mainTreeOffset;

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t\t" << "selectedCloud = " << selectedCloud << std::endl;
		cout << "\t\t" << "total = " << myClouds.at(selectedCloud).total << std::endl;
		cout << "\t\t" << "Size of myClouds = " << myClouds.size() << std::endl;

		int lolcat = 0;
		cin >> lolcat;
#endif

		rxnToExecute = myClouds.at(selectedCloud).ChooseReaction(randNum);

		if(rxnToExecute == -1)
			continue;

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t" << "Chose a reaction from Clouds." << std::endl;
#endif

#ifdef ULTRA_SPECIAL_DEBUG_MODE
		int tempRxnLocation = rxnToExecute - myClouds.at(selectedCloud).reBoost;
		cout << "tempRxnLocation = " << tempRxnLocation << std::endl;
		cout << "myClouds.at(selectedCloud).reBoost = " << myClouds.at(selectedCloud).reBoost << std::endl;
		cout << "myClouds.at(selectedCloud).offset = " << myClouds.at(selectedCloud).offset << std::endl;
		cout << "myClouds.at(selectedCloud).debugReactions.size() = " << myClouds.at(selectedCloud).debugReactions.size() << std::endl;
		int tempR1 = myClouds.at(selectedCloud).debugReactions.at(tempRxnLocation).r1;
		int tempR2 = myClouds.at(selectedCloud).debugReactions.at(tempRxnLocation).r2;
		cout << "Reaction = " << tempR1 << " + " << tempR2 << " --> ..." << std::endl;
		cout << "speciesConc[tempRxn.r1] = " << speciesConc[tempR1] << std::endl;
		cout << "speciesConc[tempRxn.r2] = " << speciesConc[tempR2] << std::endl;
#endif

		RXN_TO_EXECUTE_HAS_BEEN_DECIDED:

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t" << "RxnToExecute = " << rxnToExecute << std::endl;
		cout << "\t\t" << "Rxn = " << SystemReactions.at(rxnToExecute).r1 << " + " << SystemReactions.at(rxnToExecute).r2 << " --> " << SystemReactions.at(rxnToExecute).p1 << " + " << SystemReactions.at(rxnToExecute).p2 << ", k =  " << SystemReactions.at(rxnToExecute).rate << std::endl;
#endif

		simulationTime += -1.0/mainTree[1] * log(dsfmt_genrand_close_open(&dsfmt));

		// Now load up the tape and execute it!
		MyCloud *curCloud = NULL;

		int startingPosition = rxnTapeIndex.at(rxnToExecute);
		int numChangingChems = rxnTape.at(startingPosition);
		++startingPosition;

#ifdef VERBOSE_DURING_ALGORITHM
		cout << "\t" << "numChangingChems = " << numChangingChems << std::endl;
		cout << "\t" << "startingPosition = " << startingPosition << std::endl;
#endif
		for(int cntdown = numChangingChems; cntdown; --cntdown, startingPosition += 2) {
			int curChem = rxnTape.at(startingPosition);
			NUMBER curDelta = 1.0 * rxnTape[startingPosition + 1];

			speciesConc[curChem] += curDelta;

#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2

			cout << "\t\t" << "curChem = " << curChem << std::endl;
			cout << "\t\t" << "curDelta = " << curDelta << std::endl;
			cout << "\t\t" << "speciesConc[curChem] = " << speciesConc.at(curChem) << std::endl;
			cout << "\t\t" << "vmTapeIndex[curChem] = " << vmTapeIndex.at(curChem) << std::endl;
			cout << "\t\t" << "vmTape[vmTapeIndex[curChem]] = " << vmTape.at(vmTapeIndex.at(curChem)) << std::endl;
			if(speciesConc[curChem] - ((int)speciesConc[curChem]) != 0.0)
			{
				cout << "CUR_DELTA = " << curDelta << std::endl;
				exit(0);
			}
			if(speciesConc[curChem] < 0.0)
			{
				cout << "EPIC FAIL. speciesConc[curChem] < 0.0" << std::endl;
				exit(0);
			}
#endif

			for(int k = vmTape[vmTapeIndex[curChem]], s = vmTapeIndex[curChem] + 1; k; --k) {
				int curInstruction = vmTape[s], opCode = vmTape[s + 1];

#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
				NUMBER tempIncrement = 0.0;
#endif
				switch(curInstruction) {
					case UPDATE_SINGLE_RXN_IN_SUPERCACHE:
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						tempIncrement = -1.0 * superCacheWeights[opCode];
#endif
						SuperCacheTotalWeight -= superCacheWeights[opCode];
						superCacheWeights[opCode] = speciesConc[superCacheRxns[opCode].r1] * speciesConc[superCacheRxns[opCode].r2] * superCacheRxns[opCode].rate;
						SuperCacheTotalWeight += superCacheWeights[opCode];
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						tempIncrement += superCacheWeights[opCode];
						cout << "\t\t\t" << "UPDATE_SINGLE_RXN_IN_SUPERCACHE... " << std::endl;
						cout << "\t\t\t\t" << "superCacheRxns[opCode].r1 = " << superCacheRxns[opCode].r1 << std::endl;
						cout << "\t\t\t\t\t" << "speciesConc[superCacheRxns[opCode].r1] = " << speciesConc[superCacheRxns[opCode].r1] << std::endl;
						cout << "\t\t\t\t" << "superCacheRxns[opCode].r2 = " << superCacheRxns[opCode].r2 << std::endl;
						cout << "\t\t\t\t\t" << "speciesConc[superCacheRxns[opCode].r2] = " << speciesConc[superCacheRxns[opCode].r2] << std::endl;
						cout << "\t\t\t\t" << "superCacheRxns[opCode].rate = " << superCacheRxns[opCode].rate << std::endl;
						cout << "\t\t\t\t\t" << "superCacheWeights[opCode] = " << superCacheWeights[opCode] << std::endl;
						cout << "\t\t\t\t\t" << "Super-Cache += " << tempIncrement << std::endl;
						if(superCacheWeights[opCode] < 0.0 && superCacheWeights[opCode] < -0.000000001)
						{
							cout << "EPIC FAIL, superCacheWeights[opCode] < 0.0. " << std::endl;
							exit(0);
						}
#endif
						s += 2;
						break;
					case UPDATE_DOUBLET_RXN_IN_SUPERCACHE:
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						tempIncrement = -1.0 * superCacheWeights[opCode];
#endif
						SuperCacheTotalWeight -= superCacheWeights[opCode];
						NUMBER temp = speciesConc[curChem];
						superCacheWeights[opCode] = temp * (temp - 1.0) * superCacheRxns[opCode].rate;
						SuperCacheTotalWeight += superCacheWeights[opCode];
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						tempIncrement += superCacheWeights[opCode];
						cout << "\t\t\t" << "UPDATE_DOUBLET_RXN_IN_SUPERCACHE... " << std::endl;
						cout << "\t\t\t\t" << "Super-Cache += " << tempIncrement << std::endl;
#endif
						s += 2;
						break;
					case SET_CLOUD:
						curCloud = &(myClouds[opCode]);
						curCloud->temp = 0.0;
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SET_CLOUD" << std::endl;
						cout << "\t\t\t\t" << "curCloud = " << opCode << std::endl;
						cout << "\t\t\t\t" << "curCloud->deMultiplyingChemical = " << curCloud->deMultiplyingChemical << std::endl;
#endif
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Check..." << std::endl;
						EpsilonPreCheck(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 2;
						break;
					case UPDATE_RXN_IN_CLOUD:
						curCloud->UpdateRxnInCloud(opCode, curDelta * curCloud->rxnRates[opCode]);
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "UPDATE_RXN_IN_CLOUD" << std::endl;
#endif
#ifdef ULTRA_SPECIAL_DEBUG_MODE
						cout << "\t\t\t\t" << "curCloud->debugReactions[opCode].r1 = " << curCloud->debugReactions[opCode].r1 << std::endl;
						cout << "\t\t\t\t" << "curCloud->debugReactions[opCode].r2 = " << curCloud->debugReactions[opCode].r2 << std::endl;
						cout << "\t\t\t\t" << "curCloud->debugReactions[opCode].rate = " << curCloud->debugReactions[opCode].rate << std::endl;
#endif
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t\t" << "curCloud->rxnRates[opCode] = " << curCloud->rxnRates[opCode] << std::endl;
						cout << "\t\t\t\t" << "curCloud->temp = " << curCloud->temp << std::endl;

						if(curCloud->heap[curCloud->offset + opCode] < -0.000000001 && speciesConc[curCloud->debugReactions[opCode].r2] > 0.0) {
							cout << "EPIC FAIL: curCloud->heap[curCloud.offset + opCode] < 0.0" << std::endl;
							cout << "speciesConc[curCloud->debugReactions[opCode].r1] = " << speciesConc[curCloud->debugReactions[opCode].r1] << std::endl;
							cout << "speciesConc[curCloud->debugReactions[opCode].r2] = " << speciesConc[curCloud->debugReactions[opCode].r2] << std::endl;
							cout << "\t" << "curCloud->heap[curCloud.ofset + opCode] = " << curCloud->heap[curCloud->offset + opCode] << std::endl;
							exit(0);
						}


#endif
						s += 2;
						break;
					case UPDATE_SUBTREE_IN_CLOUD:
						curCloud->subTreeWeights[opCode] += curDelta * curCloud->subTreeTotals[opCode];
						curCloud->temp += curCloud->subTreeTotals[opCode];
						cout << "OMG LOLZ SUBTREE UPDATE!!!" << std::endl;
						exit(0);
						s += 2;
						break;
					case UPDATE_CLOUD_NO_PRIMARY:
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						tempIncrement = speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp;
						cout << "\t\t\t" << "UPDATE_CLOUD_NO_PRIMARY" << std::endl;
						cout << "\t\t\t\t" << "mainCloudTree += " << tempIncrement << std::endl;
						cout << "\t\t\t\t" << "curCloud->temp = " << curCloud->temp << std::endl;
#endif

						curCloud->total += curCloud->temp;
						UpdateMainTree(mainTree, mainTreeOffset + opCode, speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp);
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 2;
						break;
					case UPDATE_CLOUD_WITH_PRIMARY:
						// a = deMulChem + deltaDeMulChem
						// b = total + temp
						// (a + da)(b + db) - a*b = a*db + da*b + da*db = db*(da+a) * da*b
						NUMBER zebra = curCloud->temp * speciesConc[curCloud->deMultiplyingChemical] + (curDelta * curCloud->total);
						curCloud->total += curCloud->temp;

#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "UPDATE_CLOUD_WITH_PRIMARY" << std::endl;
						cout << "\t\t\t\t" << "curCloud->temp = " << curCloud->temp << std::endl;
						cout << "\t\t\t\t" << "curCloud->deMultiplyingChemical = " << curCloud->deMultiplyingChemical << std::endl;
						cout << "\t\t\t\t" << "speciesConc[curCloud->deMultiplyingChemical] = " << speciesConc[curCloud->deMultiplyingChemical] << std::endl;
						cout << "\t\t\t\t" << "curCloud->total = " << curCloud->total << std::endl;
						cout << "\t\t\t\t" << "mainCloudTree += " << zebra << std::endl;
#endif

						UpdateMainTree(mainTree, mainTreeOffset + opCode, zebra);

#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 2;
						break;
					case SET_AND_UPDATE_CLOUD_ONE_PASS:
						curCloud = &(myClouds[opCode]);
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SET_AND_UPDATE_CLOUD_ONE_PASS" << std::endl;
						cout << "\t\t\t\t" << "curCloud = " << opCode << std::endl;
						cout << "\t\t\t\t" << "mainCloudTree += " << curDelta * curCloud->total << std::endl;
#endif
						UpdateMainTree(mainTree, mainTreeOffset + opCode, curDelta * curCloud->total);

#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 2;
						break;
					case SC_URIC_UCNP:
						int opCode2_PLAIN = vmTape[s + 2];
						curCloud = &(myClouds[opCode]);
						curCloud->temp = 0.0;
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SC_URIC_UCNP" << std::endl;
						cout << "\t\t\t\t" << "opCode = " << opCode << std::endl;
						cout << "\t\t\t\t" << "opCode2_PLAIN = " << opCode2_PLAIN << std::endl;
#endif
						curCloud->UpdateRxnInCloud(opCode2_PLAIN, curDelta * curCloud->rxnRates[opCode2_PLAIN]);
						curCloud->total += curCloud->temp;
						UpdateMainTree(mainTree, mainTreeOffset + opCode, speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp);
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 3;
						break;
					case SC_URIC_UCNP_UP_ONE:
						int opCode2_UP_ONE = vmTape[s + 2];
						curCloud = &(myClouds[opCode]);
						curCloud->temp = 0.0;
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SC_URIC_UCNP_UP_ONE" << std::endl;
						cout << "\t\t\t\t" << "opCode = " << opCode << std::endl;
						cout << "\t\t\t\t" << "opCode2_UP_ONE = " << opCode2_UP_ONE << std::endl;
#endif
						curCloud->UpdateRxnInCloud_UP_ONE(opCode2_UP_ONE, curDelta * curCloud->rxnRates[opCode2_UP_ONE]);
						curCloud->total += curCloud->temp;
						UpdateMainTree(mainTree, mainTreeOffset + opCode, speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp);
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 3;
						break;
					case SC_URIC_UCNP_UP_TWO:
						int opCode2_UP_TWO = vmTape[s + 2];
						curCloud = &(myClouds[opCode]);
						curCloud->temp = 0.0;
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SC_URIC_UCNP_UP_TWO" << std::endl;
						cout << "\t\t\t\t" << "opCode = " << opCode << std::endl;
						cout << "\t\t\t\t" << "opCode2_UP_TWO = " << opCode2_UP_TWO << std::endl;
#endif
						curCloud->UpdateRxnInCloud_UP_TWO(opCode2_UP_TWO, curDelta * curCloud->rxnRates[opCode2_UP_TWO]);
						curCloud->total += curCloud->temp;
						UpdateMainTree(mainTree, mainTreeOffset + opCode, speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp);
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 3;
						break;
					case SC_URIC_UCNP_UP_THREE:
						int opCode2_UP_THREE = vmTape[s + 2];
						curCloud = &(myClouds[opCode]);
						curCloud->temp = 0.0;
#ifdef VERBOSE_DURING_ALGORITHM_LEVEL_2
						cout << "\t\t\t" << "SC_URIC_UCNP_UP_THREE" << std::endl;
						cout << "\t\t\t\t" << "opCode = " << opCode << std::endl;
						cout << "\t\t\t\t" << "opCode2_UP_THREE = " << opCode2_UP_THREE << std::endl;
#endif
						curCloud->UpdateRxnInCloud_UP_THREE(opCode2_UP_THREE, curDelta * curCloud->rxnRates[opCode2_UP_THREE]);
						curCloud->total += curCloud->temp;
						UpdateMainTree(mainTree, mainTreeOffset + opCode, speciesConc[curCloud->deMultiplyingChemical] * curCloud->temp);
#ifdef EPSILON_CLOUD_UPDATE_TEST
						cout << "\t\t\t\t" << "Doing a post Epsilon-Cloud-Update-Test..." << std::endl;
						EpsilonCheckTest(mainTree, mainTreeOffset + opCode, curCloud);
#endif
						s += 3;
						break;
				}
			}
		}

#ifdef VERBOSE_DURING_ALGORITHM
		NUMBER sagar;
		cin >> sagar;
#endif

#ifdef PRINT_OUT_SPECIES_CONCENTRATIONS
		if(i % 100)
		{
			outputSpecies << simulationTime;
			for(int horse = 0; horse < numSpecies; horse++)
				outputSpecies << "	" << speciesConc[horse];
			outputSpecies << std::endl;
		}
#endif
		// Now output whatever chemicals need to be output-ed
		++i;
	}

	gettimeofday(&timVal, NULL);
	double EndTime = timVal.tv_sec + (timVal.tv_usec/1000000.0);

	cout << "Time Taken = " << EndTime - StartTime << " seconds." << std::endl;
	
	ofstream outputTime(argv[6]);
	outputTime << (EndTime - StartTime) << " seconds." << std::endl;

	return 0;
}
