updated some pointer / loop architecture

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

@@ -165,10 +165,34 @@ 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 +212,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/src/particle.cpp

@@ -53,7 +53,7 @@ namespace Spatialpy{
     }
 
     NeighborNode::NeighborNode(Particle *data, double dist, double dWdr, double D_i_j):data(data), dist(dist), dWdr(dWdr), D_i_j(D_i_j){}
-    EventNode::EventNode(Particle *data, double tt):data(data), tt(tt){}
+    //EventNode::EventNode(Particle *data, double tt):data(data), tt(tt){}
 
     void Particle::check_particle_nan(){
         if(

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

@@ -70,8 +70,8 @@ namespace Spatialpy{
     }
     /**************************************************************************/
     void destroy_rdme(ParticleSystem*system){
-        if(debug_flag) printf("NSM: total # reacton events %lu\n",system->rdme->total_reactions);
-        if(debug_flag) printf("NSM: total # diffusion events %lu\n",system->rdme->total_diffusion);
+        if(debug_flag) printf("NSM: total # reacton events %lu\n",system->total_reactions);
+        if(debug_flag) printf("NSM: total # diffusion events %lu\n",system->total_diffusion);
     }
 
 
@@ -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]);
@@ -333,14 +346,15 @@ namespace Spatialpy{
             if(propensities[i] > 0){
                 activeChannels++
             }
-            if(p.particle_index != i){
+            if(p->particle_index != 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 = i;
+                p->particle_index = i;
             }
 
             //system->event_v.emplace_back(p, tt) ;
         }
+    **/
         //initialize_heap(rdme->rtimes,rdme->node,rdme->heap,rdme->Ncells);
         // ordered_list_sort(system->heap);
 
@@ -533,7 +547,7 @@ namespace Spatialpy{
         double totrate,cum,rdelta,rrdelta;
         int event,errcode = 0;
         long unsigned int re, spec = 0 ;
-        Particle*subvol;
+        Particle* subvol;
         size_t i,j = 0;
         //double old_rrate = 0.0,old_drate = 0.0;
         double rand1,rand2,cum2,old;
@@ -564,7 +578,7 @@ namespace Spatialpy{
 
 
         std::pair<double,int> timeRxnPair;
-        int reactionIndex;
+        int reactionIndex, subvol_index;
         /* Main loop. */
         while(tt <= end_time){
 
@@ -579,7 +593,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);}
@@ -855,7 +869,7 @@ namespace Spatialpy{
 //                }
                 //ordered_list_bubble_up_down(system->heap, dest_subvol->heap_index);
 
-                system->rdme_event_q.update(dest_subvol.particle_index, totrate, tt, rng);
+                system->rdme_event_q.update(dest_subvol->particle_index, totrate, tt, rng);
             }
 
             // re-sort the heap