Multiple runs, more rate args, end simulation when all reactions exhausted.

diff --git a/ss.cpp b/ss.cpp
index 094bba7..413e379 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -64,12 +64,10 @@ void run(RunParameters_SS &r, unsigned seed) {
     Species_parray species;
     Rxn_parray reactions;
 
-    //std::string couplingsFile = "neutral_2site.j";
     Hamiltonian H(r.couplingsInfile);
     H.set_temp(r.bh, r.bJ);
-    //H.set_temp(0.01);
-    //double mu = 6.0e-5;
-    //double mu = 1.0e-1;
+
+    for (unsigned n = 0; n < r.num_runs; n++) {
 
     // initialize virus
     VirusSpecies *s1 = new VirusSpecies("I");
@@ -95,26 +93,19 @@ void run(RunParameters_SS &r, unsigned seed) {
 
     // V -> V + V'
     VirusReaction* r1 = new VirusReaction(r.rate_s);  // s_k
-    //r1->rate_constant = 1.5;
     r1->mu = r.mu;
     r1->H = H;
     r1->V = s1;
     reactions.push_back(r1);
     // V -> 0
     VirusDeathReaction* r2 = new VirusDeathReaction(r.rate_u);  // u
-    //r2->rate_constant = 0.5;
     r2->H = H;
     r2->V = s1;
     reactions.push_back(r2);
 
 
-    //double t_end = 1000.0;
-    //long max_steps = LONG_MAX;
-    //double sample_interval = 1.0;
     simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, print_spec, r.useVerbose);
 
-    gsl_rng_free(rnd);  //Free up memory from random number generator
-
     // clean up reaction array
     for (Rxn_parray::iterator  it = reactions.begin(),
                               end = reactions.end();
@@ -130,6 +121,11 @@ void run(RunParameters_SS &r, unsigned seed) {
         delete (*it);
     }
     species.clear();
+    print_spec.clear();
+
+    }  // num_runs
+
+    gsl_rng_free(rnd);  //Free up memory from random number generator
 
 }
 
@@ -186,6 +182,7 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     double total_propensity = 0.0;
 
     print_species_header(print_spec);
+    print_species_traj(t,print_spec,verbose);   // XXX
 
     Rxn_parray::iterator rxn, end;
 
@@ -195,8 +192,8 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     }
 
 
-    // stop if time runs out or too many steps
-    for (; (max_steps == 0 || step < max_steps) && (t < t_end); ++step) {
+    // stop if time runs out or too many steps or no more reactions will occur
+    for (; (max_steps == 0 || step < max_steps) && (t < t_end) && (total_propensity > 0.0); ++step) {
 
         double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity;  // uniform on (0,R]
         double a_sum = 0.0;
@@ -466,10 +463,9 @@ void usage()
 " -E  (int)           number of steps (default 0 = no limit)\n"
 " -mu (double)        mutation rate\n"
 " -b  (double)        \"inverse temperature\" multiplier\n"
-" -esc                flag to run until escape observed (or max number of generations reached)\n"
-" -esc_all            flag to run until escape from *all* immune pressure (or max num gen)\n"
 " -v                  flag for verbose output\n"
 " -numruns (int)      number of trajectories to simulate\n"
+" -rate_X (double)    X in [s,u,p,a,r,d,dprime,b,e,w,aprime,g,h]\n"
 ;   std::cout.flush();
 }
 
@@ -519,14 +515,14 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-E")==0)  { if (++i==argc) break; else r.max_steps=(long)strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-t")==0)  { if (++i==argc) break; else r.sample_interval=strtodouble(argv[i]); }
 
-        else if (strcmp(argv[i],"-esc")==0) { r.runUntilEscape=true; }
-        else if (strcmp(argv[i],"-esc_all")==0) { r.runUntilEscape_all=true; }
         else if (strcmp(argv[i],"-v")==0)   { r.useVerbose=true;     }
 
         else if (strcmp(argv[i],"-numruns")==0)  { if (++i==argc) break; else r.num_runs=(unsigned int)strtodouble(argv[i]); }
 
         else if (strcmp(argv[i],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
 
+        else if (strcmp(argv[i],"-rate_s")==0) { if (++i==argc) break; else r.rate_s=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_u")==0) { if (++i==argc) break; else r.rate_u=strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-rate_p")==0) { if (++i==argc) break; else r.rate_p=strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-rate_a")==0) { if (++i==argc) break; else r.rate_a=strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-rate_r")==0) { if (++i==argc) break; else r.rate_r=strtodouble(argv[i]); }

diff --git a/ss.h b/ss.h
index da25a2e..5738ecd 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -163,7 +163,7 @@ public:
         directory="";
 
         n  = 100000;
-        num_runs = 1000;
+        num_runs = 1;
         mu = 6.0e-5;
         bh = 1.0;
         bJ = 1.0;