Resolve merge conflicts

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 NRMConstant_v5.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>
@@ -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/particle.cpp

@@ -17,10 +17,10 @@ See the file LICENSE.txt for details.
 #include <ANN/ANN.h>
 
 namespace Spatialpy{
-    ParticleSystem::ParticleSystem(size_t num_types, size_t num_chem_species, size_t num_chem_rxns, 
-                         size_t num_stoch_species, size_t num_stoch_rxns,size_t num_data_fn): 
-                            num_types(num_types), num_chem_species(num_chem_species), 
-                            num_chem_rxns(num_chem_rxns), num_stoch_species(num_stoch_species), 
+    ParticleSystem::ParticleSystem(size_t num_types, size_t num_chem_species, size_t num_chem_rxns,
+                         size_t num_stoch_species, size_t num_stoch_rxns,size_t num_data_fn):
+                            num_types(num_types), num_chem_species(num_chem_species),
+                            num_chem_rxns(num_chem_rxns), num_stoch_species(num_stoch_species),
                             num_stoch_rxns(num_stoch_rxns), num_data_fn(num_data_fn){
         dimension = 3;
         boundary_conditions[0] = 'n';
@@ -38,7 +38,7 @@ namespace Spatialpy{
     ParticleSystem::~ParticleSystem(){}
 
     Particle::Particle(ParticleSystem *sys, unsigned int id):sys(sys), id(id){
-    	nu = 0.01; 
+    	nu = 0.01;
     	mass = 1;
     	rho = 1;
     	solidTag = 0;
@@ -57,7 +57,7 @@ namespace Spatialpy{
 
     void Particle::check_particle_nan(){
         if(
-            isnan(x[0]) || !isfinite(x[0]) || 
+            isnan(x[0]) || !isfinite(x[0]) ||
             isnan(x[1]) || !isfinite(x[1]) ||
             isnan(x[2]) || !isfinite(x[2]) ||
             isnan(v[0]) || !isfinite(v[0]) ||
@@ -200,7 +200,7 @@ namespace Spatialpy{
     	//TODO: empty_neighbor_list(neighbors);
     	// search for points forward: (assumes the list is sorted ascending)
     	//printf("find_neighbors.search forward\n");
-    	
+
         // ANN KD Tree k fixed radius nearest neighbor search
         ANNpoint queryPt = annAllocPt(system->dimension);
         for(int i = 0; i < system->dimension; i++) {
@@ -211,7 +211,7 @@ namespace Spatialpy{
         // Number of neighbors to search for
         int k;
         if(use_exact_k) {
-            k = get_k__exact(queryPt, dist, system->kdTree);    
+            k = get_k__exact(queryPt, dist, system->kdTree);
         }else{
             k = get_k__approx(system);
         }
@@ -255,7 +255,7 @@ namespace Spatialpy{
     	// 	//TODO: Should these be breaks instead of continue?
     	//     if( (n.data.x[2] > (x[2] + system.h)) || (n.data.x[2] < (x[2] - system.h) ) ) continue;
     	//     add_to_neighbor_list(n, system);
-    	//     if(debug_flag>2){ 
+    	//     if(debug_flag>2){
     	// 	printf("find_neighbors(%i) forward found %i dist: %e    dx: %e   dy: %e   dz: %e\n",
     	// 	    id,n.data.id, particle_dist(n.data),
     	// 	    x[0] - n.data.x[0],
@@ -272,11 +272,11 @@ namespace Spatialpy{
     	//     if( (n->data.x[1] > (x[1] + system.h)) || (n->data.x[1] < (x[1] - system.h) ) ){
     	// 	continue;
     	//     }
-    	//     if( (n->data.x[2] > (x[2] + system.h)) || (n->data.x[2] < (x[2] - system.h) ) ){ 
+    	//     if( (n->data.x[2] > (x[2] + system.h)) || (n->data.x[2] < (x[2] - system.h) ) ){
     	// 	continue;
     	//     }
     	//     add_to_neighbor_list(n->data, system);
-    	//     if(debug_flag>2){ 
+    	//     if(debug_flag>2){
     	// 	printf("find_neighbors(%i) backwards found %i dist: %e    dx: %e   dy: %e   dz: %e\n",
     	// 	    id,n->data.id, particle_dist(n->data),
     	// 	    x[0] - n->data.x[0],
@@ -286,4 +286,3 @@ namespace Spatialpy{
     	// }
 
 }
-

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);
 
@@ -533,8 +547,7 @@ namespace Spatialpy{
         double totrate,cum,rdelta,rrdelta;
         int event,errcode = 0;
         long unsigned int re, spec = 0 ;
-        Particle& subvol;
-        int subvol_index;
+        Particle* subvol;
         size_t i,j = 0;
         //double old_rrate = 0.0,old_drate = 0.0;
         double rand1,rand2,cum2,old;
@@ -565,7 +578,7 @@ namespace Spatialpy{
 
 
         std::pair<double,int> timeRxnPair;
-        int reactionIndex;
+        int reactionIndex, subvol_index;
         /* Main loop. */
         while(tt <= end_time){
 
@@ -580,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);}