Fix some dependents. Modify total propensity instead of totally re-sum.

diff --git a/ss.cpp b/ss.cpp
index d92bd58..9705e55 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -125,10 +125,10 @@ bool does_affect(Reaction *newrxn, Reaction *rxn)
         std::cout << "newrxn = Killing: " << newrxn->name << '\n';
         Species *s = dynamic_cast<KillingReaction*>(newrxn)->V;
         affects = does_species_appear_in_reactants(s,rxn);
-        //return affects;
-        Species *s2 = dynamic_cast<KillingReaction*>(newrxn)->T;
-        affects |= does_species_appear_in_reactants(s2,rxn);
         return affects;
+        //Species *s2 = dynamic_cast<KillingReaction*>(newrxn)->T;
+        //affects |= does_species_appear_in_reactants(s2,rxn);
+        //return affects;
         break;
     }
 
@@ -175,11 +175,13 @@ bool does_affect(Reaction *newrxn, Reaction *rxn)
 
     case RxnType::CTLaaActivation: {
         std::cout << "newrxn = CTLaaActivation: " << newrxn->name << '\n';
-        Species *s = dynamic_cast<CTLaaActivationReaction*>(newrxn)->V;
-        affects = does_species_appear_in_reactants(s,rxn);
+        //Species *s = dynamic_cast<CTLaaActivationReaction*>(newrxn)->V;
+        //affects = does_species_appear_in_reactants(s,rxn);
         //return affects;
         Species *s2 = dynamic_cast<CTLaaActivationReaction*>(newrxn)->Tfrom;
         affects |= does_species_appear_in_reactants(s2,rxn);
+        s2 = dynamic_cast<CTLaaActivationReaction*>(newrxn)->Tto;
+        affects |= does_species_appear_in_reactants(s2,rxn);
         return affects;
         break;
     }
@@ -367,7 +369,7 @@ void run(RunParameters_SS &r, unsigned seed) {
 
 void print_species_header(const Species_parray &print_spec)
 {
-    std::cout << "time\tstep";
+    std::cout << "step\ttime";
     for (Species_parray::const_iterator spec = print_spec.begin(),
             end = print_spec.end();
             spec != end; ++spec)
@@ -448,8 +450,6 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         total_propensity += (*rxn)->propensity;
     }
 
-    std::cout << "total propensity = " << total_propensity << '\n';
-
 
     // stop if time runs out or too many steps or no more reactions will occur
     for (; (max_steps == 0 || step < max_steps) && (t < t_end) && (total_propensity > 0.0); ++step) {
@@ -475,29 +475,27 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         if (rxn == end)
             break;
 
+        // update reaction dependencies
+        double oldA = (*rxn)->propensity;
         (*rxn)->recalc();
-        std::cout << "fired: " << (*rxn)->name << '\n';
+        double newA = (*rxn)->propensity;
+        total_propensity = total_propensity - oldA + newA;
         for (rxn2 = (*rxn)->affects.begin(), end2 = (*rxn)->affects.end(); rxn2 != end2; ++rxn2) {
-            std::cout << "affects: " << (*rxn2)->name << '\n';
-            double oldA = (*rxn2)->propensity;
-            //std::cout << step << "\thit it\n";
-            //std::cout << step << "\told = " << (*rxn2)->propensity;
-            printf("%ld\told = %.4f",step,oldA);
+            oldA = (*rxn2)->propensity;
             (*rxn2)->recalc();
-            //std::cout << "\tnew = " << (*rxn2)->propensity << '\n';
-            double newA = (*rxn2)->propensity;
-            printf("\tnew = %.4f",newA);
-            if (oldA == newA)
-                printf("\tSAME");
-            printf("\n");
+            newA = (*rxn2)->propensity;
+            if (oldA == newA && newA != 0.0) {
+                std::cout << "SAME: " << (*rxn)->name << " affected " << (*rxn2)->name << " with propensity " << newA << '\n';
+            }
+            total_propensity = total_propensity - oldA + newA;
         }
 
         // recalculate rates
-        total_propensity = 0.0;
-        for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
-            //(*rxn)->recalc();     // dependencies handled above
-            total_propensity += (*rxn)->propensity;
-        }
+        //total_propensity = 0.0;
+        //for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
+        //    //(*rxn)->recalc();     // dependencies handled above
+        //    total_propensity += (*rxn)->propensity;
+        //}
 
         // print copy numbers
         if (t_next_sample <= t) {
@@ -508,7 +506,6 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         if ((step_sample_interval != 0) && (step_next_sample <= step)) {
             step_next_sample += step_sample_interval;
             print_species_traj(t,step,print_spec,verbose);
-            std::cout << "total propensity = " << total_propensity << '\n';
         }
 
     }