CTL killing reactions are smart about virus-changing reactions (by tracking last...

diff --git a/pop_ss.cpp b/pop_ss.cpp
index 69ebd18..6ef8419 100644 (file)
--- 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> &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 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -23,6 +23,7 @@ public:
     double energy;
     std::set<unsigned int> mutated_sites;
     
+    Virus();
     Virus(const Virus &v);
     Virus(const Hamiltonian &H);
     Virus(const Hamiltonian &H, const std::set<unsigned int> &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 (file)
--- 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<Reaction*>(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<Reaction*>(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
+    }
 }