CPP KD Tree Review Fixes

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>
@@ -142,7 +140,6 @@ class NRMConstant_v5 {
     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)
@@ -165,10 +162,35 @@ class NRMConstant_v5 {
         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;
+            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();
             }
@@ -188,6 +210,7 @@ class NRMConstant_v5 {
                 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;

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

@@ -13,8 +13,8 @@ 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
+
 #include "NRMConstant_v5.hpp"
 
 extern int debug_flag ;

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

@@ -315,8 +315,21 @@ namespace Spatialpy{
         double propensitySum=0.0;
         std::size_t activeChannels=0; // initial number of nonzero propensities (typically improves initial performance)
 
-
-
+        long unsigned int p_i = 0 ;
+        for(auto p = system->particles.begin(); p!=system->particles.end(); p++){
+            propensities[p_i] = p->srrate + p->sdrate;
+            propensitySum += propensities[p_i];
+            if(propensities[p_i] > 0){
+                activeChannels++ ;
+            }
+            if(p->particle_index != p_i){
+                // particles can move around in the sys->particles vector,
+                // this ensures that we know where in the vector each particle is
+                p->particle_index = p_i;
+            }
+            p_i++ ;
+        }
+        /**
         for(long unsigned int i=0; i<system->particles.size(); i++){
         //for (i = 0; i < rdme->Ncells; i++) {
             //rdme->rtimes[i] = -log(1.0-dsfmt_genrand_close_open(&dsfmt))/(rdme->srrate[i]+rdme->sdrate[i]);
@@ -341,6 +354,7 @@ namespace Spatialpy{
 
             //system->event_v.emplace_back(p, tt) ;
         }
+    **/
         //initialize_heap(rdme->rtimes,rdme->node,rdme->heap,rdme->Ncells);
         // ordered_list_sort(system->heap);
 
@@ -565,7 +579,7 @@ namespace Spatialpy{
 
 
         std::pair<double,int> timeRxnPair;
-        int reactionIndex;
+        int reactionIndex, subvol_index;
         /* Main loop. */
         while(tt <= end_time){
 
@@ -580,7 +594,7 @@ namespace Spatialpy{
             timeRxnPair = system->rdme_event_q.selectReaction();
             tt = timeRxnPair.first;
             subvol_index = timeRxnPair.second;
-            subvol = system->particles[subvol_index];
+            subvol = &system->particles[subvol_index];
             vol = (subvol->mass / subvol->rho);
 
             if(debug_flag){printf("nsm: tt=%e subvol=%i\n",tt,subvol->id);}