still trying to compile, moved some things around including template …

…functions to cpp

spatialpy/ssa_sdpd-c-simulation-engine/build/Makefile

@@ -4,7 +4,7 @@ SIMFLAGS = -l. -std=c++14 -Wall -O3
 LFLAGS = -pthread -lm
 INCDIR = $(ROOTINC)/include/ $(ROOTINC)/external/ $(ROOTINC)/external/ANN/include/
 INCDIRPARAMS = $(INCDIR:%=-I%)
-OBJ = particle.o simulate.o count_cores.o output.o simulate_rdme.o simulate_threads.o model.o pthread_barrier.o
+OBJ = NRMConstant_v5.o particle.o simulate.o count_cores.o output.o simulate_rdme.o simulate_threads.o model.o pthread_barrier.o
 ANNOBJECTS = ANN.o brute.o kd_tree.o kd_util.o kd_split.o \
 	kd_dump.o kd_search.o kd_pr_search.o kd_fix_rad_search.o \
 	bd_tree.o bd_search.o bd_pr_search.o bd_fix_rad_search.o \

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

@@ -18,8 +18,6 @@
 // removes vectors and replace with hashmap (unordered_map) to allow adding or removing "channels"
 // removed PROFILE code (look at old version for old profiling data)
 
-#include "particle.hpp"
-
 #include <iostream>
 #include <vector>
 #include <random>
@@ -31,6 +29,7 @@
 #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)
@@ -69,227 +68,27 @@ 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()) {
-//        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";
-    }
+               double simulationEndTime=std::numeric_limits<double>::max()) ;
 
     // 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()) {
-        //        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 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 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) {
-//        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;
-            }
-        }
-
-    }
+    void update(std::size_t reactionIndex, double newPropensity, double currentTime, generatorType& generator) ;
 
 	virtual ~NRMConstant_v5();
 
@@ -298,7 +97,6 @@ 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;
 
 };
 

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

@@ -13,8 +13,9 @@ See the file LICENSE.txt for details.
 #include <queue>
 #include <memory>
 #include "propensities.hpp"
-// Include ANN KD Tree
-#include <ANN/ANN.h>
+#include <ANN/ANN.h> // ANN KD Tree
+
+class NRMConstant_v5 ;
 #include "NRMConstant_v5.hpp"
 
 extern int debug_flag ;

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

@@ -9,6 +9,7 @@ See the file LICENSE.txt for details.
 #define simulate_rdme_h
 #include "particle.hpp"
 #include "propensities.hpp"
+#include "NRMConstant_v5.hpp"
 #include <random>
 
 namespace Spatialpy{