Merge pull request #660 from StochSS/develop

Release v1.6.7

gillespy2/__version__.py

@@ -24,7 +24,7 @@
 # =============================================================================
 
 
-__version__      = '1.6.6'
+__version__      = '1.6.7'
 
 __title__        = 'GillesPy2'
 __description__  = 'Python interface for Gillespie-style biochemical simulations'

gillespy2/core/assignmentrule.py

@@ -16,6 +16,8 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+import uuid
+
 from gillespy2.core.sortableobject import SortableObject
 from gillespy2.core.jsonify import Jsonify
 
@@ -36,9 +38,12 @@ class AssignmentRule(SortableObject, Jsonify):
     """
 
     def __init__(self, variable=None, formula=None, name=None):
+        if name in (None, ""):
+            self.name = f'ar{uuid.uuid4()}'.replace('-', '_')
+        else:
+            self.name = name
         self.variable = variable
         self.formula = formula
-        self.name = name
 
     def __str__(self):
         return self.variable + ': ' + self.formula

gillespy2/core/events.py

@@ -16,6 +16,8 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+import uuid
+
 from gillespy2.core.gillespyError import *
 from gillespy2.core.jsonify import Jsonify
 
@@ -39,7 +41,12 @@ class EventAssignment(Jsonify):
 
     """
 
-    def __init__(self, variable=None, expression=None):
+    def __init__(self, name=None, variable=None, expression=None):
+
+        if name in (None, ""):
+            self.name = f'evn{uuid.uuid4()}'.replace('-', '_')
+        else:
+            self.name = name
 
         self.variable = variable
         self.expression = expression

gillespy2/core/functiondefinition.py

@@ -16,6 +16,8 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+import uuid
+
 from gillespy2.core.sortableobject import SortableObject
 from gillespy2.core.jsonify import Jsonify
 
@@ -38,7 +40,10 @@ def __init__(self, name="", function=None, args=[]):
         if function is None:
             raise TypeError("Function string provided for FunctionDefinition cannot be None")
 
-        self.name = name
+        if name in (None, ""):
+            self.name = f'fd{uuid.uuid4()}'.replace('-', '_')
+        else:
+            self.name = name
         self.function_string = function
         self.args = args
 

gillespy2/core/gillespySolver.py

@@ -61,19 +61,16 @@ def get_solver_settings(self):
 
         raise SimulationError("This abstract solver class cannot be used directly")
 
-    def get_increment(self, model, increment):
+    def get_increment(self, increment):
         """
         Set the default increment value if it was not provided
 
-        :param model: The model on which the tspan can be found.
-        :type model: gillespy.Model
-
         :param increment: The current value of increment.
         :type increment: int
         """
         if increment is None:
-            return model.tspan[-1] - model.tspan[-2]
-        if model.user_set_tspan:
+            return self.model.tspan[-1] - self.model.tspan[-2]
+        if self.model.user_set_tspan:
             raise  SimulationError(
                 """
                 Failed while preparing to run the model. Both increment and timespan are set.

gillespy2/core/raterule.py

@@ -16,6 +16,8 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
+import uuid
+
 from gillespy2.core.sortableobject import SortableObject
 from gillespy2.core.jsonify import Jsonify
 
@@ -34,9 +36,12 @@ class RateRule(SortableObject, Jsonify):
     """
 
     def __init__(self, variable=None, formula='', name=None):
+        if name in (None, ""):
+            self.name = f'rr{uuid.uuid4()}'.replace('-', '_')
+        else:
+            self.name = name
         self.formula = formula
         self.variable = variable
-        self.name = name
 
     def __str__(self):
         try:

gillespy2/core/reaction.py

@@ -77,7 +77,10 @@ def __init__(self, name="", reactants={}, products={}, propensity_function=None,
         """
 
         # Metadata
-        self.name = name
+        if name in (None, ""):
+            self.name = f'rxn{uuid.uuid4()}'.replace('-', '_')
+        else:
+            self.name = name
         self.annotation = ""
 
         # We might use this flag in the future to automatically generate

gillespy2/solvers/cpp/build/build_engine.py

@@ -110,14 +110,14 @@ 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)
+            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)}")
 
         # 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)
         template_file.unlink()
-        template_gen.write_definitions(str(template_file), model.get_template(variable=variable))
+        template_gen.write_definitions(str(template_file), model.get_template())
         custom_definitions = model.get_options()
         if custom_definitions is not None:
             options_file = self.template_dir.joinpath(self.template_options_name)

gillespy2/solvers/cpp/build/template_gen.py

@@ -34,23 +34,32 @@ class SanitizedModel:
     :type model: gillespy2.Model
     """
     reserved_names = {
-        "vol": "V",
         "t": "t",
     }
 
-    def __init__(self, model: Model):
+    def __init__(self, model: Model, variable=False):
         self.model = model
+        self.variable = variable
 
         self.species: "OrderedDict[str, Species]" = OrderedDict()
         self.species_names = model.sanitized_species_names()
-        for species_name, sanitized_name in self.species_names.items():
+        self.species_id: "OrderedDict[str, int]" = OrderedDict()
+        for spec_id, spec_entry in enumerate(self.species_names.items()):
+            species_name, sanitized_name = spec_entry
             self.species[sanitized_name] = model.get_species(species_name)
+            self.species_id[species_name] = spec_id
 
         self.parameters: "OrderedDict[str, Parameter]" = OrderedDict()
-        self.parameter_names = model.sanitized_parameter_names()
+        self.parameter_names: "OrderedDict[str, str]" = OrderedDict()
+        self.parameter_names["vol"] = "P[0]" if variable else "C[0]"
+        self.parameter_id: "OrderedDict[str, int]" = OrderedDict()
+        for param_id, param_name in enumerate(model.listOfParameters.keys(), start=1):
+            if param_name not in self.parameter_names:
+                self.parameter_names[param_name] = f"P[{param_id}]" if variable else f"C[{param_id}]"
+                self.parameter_id[param_name] = param_id
         for parameter_name, sanitized_name in self.parameter_names.items():
             self.parameters[sanitized_name] = model.get_parameter(parameter_name) \
-                if parameter_name != "vol" else Parameter(name="V", expression=model.volume)
+                if parameter_name != "vol" else Parameter(name=sanitized_name, expression=str(model.volume))
 
         base_namespace = {
             # ORDER IS IMPORTANT HERE!
@@ -153,7 +162,7 @@ def use_rate_rule(self, rate_rule: "RateRule") -> "SanitizedModel":
                 log.warning(f"Could not sanitize rate rule formula expression: {rate_rule.formula}")
         return self
 
-    def get_template(self, variable=False) -> "dict[str, str]":
+    def get_template(self) -> "dict[str, str]":
         """
         Creates a dictionary of C++ macro definitions from the given model.
         The keys of the dictionary contain the name of the macro definition.
@@ -168,7 +177,7 @@ def get_template(self, variable=False) -> "dict[str, str]":
         results = dict({})
 
         # Get definitions for variables
-        parameter_definitions = template_def_variables(self, variable)
+        parameter_definitions = template_def_variables(self, self.variable)
         results.update(parameter_definitions)
 
         # Get definitions for species
@@ -332,11 +341,15 @@ def template_def_variables(model: SanitizedModel, variable=False) -> "dict[str,
     parameter_type = "VARIABLE" if variable else "CONSTANT"
     # Parameter entries, parsed and formatted
     parameter_set = []
-    for param_name, parameter in model.parameters.items():
-        parameter_set.append(f"{parameter_type}({param_name},{parameter.expression})")
+    for param_id, parameter in enumerate(model.parameters.values()):
+        parameter_set.append(f"{parameter_type}({param_id},{parameter.expression})")
 
     return {
-        "GPY_PARAMETER_VALUES": " ".join(parameter_set)
+        "GPY_PARAMETER_VALUES": " ".join(parameter_set),
+        # Currently assumes all variable or all constant.
+        # For partially variable models, modify to compute these two separately.
+        "GPY_PARAMETER_NUM_VARIABLES": str(len(parameter_set)) if variable else "0",
+        "GPY_PARAMETER_NUM_CONSTANTS": str(len(parameter_set)) if not variable else "0",
     }
 
 

gillespy2/solvers/cpp/c_base/model.cpp

@@ -21,6 +21,11 @@
 
 namespace Gillespy {
 
+	int Reaction::s_num_constants;
+	int Reaction::s_num_variables;
+	std::shared_ptr<double> Reaction::s_variables;
+	std::shared_ptr<const double> Reaction::s_constants;
+
 	template <typename PType>
 	Model<PType>::Model(
 		std::vector<std::string> species_names,

gillespy2/solvers/cpp/c_base/model.h

@@ -18,13 +18,17 @@
 
 #pragma once
 
+#include "template.h"
+
 #include <cmath>
 #include <memory>
 #include <string>
 #include <vector>
 #include <iostream>
 
-namespace Gillespy {
+namespace Gillespy
+{
+	typedef unsigned int ReactionId;
 
 	template <typename PType>
 	struct Species {
@@ -48,6 +52,59 @@ namespace Gillespy {
 
 		// List of reactions who's propensities will change when this reaction fires.
 		std::unique_ptr<int[]> species_change;
+
+		inline static double propensity(
+				ReactionId reaction_id,
+				double *state,
+				double *parameters,
+				const double *constants)
+		{
+			return map_propensity(reaction_id, state, parameters, constants);
+		}
+
+		inline static double propensity(
+				ReactionId reaction_id,
+				unsigned int *state,
+				double *parameters,
+				const double *constants)
+		{
+			return map_propensity(reaction_id, state, parameters, constants);
+		}
+
+		inline static double propensity(
+				ReactionId reaction_id,
+				int *state,
+				double *parameters,
+				const double *constants)
+		{
+			return map_propensity(reaction_id, state, parameters, constants);
+		}
+
+		inline static double propensity(ReactionId reaction_id, double *state)
+		{
+			return map_propensity(reaction_id, state, s_variables.get(), s_constants.get());
+		}
+
+		inline static double propensity(ReactionId reaction_id, int *state)
+		{
+			return map_propensity(reaction_id, state, s_variables.get(), s_constants.get());
+		}
+
+		inline static double propensity(ReactionId reaction_id, unsigned int *state)
+		{
+			return map_propensity(reaction_id, state, s_variables.get(), s_constants.get());
+		}
+
+		inline static void load_parameters()
+		{
+			s_variables = std::shared_ptr<double>(get_variables(&s_num_variables));
+			s_constants = std::shared_ptr<const double>(get_constants(&s_num_constants));
+		}
+
+		static int s_num_variables;
+		static int s_num_constants;
+		static std::shared_ptr<double> s_variables;
+		static std::shared_ptr<const double> s_constants;
 	};
 
 	template <typename PType>
@@ -67,15 +124,6 @@ namespace Gillespy {
 		);
 	};
 
-	class IPropensityFunction {
-	public:
-		virtual double evaluate(unsigned int reaction_number, unsigned int *state) = 0;
-		virtual double TauEvaluate(unsigned int reaction_number, const int *S) = 0;
-		virtual double ODEEvaluate(int reaction_number, const std::vector<double> &S) = 0;
-
-		virtual ~IPropensityFunction() {};
-	};
-
 	template <typename PType>
 	struct Simulation {
 		int random_seed;
@@ -95,8 +143,6 @@ namespace Gillespy {
 
 		Model<PType> *model;
 
-		IPropensityFunction *propensity_function;
-
 		template <class T> friend std::ostream &operator << (std::ostream &os, const Simulation<T> &simulation);
 
 		// output_results_buffer: Writes the contents of the entire simulation trajectory.

gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp

@@ -39,25 +39,6 @@ unsigned int number_timesteps = 0;
 double end_time = 100.0;
 double increment = 0;
 
-class PropensityFunction : public IPropensityFunction
-{
-public:
-	double evaluate(unsigned int reaction_number, unsigned int *S)
-	{
-		return 1.0;
-	}
-
-	double TauEvaluate(unsigned int reaction_number, const int *S)
-	{
-		return 1.0;
-	}
-
-	double ODEEvaluate(int reaction_number, const std::vector<double> &S)
-	{
-		return map_ode_propensity(reaction_number, S);
-	}
-};
-
 int main(int argc, char *argv[]) {
 	ArgParser parser = ArgParser(argc, argv);
 
@@ -72,6 +53,7 @@ int main(int argc, char *argv[]) {
 	number_timesteps = parser.timesteps;
 	increment = parser.increment;
 
+	Reaction::load_parameters();
 	Model<double> model(species_names, species_populations, reaction_names);
 	add_reactions(model);
 
@@ -80,15 +62,13 @@ int main(int argc, char *argv[]) {
 		random_seed = time(NULL);
 	}
 
-	IPropensityFunction *propensity_function = new PropensityFunction();
 	Simulation<double> simulation;
 
 	simulation.model = &model;
 	simulation.end_time = end_time;
 	simulation.random_seed = random_seed;
 	simulation.number_timesteps = number_timesteps;
 	simulation.number_trajectories = number_trajectories;
-	simulation.propensity_function = propensity_function;
 	simulation.current_time = 0.0;
 	simulation.output_interval = parser.output_interval;
 
@@ -98,7 +78,5 @@ int main(int argc, char *argv[]) {
 	ODESolver(&simulation, increment);
 	simulation.output_buffer_final(std::cout);
 
-	delete propensity_function;
-
 	return 0;
 }