Fix turing pattern example

spatialpy/Model.py

@@ -736,8 +736,12 @@ def create_mass_action(self):
         for r in self.reactants:
             total_stoch += self.reactants[r]
         if total_stoch > 2:
-            raise ReactionError("Reaction: A mass-action reaction cannot involve more than two of one species or one "
-                                "of two species.")
+            raise ReactionError(
+                      """Reaction: A mass-action reaction cannot involve more than two of one species or one "
+                         of two species (no more than 2 total reactants).
+                         SpatialPy support zeroth, first and second order propensities only.
+                         There is no theoretical justification for higher order propensities.
+                         Users can still create such propensities using a 'custom propensity'.""")
         # Case EmptySet -> Y
 
         if isinstance(self.marate, str):

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

@@ -93,11 +93,6 @@ namespace Spatialpy{
                    std::mt19937_64& rng, double timeOffset=0.0,
                    double simulationEndTime=std::numeric_limits<double>::max()) ;
 
-//        void build(std::vector<Spatialpy::Particle>& particles,
-//                   double propensitySum, std::size_t activeChannels,
-//                   std::mt19937_64& rng, 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
@@ -110,7 +105,7 @@ namespace Spatialpy{
         int computeBinIndex(double firingTime);
         //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);
+        bool setBinNumberAndBounds(double newLowerBound, double propensitySum, int activeChannels);
         std::size_t activeChannelCounter;
 
     };

spatialpy/ssa_sdpd-c-simulation-engine/src/NRMConstant_v5.cpp

@@ -47,7 +47,10 @@ namespace Spatialpy{
             return false;
         }
 
-        setBinNumberAndBounds(timeOffset,propensitySum,activeChannels); //choose lowerBound, upperBound, and numberOfBins
+        if(!setBinNumberAndBounds(timeOffset,propensitySum,activeChannels)){ //choose lowerBound, upperBound, and numberOfBins
+            std::cerr << "ERROR: setBinNumberAndBound failed in build()" << std::endl;
+            return false;
+        }
         minBin=0;
 
         std::vector<std::pair<double, std::size_t> > emptyVector;
@@ -90,83 +93,14 @@ namespace Spatialpy{
 
         if (active_channel_count!=activeChannelCounter) {
             std::cout << "ERROR: active channel count is inconsistent.\n";
-            exit(1);
+            //exit(1);
+            return false;
         }
 
     //  std::cout << "exiting build...\n";
         return true;
     }
 
-//    void NRMConstant_v5::build(std::vector<Spatialpy::Particle>& particles,
-//           double propensitySum, std::size_t activeChannels,
-//           std::mt19937_64& rng, double timeOffset,
-//           double simulationEndTime) {
-//    //  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(rng)/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";
-//
-//    }
 
     // reactionIndex is particle id
     void NRMConstant_v5::update(std::size_t reactionIndex, double newPropensity, double currentTime, std::mt19937_64& rng) {
@@ -240,8 +174,9 @@ namespace Spatialpy{
         while (theHashTable[minBin].size()==0) {
             ++minBin;
             if (minBin==theHashTable.size()) {
-            //  std::cout << "rebuilding...\n";
+                //std::cout << "rebuilding...\n";
                 if (!rebuild()) {
+                    //std::cerr << "in selectReaction, rebuild() returned false" << std::endl;
                     return std::make_pair(std::numeric_limits<double>::max(),-1);
                 }
             }
@@ -303,7 +238,11 @@ namespace Spatialpy{
             //  std::cout << "previousBinWidth was " << previousBinWidth <<"\n";
         }
 
-        setBinNumberAndBounds(currentUpperBound,propensitySumEstimate,activeChannelCounter);//
+        if(!setBinNumberAndBounds(currentUpperBound,propensitySumEstimate,activeChannelCounter)){//
+            //std::cerr << "ERROR: setBinNumberAndBound failed in rebuild()" << std::endl;
+            minBin=0;
+            return false;
+        }
         minBin=0;
 
         std::vector<std::pair<double, std::size_t> > emptyVector;
@@ -348,17 +287,19 @@ namespace Spatialpy{
 
     //fixed number of bins
     //currentUpperBound=
-    void NRMConstant_v5::setBinNumberAndBounds(double newLowerBound, double propensitySum, int activeChannels) {
+    bool NRMConstant_v5::setBinNumberAndBounds(double newLowerBound, double propensitySum, int activeChannels) {
         if (activeChannels==0) {
             std::cout << "ERROR: setBinNumberAndBounds not set up to handle activeChannels=0\n";
-            exit(1);
+            //exit(1);
+            return false;
         }
 
         //  std::cout << "...in setBinNumberAndBounds...newLowerBound=" << newLowerBound << std::endl;
         currentLowerBound=newLowerBound;
         if (currentLowerBound>endTime) {
-            std::cerr << "ERROR: calling rebuild when simulation end time exceeded. Terminating. newLowerBound="<<newLowerBound<<", endTime="<<endTime<<"\n";
-            exit(1);
+            //std::cerr << "ERROR: calling rebuild when simulation end time exceeded. Terminating. newLowerBound="<<newLowerBound<<", endTime="<<endTime<<"\n";
+            //exit(1);
+            return false;
         }
 
         //  std::cout << "propensitySum is " << propensitySum << ", so step size is " << 1.0/propensitySum << "\n";
@@ -398,5 +339,6 @@ namespace Spatialpy{
     //  std::cout << "currentLowerBound=" << currentLowerBound << ", currentUpperBound=" << currentUpperBound << std::endl;
     //  std::cout << "binWidth=" << binWidth << std::endl;
     //  std::cout << "numberOfBins=" << numberOfBins << std::endl;
+        return true;
     }
 }

spatialpy/ssa_sdpd-c-simulation-engine/src/simulate_rdme.cpp

@@ -223,16 +223,35 @@ namespace Spatialpy{
         Particle *dest_subvol = NULL;
         NeighborNode *nn = NULL;
 
+        if(debug_flag>1){
+            printf("take_step() tt=%e end_time=%e\n",tt,end_time);
+            fflush(stdout);
+        }
+
         std::pair<double,int> timeRxnPair;
         int reactionIndex, subvol_index;
         /* Main loop. */
         while(tt <= end_time){
             /* Get the subvolume in which the next event occurred.
              This subvolume is on top of the heap. */
 
+            if(debug_flag>1){
+                printf("selectReaction()\n");
+                fflush(stdout);
+            }
             timeRxnPair = system->rdme_event_q.selectReaction();
             tt = timeRxnPair.first;
             subvol_index = timeRxnPair.second;
+            if(debug_flag>1){
+                printf("selectReaction: subvol_index=%d tt=%e\n",subvol_index,tt);
+                fflush(stdout);
+            }
+            if(subvol_index == -1){ // catch special case of an empty heap
+                //printf("ending take_step, subvol_index=%d tt=%e\n",subvol_index,tt);
+                //fflush(stdout);
+                //debug_flag = 2;
+                return;
+            }
             subvol = &system->particles[subvol_index];
             vol = (subvol->mass / subvol->rho);
             if(debug_flag){printf("nsm: tt=%e subvol=%i\n",tt,subvol->id);}