Option to print trajectory every N steps.

diff --git a/ss.cpp b/ss.cpp
index e6f4138..e10176d 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -158,7 +158,7 @@ void run(RunParameters_SS &r, unsigned seed) {
     else
         initialize_Ising(r,species,reactions,print_spec);
 
-    simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, print_spec, r.useVerbose);
+    simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, r.step_sample_interval, print_spec, r.useVerbose);
 
     // clean up reaction array
     for (Rxn_parray::iterator  it = reactions.begin(),
@@ -186,7 +186,7 @@ void run(RunParameters_SS &r, unsigned seed) {
 
 void print_species_header(const Species_parray &print_spec)
 {
-    std::cout << "time";
+    std::cout << "time\tstep";
     for (Species_parray::const_iterator spec = print_spec.begin(),
             end = print_spec.end();
             spec != end; ++spec)
@@ -197,9 +197,9 @@ void print_species_header(const Species_parray &print_spec)
     std::cout << '\n';
 }
 
-void print_species_traj(double t, const Species_parray &print_spec, bool verbose)
+void print_species_traj(double t, long step, const Species_parray &print_spec, bool verbose)
 {
-    std::cout << t;
+    std::cout << step << '\t' << t;
     for (Species_parray::const_iterator spec = print_spec.begin(),
                                          end = print_spec.end();
             spec != end; ++spec)
@@ -244,7 +244,7 @@ void print_species_traj(double t, const Species_parray &print_spec, bool verbose
 }
 
 
-void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec, bool verbose)
+void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, unsigned step_sample_interval, const Species_parray &print_spec, bool verbose)
 {
     long step = 0;
 
@@ -252,11 +252,12 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     double dt = 0.0;
 
     double t_next_sample = t + sample_interval;
+    long step_next_sample = step + step_sample_interval;
 
     double total_propensity = 0.0;
 
     print_species_header(print_spec);
-    print_species_traj(t,print_spec,verbose);   // XXX
+    print_species_traj(t,step,print_spec,verbose);   // XXX
 
     Rxn_parray::iterator rxn, end;
 
@@ -300,13 +301,18 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         // print copy numbers
         if (t_next_sample <= t) {
             t_next_sample += sample_interval;
-            print_species_traj(t,print_spec,verbose);
+            print_species_traj(t,step,print_spec,verbose);
+        }
+
+        if ((step_sample_interval != 0) && (step_next_sample <= step)) {
+            step_next_sample += step_sample_interval;
+            print_species_traj(t,step,print_spec,verbose);
         }
 
     }
 
     if (t_end == HUGE_VAL)
-        print_species_traj(t,print_spec,verbose);
+        print_species_traj(t,step,print_spec,verbose);
 
 }
 
@@ -783,6 +789,7 @@ void usage()
 " -e  (int/double)    time to end simulation\n"
 " -t  (int/double)    sampling time interval\n"
 " -E  (int)           number of steps (default 0 = no limit)\n"
+" -T  (int)           sampling step interval (default 0 = print on time only)\n"
 " -mu (double)        mutation rate\n"
 " -b  (double)        \"inverse temperature\" multiplier\n"
 " -v                  flag for verbose output\n"
@@ -839,6 +846,7 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-e")==0)  { if (++i==argc) break; else r.t_end=strtodouble(argv[i]); }
         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],"-T")==0)  { if (++i==argc) break; else r.step_sample_interval=strtouint(argv[i]); }
 
         else if (strcmp(argv[i],"-v")==0)   { r.useVerbose=true;     }
 

diff --git a/ss.h b/ss.h
index d971f19..3a8f93c 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -151,6 +151,7 @@ public:
     double t_end;
     double sample_interval;
     long max_steps;
+    unsigned step_sample_interval;
 
 
     std::vector<std::string> epitopefiles;
@@ -197,6 +198,7 @@ public:
 
         t_end = HUGE_VAL;
         sample_interval = HUGE_VAL;
+        step_sample_interval = 0;
         max_steps = 0;
 
         num_CTL_clones = 0;
@@ -270,7 +272,7 @@ void importState(RunParameters_SS &r);
 void importState_Potts(RunParameters_SS &r);
 void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec);
 void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, NTVirusSpecies *V, Species_parray &print_spec);
-void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec, bool verbose);
+void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, unsigned step_sample_interval, const Species_parray &print_spec, bool verbose);
 void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);
 void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);