Merge pull request #853 from StochSS/fix_hybrid_switching

Fix hybrid switching

gillespy2/solvers/cpp/build/build_engine.py

@@ -115,8 +115,8 @@ def prepare(self, model: "Union[Model, template_gen.SanitizedModel]", variable=F
         # If a raw GillesPy2 model was provided, convert it to a sanitized model.
         if isinstance(model, gillespy2.Model):
             model = template_gen.SanitizedModel(model, variable=variable)
-        elif not isinstance(model, template_gen.SanitizedModel):
-            raise TypeError(f"Build engine expected gillespy2.Model or SanitizedModel type: received {type(model)}")
+        elif not isinstance(model, template_gen.SanitizedModel) and type(model).__name__ == "SanitizedModel":
+            raise TypeError(f"Build engine expected gillespy2.Model or SanitizedModel type: received {type(model)} , __name__={type(model).__name__}")
 
         # Build the template and write it to the temp directory and remove the sample template_definitions header.
         template_file = self.template_dir.joinpath(self.template_definitions_name)

gillespy2/solvers/cpp/build/template_gen.py

@@ -87,6 +87,8 @@ def __init__(self, model: Model, variable=False):
         # Reactions: maps reaction names to their stoichiometry matrix.
         # Stoichiometry matrix maps a sanitized species name to its stoichiometry.
         self.reactions: "OrderedDict[str, dict[str, int]]" = OrderedDict()
+        self.reaction_reactants: "OrderedDict[str, dict[str, int]]" = OrderedDict()
+        self.reaction_products: "OrderedDict[str, dict[str, int]]" = OrderedDict()
         # Rate Rules: maps sanitized species names to their corresponding rate rule expression.
         self.rate_rules: "OrderedDict[str, str]" = OrderedDict()
         # Options: custom definitions that can be supplied by the solver, maps macros to their definitions.
@@ -122,15 +124,17 @@ def use_reaction(self, reaction: "Reaction") -> "SanitizedModel":
         :type reaction: gillespy2.Reaction
         """
         self.reactions[reaction.name] = {spec: int(0) for spec in self.species_names.values()}
+        self.reaction_reactants[reaction.name] = {spec: int(0) for spec in self.species_names.values()}
+        self.reaction_products[reaction.name] = {spec: int(0) for spec in self.species_names.values()}
         for reactant, stoich_value in reaction.reactants.items():
-            if isinstance(reactant, Species):
-                reactant = self.species_names[reactant.name]
+            reactant = self.species_names[reactant.name]
             self.reactions[reaction.name][reactant] -= int(stoich_value)
+            self.reaction_reactants[reaction.name][reactant] = int(stoich_value)
 
         for product, stoich_value in reaction.products.items():
-            if isinstance(product, Species):
-                product = self.species_names[product.name]
+            product = self.species_names[product.name]
             self.reactions[reaction.name][product] += int(stoich_value)
+            self.reaction_products[reaction.name][product] = int(stoich_value)
 
         return self
 
@@ -401,18 +405,30 @@ def template_def_reactions(model: SanitizedModel, ode=False) -> "dict[str, str]"
     """
     num_reactions = str(len(model.reactions))
     reaction_set = OrderedDict()
+    reactants_set = OrderedDict()
+    products_set = OrderedDict()
 
     for rxn_name, reaction in model.reactions.items():
         stoich = [str(int(reaction[species])) for species in model.species_names.values()]
         reaction_set[rxn_name] = f"{{{','.join(stoich)}}}"
+    for rxn_name, reaction_reactants in model.reaction_reactants.items():
+        reactants_count = [str(int(reaction_reactants[species])) for species in model.species_names.values()]
+        reactants_set[rxn_name] = f"{{{','.join(reactants_count)}}}"
+    for rxn_name, reaction_products in model.reaction_products.items():
+        products_count = [str(int(reaction_products[species])) for species in model.species_names.values()]
+        products_set[rxn_name] = f"{{{','.join(products_count)}}}"
 
     reaction_names = " ".join([f"REACTION_NAME({rxn})" for rxn in reaction_set.keys()])
     reaction_set = f"{{{','.join(reaction_set.values())}}}"
+    reactants_set = f"{{{','.join(reactants_set.values())}}}"
+    products_set = f"{{{','.join(products_set.values())}}}"
 
     return {
         "GPY_NUM_REACTIONS": num_reactions,
-        "GPY_REACTIONS": reaction_set,
         "GPY_REACTION_NAMES": reaction_names,
+        "GPY_REACTIONS": reaction_set,
+        "GPY_REACTION_REACTANTS": reactants_set,
+        "GPY_REACTION_PRODUCTS": products_set,
     }
 
 

gillespy2/solvers/cpp/c_base/Tau/tau.cpp

@@ -42,12 +42,12 @@ namespace Gillespy
 
             for (int spec = 0; spec < model.number_species; spec++)
             {
-                if (model.reactions[r].species_change[spec] > 0)
+                if (model.reactions[r].products_change[spec] > 0)
                 {
                     tau_args.products[r].push_back(spec);
                 }
 
-                else if (model.reactions[r].species_change[spec] < 0)
+                else if (model.reactions[r].reactants_change[spec] > 0)
                 {
                     rxn_order += 1;
                     tau_args.reactions_reactants[r].push_back(spec);

gillespy2/solvers/cpp/c_base/model.cpp

@@ -46,9 +46,13 @@ namespace Gillespy {
         for (unsigned int reaction = 0; reaction < number_reactions; reaction++) {
             reactions[reaction].name = reaction_names[reaction];
             reactions[reaction].species_change = std::make_unique<int[]>(number_species);
+            reactions[reaction].reactants_change = std::make_unique<int[]>(number_species);
+            reactions[reaction].products_change = std::make_unique<int[]>(number_species);
 
             for (unsigned int species = 0; species < number_species; species++) {
                 reactions[reaction].species_change[species] = 0;
+                reactions[reaction].reactants_change[species] = 0;
+                reactions[reaction].products_change[species] = 0;
             }
 
             reactions[reaction].affected_reactions = std::vector<unsigned int>();
@@ -66,7 +70,8 @@ namespace Gillespy {
         for (unsigned int r1 = 0; r1 < number_reactions; r1++) {
             for (unsigned int r2 = 0; r2 < number_reactions; r2++) {
                 for (unsigned int s = 0; s < number_species; s++) {
-                    if (reactions[r2].species_change[s] != 0) {
+                     if(reactions[r1].species_change[s]  != 0  &&
+                        reactions[r2].reactants_change[s] > 0 ){
                         reactions[r1].affected_reactions.push_back(r2);
                     }
                 }

gillespy2/solvers/cpp/c_base/model.h

@@ -72,6 +72,8 @@ namespace Gillespy
 
         // List of reactions who's propensities will change when this reaction fires.
         std::unique_ptr<int[]> species_change;
+        std::unique_ptr<int[]> reactants_change;
+        std::unique_ptr<int[]> products_change;
 
         inline static double propensity(
                 ReactionId reaction_id,

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp

@@ -304,37 +304,35 @@ namespace Gillespy
             int num_species = species.size();
             std::set<int> det_rxns;
 
-            for (int rxn_i = 0; rxn_i < reactions.size(); ++rxn_i)
+            for (int rxn_i = 0; rxn_i < num_reactions; ++rxn_i)
             {
                 // start with the assumption that reaction is determinstic
                 HybridReaction &rxn = reactions[rxn_i];
                 rxn.mode = SimulationState::CONTINUOUS;
 
                 // iterate through the dependent species of this reaction
-                // Loop breaks if we've already determined that it is to be marked as discrete.
                 for (int spec_i = 0; spec_i < num_species && rxn.mode == SimulationState::CONTINUOUS; ++spec_i)
                 {
-                    // Reaction has a dependency on a species if its dx is positive or negative.
-                    // Any species with "dependency" change of 0 is by definition not a dependency.
-                    if (rxn.get_base_reaction()->species_change[spec_i] == 0)
-                    {
+                    if (rxn.get_base_reaction()->reactants_change[spec_i] == 0 && 
+                        rxn.get_base_reaction()->products_change[spec_i] == 0) {
                         continue;
                     }
 
                     // if any of the dependencies are set by the user as discrete OR
                     // have been set as dynamic and has not been flagged as deterministic,
                     // allow it to be modelled discretely
-                    if (species[spec_i].user_mode == SimulationState::DYNAMIC)
-                    {
+                    if (species[spec_i].user_mode == SimulationState::DYNAMIC) {
                         rxn.mode = species[spec_i].partition_mode;
-                    } else
-                    {
+                    } else {
                         rxn.mode = species[spec_i].user_mode;
                     }
+                    // Loop breaks if we've already determined that it is to be marked as discrete.
+                    if(rxn.mode == SimulationState::DISCRETE){
+                        break;
+                    }
                 }
 
-                if (rxn.mode == SimulationState::CONTINUOUS)
-                {
+                if (rxn.mode == SimulationState::CONTINUOUS) {
                     det_rxns.insert(rxn_i);
                 }
             }
@@ -343,90 +341,120 @@ namespace Gillespy
         }
 
         void partition_species(
+                double current_time,
                 std::vector<HybridReaction> &reactions,
                 std::vector<HybridSpecies> &species,
                 const std::vector<double> &propensity_values,
                 std::vector<double> &curr_state,
                 double tau_step,
                 const TauArgs<double> &tauArgs)
         {
+
+
             // coefficient of variance- key:species id, value: cv
-            std::map<int, double> cv;
+            std::map<int, double> CV;
             // means
             std::map<int, double> means;
             // standard deviation
             std::map<int, double> sd;
 
             // Initialize means and sd's
-            for (int spec_i = 0; spec_i < species.size(); ++spec_i)
-            {
+            for (int spec_i = 0; spec_i < species.size(); ++spec_i) {
                 HybridSpecies &spec = species[spec_i];
 
-                if (spec.user_mode == SimulationState::DYNAMIC)
-                {
+                if (spec.user_mode == SimulationState::DYNAMIC) {
                     means.insert({spec_i, curr_state[spec_i]});
                     sd.insert({spec_i, 0});
                 }
             }
 
             // calculate means and standard deviations for dynamic-mode species involved in reactions
-            for (int rxn_i = 0; rxn_i < reactions.size(); ++rxn_i)
-            {
+            for (int rxn_i = 0; rxn_i < reactions.size(); ++rxn_i) {
                 HybridReaction &rxn = reactions[rxn_i];
-
-                for (int spec_i = 0; spec_i < species.size(); ++spec_i)
-                {
-                    // Only dynamic species whose mean/SD is requested are to be considered.
-                    if (means.count(spec_i) <= 0)
-                    {
+                for (int spec_i = 0; spec_i < species.size(); ++spec_i) {
+                    HybridSpecies &spec = species[spec_i];
+                    // Only dynamic species are to be considered.
+                    if (spec.user_mode != SimulationState::DYNAMIC) {
                         continue;
                     }
-                    // Selected species is either a reactant or a product, depending on whether
-                    //   dx is positive or negative.
-                    // 0-dx species are not dependencies of this reaction, so dx == 0 is ignored.
-                    int spec_dx = rxn.get_base_reaction()->species_change[spec_i];
-                    if (spec_dx < 0)
-                    {
+                    // Selected species is either a reactant or a product, 
+                    if(rxn.get_base_reaction()->reactants_change[spec_i] > 0) {
                         // Selected species is a reactant.
-                        means[spec_i] -= (tau_step * propensity_values[rxn_i] * spec_dx);
-                        sd[spec_i] += (tau_step * propensity_values[rxn_i] * std::pow(spec_dx, 2));
-                    } else if (spec_dx > 0)
-                    {
+                        means[spec_i] -= (propensity_values[rxn_i] * rxn.get_base_reaction()->reactants_change[spec_i]);
+                        sd[spec_i] += (propensity_values[rxn_i] * std::pow(rxn.get_base_reaction()->reactants_change[spec_i], 2));
+                    } 
+                    if(rxn.get_base_reaction()->products_change[spec_i] > 0) {
                         // Selected species is a product.
                         HybridSpecies &product = species[spec_i];
-                        means[spec_i] += (tau_step * propensity_values[rxn_i] * spec_dx);
-                        sd[spec_i] += (tau_step * propensity_values[rxn_i] * std::pow(spec_dx, 2));
+                        means[spec_i] += (propensity_values[rxn_i] * rxn.get_base_reaction()->products_change[spec_i]);
+                        sd[spec_i] += (propensity_values[rxn_i] * std::pow(rxn.get_base_reaction()->products_change[spec_i], 2));
                     }
                 }
             }
-
             // calculate coefficient of variation using means and sd
+            for (int spec_i = 0; spec_i < species.size(); ++spec_i) {
+                HybridSpecies &spec = species[spec_i];
+                // Only dynamic species are to be considered.
+                if (spec.user_mode != SimulationState::DYNAMIC) {
+                    continue;
+                }
+                if (means[spec_i] > 0 && sd[spec_i] > 0) {
+                    CV[spec_i] = (sqrt(sd[spec_i]) / means[spec_i]);
+                } else {
+                    CV[spec_i] = 1;
+                }
+            }
+
+            // keep a history of the past CV values, and calculate a time-average values
+            std::map<int, double> CV_a;
+            static std::map<int, std::queue<double>> cv_history;
+            static std::map<int, double> cv_history_sum;
+            int history_length = 12;
+            // 
+            if( current_time == 0.0 ){  // reset cv_history at start of trajectory
+                cv_history.clear();
+                cv_history_sum.clear();
+            }
+
             for (int spec_i = 0; spec_i < species.size(); ++spec_i)
             {
-                if (means.count(spec_i) <= 0)
-                {
+
+                HybridSpecies &spec = species[spec_i];
+                // Only dynamic species are to be considered.
+                if (spec.user_mode != SimulationState::DYNAMIC){
                     continue;
                 }
-
+                if(cv_history.count(spec_i) == 0){
+                    cv_history[spec_i] = std::queue<double>();
+                }
+                if(cv_history_sum.count(spec_i) == 0){
+                    cv_history_sum[spec_i] = 0.0;
+                }
+                cv_history[spec_i].push(CV[spec_i]);
+                cv_history_sum[spec_i] += CV[spec_i];
+                if(cv_history[spec_i].size() > history_length){
+                    double removed = cv_history[spec_i].front();
+                    cv_history[spec_i].pop();
+                    cv_history_sum[spec_i] -= removed;
+                }
+                CV_a[spec_i] = cv_history_sum[spec_i]/cv_history[spec_i].size();
+            }
 
+            // Select DISCRETE or CONTINOUS mode for each species
+            for (int spec_i = 0; spec_i < species.size(); ++spec_i) {
                 HybridSpecies &spec = species[spec_i];
-                if (spec.switch_min == 0)
-                {
-                    // (default value means switch  min not set, use switch tol)
-                    if (means[spec_i] > 0)
-                    {
-                        cv[spec_i] = (sd[spec_i] / means[spec_i]);
-                    } else
-                    {
-                        cv[spec_i] = 1;
-                    }
-                    //std::cerr<<"\t\tspec"<<spec_i<<" cv="<<cv[spec_i]<<"=("<<sd[spec_i]<<" / "<<means[spec_i]<<") switch_tol="<<spec.switch_tol<<"\n";
 
-                    spec.partition_mode = cv[spec_i] < spec.switch_tol
+                // Only dynamic species are to be considered.
+                if (spec.user_mode != SimulationState::DYNAMIC) {
+                    continue;
+                }
+                unsigned int prev_partition = spec.partition_mode;
+                if (spec.switch_min == 0) {
+                    // (default value means switch  min not set, use switch tol)
+                    spec.partition_mode = CV_a[spec_i] < spec.switch_tol
                                           ? SimulationState::CONTINUOUS
                                           : SimulationState::DISCRETE;
                 } else {
-                    //std::cerr<<"\t\tspec"<<spec_i<<" mean="<<means[spec_i]<<" switch_min="<<spec.switch_min<<"\n";
                     spec.partition_mode = means[spec_i] > spec.switch_min
                                           ? SimulationState::CONTINUOUS
                                           : SimulationState::DISCRETE;

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h

@@ -308,6 +308,7 @@ namespace Gillespy
                 std::vector<HybridSpecies> &species);
 
         void partition_species(
+                double current_time,
                 std::vector<HybridReaction> &reactions,
                 std::vector<HybridSpecies> &species,
                 const std::vector<double> &propensity_values,