From: Dariusz Murakowski Date: Thu, 14 May 2015 21:52:56 +0000 (-0400) Subject: CTL killing reactions are smart about virus-changing reactions (by tracking last... X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=bc04fad4ea93af8967758c84de0b538de50475c3;p=VirEvoDyn.git CTL killing reactions are smart about virus-changing reactions (by tracking last-changed species). --- diff --git a/pop_ss.cpp b/pop_ss.cpp index 69ebd18..6ef8419 100644 --- a/pop_ss.cpp +++ b/pop_ss.cpp @@ -11,6 +11,12 @@ #include "seqTools.h" +Virus::Virus() +//: energy(0.0) +//, mutated_sites() +{ +} + // copy existing virus Virus::Virus(const Virus &v) : energy(v.energy) @@ -67,6 +73,16 @@ std::string aaseq2str(const std::vector &aa_seq) return out; } +NTVirus::NTVirus() +//: energy(0.0) +//, nt_seq() +//, aa_seq() +//, config() +//, L_nt(0) +//, L_aa(0) +{ +} + // copy existing virus NTVirus::NTVirus(const NTVirus &v) : energy(v.energy) diff --git a/pop_ss.h b/pop_ss.h index 4312b48..45154b2 100644 --- a/pop_ss.h +++ b/pop_ss.h @@ -23,6 +23,7 @@ public: double energy; std::set mutated_sites; + Virus(); Virus(const Virus &v); Virus(const Hamiltonian &H); Virus(const Hamiltonian &H, const std::set &mutated_sites); @@ -58,6 +59,7 @@ public: unsigned int L_nt; unsigned int L_aa; + NTVirus(); NTVirus(const NTVirus &v); NTVirus(const std::string &NT); NTVirus(const PottsHamiltonian &H, const std::string &NT); @@ -121,6 +123,8 @@ public: VirusSpecies(std::string x) : Species(x) {} VirusSpecies(std::string x, long c) : Species(x,c) {} virus_map pop; + virus_map::iterator lastModified; + Virus lastDeleted; }; class NTVirusSpecies : public Species { @@ -130,6 +134,8 @@ public: NTVirusSpecies(std::string x) : Species(x) {} NTVirusSpecies(std::string x, long c) : Species(x,c) {} NTVirus_map pop; + NTVirus_map::iterator lastModified; + NTVirus lastDeleted; }; class CTLSpecies : public Species { diff --git a/reaction.cpp b/reaction.cpp index 2447661..84bd466 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -138,6 +138,7 @@ void VirusReaction::fire() if (V->pop.count(v)==0) V->pop[v] = 1; // new variant else V->pop[v] += 1; // add to existing pool + V->lastModified = V->pop.find(v); } @@ -179,9 +180,15 @@ void VirusDeathReaction::fire() // update copy number V->pop[iter->first] -= 1; + V->lastModified = iter; + // delete zero-population strains - if (V->pop[iter->first] == 0) - V->pop.erase(iter); + if (V->pop[iter->first] == 0) { + V->lastDeleted = iter->first; + V->lastModified = V->pop.end(); + V->pop.erase(iter); // does not modify iter in C++11 + } + } @@ -202,9 +209,52 @@ double KillingReaction::recalc() return propensity; } -double KillingReaction::recalc(Reaction *rxn) +double KillingReaction::recalc(Reaction *rxn_in) { - return this->recalc(); + RxnType type = get_rxn_type(rxn_in); + + switch (type) { + + case RxnType::Virus: { // increased V->count + //VirusReaction *rxn = dynamic_cast(rxn_in); + virus_map::iterator iter = V->lastModified; + propensity += rate_constant * T->count * T->recognized(iter->first); + break; + } + + case RxnType::VirusDeath: + case RxnType::Killing: { // decreased V->count + virus_map::iterator iter = V->lastModified; + if (iter == V->pop.end()) { + propensity -= rate_constant * T->count * T->recognized(V->lastDeleted); + } + else { + propensity -= rate_constant * T->count * T->recognized(iter->first); + } + + break; + } + + case RxnType::Elementary: + case RxnType::CTLActivation: + case RxnType::NTVirus: + case RxnType::NTVirusDeath: + case RxnType::AAKilling: + case RxnType::CTLaaActivation: { // did not change V->count + propensity = this->recalc(); + // TODO + //std::cerr << "shouldn't have hit here in KillingReaction::recalc(Reaction*)....\n"; + break; + } + + } // switch (type) + + if (V->lastModified == V->pop.end()) { // deleted such that none is left + // TODO + } + + //return this->recalc(); + return propensity; } void KillingReaction::fire() @@ -222,9 +272,13 @@ void KillingReaction::fire() // update copy number V->pop[iter->first] -= 1; + V->lastModified = iter; // delete zero-population strains - if (V->pop[iter->first] == 0) - V->pop.erase(iter); + if (V->pop[iter->first] == 0) { + V->lastDeleted = iter->first; + V->lastModified = V->pop.end(); + V->pop.erase(iter); // does not modify iter in C++11 + } } @@ -353,6 +407,8 @@ void NTVirusReaction::fire() if (V->pop.count(v)==0) V->pop[v] = 1; // new variant else V->pop[v] += 1; // add to existing pool + V->lastModified = V->pop.find(v); + } @@ -394,9 +450,14 @@ void NTVirusDeathReaction::fire() // update copy number V->pop[iter->first] -= 1; + V->lastModified = iter; + // delete zero-population strains - if (V->pop[iter->first] == 0) - V->pop.erase(iter); + if (V->pop[iter->first] == 0) { + V->lastDeleted = iter->first; + V->lastModified = V->pop.end(); + V->pop.erase(iter); // does not modify iter in C++11 + } } @@ -417,9 +478,52 @@ double AAKillingReaction::recalc() return propensity; } -double AAKillingReaction::recalc(Reaction *rxn) +double AAKillingReaction::recalc(Reaction *rxn_in) { - return this->recalc(); + RxnType type = get_rxn_type(rxn_in); + + switch (type) { + + case RxnType::NTVirus: { // increased V->count + //VirusReaction *rxn = dynamic_cast(rxn_in); + NTVirus_map::iterator iter = V->lastModified; + propensity += rate_constant * T->count * T->recognized(iter->first); + break; + } + + case RxnType::NTVirusDeath: + case RxnType::AAKilling: { // decreased V->count + NTVirus_map::iterator iter = V->lastModified; + if (iter == V->pop.end()) { + propensity -= rate_constant * T->count * T->recognized(V->lastDeleted); + } + else { + propensity -= rate_constant * T->count * T->recognized(iter->first); + } + + break; + } + + case RxnType::Elementary: + case RxnType::CTLActivation: + case RxnType::Virus: + case RxnType::VirusDeath: + case RxnType::Killing: + case RxnType::CTLaaActivation: { // did not change V->count + propensity = this->recalc(); + // TODO + //std::cerr << "shouldn't have hit here in KillingReaction::recalc(Reaction*)....\n"; + break; + } + + } // switch (type) + + if (V->lastModified == V->pop.end()) { // deleted such that none is left + // TODO + } + + //return this->recalc(); + return propensity; } void AAKillingReaction::fire() @@ -437,9 +541,13 @@ void AAKillingReaction::fire() // update copy number V->pop[iter->first] -= 1; + V->lastModified = iter; // delete zero-population strains - if (V->pop[iter->first] == 0) - V->pop.erase(iter); + if (V->pop[iter->first] == 0) { + V->lastDeleted = iter->first; + V->lastModified = V->pop.end(); + V->pop.erase(iter); // does not modify iter in C++11 + } }