From 81d6f2da9b71640a40659ba401cd0fc4af4dc887 Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Fri, 17 Apr 2015 20:32:55 -0400 Subject: [PATCH] Output at intervals, remove more WF. --- pop_ss.h | 8 +-- ss.cpp | 212 +++++++++++++++++++++++++++++++++++---------------------------- ss.h | 56 ++++++++++++----- 3 files changed, 164 insertions(+), 112 deletions(-) diff --git a/pop_ss.h b/pop_ss.h index 8c264b5..83cd367 100644 --- 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 epitopes; std::vector affinity; diff --git a/ss.cpp b/ss.cpp index 085a74e..dbeac56 100644 --- 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); // .st + for (size_t i=0; i.st if (r.useEpitope) importEpitope(r,species,reactions,&s1,print_spec); // .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 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 --- a/ss.h +++ b/ss.h @@ -9,6 +9,8 @@ #include #include #include +#include // LONG_MAX +#include // 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 > initPop; // Initial population sequences std::vector initFrac; // Initial population fraction - std::vector > eWT; // Sites that are consensus (WT) in the targeted epitope - std::vector > 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 epitopefiles; std::vector 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() { -- 2.7.4