undo bogus attempt to move template function definition to cpp

spatialpy/ssa_sdpd-c-simulation-engine/include/NRMConstant_v5.hpp

@@ -29,7 +29,6 @@
 #include <chrono>
 #include <limits>
 
-using namespace Spatialpy ;
 class NRMConstant_v5 {
     typedef double firingTimeT;
     typedef std::size_t reactionIndexT;//assumes channels have unique nonnegative ID (within range of size_t)
@@ -68,27 +67,227 @@ class NRMConstant_v5 {
 	NRMConstant_v5();
 
     // activeChannels is count of nonzero propensities
-    std::size_t activeChannelCounter;
     template <typename generatorType>
+
     void build(std::vector<double>& propensities, generatorType& generator,
                double propensitySum, std::size_t activeChannels,
                double timeOffset=0.0, 
-               double simulationEndTime=std::numeric_limits<double>::max()) ;
+               double simulationEndTime=std::numeric_limits<double>::max()) {
+//        std::cout << "in build..." << std::endl;
+
+        activeChannelCounter=activeChannels;
+        endTime=simulationEndTime;
+        previousFiringTime=timeOffset;//initialize (could set to nan or something instead)
+
+//        nextFiringTime.resize(propensities.size());
+//        binIndexAndPositionInBin.resize(propensities.size());
+
+        if (propensitySum==0.0) {
+            std::cerr << "ERROR: propensity sum is zero on initial build. Terminating." << std::endl;
+            exit(1);
+        }
+
+        setBinNumberAndBounds(timeOffset,propensitySum,activeChannels);//choose lowerBound, upperBound, and numberOfBins
+        minBin=0;
+
+        std::vector<std::pair<double, std::size_t> > emptyVector;
+
+        theHashTable.resize(numberOfBins,emptyVector);
+
+        std::size_t active_channel_count=0; // for verifying build input
+
+        for (std::size_t i=0; i!=propensities.size(); ++i) {
+            double firingTime;
+            if (propensities[i]==0.0) {
+                firingTime=std::numeric_limits<double>::infinity();
+            }
+            else {
+                firingTime=exponential(generator)/propensities[i]+timeOffset;
+                ++active_channel_count;
+            }
+
+            nextFiringTime[i]=firingTime;//.insert(std::make_pair<reactionIndexT,firingTimeT>(i,firingTime));
+//            std::cout << "nextFiringTime["<<i<<"]="<<nextFiringTime[i]<<"\n";
+            //insert into hashTable
+            int bin=computeBinIndex(firingTime);
+            if (bin>=0) {
+                theHashTable[bin].push_back(std::make_pair(firingTime,i));//place this rxn at back of bin
+                binIndexAndPositionInBin[i]=std::make_pair(bin,theHashTable[bin].size()-1);
+            } else {
+                binIndexAndPositionInBin[i]=std::make_pair<int,int>(-1,-1);// bin (and index within bin) is -1 if not in the hash table
+            }
+        }
+
+        //set rxn counter to 0
+        rxnCountThisBuildOrRebuild=0;
+
+        //record info from build
+        rxnsBetweenBuildOrRebuild.push_back(0);
+        lowerBoundsBuildOrRebuild.push_back(currentLowerBound);//keep entry for each build/rebuild; lowerBound.back() corresponds to currentLowerBound
+        upperBoundsBuildOrRebuild.push_back(currentUpperBound);
+
+//        printHashTable();
+
+        if (active_channel_count!=activeChannelCounter) {
+            std::cout << "ERROR: active channel count is inconsistent.\n";
+            exit(1);
+        }
+
+//        std::cout << "exiting build...\n";
+    }
 
     // activeChannels is count of nonzero propensities
     template <typename generatorType>
-    void build(std::vector<Spatialpy::Particle> particles, generatorType& generator, 
-               double propensitySum, std::size_t activeChannels, 
-               double timeOffset=0.0, 
-               double simulationEndTime=std::numeric_limits<double>::max()) ;
+    void build(std::vector<Spatialpy::Particle> particles, generatorType& generator, double propensitySum, std::size_t activeChannels, double timeOffset=0.0, double simulationEndTime=std::numeric_limits<double>::max()) {
+        //        std::cout << "in build..." << std::endl;
+        activeChannelCounter=activeChannels;
+        endTime=simulationEndTime;
+        previousFiringTime=timeOffset;//initialize (could set to nan or something instead)
+
+        //        nextFiringTime.resize(propensities.size());
+        //        binIndexAndPositionInBin.resize(propensities.size());
+
+        if (propensitySum==0.0) {
+            std::cerr << "ERROR: propensity sum is zero on initial build. Terminating." << std::endl;
+            exit(1);
+        }
+
+        setBinNumberAndBounds(timeOffset,propensitySum,activeChannels);//choose lowerBound, upperBound, and numberOfBins
+        minBin=0;
+
+        std::vector<std::pair<double, std::size_t> > emptyVector;
+
+        theHashTable.resize(numberOfBins,emptyVector);
+
+        std::size_t active_channel_count=0; // for verifying build input
+
+        std::cout << "?? is index into particles vector a valid unique id?\n";
+
+        for(auto p = particles.begin(); p!=particles.end(); p++){
+            double firingTime;
+            double this_propensity = p->srrate+p->sdrate;
+            unsigned this_id = p->id;
+            if (this_propensity==0.0) {
+                firingTime=std::numeric_limits<double>::infinity();
+            }
+            else {
+                firingTime=exponential(generator)/this_propensity+timeOffset;
+                ++active_channel_count;
+            }
+
+            nextFiringTime[this_id]=firingTime;//.insert(std::make_pair<reactionIndexT,firingTimeT>(i,firingTime));
+            //            std::cout << "nextFiringTime["<<i<<"]="<<nextFiringTime[i]<<"\n";
+            //insert into hashTable
+            int bin=computeBinIndex(firingTime);
+            if (bin>=0) {
+                theHashTable[bin].push_back(std::make_pair(firingTime,this_id));//place this rxn at back of bin
+                binIndexAndPositionInBin[this_id]=std::make_pair(bin,theHashTable[bin].size()-1);
+            } else {
+                binIndexAndPositionInBin[this_id]=std::make_pair<int,int>(-1,-1);// bin (and index within bin) is -1 if not in the hash table
+            }
+        }
+        /**
+        for (std::size_t i=0; i!=particles.size(); ++i) {
+            double firingTime;
+            double this_propensity = particles[i]->srrate+particles[i]->sdrate;
+            unsigned this_id = particles[i]->id;
+            if (this_propensity==0.0) {
+                firingTime=std::numeric_limits<double>::infinity();
+            }
+            else {
+                firingTime=exponential(generator)/this_propensity+timeOffset;
+                ++active_channel_count;
+            }
+
+            nextFiringTime[this_id]=firingTime;//.insert(std::make_pair<reactionIndexT,firingTimeT>(i,firingTime));
+            //            std::cout << "nextFiringTime["<<i<<"]="<<nextFiringTime[i]<<"\n";
+            //insert into hashTable
+            int bin=computeBinIndex(firingTime);
+            if (bin>=0) {
+                theHashTable[bin].push_back(std::make_pair(firingTime,this_id));//place this rxn at back of bin
+                binIndexAndPositionInBin[this_id]=std::make_pair(bin,theHashTable[bin].size()-1);
+            } else {
+                binIndexAndPositionInBin[this_id]=std::make_pair<int,int>(-1,-1);// bin (and index within bin) is -1 if not in the hash table
+            }
+        }
+        **/
+
+        //set rxn counter to 0
+        rxnCountThisBuildOrRebuild=0;
+
+        //record info from build
+        rxnsBetweenBuildOrRebuild.push_back(0);
+        lowerBoundsBuildOrRebuild.push_back(currentLowerBound);//keep entry for each build/rebuild; lowerBound.back() corresponds to currentLowerBound
+        upperBoundsBuildOrRebuild.push_back(currentUpperBound);
+
+        //        printHashTable();
+
+        if (active_channel_count!=activeChannelCounter) {
+            std::cout << "ERROR: active channel count is inconsistent.\n";
+            exit(1);
+        }
+
+        //        std::cout << "exiting build...\n";
+
+    }
 
     void printHashTable();//for debug/dev
 
     timeIndexPair selectReaction(); // returns pair of firing time, rxn index
 
     // reactionIndex is particle id
     template <typename generatorType>
-    void update(std::size_t reactionIndex, double newPropensity, double currentTime, generatorType& generator) ;
+    void update(std::size_t reactionIndex, double newPropensity, double currentTime, generatorType& generator) {
+//        std::cout << "updating reactionIndex=" << reactionIndex << ", newPropensity=" << newPropensity << ", currentTime=" << currentTime << std::endl;
+
+        double firingTime;
+        if (newPropensity==0) {
+            firingTime=std::numeric_limits<double>::infinity();
+            if (nextFiringTime.at(reactionIndex)!=std::numeric_limits<double>::infinity()) {
+                activeChannelCounter--;
+            }
+        }
+        else {
+            firingTime=exponential(generator)/newPropensity+currentTime;
+            if (nextFiringTime.at(reactionIndex)==std::numeric_limits<double>::infinity()) {
+                activeChannelCounter++;
+            }
+        }
+//        std::cout << "new firingTime=" << firingTime << std::endl;
+        nextFiringTime[reactionIndex]=firingTime;
+        int newBin=computeBinIndex(firingTime);
+
+//        std::cout << "newBin=" << newBin << std::endl;
+        int oldBin=binIndexAndPositionInBin[reactionIndex].first;
+//        std::cout << "oldBin=" << oldBin << std::endl;
+        int oldPositionInBin=binIndexAndPositionInBin[reactionIndex].second;
+        if (newBin!=oldBin) {
+            if (oldBin>=0) {
+                //remove from old bin
+                if (theHashTable[oldBin].size()>1) {
+
+                    //take last element in old bin and place in this element's spot
+                    theHashTable[oldBin][oldPositionInBin]=theHashTable[oldBin].back();
+                    std::size_t movedElementIndex=theHashTable[oldBin][oldPositionInBin].second;
+                    //update old last element's ...within bin index
+                    binIndexAndPositionInBin[movedElementIndex].second=oldPositionInBin;
+                }
+                theHashTable[oldBin].pop_back();
+            }
+            binIndexAndPositionInBin[reactionIndex].first=newBin;
+            if (newBin>=0) {
+                theHashTable[newBin].push_back(std::pair<double,std::size_t>(firingTime,reactionIndex));
+                binIndexAndPositionInBin[reactionIndex].second=theHashTable[newBin].size()-1;
+            }
+        }
+        else {
+            //just update firing time
+            if (newBin>=0) {
+                theHashTable[newBin][binIndexAndPositionInBin[reactionIndex].second].first=firingTime;
+            }
+        }
+
+    }
 
 	virtual ~NRMConstant_v5();
 
@@ -97,6 +296,7 @@ class NRMConstant_v5 {
     //std::size_t insertIntoBin(std::size_t binIndex, double firingTime, std::size_t reactionIndex);
     bool rebuild();//returns false if all propensities are 0
     void setBinNumberAndBounds(double newLowerBound, double propensitySum, int activeChannels);
+    std::size_t activeChannelCounter;
 
 };