Fixes #825 for C solvers, changes argument to constructor

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp

@@ -35,6 +35,9 @@
 #include "integrator.h"
 #include "tau.h"
 
+#include "template_defaults.h"
+
+
 static void silent_error_handler(int error_code, const char *module, const char *function_name,
                           char *message, void *eh_data);
 
@@ -313,15 +316,19 @@ namespace Gillespy
                     }
 
                     // Expected tau step is determined.
-                    tau_step = select<double, double>(
-                            model,
-                            tau_args,
-                            tau_tol,
-                            simulation->current_time,
-                            save_time,
-                            sol.data.propensities,
-                            current_state
-                    );
+                    if(GPY_CONSTANT_TAU_STEPSIZE > 0){
+                        tau_step = GPY_CONSTANT_TAU_STEPSIZE;
+                    }else{
+                        tau_step = select<double, double>(
+                                model,
+                                tau_args,
+                                tau_tol,
+                                simulation->current_time,
+                                save_time,
+                                sol.data.propensities,
+                                current_state
+                        );
+                    }
                     partition_species(
                             simulation->current_time,
                             simulation->reaction_state,

gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSolver.cpp

@@ -31,6 +31,8 @@
 #include "TauLeapingSolver.h"
 #include "tau.h"
 
+#include "template_defaults.h"
+
 namespace Gillespy
 {
     static volatile bool interrupted = false;
@@ -149,8 +151,11 @@ namespace Gillespy
                     {
                         propensity_values[reaction_number] = Reaction::propensity(reaction_number, current_state.data());
                     }
-
-                    tau_step = select(*(simulation->model), tau_args, tau_tol, simulation->current_time, save_time, propensity_values, current_state);
+                    if(GPY_CONSTANT_TAU_STEPSIZE > 0){
+                        tau_step = GPY_CONSTANT_TAU_STEPSIZE;
+                    }else{
+                        tau_step = select(*(simulation->model), tau_args, tau_tol, simulation->current_time, save_time, propensity_values, current_state);
+                    }
                     prev_curr_state = current_state;
                     double prev_curr_time = simulation->current_time;
                     int loop_cnt = 0;

gillespy2/solvers/cpp/c_base/template/template.cpp

@@ -52,6 +52,7 @@ namespace Gillespy
     int reaction_reactants[GPY_NUM_REACTIONS][GPY_NUM_SPECIES] = GPY_REACTION_REACTANTS;
     int reaction_products[GPY_NUM_REACTIONS][GPY_NUM_SPECIES] = GPY_REACTION_PRODUCTS;
     int reaction_props_deps[GPY_NUM_REACTIONS][GPY_NUM_SPECIES] = GPY_REACTION_PROPS_DEPS;
+    //double constant_tau_stepsize = GPY_CONSTANT_TAU_STEPSIZE;
     std::string r_names[GPY_NUM_REACTIONS] = 
     {
         #define REACTION_NAME(name) #name,

gillespy2/solvers/cpp/c_base/template/template_defaults.h

@@ -63,6 +63,7 @@
 #define GPY_REACTION_NAMES
 #endif
 
+
 // ===============================================================
 // ================ HYBRID SOLVER OPTION DEFAULTS ================
 // ===============================================================
@@ -87,5 +88,9 @@
 #define GPY_RATE_RULES
 #endif
 
+#ifndef GPY_CONSTANT_TAU_STEPSIZE
+#define GPY_CONSTANT_TAU_STEPSIZE
+#endif
+
 #endif
 

gillespy2/solvers/cpp/tau_hybrid_c_solver.py

@@ -22,6 +22,12 @@ class TauHybridCSolver(GillesPySolver, CSolver):
     name = "TauHybridCSolver"
     target = "hybrid"
 
+    def __init__(self, model = None, output_directory = None, delete_directory = True, resume=None, variable = False,  constant_tau_stepsize=None):
+
+        self.constant_tau_stepsize = constant_tau_stepsize
+        super().__init__(model=model, output_directory=output_directory, 
+                         delete_directory=delete_directory, resume=resume, variable=variable)
+
     class ErrorStatus(IntEnum):
         UNKNOWN = 1
         LOOP_OVER_INTEGRATE = 2
@@ -30,8 +36,7 @@ class ErrorStatus(IntEnum):
         NEGATIVE_STATE_NO_SSA_REACTION = 5
         NEGATIVE_STATE_AT_BEGINING_OF_STEP = 6
 
-    @classmethod
-    def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel":
+    def __create_options(self, sanitized_model: "SanitizedModel") -> "SanitizedModel":
         """
         Populate the given list of species modes into a set of template macro definitions.
         Generated options are specific to the Tau Hybrid solver,
@@ -169,14 +174,25 @@ def get_supported_features(cls):
     def get_supported_integrator_options(cls) -> "Set[str]":
         return { "rtol", "atol", "max_step" }
 
+
+
     def _build(self, model: "Union[Model, SanitizedModel]", simulation_name: str, variable: bool, debug: bool = False,
                custom_definitions=None) -> str:
         variable = variable or len(model.listOfEvents) > 0
-        sanitized_model = TauHybridCSolver.__create_options(SanitizedModel(model, variable=variable))
+        sanitized_model = self.__create_options(SanitizedModel(model, variable=variable))
         for rate_rule in model.listOfRateRules.values():
             sanitized_model.use_rate_rule(rate_rule)
+
+        # determine if a constant stepsize has been requested
+        if self.constant_tau_stepsize is not None:
+            sanitized_model.options['GPY_CONSTANT_TAU_STEPSIZE'] = str(float(self.constant_tau_stepsize))
+        else:
+            sanitized_model.options['GPY_CONSTANT_TAU_STEPSIZE'] = '0'
         return super()._build(sanitized_model, simulation_name, variable, debug)
 
+
+
+
     def _handle_return_code(self, return_code: "int") -> "int":
         if return_code == TauHybridCSolver.ErrorStatus.UNKNOWN:
             raise SimulationError("C++ solver failed (no error code given).")
@@ -270,10 +286,12 @@ def run(self=None, model: Model = None, t: int = None, number_of_trajectories: i
                 raise SimulationError("A model is required to run the simulation.")
             self._set_model(model=model)
 
+
         self.model.compile_prep()
         self.validate_model(self.model, model)
         self.validate_sbml_features(model=self.model)
 
+
         self.validate_tspan(increment=increment, t=t)
         if increment is None:
             increment = self.model.tspan[-1] - self.model.tspan[-2]

gillespy2/solvers/cpp/tau_leaping_c_solver.py

@@ -23,6 +23,9 @@
 
 from .c_solver import CSolver, SimulationReturnCode
 
+from gillespy2.solvers.cpp.build.template_gen import SanitizedModel
+
+
 class TauLeapingCSolver(GillesPySolver, CSolver):
 
     """
@@ -35,6 +38,23 @@ class TauLeapingCSolver(GillesPySolver, CSolver):
     name = "TauLeapingCSolver"
     target = "tau_leap"
 
+    def __init__(self, model = None, output_directory = None, delete_directory = True, resume=None, variable = False,  constant_tau_stepsize=None):
+
+        self.constant_tau_stepsize = constant_tau_stepsize
+        super().__init__(model=model, output_directory=output_directory, 
+                         delete_directory=delete_directory, resume=resume, variable=variable) 
+
+    def _build(self, model, simulation_name, variable, debug = False,
+               custom_definitions=None):
+        sanitized_model = SanitizedModel(model, variable=variable)
+        # determine if a constant stepsize has been requested
+        if self.constant_tau_stepsize is not None:
+            sanitized_model.options['GPY_CONSTANT_TAU_STEPSIZE'] = str(float(self.constant_tau_stepsize))
+        else:
+            sanitized_model.options['GPY_CONSTANT_TAU_STEPSIZE'] = '0'
+        return super()._build(sanitized_model, simulation_name, variable, debug)
+
+
     @classmethod
     def get_solver_settings(cls):
         """
@@ -46,7 +66,7 @@ def get_solver_settings(cls):
 
     def run(self=None, model: Model = None, t: int = None, number_of_trajectories: int = 1, timeout: int = 0,
             increment: int = None, seed: int = None, debug: bool = False, profile: bool = False, variables={},
-            resume=None, live_output: str = None, live_output_options: dict = {}, tau_tol=0.03, **kwargs):
+            resume=None, live_output: str = None, live_output_options: dict = {}, tau_tol=0.03, constant_tau_stepsize=None, **kwargs):
 
         """
         :param model: The model on which the solver will operate. (Deprecated)
@@ -85,6 +105,8 @@ def run(self=None, model: Model = None, t: int = None, number_of_trajectories: i
         result in larger tau steps. Default value is 0.03.
         :type tau_tol: float
 
+        :param constant_tau_stepsize: If set, overrides the automatic stepsize selection and uses the given
+                    value as the stepsize on each step.
         :returns: A result object containing the results of the simulation
         :rtype: gillespy2.Results
         """
@@ -111,6 +133,8 @@ def run(self=None, model: Model = None, t: int = None, number_of_trajectories: i
                 raise SimulationError("A model is required to run the simulation.")
             self._set_model(model=model)
 
+        self.constant_tau_stepsize = constant_tau_stepsize
+
         self.model.compile_prep()
         self.validate_model(self.model, model)
         self.validate_sbml_features(model=self.model)

gillespy2/solvers/numpy/tau_hybrid_solver.py

@@ -73,7 +73,7 @@ class TauHybridSolver(GillesPySolver):
     result = None
     stop_event = None
 
-    def __init__(self, model=None, profile_reactions=False):
+    def __init__(self, model=None, profile_reactions=False, constant_tau_stepsize=None):
         if model is None:
             raise SimulationError("A model is required to run the simulation.")
 
@@ -93,6 +93,7 @@ def __init__(self, model=None, profile_reactions=False):
                 self.non_negative_species.add(key.name)
             for key, value in model.listOfReactions[reaction].products.items():
                 self.non_negative_species.add(key.name)
+        self.constant_tau_stepsize = constant_tau_stepsize
 
     def __save_state_to_output(self, curr_time, save_index, curr_state, species, 
                                 trajectory, save_times):
@@ -931,7 +932,7 @@ def get_supported_features(cls):
 
     def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None, seed=None,
             debug=False, profile=False, tau_tol=0.03, event_sensitivity=100,integrator_options={},
-            live_output=None, live_output_options={}, timeout=None, constant_tau_stepsize=None, **kwargs):
+            live_output=None, live_output_options={}, timeout=None, **kwargs):
         """
         Function calling simulation of the model. This is typically called by the run function in GillesPy2 model
         objects and will inherit those parameters which are passed with the model as the arguments this run function.
@@ -981,9 +982,6 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None,
 
         :param timeout: If set, if simulation takes longer than timeout, will exit.
         :type timeout: int
-
-        :param constant_tau_stepsize: If set, overrides the automatic stepsize selection and uses the given
-            value as the stepsize on each step.
        
         :returns: A result object containing the results of the simulation.
         :rtype: gillespy2.Results
@@ -1016,7 +1014,6 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None,
         self.validate_model(self.model, model)
         self.validate_sbml_features(model=self.model)
 
-        self.constant_tau_stepsize = constant_tau_stepsize
 
         self.validate_tspan(increment=increment, t=t)
         if increment is None:

gillespy2/solvers/numpy/tau_leaping_solver.py

@@ -43,7 +43,7 @@ class TauLeapingSolver(GillesPySolver):
     pause_event = None
     result = None
 
-    def __init__(self, model=None, debug=False):
+    def __init__(self, model=None, debug=False, constant_tau_stepsize=None):
         if model is None:
             raise SimulationError("A model is required to run the simulation.")
 
@@ -55,6 +55,7 @@ def __init__(self, model=None, debug=False):
         self.model = copy.deepcopy(model)
         self.debug = debug
         self.is_instantiated = True
+        self.constant_tau_stepsize = True
 
     def __get_reactions(self, step, curr_state, curr_time, save_time, propensities, reactions):
         """
@@ -97,7 +98,7 @@ def get_solver_settings(cls):
 
     def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None, seed=None,
             debug=False, profile=False,  live_output=None, live_output_options={},
-            timeout=None, resume=None, tau_tol=0.03, constant_tau_stepsize=None, **kwargs):
+            timeout=None, resume=None, tau_tol=0.03, **kwargs):
         """
         Function calling simulation of the model.
         This is typically called by the run function in GillesPy2 model objects
@@ -143,9 +144,6 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None,
             result in larger tau steps. Default value is 0.03.
         :type tau_tol: float
 
-        :param constant_tau_stepsize: If set, overrides the automatic stepsize selection and uses the given
-            value as the stepsize on each step.
-        
         :returns: A result object containing the results of the simulation.
         :rtype: gillespy2.Results
         """
@@ -165,7 +163,6 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None,
             )
             self = TauLeapingSolver(model=model, debug=debug, profile=profile)
 
-        self.constant_tau_stepsize = constant_tau_stepsize
 
         if model is not None:
             log.warning('model = gillespy2.model is deprecated. Future releases '