Output at intervals, remove more WF.

diff --git a/pop_ss.h b/pop_ss.h
index 8c264b5..83cd367 100644 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -28,8 +28,8 @@ class Species {
 public:
     long count;
     std::string name;
-    Species() : count(0) {}
-    Species(long c) : count(c) {}
+    Species() : count(0), name("") {}
+    Species(long c) : count(c), name("") {}
 };
 
 class SimpleSpecies : public Species {
@@ -46,8 +46,8 @@ public:
 
 class CTLSpecies : public Species {
 public:
-    CTLSpecies() : Species() {}
-    CTLSpecies(long c) : Species(c) {}
+    CTLSpecies() : Species(), num_ep(0) {}
+    CTLSpecies(long c) : Species(c), num_ep(0) {}
 
     std::vector<EpitopeRecognizer> epitopes;
     std::vector<double> affinity;

diff --git a/ss.cpp b/ss.cpp
index 085a74e..dbeac56 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -51,6 +51,13 @@ void run(RunParameters_SS &r, unsigned seed) {
     //gsl_rng_set(rnd,rand());
     gsl_rng_set(rnd,seed);
 
+    if (r.t_end == HUGE_VAL && r.max_steps == 0 && r.sample_interval == HUGE_VAL) {
+        std::cerr << "ERROR: must specify one of -e or -t or -E; see -h for usage" << std::endl;
+        exit(1);
+    }
+    if (r.sample_interval == HUGE_VAL)
+        r.sample_interval = r.t_end;  // ensure last time will be printed
+
     r.setFiles();
 
     // master arrays of (pointers to) species and reactions
@@ -65,51 +72,90 @@ void run(RunParameters_SS &r, unsigned seed) {
     //double mu = 1.0e-1;
 
     // initialize virus
-    VirusSpecies s1; s1.count = r.n;
-    s1.pop[Virus(H,r.mu)] = s1.count;   // default mu = 6.0e-5
+    VirusSpecies s1;
     species.push_back(&s1);
+    if (r.importState) {
+        importState(r);      // <infile>.st
+        for (size_t i=0; i<r.initPop.size(); i++) {
+            Virus V(H, r.mu, r.initPop[i]);
+            unsigned int N = (unsigned int) r.n * r.initFrac[i];
+            s1.pop[V] = N;
+            s1.count += N;
+        }
+    }
+    else {
+        s1.count = r.n;
+        s1.pop[Virus(H,r.mu)] = s1.count;   // default mu = 6.0e-5
+    }
 
     Species_parray print_spec;
     print_spec.push_back(&s1);
 
-    if (r.importState) importState(r);      // <infile>.st
     if (r.useEpitope)  importEpitope(r,species,reactions,&s1,print_spec);    // <infile>.ep
 
     // V -> V + V'
-    VirusReaction* r1 = new VirusReaction(1.5);  // s_k
+    VirusReaction* r1 = new VirusReaction(r.rate_s);  // s_k
     //r1->rate_constant = 1.5;
     r1->H = H;
     r1->V = &s1;
     reactions.push_back(r1);
     // V -> 0
-    VirusDeathReaction* r2 = new VirusDeathReaction(0.5);  // u
+    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, t_end, max_steps, sample_interval, print_spec);
+    //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);
 
     gsl_rng_free(rnd);  //Free up memory from random number generator
 }
 
 
+void print_species_header(const Species_parray &print_spec)
+{
+    std::cout << "time";
+    for (Species_parray::const_iterator spec = print_spec.begin(),
+            end = print_spec.end();
+            spec != end; ++spec)
+    {
+        // TODO: prettier representation
+        std::cout << '\t' << (*spec)->name;
+    }
+    std::cout << '\n';
+}
+
+void print_species_traj(double t, const Species_parray &print_spec)
+{
+    std::cout << t;
+    for (Species_parray::const_iterator spec = print_spec.begin(),
+            end = print_spec.end();
+            spec != end; ++spec)
+    {
+        // TODO: prettier representation
+        std::cout << '\t' << (*spec)->count;
+    }
+    std::cout << '\n';
+}
+
+
 void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec)
 {
     long step = 0;
 
     double t = 0.0;
-    double t_remain = t_end;
     double dt = 0.0;
 
     double t_next_sample = t + sample_interval;
 
     double total_propensity = 0.0;
 
+    print_species_header(print_spec);
+
     Rxn_parray::iterator rxn, end;
 
     for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
@@ -118,7 +164,8 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     }
 
 
-    while (1) {
+    // stop if time runs out or too many steps
+    for (; (max_steps == 0 || step < max_steps) && (t < t_end); ++step) {
 
         double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity;  // uniform on (0,R]
         double a_sum = 0.0;
@@ -126,11 +173,6 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         // time to reaction
         dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
         t += dt;
-        t_remain -= dt;
-
-        // stop if time ran out
-        if (t_remain < 0.0)
-            break;
 
         // determine next reaction
         for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
@@ -153,31 +195,16 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
             total_propensity += (*rxn)->propensity;
         }
 
-        if (step >= max_steps)
-            break;
-
-        ++step;
-
-
         // print copy numbers
         if (t_next_sample <= t) {
             t_next_sample += sample_interval;
-
-            std::cout << t << '\t';
-            for (Species_parray::const_iterator spec = print_spec.begin(),
-                                           end = print_spec.end();
-              spec != end; ++spec)
-            {
-                // TODO: prettier representation
-                std::cout << (*spec)->count << '\t';
-            }
-            std::cout << '\n';
+            print_species_traj(t,print_spec);
         }
 
-
     }
 
-
+    if (t_end == HUGE_VAL)
+        print_species_traj(t,print_spec);
 
 }
 
@@ -237,7 +264,7 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
 
     if ((r.num_CTL_clones != r.epitopefiles.size()) || 
         (r.num_CTL_clones != r.init_CTL_num.size())) {
-        perror((std::string("ERROR in importEpitope: invalid number of CTL clones").c_str()));
+        std::cerr << "ERROR in importEpitope: invalid number of CTL clones" << std::endl;
         exit(1);
     }
 
@@ -320,61 +347,61 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
         print_spec.push_back(M);
 
         {   // p
-            KillingReaction* r = new KillingReaction(1e-5);
-            r->T = T;  r->V = V;
-            reactions.push_back(r);
+            KillingReaction* rx = new KillingReaction(r.rate_p);
+            rx->T = T;  rx->V = V;
+            reactions.push_back(rx);
         }
         {   // a
-            CTLActivationReaction* r = new CTLActivationReaction(1e-7);
-            r->Tfrom = N;  r->Tto = T;  r->V = V;
-            reactions.push_back(r);
+            CTLActivationReaction* rx = new CTLActivationReaction(r.rate_a);
+            rx->Tfrom = N;  rx->Tto = T;  rx->V = V;
+            reactions.push_back(rx);
         }
         {   // r
-            ElementaryReaction* r = new ElementaryReaction(4.0);
-            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
-            r->products.push_back(T); r->product_stoich.push_back(1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_r);
+            rx->reactants.push_back(T); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(T); rx->product_stoich.push_back(1);
+            reactions.push_back(rx);
         }
         {   // d = 0.5  // note d' = 3.0
-            ElementaryReaction* r = new ElementaryReaction(3.0);
-            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
-            r->products.push_back(T); r->product_stoich.push_back(-1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_dprime);
+            rx->reactants.push_back(T); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(T); rx->product_stoich.push_back(-1);
+            reactions.push_back(rx);
         }
         {   // b
-            ElementaryReaction* r = new ElementaryReaction(1e-4);
-            r->reactants.push_back(N); r->reactant_stoich.push_back(1);
-            r->products.push_back(N); r->product_stoich.push_back(1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_b);
+            rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(N); rx->product_stoich.push_back(1);
+            reactions.push_back(rx);
         }
         {   // e
-            ElementaryReaction* r = new ElementaryReaction(3e-4);
-            r->reactants.push_back(N); r->reactant_stoich.push_back(1);
-            r->products.push_back(N); r->product_stoich.push_back(-1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_e);
+            rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(N); rx->product_stoich.push_back(-1);
+            reactions.push_back(rx);
         }
         {   // w
-            ElementaryReaction* r = new ElementaryReaction(7.5);
-            r->products.push_back(N); r->product_stoich.push_back(1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_w);
+            rx->products.push_back(N); rx->product_stoich.push_back(1);
+            reactions.push_back(rx);
         }
         {   // a'
-            CTLActivationReaction* r = new CTLActivationReaction(1e-6);
-            r->Tfrom = M;  r->Tto = T;  r->V = V;
-            reactions.push_back(r);
+            CTLActivationReaction* rx = new CTLActivationReaction(r.rate_aprime);
+            rx->Tfrom = M;  rx->Tto = T;  rx->V = V;
+            reactions.push_back(rx);
         }
         {   // g
-            ElementaryReaction* r = new ElementaryReaction(0.03);
-            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
-            r->products.push_back(T); r->product_stoich.push_back(-1);
-            r->products.push_back(M); r->product_stoich.push_back(1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_g);
+            rx->reactants.push_back(T); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(T); rx->product_stoich.push_back(-1);
+            rx->products.push_back(M); rx->product_stoich.push_back(1);
+            reactions.push_back(rx);
         }
         {   // h
-            ElementaryReaction* r = new ElementaryReaction(0.03);
-            r->reactants.push_back(M); r->reactant_stoich.push_back(1);
-            r->products.push_back(M); r->product_stoich.push_back(-1);
-            reactions.push_back(r);
+            ElementaryReaction* rx = new ElementaryReaction(r.rate_h);
+            rx->reactants.push_back(M); rx->reactant_stoich.push_back(1);
+            rx->products.push_back(M); rx->product_stoich.push_back(-1);
+            reactions.push_back(rx);
         }
 
         for (size_t i = 0; i < T->num_ep; ++i) {  // XXX
@@ -388,20 +415,6 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
 }
 
 
-// insert into initial state of population
-
-void add_to_two_site_pop(RunParameters_SS &r, bool s1, bool s2, double frac) {
-
-    r.initFrac.push_back(frac);
-
-    std::set<unsigned int> tempSet;
-    if (s1) tempSet.insert(0);
-    if (s2) tempSet.insert(1);
-    r.initPop.push_back(tempSet);
-
-}
-
-
 void usage()
 {
     std::cout <<
@@ -414,7 +427,9 @@ void usage()
 " -ep (string)        input file containing targeted epitope\n"
 " -en (int)           corresponding number of that CTL (ordered!)\n"
 " -n  (int/double)    population size\n"
-" -g  (int/double)    number of generations\n"
+" -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"
 " -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"
@@ -439,7 +454,6 @@ void usage()
  *  -esc                flag to run until escape observed (or max number of generations reached)
  *  -v                  flag for verbose output
  *  -numruns (int)      number of trajectories to simulate
- *  -2site              flag for two-site two-allele system
  */
 
 int main(int argc, char *argv[]) {
@@ -460,15 +474,15 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-en")==0) { if (++i==argc) break; else { r.init_CTL_num.push_back((long)strtodouble(argv[i])); r.useEpitope=true; } }
         
         else if (strcmp(argv[i],"-n")==0)  { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
-        else if (strcmp(argv[i],"-g")==0)  { if (++i==argc) break; else r.g=(unsigned int)strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-mu")==0) { if (++i==argc) break; else r.mu=strtodouble(argv[i]);              }
         else if (strcmp(argv[i],"-b")==0)  { if (++i==argc) break; else { r.bh=strtodouble(argv[i]); r.bJ=r.bh; } }
         else if (strcmp(argv[i],"-bh")==0) { if (++i==argc) break; else r.bh=strtodouble(argv[i]);              }
         else if (strcmp(argv[i],"-bJ")==0) { if (++i==argc) break; else r.bJ=strtodouble(argv[i]);              }
 
-        else if (strcmp(argv[i],"-penaltyEach")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters_SS::PenaltyEACH; } }
-        else if (strcmp(argv[i],"-penaltyTotal")==0) { if (++i==argc) break; else { r.penalty=strtodouble(argv[i]); r.penaltyType=RunParameters_SS::PenaltyTOTAL; } }
-        
+        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],"-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;     }
@@ -477,7 +491,17 @@ int main(int argc, char *argv[]) {
 
         else if (strcmp(argv[i],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
 
-        else if (strcmp(argv[i],"-write_mod")==0) { if (++i==argc) break; else r.write_mod=(unsigned int)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]); }
+        else if (strcmp(argv[i],"-rate_d")==0) { if (++i==argc) break; else r.rate_d=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_dprime")==0) { if (++i==argc) break; else r.rate_dprime=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_b")==0) { if (++i==argc) break; else r.rate_b=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_e")==0) { if (++i==argc) break; else r.rate_e=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_w")==0) { if (++i==argc) break; else r.rate_w=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_aprime")==0) { if (++i==argc) break; else r.rate_aprime=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_g")==0) { if (++i==argc) break; else r.rate_g=strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-rate_h")==0) { if (++i==argc) break; else r.rate_h=strtodouble(argv[i]); }
 
         else if (strcmp(argv[i],"-h")==0 || strcmp(argv[i],"--help")==0)   { usage(); return 0;    }
         
@@ -491,4 +515,4 @@ int main(int argc, char *argv[]) {
     return 0;
     
 }      
-       
+

diff --git a/ss.h b/ss.h
index 87d4f3f..e3b5b63 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -9,6 +9,8 @@
 #include <cstring>
 #include <iostream>
 #include <stdio.h>
+#include <climits>  // LONG_MAX
+#include <cfloat>   // DBL_MAX
 
 #include "reaction.h"
 #include "pop_ss.h"
@@ -109,7 +111,6 @@ public:
     std::string epitopefile;    // Input file name for targeted epitopes
     
     unsigned int n;             // Population size
-    unsigned int g;             // Number of generations
     unsigned int num_runs;      // Number of times to run the simulation for gathering overall statistics
     double mu;                  // Mutation rate
     double bh;                  // Field "inverse temperature" multiplier
@@ -120,13 +121,9 @@ public:
     bool useEpitope;            // Include selective pressure on an epitope 
     bool useVerbose;            // If true, print extra information while program is running
 
-    unsigned int write_mod;     // Write state of the population every __ generations
     
     std::vector<std::set<unsigned int> > initPop;   // Initial population sequences
     std::vector<double>                  initFrac;  // Initial population fraction
-    std::vector<std::vector<unsigned int> >           eWT;       // Sites that are consensus (WT) in the targeted epitope
-    std::vector<std::vector<unsigned int> >           eMut;      // Sites that are mutant in the targeted epitope
-    double                               penalty;   // Energy penalty if sequence contains the targeted epitope
     
     std::string couplingsInfile;
     std::string trajectoryOutfile;
@@ -134,22 +131,36 @@ public:
     std::string stateInfile;
     std::string epitopeInfile;
 
-    enum PenaltyTYPE { PenaltyEACH, PenaltyTOTAL };
-
-    PenaltyTYPE penaltyType;
+    double t_end;
+    double sample_interval;
+    long max_steps;
 
 
     std::vector<std::string> epitopefiles;
     std::vector<long> init_CTL_num;
     unsigned int num_CTL_clones;
-    
-    
+
+    double rate_s;
+    double rate_u;
+
+    double rate_p;
+    double rate_a;
+    double rate_r;
+    double rate_d;
+    double rate_dprime;
+    double rate_b;
+    double rate_e;
+    double rate_w;
+    double rate_aprime;
+    double rate_g;
+    double rate_h;
+
+
     RunParameters_SS() {
-    
+
         directory="";
         
         n  = 100000;
-        g  = 150;
         num_runs = 1000;
         mu = 6.0e-5;
         bh = 1.0;
@@ -158,10 +169,27 @@ public:
         importState=false;
         useVerbose=false;
 
-        write_mod = 1;
+        t_end = HUGE_VAL;
+        sample_interval = HUGE_VAL;
+        max_steps = 0;
 
         num_CTL_clones = 0;
-        
+
+        rate_s = 1.5;
+        rate_u = 0.5;
+
+        rate_p = 1e-5;
+        rate_a = 1e-7;
+        rate_r = 4.0;
+        rate_d = 0.5;
+        rate_dprime = 3.0;
+        rate_b = 1e-4;
+        rate_e = 3e-4;
+        rate_w = 7.5;
+        rate_aprime = 1e-6;
+        rate_g = 0.03;
+        rate_h = 0.03;
+
     }
 
     void setFiles() {