From 062fd84e9223851b6e7755c025d73a1e7800b33a Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Thu, 9 Apr 2015 09:11:09 -0400 Subject: [PATCH] Separate function for running the dynamics. --- reaction.cpp | 24 +++++++++++ reaction.h | 3 +- ss.cpp | 133 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- ss.h | 19 +++++++++ 4 files changed, 174 insertions(+), 5 deletions(-) diff --git a/reaction.cpp b/reaction.cpp index 3eca7ee..20f55e3 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -202,3 +202,27 @@ void KillingReaction::fire() V->pop.erase(iter); } + +double CTLActivationReaction::recalc() +{ + propensity = rate_constant * Tfrom->count; + + double partial_sum = 0.0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + partial_sum += iter->second * Tfrom->recognized(iter->first); + } + + propensity *= partial_sum; + return propensity; + +} + +void CTLActivationReaction::fire() +{ + Tfrom->count -= 1; + Tto->count += 1; +} + + diff --git a/reaction.h b/reaction.h index 7501941..9865256 100644 --- a/reaction.h +++ b/reaction.h @@ -59,7 +59,8 @@ public: void fire(); }; -class CTLTransitionReaction : public Reaction { +// V_s + T_k -> V_s + T_k' +class CTLActivationReaction : public Reaction { public: VirusSpecies* V; CTLSpecies* Tfrom; diff --git a/ss.cpp b/ss.cpp index d7b09b8..7f08189 100644 --- a/ss.cpp +++ b/ss.cpp @@ -5,6 +5,7 @@ #include #include #include // for partial_sum +#include #include @@ -14,8 +15,6 @@ #include "ss.h" #include "reaction.h" -typedef std::vector Rxn_parray; - gsl_rng *rnd; //pointer to random number generator state // generate a random number generator seed from /dev/urandom @@ -52,7 +51,10 @@ void run(RunParameters_SS &r, unsigned seed) { //gsl_rng_set(rnd,rand()); gsl_rng_set(rnd,seed); -#if 1 + r.setFiles(); + if (r.importState) importState(r); // .st + if (r.useEpitope) importEpitope(r); // .ep + std::string couplingsFile = "neutral_2site.j"; Hamiltonian H(couplingsFile); //double mu = 6.0e-5; @@ -77,6 +79,16 @@ void run(RunParameters_SS &r, unsigned seed) { r2->V = &s1; reactions.push_back(r2); + Species_parray print_spec; + print_spec.push_back(&s1); + + double t_end = 1.0; + long max_steps = LONG_MAX; + simulate(reactions, species, t_end, max_steps, print_spec); + + +#if 0 + /* double total_propensity = 0.0; for (std::vector::iterator rxn = reactions.begin(), @@ -485,6 +497,111 @@ void run(RunParameters_SS &r, unsigned seed) { } +void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const Species_parray &print_spec) +{ + long step = 0; + + double t = 0.0; + double t_remain = t_end; + double dt = 0.0; + + double total_propensity = 0.0; + + Rxn_parray::iterator rxn, end; + + for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { + (*rxn)->recalc(); + total_propensity += (*rxn)->propensity; + } + + + while (1) { + + double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity; // uniform on (0,R] + double a_sum = 0.0; + + // 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) { + a_sum += (*rxn)->propensity; + if (rand <= a_sum) { + (*rxn)->fire(); + break; + } + } + + // TODO: better handle case with total propensity == 0.0 + + if (rxn == end) + break; + + // recalculate rates + total_propensity = 0.0; + for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { + (*rxn)->recalc(); + total_propensity += (*rxn)->propensity; + } + + if (step >= max_steps) + break; + + ++step; + + 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'; + +#if 0 + // print copy numbers + std::cout << t << '\t'; + for (std::vector::iterator spec = species.begin(), + end = species.end(); + spec != end; ++spec) { + std::cout << (*spec)->count << '\t'; + } + + MutatedSiteSequence ms; + FILE* output = stdout; + for (virus_map::iterator iter = s1.pop.begin(), + end = s1.pop.end(); + iter != end; ++iter) { + std::cout << iter->second << '('; + ms = iter->first.mutated_sites; + // print sequence as CSV (comma-separated) + MutatedSiteSequence::iterator ms_iter=ms.begin(); + if (ms_iter != ms.end()) { + fprintf(output,"%d",*ms_iter); + ++ms_iter; + } + for (; ms_iter!=ms.end(); ++ms_iter) + fprintf(output,",%d",*ms_iter); + std::cout << ')'; + std::cout << '\t'; + } + std::cout << '\n'; +#endif + + } + + + +} + + // Import initial state from a state file void importState(RunParameters_SS &r) { @@ -545,6 +662,7 @@ void importEpitope(RunParameters_SS &r) { std::string word; unsigned int site; + double aff; std::vector tmpEp; size_t pos0 = 0; size_t posBar = readStr.find('|',pos0); @@ -553,6 +671,13 @@ void importEpitope(RunParameters_SS &r) { // could use std::noskipws to keep tab ('\t') characters? std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0)); + + // first entry is affinity + readStrStrm >> word; + aff = strtodouble(word); + printf("%.6e ",aff); + + // next entries are WT sites while (readStrStrm >> word) { std::cout << word << " "; // XXX std::istringstream i(word); @@ -669,7 +794,7 @@ int main(int argc, char *argv[]) { else if (strcmp(argv[i],"-i")==0) { if (++i==argc) break; else r.infile=argv[i]; } else if (strcmp(argv[i],"-o")==0) { if (++i==argc) break; else r.outfile=argv[i]; } else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } } - else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefile=argv[i]; r.useEpitope=true; } } + else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(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]); } diff --git a/ss.h b/ss.h index ad52952..747094a 100644 --- a/ss.h +++ b/ss.h @@ -10,6 +10,12 @@ #include #include +#include "reaction.h" +#include "pop_ss.h" + +typedef std::vector Rxn_parray; +typedef std::vector Species_parray; + // STRING MANIPULATION @@ -118,6 +124,10 @@ public: enum PenaltyTYPE { PenaltyEACH, PenaltyTOTAL }; PenaltyTYPE penaltyType; + + + + std::vector epitopefiles; RunParameters_SS() { @@ -163,6 +173,10 @@ public: supplementaryOutfile=directory+"/"+outfile+".sum"; stateInfile=directory+"/"+statefile+".st"; epitopeInfile=directory+"/"+epitopefile+".ep"; + + for (size_t i = 0; i != epitopefiles.size(); ++i) { + epitopefiles[i] = directory+"/"+epitopefiles[i]+".ep"; + } } @@ -173,6 +187,10 @@ public: supplementaryOutfile=outfile+".sum"; stateInfile=statefile+".st"; epitopeInfile=epitopefile+".ep"; + + for (size_t i = 0; i != epitopefiles.size(); ++i) { + epitopefiles[i] = epitopefiles[i]+".ep"; + } } @@ -187,6 +205,7 @@ public: void run(RunParameters_SS &r); void importState(RunParameters_SS &r); void importEpitope(RunParameters_SS &r); +void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const Species_parray &print_spec); #endif -- 2.7.4