Recalculate differently based on what reaction caused propensity "change" (currently...

diff --git a/reaction.cpp b/reaction.cpp
index f77aecc..2447661 100644 (file)
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -31,6 +31,11 @@ double ElementaryReaction::recalc()
     return this->propensity;
 }
 
+double ElementaryReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void ElementaryReaction::fire()
 {
     std::vector<Species*>::iterator reac = this->reactants.begin(),
@@ -84,6 +89,11 @@ double VirusReaction::recalc()
     return this->propensity;
 }
 
+double VirusReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void VirusReaction::fire()
 {
     // pick which virus reproduces
@@ -148,6 +158,11 @@ double VirusDeathReaction::recalc()
     return propensity;
 }
 
+double VirusDeathReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void VirusDeathReaction::fire()
 {
     // pick which virus dies
@@ -187,6 +202,11 @@ double KillingReaction::recalc()
     return propensity;
 }
 
+double KillingReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void KillingReaction::fire()
 {
     // pick which virus to kill
@@ -224,6 +244,11 @@ double CTLActivationReaction::recalc()
 
 }
 
+double CTLActivationReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void CTLActivationReaction::fire()
 {
     Tfrom->count -= 1;
@@ -254,6 +279,11 @@ double NTVirusReaction::recalc()
     return this->propensity;
 }
 
+double NTVirusReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void NTVirusReaction::fire()
 {
     // pick which virus reproduces
@@ -343,6 +373,11 @@ double NTVirusDeathReaction::recalc()
     return propensity;
 }
 
+double NTVirusDeathReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void NTVirusDeathReaction::fire()
 {
     // pick which virus dies
@@ -382,6 +417,11 @@ double AAKillingReaction::recalc()
     return propensity;
 }
 
+double AAKillingReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void AAKillingReaction::fire()
 {
     // pick which virus to kill
@@ -419,6 +459,11 @@ double CTLaaActivationReaction::recalc()
 
 }
 
+double CTLaaActivationReaction::recalc(Reaction *rxn)
+{
+    return this->recalc();
+}
+
 void CTLaaActivationReaction::fire()
 {
     Tfrom->count -= 1;

diff --git a/reaction.h b/reaction.h
index d5eb1d8..f8a29b7 100644 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -19,6 +19,7 @@ public:
     std::string name;
     std::vector<Reaction*> affects; // reactions whose propensity must be updated when this one fires
     virtual double recalc() = 0;
+    virtual double recalc(Reaction *rxn) = 0;   // what other reaction fired to cause this one's propensity to be recalculated?
     virtual void fire() = 0;
 };
 
@@ -35,6 +36,7 @@ public:
     std::vector<int> reactant_stoich;   // reaction order
     std::vector<int> product_stoich;  // must be negative for reactants that are consumed!
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -51,6 +53,7 @@ public:
     Hamiltonian H;
 
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -65,6 +68,7 @@ public:
     VirusSpecies* V;
     Hamiltonian H;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -78,6 +82,7 @@ public:
     VirusSpecies* V;
     CTLSpecies* T;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -91,6 +96,7 @@ public:
     CTLSpecies* Tfrom;
     CTLSpecies* Tto;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -111,6 +117,7 @@ public:
     PottsHamiltonian H;
 
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -125,6 +132,7 @@ public:
     NTVirusSpecies* V;
     PottsHamiltonian H;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -138,6 +146,7 @@ public:
     NTVirusSpecies* V;
     CTLaaSpecies* T;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 
@@ -151,6 +160,7 @@ public:
     CTLaaSpecies* Tfrom;
     CTLaaSpecies* Tto;
     double recalc();
+    double recalc(Reaction *rxn);
     void fire();
 };
 

diff --git a/ss.cpp b/ss.cpp
index 3356b2e..4a8e438 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -456,7 +456,7 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         total_propensity = total_propensity - oldA + newA;
         for (rxn2 = (*rxn)->affects.begin(), end2 = (*rxn)->affects.end(); rxn2 != end2; ++rxn2) {
             oldA = (*rxn2)->propensity;
-            (*rxn2)->recalc();
+            (*rxn2)->recalc(*rxn);
             newA = (*rxn2)->propensity;
             //if (oldA == newA && newA != 0.0) {
             //    std::cout << "SAME: " << (*rxn)->name << " affected " << (*rxn2)->name << " with propensity " << newA << '\n';