From f12f64122cf8ffc598dd456f2c15befbc7c8dd3c Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 10 Jan 2022 12:22:34 -0800 Subject: [PATCH 01/54] Adding test to check performance across versions of GillesPy2 --- .../run_version_test.py | 47 +++++++ .../tyson_oscillator.py | 129 ++++++++++++++++++ 2 files changed, 176 insertions(+) create mode 100755 test/version_performance_tests/run_version_test.py create mode 100755 test/version_performance_tests/tyson_oscillator.py diff --git a/test/version_performance_tests/run_version_test.py b/test/version_performance_tests/run_version_test.py new file mode 100755 index 000000000..5022bdd0a --- /dev/null +++ b/test/version_performance_tests/run_version_test.py @@ -0,0 +1,47 @@ +#!/usr/bin/env python3 + +import subprocess +import os +import re + +versions = [ +] + +cmd3 = f"git clone https://github.com/StochSS/GillesPy2" +print(cmd3) +output3 = subprocess.run(cmd3, shell=True, capture_output=True) + +cmd4 = f"cd GillesPy2; git tag" +print(cmd4) +output4 = subprocess.run(cmd4, shell=True, capture_output=True) +versions = output4.stdout.decode("utf-8").split("\n") + + +def SortFn(ver): + x = re.findall("\d+",ver) + return int(x[0])*10000+int(x[1])*100+int(x[2]) + +while("" in versions) : + versions.remove("") + +versions.sort(key=SortFn) + +#print(f"versions={versions}") + + +for ver in versions: + print(f"{ver}: ",end='') + cmd2 = f"cd GillesPy2; git checkout {ver}" + output2 = subprocess.run(cmd2, shell=True, capture_output=True) + cmd = f'export PYTHONPATH="{os.getcwd()}/GillesPy2"; python3 ./tyson_oscillator.py' + #print(cmd) + + output = subprocess.run(cmd, shell=True, capture_output=True) + out = output.stdout.decode("utf-8").strip() + try: + print(float(out)) + except: + print("error"); + #print(output) + + diff --git a/test/version_performance_tests/tyson_oscillator.py b/test/version_performance_tests/tyson_oscillator.py new file mode 100755 index 000000000..507ece707 --- /dev/null +++ b/test/version_performance_tests/tyson_oscillator.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +# ============================================================================= +# File: tyson_example.py +# Summary: two-state oscillator from a paper by Novak & Tyson, 2008 +# +# ----------------------- +# How to run this example +# ----------------------- +# +# You can run this program from Python by starting a terminal shell on your +# computer, then changing the working directory in your shell to the +# directory where this file is located, and then running the following command: +# +# python3 -m tyson_oscillator.py +# +# ------------------ +# Author and history +# ------------------ +# 2019-06-13 Sean Matthew +# ============================================================================= + +import numpy as np +import sys + + +try: + import gillespy2 +except: + print(f'Could not import the gillespy2 package. path={sys.path}') + sys.exit() + + +# Model definition. +# ............................................................................. +# In GillesPy2, a model to be simulated is expressed as an object having the +# parent class "Model". Components of the model to be simulated, such as the +# reactions, molecular species, and the time span for simualtion, are all +# defined within the class definition. + +class Tyson2StateOscillator(gillespy2.Model): + """ + Here, as a test case, we run a simple two-state oscillator (Novak & Tyson + 2008) as an example of a stochastic reaction system. + """ + def __init__(self, parameter_values=None): + """ + """ + system_volume = 300 #system volume + gillespy2.Model.__init__(self, name="tyson-2-state", volume=system_volume) + self.timespan(np.linspace(0,100,101)) + + # ============================================= + # Define model species, initial values, parameters, and volume + # ============================================= + + # Parameter values for this biochemical system are given in + # concentration units. However, stochastic systems must use population + # values. For example, a concentration unit of 0.5mol/(L*s) + # is multiplied by a volume unit, to get a population/s rate + # constant. Thus, for our non-mass action reactions, we include the + # parameter "vol" in order to convert population units to concentration + # units. Volume here = 300. + + P = gillespy2.Parameter(name='P', expression=2.0) + kt = gillespy2.Parameter(name='kt', expression=20.0) + kd = gillespy2.Parameter(name='kd', expression=1.0) + a0 = gillespy2.Parameter(name='a0', expression=0.005) + a1 = gillespy2.Parameter(name='a1', expression=0.05) + a2 = gillespy2.Parameter(name='a2', expression=0.1) + kdx = gillespy2.Parameter(name='kdx', expression=1.0) + self.add_parameter([P, kt, kd, a0, a1, a2, kdx]) + + # Species + # Initial values of each species (concentration converted to pop.) + X = gillespy2.Species(name='X', initial_value=int(0.65609071*system_volume)) + Y = gillespy2.Species(name='Y', initial_value=int(0.85088331*system_volume)) + self.add_species([X, Y]) + + # ============================================= + # Define the reactions within the model + # ============================================= + + # creation of X: + rxn1 = gillespy2.Reaction(name = 'X production', + reactants = {}, + products = {X:1}, + propensity_function = 'vol*1/(1+(Y*Y/((vol*vol))))') + + # degradadation of X: + rxn2 = gillespy2.Reaction(name = 'X degradation', + reactants = {X:1}, + products = {}, + rate = kdx) + + # creation of Y: + rxn3 = gillespy2.Reaction(name = 'Y production', + reactants = {X:1}, + products = {X:1, Y:1}, + rate = kt) + + # degradation of Y: + rxn4 = gillespy2.Reaction(name = 'Y degradation', + reactants = {Y:1}, + products = {}, + rate = kd) + + # nonlinear Y term: + rxn5 = gillespy2.Reaction(name = 'Y nonlin', + reactants = {Y:1}, + products = {}, + propensity_function = 'Y/(a0 + a1*(Y/vol)+a2*Y*Y/(vol*vol))') + + self.add_reaction([rxn1,rxn2,rxn3,rxn4,rxn5]) + + +# Model simulation. +# ............................................................................. + +if __name__ == '__main__': + # A simulation in GillesPy2 is performed by first instantiating the model + # to be simulated, and then invoking the "run" method on that object. + # The results of the simulation are the output of "run". + import time + from gillespy2 import SSACSolver + tic = time.time() + for _ in range(10): + tyson_model = Tyson2StateOscillator() + trajectories = tyson_model.run(solver=SSACSolver) + print((time.time() - tic)/10) From d89a4ac4df30fa465fd3e1e977ee263feae17818 Mon Sep 17 00:00:00 2001 From: BryanRumsey <44621966+BryanRumsey@users.noreply.github.com> Date: Tue, 11 Jan 2022 15:44:34 -0500 Subject: [PATCH 02/54] Changed version to v1.6.8 --- gillespy2/__version__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/__version__.py b/gillespy2/__version__.py index 38a796ba0..8e42af11f 100644 --- a/gillespy2/__version__.py +++ b/gillespy2/__version__.py @@ -24,7 +24,7 @@ # ============================================================================= -__version__ = '1.6.7' +__version__ = '1.6.8' __title__ = 'GillesPy2' __description__ = 'Python interface for Gillespie-style biochemical simulations' From e9554328718af0b452dfa1d7ff11bb7e018d49fd Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 13 Jan 2022 14:26:05 -0500 Subject: [PATCH 03/54] Fixed 'SimulationError' import issue. --- gillespy2/solvers/cpp/ode_c_solver.py | 3 ++- gillespy2/solvers/cpp/ssa_c_solver.py | 3 ++- gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 3 ++- gillespy2/solvers/cpp/tau_leaping_c_solver.py | 3 ++- gillespy2/solvers/numpy/ode_solver.py | 2 +- gillespy2/solvers/numpy/ssa_solver.py | 3 ++- gillespy2/solvers/numpy/tau_leaping_solver.py | 2 +- 7 files changed, 12 insertions(+), 7 deletions(-) diff --git a/gillespy2/solvers/cpp/ode_c_solver.py b/gillespy2/solvers/cpp/ode_c_solver.py index 43a59a776..4bde7118b 100644 --- a/gillespy2/solvers/cpp/ode_c_solver.py +++ b/gillespy2/solvers/cpp/ode_c_solver.py @@ -19,7 +19,8 @@ from gillespy2.solvers.cpp.c_decoder import BasicSimDecoder from gillespy2.solvers.utilities import solverutils as cutils -from gillespy2.core import GillesPySolver, gillespyError, Model +from gillespy2.core import GillesPySolver, Model +from gillespy2.core.gillespyError import * from gillespy2.core import Results from .c_solver import CSolver, SimulationReturnCode diff --git a/gillespy2/solvers/cpp/ssa_c_solver.py b/gillespy2/solvers/cpp/ssa_c_solver.py index 34867e40e..7aa81d2e7 100644 --- a/gillespy2/solvers/cpp/ssa_c_solver.py +++ b/gillespy2/solvers/cpp/ssa_c_solver.py @@ -20,7 +20,8 @@ from gillespy2.solvers.cpp.c_decoder import IterativeSimDecoder from gillespy2.solvers.utilities import solverutils as cutils -from gillespy2.core import GillesPySolver, gillespyError, Model +from gillespy2.core import GillesPySolver, Model +from gillespy2.core.gillespyError import * from gillespy2.core import Results from .c_solver import CSolver, SimulationReturnCode diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index f3d41a02f..a18a9b850 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -1,7 +1,8 @@ import gillespy2 from gillespy2.solvers.cpp.c_decoder import IterativeSimDecoder from gillespy2.solvers.utilities import solverutils as cutils -from gillespy2.core import GillesPySolver, gillespyError, Model +from gillespy2.core import GillesPySolver, Model +from gillespy2.core.gillespyError import * from typing import Union from gillespy2.core import Results diff --git a/gillespy2/solvers/cpp/tau_leaping_c_solver.py b/gillespy2/solvers/cpp/tau_leaping_c_solver.py index c56a98ea9..36d1d34aa 100644 --- a/gillespy2/solvers/cpp/tau_leaping_c_solver.py +++ b/gillespy2/solvers/cpp/tau_leaping_c_solver.py @@ -19,7 +19,8 @@ from gillespy2.solvers.cpp.c_decoder import IterativeSimDecoder from gillespy2.solvers.utilities import solverutils as cutils -from gillespy2.core import GillesPySolver, gillespyError, Model +from gillespy2.core import GillesPySolver, Model +from gillespy2.core.gillespyError import * from gillespy2.core import Results from .c_solver import CSolver, SimulationReturnCode diff --git a/gillespy2/solvers/numpy/ode_solver.py b/gillespy2/solvers/numpy/ode_solver.py index acd5efc58..957346903 100644 --- a/gillespy2/solvers/numpy/ode_solver.py +++ b/gillespy2/solvers/numpy/ode_solver.py @@ -22,7 +22,7 @@ from scipy.integrate import ode from collections import OrderedDict import numpy as np -from gillespy2.core import GillesPySolver, log, gillespyError +from gillespy2.core import GillesPySolver, log, gillespyError, SimulationError from gillespy2.solvers.utilities import solverutils as nputils from gillespy2.core.results import Results diff --git a/gillespy2/solvers/numpy/ssa_solver.py b/gillespy2/solvers/numpy/ssa_solver.py index 5056b6fba..46c1984cf 100644 --- a/gillespy2/solvers/numpy/ssa_solver.py +++ b/gillespy2/solvers/numpy/ssa_solver.py @@ -18,7 +18,8 @@ from threading import Thread, Event from gillespy2.core.results import Results -from gillespy2.core import GillesPySolver, log, gillespyError +from gillespy2.core import GillesPySolver, log +from gillespy2.core.gillespyError import * from gillespy2.solvers.utilities import solverutils as nputils import random import math diff --git a/gillespy2/solvers/numpy/tau_leaping_solver.py b/gillespy2/solvers/numpy/tau_leaping_solver.py index 09e510af7..a54766a47 100644 --- a/gillespy2/solvers/numpy/tau_leaping_solver.py +++ b/gillespy2/solvers/numpy/tau_leaping_solver.py @@ -24,7 +24,7 @@ from gillespy2.solvers.utilities import Tau from gillespy2.solvers.utilities import solverutils as nputils from gillespy2.core import GillesPySolver, log, liveGraphing -from gillespy2.core import ModelError, ExecutionError +from gillespy2.core.gillespyError import * from gillespy2.core.results import Results class TauLeapingSolver(GillesPySolver): From 20c4b6d59442c285bd372491d8b8e665c2c867ff Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 13 Jan 2022 11:33:51 -0800 Subject: [PATCH 04/54] Fixes #679 --- gillespy2/core/model.py | 1 + gillespy2/solvers/cpp/build/template_gen.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 965931401..00d05135c 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -26,6 +26,7 @@ from gillespy2.core.results import Trajectory,Results from collections import OrderedDict from gillespy2.core.gillespyError import * +from .gillespyError import SimulationError try: import lxml.etree as eTree diff --git a/gillespy2/solvers/cpp/build/template_gen.py b/gillespy2/solvers/cpp/build/template_gen.py index 164977ee9..4e73cec1d 100644 --- a/gillespy2/solvers/cpp/build/template_gen.py +++ b/gillespy2/solvers/cpp/build/template_gen.py @@ -21,6 +21,8 @@ from gillespy2.core import Species, Reaction, Parameter, Model, RateRule from gillespy2.solvers.cpp.build.expression import Expression from gillespy2.core import log +from gillespy2.core.gillespyError import SimulationError + import math From f4952e37014cc1194c55113a22f078245348c9fa Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Thu, 13 Jan 2022 17:39:41 -0500 Subject: [PATCH 05/54] Bug fixes for models with >2 events - Fixed return value for `get_event_id()` - Removed duplicate trigger state update - Added validation to `Integrator#integrate` wrapper method --- .../cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp | 6 ------ .../solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h | 2 +- .../solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp | 4 ++++ 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp index 606598013..5295c9858 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp @@ -604,12 +604,6 @@ namespace Gillespy m_trigger_pool.erase(event.get_event_id()); } - // Step 5: Update any trigger states to reflect the new trigger value. - for (auto &event : m_events) - { - m_trigger_state.find(event.get_event_id())->second = event.trigger(t, event_state); - } - return has_active_events(); } } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h index b21fd4f45..cc8223108 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h @@ -88,7 +88,7 @@ namespace Gillespy } inline bool is_persistent() const { return m_use_persist; } - inline bool get_event_id() const { return m_event_id; } + inline int get_event_id() const { return m_event_id; } EventExecution get_execution(double t, const double *state, int num_state) const; diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index 615ec0b1d..ed641e6a3 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -143,6 +143,10 @@ IntegrationResults Integrator::integrate(double *t) IntegrationResults Integrator::integrate(double *t, std::set &event_roots, std::set &reaction_roots) { IntegrationResults results = integrate(t); + if (status != IntegrationStatus::OK) { + return results; + } + unsigned long long num_triggers = data.active_triggers.size(); unsigned long long num_rxn_roots = data.active_reaction_ids.size(); unsigned long long root_size = data.active_triggers.size() + data.active_reaction_ids.size(); From 1ef7594f3d0e4bbdfc90d4dc6904fc72abd5648e Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 14 Jan 2022 15:11:27 -0500 Subject: [PATCH 06/54] Bug fixes for pause-resume - Changed `--init_pop` to accept `float` rather than `int` to remove rounding on startup - Use `getpid()` as a temporary fix for RNG seed --- gillespy2/solvers/cpp/build/template_gen.py | 2 +- gillespy2/solvers/cpp/c_base/model.cpp | 2 +- .../tau_hybrid_cpp_solver/TauHybridSimulation.cpp | 2 +- .../cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp | 4 ---- .../cpp/c_base/tau_hybrid_cpp_solver/integrator.h | 2 +- gillespy2/solvers/utilities/solverutils.py | 12 ++++++------ 6 files changed, 10 insertions(+), 14 deletions(-) diff --git a/gillespy2/solvers/cpp/build/template_gen.py b/gillespy2/solvers/cpp/build/template_gen.py index 164977ee9..d89d9e3e9 100644 --- a/gillespy2/solvers/cpp/build/template_gen.py +++ b/gillespy2/solvers/cpp/build/template_gen.py @@ -368,7 +368,7 @@ def template_def_species(model: SanitizedModel) -> "dict[str, str]": populations = OrderedDict() for spec_name, spec in model.species.items(): - populations[spec_name] = str(spec.initial_value) + populations[spec_name] = str(float(spec.initial_value)) # Species names, parsed and formatted sanitized_names = [f"SPECIES_NAME({name})" for name in populations.keys()] populations = f"{{{','.join(populations.values())}}}" diff --git a/gillespy2/solvers/cpp/c_base/model.cpp b/gillespy2/solvers/cpp/c_base/model.cpp index fd6a1418d..41132f262 100644 --- a/gillespy2/solvers/cpp/c_base/model.cpp +++ b/gillespy2/solvers/cpp/c_base/model.cpp @@ -147,7 +147,7 @@ namespace Gillespy { os << timeline[timestep] << ','; for (int species = 0; species < model->number_species; species++) { - os << trajectories[trajectory][timestep][species] << ','; + os << (double) trajectories[trajectory][timestep][species] << ','; } } } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp index f3ba18c28..37de29e8b 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp @@ -53,7 +53,7 @@ int main(int argc, char* argv[]) add_reactions(model); if(seed_time){ - random_seed = time(NULL); + random_seed = time(nullptr) % getpid(); } //Simulation INIT TauHybrid::HybridSimulation simulation(model); diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index 615ec0b1d..17d75d439 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -217,10 +217,6 @@ bool Integrator::disable_root_finder() return validate(this, CVodeRootInit(cvode_mem, 0, NULL)); } - -URNGenerator::URNGenerator() - : uniform(0, 1) {} - URNGenerator::URNGenerator(unsigned long long seed) : uniform(0, 1), rng(seed) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h index a0ce3884c..f0a3a0f88 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h @@ -170,7 +170,7 @@ namespace Gillespy unsigned long long seed; public: double next(); - URNGenerator(); + URNGenerator() = delete; explicit URNGenerator(unsigned long long seed); }; diff --git a/gillespy2/solvers/utilities/solverutils.py b/gillespy2/solvers/utilities/solverutils.py index 3dd958d69..f5ea54fbe 100644 --- a/gillespy2/solvers/utilities/solverutils.py +++ b/gillespy2/solvers/utilities/solverutils.py @@ -91,19 +91,19 @@ def update_species_init_values(listOfSpecies, species, variables, resume = None) populations = '' for i in range(len(species) - 1): if species[i] in variables: - populations += '{} '.format(int(variables[species[i]])) + populations += '{} '.format(float(variables[species[i]])) else: if resume is not None: - populations += '{} '.format(int(resume[species[i]][-1])) + populations += '{} '.format(float(resume[species[i]][-1])) else: - populations += '{} '.format(int(listOfSpecies[species[i]].initial_value)) + populations += '{} '.format(float(listOfSpecies[species[i]].initial_value)) if species[-1] in variables: - populations += '{}'.format(int(variables[species[-1]])) + populations += '{}'.format(float(variables[species[-1]])) else: if resume is not None: - populations += '{} '.format(int(resume[species[-1]][-1])) + populations += '{} '.format(float(resume[species[-1]][-1])) else: - populations += '{}'.format(int(listOfSpecies[species[-1]].initial_value)) + populations += '{}'.format(float(listOfSpecies[species[-1]].initial_value)) return populations def change_param_values(listOfParameters, parameters, volume, variables): From 0fe86fcbf5a48291343a35503be00d12399fefcd Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 14 Jan 2022 16:54:48 -0500 Subject: [PATCH 07/54] Add C++ includes/macros for `getpid()`, port `getpid()` seed to other solvers --- gillespy2/solvers/cpp/c_base/model.h | 8 ++++++++ .../solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp | 2 +- .../solvers/cpp/c_base/ssa_cpp_solver/SSASimulation.cpp | 2 +- .../c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp | 2 +- .../tau_leaping_cpp_solver/TauLeapingSimulation.cpp | 2 +- 5 files changed, 12 insertions(+), 4 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/model.h b/gillespy2/solvers/cpp/c_base/model.h index f28e4196c..3837e96ca 100644 --- a/gillespy2/solvers/cpp/c_base/model.h +++ b/gillespy2/solvers/cpp/c_base/model.h @@ -26,6 +26,14 @@ #include #include +#if defined(WIN32) || defined(WIN32) || defined(__WIN32) || defined(__WIN32__) +#include +#define GPY_PID_GET() ((int) GetCurrentProcessId()) +#else +#include +#define GPY_PID_GET() (getpid()) +#endif + namespace Gillespy { typedef unsigned int ReactionId; diff --git a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp index 499510c32..b756c4ff2 100644 --- a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESimulation.cpp @@ -59,7 +59,7 @@ int main(int argc, char *argv[]) { if (seed_time) { - random_seed = time(NULL); + random_seed = time(nullptr) % GPY_PID_GET(); } Simulation simulation; diff --git a/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASimulation.cpp b/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASimulation.cpp index 1150714e6..29eb1df23 100644 --- a/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASimulation.cpp @@ -59,7 +59,7 @@ int main(int argc, char *argv[]) if (seed_time) { - random_seed = time(NULL); + random_seed = time(nullptr) % GPY_PID_GET(); } Simulation simulation; diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp index 37de29e8b..d13878b09 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp @@ -53,7 +53,7 @@ int main(int argc, char* argv[]) add_reactions(model); if(seed_time){ - random_seed = time(nullptr) % getpid(); + random_seed = time(nullptr) % GPY_PID_GET(); } //Simulation INIT TauHybrid::HybridSimulation simulation(model); diff --git a/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSimulation.cpp b/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSimulation.cpp index 12bbeaf42..f77ac6188 100644 --- a/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSimulation.cpp @@ -59,7 +59,7 @@ int main(int argc, char* argv[]){ add_reactions(model); if(seed_time) { - random_seed = time(NULL); + random_seed = time(nullptr) % GPY_PID_GET(); } Simulation simulation; From f2cd38bbe9ee1d9f9eb4866ce604d5e1b65bb1d3 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Tue, 18 Jan 2022 21:14:41 -0500 Subject: [PATCH 08/54] Remove `population` references in Tau Hybrid RHS function --- .../tau_hybrid_cpp_solver/HybridModel.cpp | 31 +++---------------- .../tau_hybrid_cpp_solver/HybridModel.h | 4 +-- .../tau_hybrid_cpp_solver/integrator.cpp | 16 ++-------- 3 files changed, 10 insertions(+), 41 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp index 5295c9858..81ff83799 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp @@ -242,8 +242,7 @@ namespace Gillespy double DifferentialEquation::evaluate( const double t, - double *ode_state, - int *ssa_state) + double *ode_state) { double sum = 0.0; @@ -254,7 +253,7 @@ namespace Gillespy for (auto &formula : formulas) { - sum += formula(ode_state, ssa_state); + sum += formula(ode_state); } return sum; @@ -293,29 +292,9 @@ namespace Gillespy auto &formula_set = spec.diff_equation.formulas; int spec_diff = rxn.base_reaction->species_change[spec_i]; - switch (spec.partition_mode) - { - case SimulationState::CONTINUOUS: - formula_set.push_back([rxn_i, spec_diff]( - double *ode_state, - int *ssa_state) - { - return spec_diff * Reaction::propensity(rxn_i, ode_state); - }); - break; - - case SimulationState::DISCRETE: - formula_set.push_back([rxn_i, spec_diff]( - double *ode_state, - int *ssa_state) - { - return spec_diff * Reaction::propensity(rxn_i, ssa_state); - }); - break; - - default: - break; - } + formula_set.push_back([rxn_i, spec_diff](double *state) { + return spec_diff * Reaction::propensity(rxn_i, state); + }); } } } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h index cc8223108..da680d44e 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h @@ -177,10 +177,10 @@ namespace Gillespy struct DifferentialEquation { public: - std::vector> formulas; + std::vector> formulas; std::vector> rate_rules; - double evaluate(double t, double *ode_state, int *ssa_state); + double evaluate(double t, double *ode_state); }; enum SimulationState : unsigned int diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index acc2d4cf0..a69974429 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -293,9 +293,6 @@ int Gillespy::TauHybrid::rhs(realtype t, N_Vector y, N_Vector ydot, void *user_d std::vector *species = data->species_state; std::vector *reactions = data->reaction_state; std::vector &propensities = data->propensities; - // Concentrations and reactions are both used for their respective propensity evaulations. - // They both should, roughly, reflect the same data, but tau selection requires both. - std::vector &populations = data->populations; unsigned int num_species = sim->model->number_species; unsigned int num_reactions = sim->model->number_reactions; @@ -303,17 +300,10 @@ int Gillespy::TauHybrid::rhs(realtype t, N_Vector y, N_Vector ydot, void *user_d // First half is for concentrations, second half is for reaction offsets. realtype *dydt_offsets = &dydt[num_species]; - // Populate the current ODE state into the concentrations vector. - // dy/dt results are initialized to zero, and become the change in propensity. - unsigned int spec_i; - for (spec_i = 0; spec_i < num_species; ++spec_i) - { - populations[spec_i] = static_cast(Y[spec_i]); - } - // Deterministic reactions generally are "evaluated" by generating dy/dt functions // for each of their dependent species. // To handle these, we will go ahead and evaluate each species' differential equations. + unsigned int spec_i; for (spec_i = 0; spec_i < num_species; ++spec_i) { if ((*species)[spec_i].boundary_condition) { @@ -322,7 +312,7 @@ int Gillespy::TauHybrid::rhs(realtype t, N_Vector y, N_Vector ydot, void *user_d } else { - dydt[spec_i] = (*species)[spec_i].diff_equation.evaluate(t, Y, populations.data()); + dydt[spec_i] = (*species)[spec_i].diff_equation.evaluate(t, Y); } } @@ -333,7 +323,7 @@ int Gillespy::TauHybrid::rhs(realtype t, N_Vector y, N_Vector ydot, void *user_d switch ((*reactions)[rxn_i].mode) { case SimulationState::DISCRETE: // Process stochastic reaction state by updating the root offset for each reaction. - propensity = Reaction::propensity(rxn_i, populations.data()); + propensity = Reaction::propensity(rxn_i, Y); dydt_offsets[rxn_i] = propensity; propensities[rxn_i] = propensity; break; From c5953957b61d07077707b5234b24285923e7427f Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Thu, 20 Jan 2022 13:17:34 -0500 Subject: [PATCH 09/54] Raise `SyntaxError` instead of returning `None` in `Expression` class --- gillespy2/solvers/cpp/build/expression.py | 32 +++++++++++++++-------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/gillespy2/solvers/cpp/build/expression.py b/gillespy2/solvers/cpp/build/expression.py index a12dda28f..63eb6f777 100644 --- a/gillespy2/solvers/cpp/build/expression.py +++ b/gillespy2/solvers/cpp/build/expression.py @@ -50,17 +50,17 @@ def __init__(self, self.namespace = dict({}) if namespace is None else namespace self.blacklist = dict({}) if blacklist is None else blacklist self.sanitize = sanitize - self.invalid_names = [] - self.invalid_operators = [] + self.invalid_names = set() + self.invalid_operators = set() def check_blacklist(self, operator): operator = type(operator) if operator in self.blacklist: - self.invalid_operators.append(str(self.blacklist.get(operator))) + self.invalid_operators.add(str(self.blacklist.get(operator))) def visit_Name(self, node: "ast.Name"): if node.id not in self.namespace: - self.invalid_names.append(node.id) + self.invalid_names.add(node.id) elif self.sanitize: node.id = self.namespace.get(node.id) self.generic_visit(node) @@ -68,7 +68,7 @@ def visit_Name(self, node: "ast.Name"): def visit_Call(self, node: "ast.Call"): if node.func.id not in self.namespace: - self.invalid_names.append(node.func.id) + self.invalid_names.add(node.func.id) elif self.sanitize: node.func.id = self.namespace.get(node.func.id) self.generic_visit(node) @@ -199,15 +199,25 @@ def validate(self, statement: "str") -> "ExpressionResults": return ExpressionResults(invalid_names=validator.invalid_names, invalid_operators=validator.invalid_operators) - def __get_expr(self, converter: "ExpressionConverter") -> "Optional[str]": + def __get_expr(self, statement: "str", converter: "ExpressionConverter") -> "Optional[str]": validator = Expression.ValidationVisitor(self.namespace, self.blacklist, self.sanitize) validator.visit(converter.tree) + failures_found = [] + if validator.invalid_operators: - return None + base_msg = "Blacklisted operator" + base_msg = f"{base_msg}s" if len(validator.invalid_operators) > 1 else base_msg + failures_found.append(f"{base_msg}: {','.join(validator.invalid_operators)}") if validator.invalid_names: - return None + base_msg = "Cannot resolve species name" + base_msg = f"{base_msg}s" if len(validator.invalid_names) > 1 else base_msg + failures_found.append(f"{base_msg}: {','.join(validator.invalid_names)}") + + if len(failures_found) > 0: + raise SyntaxError(f"Invalid GillesPy2 expression \"{statement}\"\n" + + "\n".join([f"* {msg}" for msg in failures_found])) return converter.get_str() @@ -220,7 +230,7 @@ def getexpr_python(self, statement: "str") -> "Optional[str]": :returns: Python expression string, if valid. Returns None if validation fails. """ expr = ast.parse(statement) - return self.__get_expr(PythonConverter(expr)) + return self.__get_expr(statement, PythonConverter(expr)) def getexpr_cpp(self, statement: "str") -> "Optional[str]": """ @@ -232,7 +242,7 @@ def getexpr_cpp(self, statement: "str") -> "Optional[str]": """ statement = ExpressionConverter.convert_str(statement) expr = ast.parse(statement) - return self.__get_expr(CppConverter(expr)) + return self.__get_expr(statement, CppConverter(expr)) class ExpressionResults: @@ -241,7 +251,7 @@ class ExpressionResults: Any expression items which indicate an invalid expression are listed on an ExpressionResults instance. Empty lists indicate that the expression is valid. """ - def __init__(self, invalid_names: "list[str]" = None, invalid_operators: "list[str]" = None, is_valid=True): + def __init__(self, invalid_names: "set[str]" = None, invalid_operators: "set[str]" = None, is_valid=True): """ Container struct for returning the results of expression validation. From 23c3e662435af75c7c8856d71d0e5528dc636636 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Thu, 20 Jan 2022 13:26:28 -0500 Subject: [PATCH 10/54] Add configurable list of "known function names" to `Expression` class --- gillespy2/solvers/cpp/build/template_gen.py | 7 +++++++ test/test_c_solvers.py | 11 ++++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/build/template_gen.py b/gillespy2/solvers/cpp/build/template_gen.py index e7e94e371..ad71c379d 100644 --- a/gillespy2/solvers/cpp/build/template_gen.py +++ b/gillespy2/solvers/cpp/build/template_gen.py @@ -39,6 +39,12 @@ class SanitizedModel: "t": "t", } + # Global functions that aren't present in the `math` package, + # as well as functions in Python that have a different name in C++. + function_map = { + "abs": "abs", + } + def __init__(self, model: Model, variable=False): self.model = model self.variable = variable @@ -68,6 +74,7 @@ def __init__(self, model: Model, variable=False): # All "system" namespace entries should always be first. # Otherwise, user-defined identifiers (like, for example, "gamma") might get overwritten. **{name: name for name in math.__dict__.keys()}, + **self.function_map, **self.species_names, **self.parameter_names, **self.reserved_names, diff --git a/test/test_c_solvers.py b/test/test_c_solvers.py index 236491cb4..c9b0d80b2 100644 --- a/test/test_c_solvers.py +++ b/test/test_c_solvers.py @@ -27,6 +27,7 @@ from gillespy2.solvers.cpp import SSACSolver, ODECSolver, TauLeapingCSolver from gillespy2.solvers.cpp import TauHybridCSolver from gillespy2.solvers.cpp.build.expression import Expression, ExpressionConverter +from gillespy2.solvers.cpp.build.template_gen import SanitizedModel class ExpressionTestCase: @@ -85,6 +86,10 @@ class TestCSolvers(unittest.TestCase): ExpressionTestCase({"x": "x", "y": "y", "z": "z"}, "(x^2/y^2/z^2)/x^2/y^2/z^2**1/x**1/y**1/z", [ [5.1, 0.1, 2.0], [0.1, 5.1, 2.0], [2.0, 0.1, 5.1], [2.0, 5.1, 0.1], ]), + # Known, builtin math expression functions work. + ExpressionTestCase({"x": "x"}, "abs(x)", [ + [100.0], [100], [-100.0], [-100], [0], + ]), ] comparisons = [ # Asserts that single comparison expressions work. @@ -189,7 +194,11 @@ def run(args: "list[str]") -> str: def test_expressions(expressions: "list[ExpressionTestCase]", use_bool=False): for entry in expressions: expression = ExpressionConverter.convert_str(entry.expression) - expr = Expression(namespace=entry.args) + expr = Expression(namespace={ + **SanitizedModel.reserved_names, + **SanitizedModel.function_map, + **entry.args, + }) cpp_expr = expr.getexpr_cpp(expression) with self.subTest(msg="Evaluating converted C expressions", expression=entry.expression, From 84cf4e7d03e0e5f2a9a93009726a01b4572593d7 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 17:07:43 -0800 Subject: [PATCH 11/54] catch errors if events are not setup correctly --- gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index a18a9b850..f7febcb88 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -108,6 +108,17 @@ def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel" assignments.append(str(assign_id)) event_assignment_list.append(assign_str) assign_id += 1 + # Check for "None"s + for a in assignments: + if a is None: raise Exception(f"a is Nonein event={event}") + if event_id is None: raise Exception(f"event_id is None in event={event}") + if trigger is None: raise Exception(f"trigger is None in event={event}") + if delay is None: raise Exception(f"delay is None in event={event}") + if priority is None: raise Exception(f"priority is None in event={event}") + if use_trigger is None: raise Exception(f"use_trigger is None in event={event}") + if use_persist is None: raise Exception(f"use_persist is None in event={event}") + if initial_value is None: raise Exception(f"initial_value is None in event={event}") + assignments: "str" = " AND ".join(assignments) event_list.append( f"EVENT(" From ac008344e33fa148839b4027e81649a6dadb12a5 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 17:18:00 -0800 Subject: [PATCH 12/54] cleanup --- gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index f7febcb88..236b58e9c 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -110,7 +110,7 @@ def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel" assign_id += 1 # Check for "None"s for a in assignments: - if a is None: raise Exception(f"a is Nonein event={event}") + if a is None: raise Exception(f"assignment={a} is None in event={event}") if event_id is None: raise Exception(f"event_id is None in event={event}") if trigger is None: raise Exception(f"trigger is None in event={event}") if delay is None: raise Exception(f"delay is None in event={event}") From 821353cc2b2782b7453582fd2f85910744c6365b Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 17:57:31 -0800 Subject: [PATCH 13/54] fixed error message --- gillespy2/core/events.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/core/events.py b/gillespy2/core/events.py index 9954a510d..8c385e040 100644 --- a/gillespy2/core/events.py +++ b/gillespy2/core/events.py @@ -70,7 +70,7 @@ def __init__(self, name=None, variable=None, expression=None): 'valid string expression') def __str__(self): - return self.variable.name + ': ' + self.expression + return f"{self.variable}: {self.expression}" class EventTrigger(Jsonify): """ @@ -261,4 +261,4 @@ def add_assignment(self, assignment): else: raise ModelError("Unexpected parameter for add_assignment. Parameter must be EventAssignment or list of " "EventAssignments") - return assignment \ No newline at end of file + return assignment From 28904865cc2f2091fdd34e73ca4bbada16ac64c6 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 18:01:46 -0800 Subject: [PATCH 14/54] improved error message --- gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index 236b58e9c..04d6372d4 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -93,7 +93,8 @@ def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel" elif variable in sanitized_model.model.listOfParameters: variable = sanitized_model.model.listOfParameters.get(variable) else: - raise ValueError(f"Invalid event assignment {assign}: received name {variable} " + raise ValueError(f"Error in event={event} " + f"Invalid event assignment {assign}: received name {variable} " f"Must match the name of a valid Species or Parameter.") if isinstance(variable, gillespy2.Species): From 4458adf20131c7392deb3b78d23a00a0d3016418 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 18:09:34 -0800 Subject: [PATCH 15/54] Raise error if element not found, otherwise return error string treated as Species obj --- gillespy2/core/model.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 00d05135c..fd04c5c3e 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -841,7 +841,8 @@ def get_element(self, ename): return self.get_assignment_rule(ename) if ename in self.listOfFunctionDefinitions: return self.get_function_definition(ename) - return 'Element not found!' + #return 'Element not found!' + raise ModelError(f"model.get_element(): element={ename} not found") def get_best_solver(self): From 6a99f9e6c20abe541cfaade864f91f2225ae21fa Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 17 Jan 2022 23:30:56 -0800 Subject: [PATCH 16/54] Fixing bug --- gillespy2/solvers/utilities/solverutils.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/gillespy2/solvers/utilities/solverutils.py b/gillespy2/solvers/utilities/solverutils.py index f5ea54fbe..668f40e42 100644 --- a/gillespy2/solvers/utilities/solverutils.py +++ b/gillespy2/solvers/utilities/solverutils.py @@ -19,6 +19,7 @@ import ast # for dependency graphing import numpy as np from gillespy2.core import log, Species +from gillespy2.core import ModelError """ NUMPY SOLVER UTILITIES BELOW @@ -143,8 +144,11 @@ def species_parse(model, custom_prop_fun): class SpeciesParser(ast.NodeTransformer): def visit_Name(self, node): - if isinstance(model.get_element(node.id), Species): - parsed_species.append(model.get_element(node.id)) + try: + if isinstance(model.get_element(node.id), Species): + parsed_species.append(model.get_element(node.id)) + except ModelError: + pass expr = custom_prop_fun expr = ast.parse(expr, mode='eval') From 70cd5c0d60c17bf0daa84a1c483c7efa4a46a16c Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Wed, 26 Jan 2022 11:30:33 -0500 Subject: [PATCH 17/54] Remove dead code from comment --- gillespy2/core/model.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index fd04c5c3e..0c881d77e 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -841,7 +841,6 @@ def get_element(self, ename): return self.get_assignment_rule(ename) if ename in self.listOfFunctionDefinitions: return self.get_function_definition(ename) - #return 'Element not found!' raise ModelError(f"model.get_element(): element={ename} not found") From 062642b1e35b2e270b46bea79cfbdb65a5ef9381 Mon Sep 17 00:00:00 2001 From: seanebum Date: Thu, 27 Jan 2022 16:16:01 -0500 Subject: [PATCH 18/54] added TauHybridCSolver into get_best_solver model method and gave it priority over the numpy version --- gillespy2/core/model.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 0c881d77e..eb36426fb 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -867,18 +867,21 @@ def get_best_solver(self): if tempMode == 'dynamic' or tempMode == 'continuous': hybrid_check = True break - if can_use_numpy and hybrid_check: + + from gillespy2.solvers.cpp.build.build_engine import BuildEngine + can_use_cpp = not len(BuildEngine.get_missing_dependencies()) + + if can_use_cpp and hybrid_check: + from gillespy2 import TauHybridCSolver + return TauHybridCSolver + elif can_use_numpy and hybrid_check: from gillespy2 import TauHybridSolver return TauHybridSolver - elif not can_use_numpy and hybrid_check: raise ModelError('TauHybridSolver is the only solver currently that supports ' 'AssignmentRules, RateRules, FunctionDefinitions, or Events. ' 'Please install Numpy.') - from gillespy2.solvers.cpp.build.build_engine import BuildEngine - can_use_cpp = not len(BuildEngine.get_missing_dependencies()) - if can_use_cpp is False and can_use_numpy and not hybrid_check: from gillespy2 import NumPySSASolver return NumPySSASolver From f4410a7e984c82993b439e9cc04a5e24e77db6f4 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 10:31:31 -0500 Subject: [PATCH 19/54] Add solver name to error message in `_validate_sbml_features()` --- gillespy2/solvers/cpp/c_solver.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/c_solver.py b/gillespy2/solvers/cpp/c_solver.py index cb4595c83..0a403de27 100644 --- a/gillespy2/solvers/cpp/c_solver.py +++ b/gillespy2/solvers/cpp/c_solver.py @@ -372,7 +372,9 @@ def _validate_sbml_features(self, unsupported_features: "dict[str, str]"): detected.append(feature_name) if len(detected): - raise gillespyError.ModelError(f"Could not run Model. SBML Feature: {detected} not supported by SSACSolver.") + raise gillespyError.ModelError(f"Could not run Model, " + f"SBML Features not supported by {self.name}: " + + ", ".join(detected)) def _validate_seed(self, seed: int): if seed is None: From 5183bd4a3fdb6b5beba3b0a5f6011ef9f2ad4027 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 11:54:36 -0500 Subject: [PATCH 20/54] Add all C++ solvers to behavioral tests - Timeout not working on all solvers - Wrapped each solver in a subtest --- test/test_all_solvers.py | 49 +++++++++++++++++++++++++--------------- 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 8e13f7228..4523d164c 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -27,11 +27,23 @@ from gillespy2 import NumPySSASolver from gillespy2 import TauLeapingSolver from gillespy2 import TauHybridSolver +from gillespy2 import ODECSolver +from gillespy2 import TauLeapingCSolver +from gillespy2 import TauHybridCSolver class TestAllSolvers(unittest.TestCase): - solvers = [SSACSolver, ODESolver, NumPySSASolver, TauLeapingSolver, TauHybridSolver] + solvers = [ + SSACSolver, + ODESolver, + NumPySSASolver, + TauLeapingSolver, + TauHybridSolver, + ODECSolver, + TauLeapingCSolver, + TauHybridCSolver, + ] model = Example() for sp in model.listOfSpecies.values(): @@ -65,15 +77,15 @@ def test_return_type_show_labels(self): self.assertTrue(isinstance(self.labeled_results_more_trajectories[solver][0]['Sp'], np.ndarray)) self.assertTrue(isinstance(self.labeled_results_more_trajectories[solver][0]['Sp'][0], np.float)) - def test_random_seed(self): for solver in self.solvers: - same_results = self.model.run(solver=solver, seed=1) - compare_results = self.model.run(solver=solver,seed=1) - self.assertTrue(np.array_equal(same_results.to_array(), compare_results.to_array())) - if solver.name == 'ODESolver': continue - diff_results = self.model.run(solver=solver, seed=2) - self.assertFalse(np.array_equal(diff_results.to_array(),same_results.to_array())) + with self.subTest(solver=solver.name): + same_results = self.model.run(solver=solver, seed=1) + compare_results = self.model.run(solver=solver,seed=1) + self.assertTrue(np.array_equal(same_results.to_array(), compare_results.to_array())) + if solver.name in ["ODESolver", "ODECSolver"]: continue + diff_results = self.model.run(solver=solver, seed=2) + self.assertFalse(np.array_equal(diff_results.to_array(), same_results.to_array())) def test_random_seed_unnamed_reactions(self): model = self.model @@ -82,25 +94,26 @@ def test_random_seed_unnamed_reactions(self): unnamed_rxn = gillespy2.Reaction(reactants={}, products={'Sp':1}, rate=k2) model.add_reaction(unnamed_rxn) for solver in self.solvers: - same_results = self.model.run(solver=solver, seed=1) - compare_results = self.model.run(solver=solver,seed=1) - self.assertTrue(np.array_equal(same_results.to_array(), compare_results.to_array())) - if solver.name == 'ODESolver': continue - diff_results = self.model.run(solver=solver, seed=2) - self.assertFalse(np.array_equal(diff_results.to_array(),same_results.to_array())) + with self.subTest(solver=solver.name): + same_results = self.model.run(solver=solver, seed=1) + compare_results = self.model.run(solver=solver,seed=1) + self.assertTrue(np.array_equal(same_results.to_array(), compare_results.to_array())) + if solver.name in ["ODESolver", "ODECSolver"]: continue + diff_results = self.model.run(solver=solver, seed=2) + self.assertFalse(np.array_equal(diff_results.to_array(), same_results.to_array())) def test_extraneous_args(self): for solver in self.solvers: - with self.assertLogs(level='WARN'): + with self.subTest(solver=solver.name), self.assertLogs(level='WARN'): model = Example() - results = model.run(solver=solver, nonsense='ABC') + model.run(solver=solver, nonsense='ABC') def test_timeout(self): for solver in self.solvers: - with self.assertLogs(level='WARN'): + with self.subTest(solver=solver.name), self.assertLogs(level='WARN'): model = Oregonator() model.timespan(np.linspace(0, 1000000, 101)) - results = model.run(solver=solver, timeout=1) + model.run(solver=solver, timeout=1) def test_basic_solver_import(self): from gillespy2.solvers.numpy.basic_tau_leaping_solver import BasicTauLeapingSolver From d09e450f0c9f903e1eeb76956fcbe9a8b71d4971 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 12:19:09 -0500 Subject: [PATCH 21/54] Macro-defined signal handlers - Added `GPY_INTERRUPT_*` macros for cross-platform signal handling - Added macro calls to each solver --- gillespy2/solvers/cpp/c_base/model.h | 14 +++++++++- .../cpp/c_base/ode_cpp_solver/ODESolver.cpp | 6 +++++ .../cpp/c_base/ssa_cpp_solver/SSASolver.cpp | 27 +++---------------- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 15 +++++------ .../TauLeapingSolver.cpp | 9 +++---- 5 files changed, 33 insertions(+), 38 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/model.h b/gillespy2/solvers/cpp/c_base/model.h index 3837e96ca..a5b909179 100644 --- a/gillespy2/solvers/cpp/c_base/model.h +++ b/gillespy2/solvers/cpp/c_base/model.h @@ -26,12 +26,24 @@ #include #include -#if defined(WIN32) || defined(WIN32) || defined(__WIN32) || defined(__WIN32__) +#if defined(WIN32) || defined(_WIN32) || defined(__WIN32) || defined(__WIN32__) #include #define GPY_PID_GET() ((int) GetCurrentProcessId()) +#define GPY_INTERRUPT_HANDLER(handler_name, handler_code) \ +static BOOL WINAPI handler_name(DWORD signum) { \ + do handler_code while(0); \ + return TRUE; \ +} +#define GPY_INTERRUPT_INSTALL_HANDLER(handler) SetConsoleCtrlHandler(handler, TRUE) #else #include +#include #define GPY_PID_GET() (getpid()) +#define GPY_INTERRUPT_HANDLER(handler_name, handler_code) \ +static void handler_name(int signum) { \ + do handler_code while(0); \ +} +#define GPY_INTERRUPT_INSTALL_HANDLER(handler) signal(SIGINT, signalHandler) #endif namespace Gillespy diff --git a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp index 33016c3c6..1e72fce9b 100644 --- a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp +++ b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp @@ -32,6 +32,12 @@ namespace Gillespy { + static volatile bool interrupted = false; + + GPY_INTERRUPT_HANDLER(signal_handler, { + interrupted = true; + }) + static int f(realtype t, N_Vector y, N_Vector y_dot, void *user_data); struct UserData diff --git a/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASolver.cpp b/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASolver.cpp index bdf4a8386..1fe068b5a 100644 --- a/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASolver.cpp +++ b/gillespy2/solvers/cpp/c_base/ssa_cpp_solver/SSASolver.cpp @@ -18,40 +18,19 @@ #include #include -#include -#include - -#ifdef _WIN32 -#include -#undef max -#endif - #include "SSASolver.h" namespace Gillespy { volatile bool interrupted = false; -#ifdef _WIN32 - BOOL WINAPI eventHandler(DWORD CtrlType) - { - interrupted = true; - return TRUE; - } -#endif - - void signalHandler(int signum) - { + GPY_INTERRUPT_HANDLER(signal_handler, { interrupted = true; - } + }) void ssa_direct(Simulation *simulation) { -#ifdef _WIN32 - SetConsoleCtrlHandler(eventHandler, TRUE); -#else - signal(SIGINT, signalHandler); -#endif + GPY_INTERRUPT_INSTALL_HANDLER(signal_handler); if (simulation) { diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index e28b9f320..e8667221f 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -34,22 +34,21 @@ namespace Gillespy { - namespace TauHybrid - { + static volatile bool interrupted = false; - bool interrupted = false; - - void signalHandler(int signum) - { - interrupted = true; - } + GPY_INTERRUPT_HANDLER(signal_handler, { + interrupted = true; + }) + namespace TauHybrid + { void TauHybridCSolver(HybridSimulation *simulation, std::vector &events, const double tau_tol) { if (simulation == NULL) { return; } + GPY_INTERRUPT_INSTALL_HANDLER(signal_handler); Model &model = *(simulation->model); int num_species = model.number_species; diff --git a/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSolver.cpp index 2a150c3b2..0d7813a5b 100644 --- a/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_leaping_cpp_solver/TauLeapingSolver.cpp @@ -33,13 +33,12 @@ namespace Gillespy { - bool interrupted = false; + static volatile bool interrupted = false; std::mt19937_64 generator; - void signalHandler(int signum) - { + GPY_INTERRUPT_HANDLER(signal_handler, { interrupted = true; - } + }) std::pair, double> get_reactions( const Gillespy::Model *model, @@ -78,7 +77,7 @@ namespace Gillespy void tau_leaper(Gillespy::Simulation *simulation, const double tau_tol) { - signal(SIGINT, signalHandler); + GPY_INTERRUPT_INSTALL_HANDLER(signal_handler); if (!simulation) { From 9aebde67a5ed3edf22bb3ea7aabcc2ff752429e7 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 13:21:34 -0500 Subject: [PATCH 22/54] Add interrupt handlers to all solvers - Reduced timeout in unit test to 0.1s - `ODECSolver` completes in <1s - Interrupt checks added to `TauHybridCSolver` --- .../solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp | 4 +++- .../c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp | 11 ++++++++--- test/test_all_solvers.py | 4 ++-- 3 files changed, 13 insertions(+), 6 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp index 1e72fce9b..127b3e63c 100644 --- a/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp +++ b/gillespy2/solvers/cpp/c_base/ode_cpp_solver/ODESolver.cpp @@ -51,6 +51,8 @@ namespace Gillespy void ODESolver(Simulation *simulation, double increment) { + GPY_INTERRUPT_INSTALL_HANDLER(signal_handler); + // CVODE constants are returned on every success or failure. // CV_SUCCESS: Operation was successful. // CV_MEM_NULL: CVODE memory block was not initialized with CVodeCreate. @@ -120,7 +122,7 @@ namespace Gillespy realtype tret = 0; int current_time = 0; - for (tout = step_length; tout < end_time || cmpf(tout, end_time); tout += step_length) + for (tout = step_length; !interrupted && tout < end_time || cmpf(tout, end_time); tout += step_length) { // CV_NORMAL causes the solver to take internal steps until it has reached or just passed the `tout` // parameter. The solver interpolates in order to return an approximate value of `y(tout)`. diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index e8667221f..1de5fde87 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -76,7 +76,7 @@ namespace Gillespy TauArgs tau_args = initialize(model, tau_tol); // Simulate for each trajectory - for (int traj = 0; traj < num_trajectories; traj++) + for (int traj = 0; !interrupted && traj < num_trajectories; traj++) { if (traj > 0) { @@ -135,7 +135,7 @@ namespace Gillespy // For now, a "guard" is put in place to prevent potentially infinite loops from occurring. unsigned int integration_guard = 1000; - while (integration_guard > 0 && simulation->current_time < simulation->end_time) + while (!interrupted && integration_guard > 0 && simulation->current_time < simulation->end_time) { // Compute current propensity values based on existing state. for (unsigned int rxn_i = 0; rxn_i < num_reactions; ++rxn_i) @@ -155,6 +155,8 @@ namespace Gillespy } sol.data.propensities[rxn_i] = propensity; } + if (interrupted) + break; // Expected tau step is determined. tau_step = select( @@ -303,7 +305,10 @@ namespace Gillespy current_populations[p_i] = (int) current_state[p_i]; } } - } while (invalid_state); + } while (invalid_state && !interrupted); + + if (interrupted) + break; // Invalid state after the do-while loop implies that an unrecoverable error has occurred. // While prior results are considered usable, the current integration results are not. diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 4523d164c..b3b26339b 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -112,8 +112,8 @@ def test_timeout(self): for solver in self.solvers: with self.subTest(solver=solver.name), self.assertLogs(level='WARN'): model = Oregonator() - model.timespan(np.linspace(0, 1000000, 101)) - model.run(solver=solver, timeout=1) + model.timespan(np.linspace(0, 1000000, 1001)) + model.run(solver=solver, timeout=0.1) def test_basic_solver_import(self): from gillespy2.solvers.numpy.basic_tau_leaping_solver import BasicTauLeapingSolver From 3b36c499f1762e1d74a24550b57a697a617aab63 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 13:31:51 -0500 Subject: [PATCH 23/54] Fix typo in interrupt handler macro --- gillespy2/solvers/cpp/c_base/model.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/c_base/model.h b/gillespy2/solvers/cpp/c_base/model.h index a5b909179..f9c1878c1 100644 --- a/gillespy2/solvers/cpp/c_base/model.h +++ b/gillespy2/solvers/cpp/c_base/model.h @@ -43,7 +43,7 @@ static BOOL WINAPI handler_name(DWORD signum) { \ static void handler_name(int signum) { \ do handler_code while(0); \ } -#define GPY_INTERRUPT_INSTALL_HANDLER(handler) signal(SIGINT, signalHandler) +#define GPY_INTERRUPT_INSTALL_HANDLER(handler) signal(SIGINT, handler) #endif namespace Gillespy From 1a3530a9c397db9a48af42d82129a9cad852d3ec Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 28 Jan 2022 15:21:42 -0500 Subject: [PATCH 24/54] Update output buffering in `TauHybridCSolver` - Removed all references to deprecated `trajectories` array in C++ solver - Replaced writes to `trajectories` with `output_buffer()` methods to `TauHybridCSolver` --- gillespy2/solvers/cpp/c_base/model.cpp | 23 ++----------------- gillespy2/solvers/cpp/c_base/model.h | 2 -- .../TauHybridSimulation.cpp | 2 +- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 9 +++++--- 4 files changed, 9 insertions(+), 27 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/model.cpp b/gillespy2/solvers/cpp/c_base/model.cpp index 41132f262..f2ece4da1 100644 --- a/gillespy2/solvers/cpp/c_base/model.cpp +++ b/gillespy2/solvers/cpp/c_base/model.cpp @@ -88,19 +88,6 @@ namespace Gillespy { void init_simulation(Model *model, Simulation &simulation) { init_timeline(model, simulation); - unsigned int trajectory_size = simulation.number_timesteps * (model->number_species); - simulation.trajectories_1D = new TNum[simulation.number_trajectories * trajectory_size]; - simulation.trajectories = new TNum * *[simulation.number_trajectories]; - - for (unsigned int trajectory = 0; trajectory < simulation.number_trajectories; trajectory++) { - simulation.trajectories[trajectory] = new TNum * [simulation.number_timesteps]; - - for (unsigned int timestep = 0; timestep < simulation.number_timesteps; timestep++) { - simulation.trajectories[trajectory][timestep] = - &(simulation.trajectories_1D[trajectory * trajectory_size + timestep * (model->number_species)]); - } - } - simulation.current_state = new TNum[model->number_species]; // Output interval must lie within the range (0, num_timesteps]. // An output interval of 0 signifies to output entire trajectories. @@ -113,13 +100,7 @@ namespace Gillespy { template Simulation::~Simulation() { delete[] timeline; - delete[] trajectories_1D; - - for (unsigned int trajectory = 0; trajectory < number_trajectories; trajectory++) { - delete[] trajectories[trajectory]; - } - - delete[] trajectories; + delete[] current_state; } template @@ -147,7 +128,7 @@ namespace Gillespy { os << timeline[timestep] << ','; for (int species = 0; species < model->number_species; species++) { - os << (double) trajectories[trajectory][timestep][species] << ','; + os << (double) current_state[species] << ','; } } } diff --git a/gillespy2/solvers/cpp/c_base/model.h b/gillespy2/solvers/cpp/c_base/model.h index 3837e96ca..e9a5c7463 100644 --- a/gillespy2/solvers/cpp/c_base/model.h +++ b/gillespy2/solvers/cpp/c_base/model.h @@ -145,8 +145,6 @@ namespace Gillespy double end_time; double *timeline; - PType *trajectories_1D; - PType ***trajectories; PType *current_state; Model *model; diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp index d13878b09..fde42c463 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp @@ -71,6 +71,6 @@ int main(int argc, char* argv[]) Gillespy::TauHybrid::Event::use_events(events); TauHybrid::TauHybridCSolver(&simulation, events, tau_tol); - simulation.output_results_buffer(std::cout); + simulation.output_buffer_final(std::cout); return 0; } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index e28b9f320..c5713ab89 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -94,8 +94,11 @@ namespace Gillespy for (int spec_i = 0; spec_i < num_species; ++spec_i) { current_state[spec_i] = species[spec_i].initial_population; + simulation->current_state[spec_i] = current_state[spec_i]; current_populations[spec_i] = species[spec_i].initial_population; } + simulation->reset_output_buffer(traj); + simulation->output_buffer_range(std::cout); // Check for initial event triggers at t=0 (based on initial_value of trigger) std::set event_roots; @@ -115,7 +118,6 @@ namespace Gillespy spec->partition_mode = spec->user_mode == SimulationState::DYNAMIC ? SimulationState::DISCRETE : spec->user_mode; - simulation->trajectories[traj][0][spec_i] = current_state[spec_i]; } // SIMULATION STEP LOOP @@ -344,9 +346,10 @@ namespace Gillespy { for (int spec_i = 0; spec_i < num_species; ++spec_i) { - simulation->trajectories[traj][save_idx][spec_i] = current_state[spec_i]; + simulation->current_state[spec_i] = current_state[spec_i]; } - save_time = simulation->timeline[++save_idx]; + simulation->output_buffer_range(std::cout, save_idx++); + save_time = simulation->timeline[save_idx]; } } From 86f3b3e26458f43c6f267bd8f993ed506f55926d Mon Sep 17 00:00:00 2001 From: seanebum Date: Mon, 31 Jan 2022 09:42:06 -0500 Subject: [PATCH 25/54] fixed solver check to select python hybrid solver in presence of assignment rules or function definitions. Fixed conflicts in unit tests where selection appropriately changed from TauHybridSolver to TauHybridCSolver --- gillespy2/core/model.py | 9 ++++++--- test/test_hybrid_solver.py | 10 +++++----- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index eb36426fb..d2407e11c 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -857,9 +857,12 @@ def get_best_solver(self): """ from gillespy2.solvers.numpy import can_use_numpy hybrid_check = False - if len(self.get_all_assignment_rules()) or len(self.get_all_rate_rules()) \ - or len(self.get_all_function_definitions()) or len(self.get_all_events()): + chybrid_check = True + if len(self.get_all_rate_rules()) or len(self.get_all_events()): hybrid_check = True + if len(self.get_all_assignment_rules()) or len(self.get_all_function_definitions()): + hybrid_check = True + chybrid_check = False if len(self.get_all_species()) and hybrid_check == False: for i in self.get_all_species(): @@ -871,7 +874,7 @@ def get_best_solver(self): from gillespy2.solvers.cpp.build.build_engine import BuildEngine can_use_cpp = not len(BuildEngine.get_missing_dependencies()) - if can_use_cpp and hybrid_check: + if can_use_cpp and hybrid_check and chybrid_check: from gillespy2 import TauHybridCSolver return TauHybridCSolver elif can_use_numpy and hybrid_check: diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index a3c391433..e583ecde7 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -41,7 +41,7 @@ def test_add_rate_rule(self): model.add_species([species]) model.add_rate_rule([rule]) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridSolver') + self.assertEqual(results[0].solver_name, 'TauHybridCSolver') def test_add_rate_rule_dict(self): model = Example() @@ -88,7 +88,7 @@ def test_add_continuous_species_dependent_event(self): event1.add_assignment([ea1, ea2]) model.add_event(event1) results = model.run() - self.assertEqual(results[0].solver_name,'TauHybridSolver') + self.assertEqual(results[0].solver_name, 'TauHybridCSolver') self.assertEqual(results['Sp'][-1], 1000) def test_add_stochastic_species_dependent_event(self): @@ -174,21 +174,21 @@ def test_ensure_hybrid_dynamic_species(self): species1 = gillespy2.Species('test_species1', initial_value=1,mode='dynamic') model.add_species(species1) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridSolver') + self.assertEqual(results[0].solver_name, 'TauHybridCSolver') def test_ensure_hybrid_continuous_species(self): model = Example() species1 = gillespy2.Species('test_species1', initial_value=1,mode='continuous') model.add_species(species1) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridSolver') + self.assertEqual(results[0].solver_name, 'TauHybridCSolver') def test_ensure_continuous_dynamic_timeout_warning(self): model = Example() species1 = gillespy2.Species('test_species1', initial_value=1, mode='dynamic') model.add_species(species1) with self.assertLogs(level='WARN'): - results = model.run(timeout=1) + results = model.run(solver=TauHybridSolver, timeout=1) def test_run_example__with_increment_only(self): model = ExampleNoTspan() From 27adc4aba82bd97f3e9907e4289b76d7fd82bcf0 Mon Sep 17 00:00:00 2001 From: seanebum Date: Mon, 31 Jan 2022 14:25:23 -0500 Subject: [PATCH 26/54] updating get_best_algo code --- gillespy2/core/model.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index d2407e11c..6be30e0f6 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -900,6 +900,9 @@ def get_best_solver_algo(self, algorithm): from gillespy2.solvers.numpy import can_use_numpy from gillespy2.solvers.cpp.build.build_engine import BuildEngine can_use_cpp = not len(BuildEngine.get_missing_dependencies()) + chybrid_check = True + if len(self.get_all_assignment_rules()) or len(self.get_all_function_definitions()): + chybrid_check = False if not can_use_cpp and can_use_numpy: raise ModelError("Please install C++ or Numpy to use GillesPy2 solvers.") @@ -927,6 +930,15 @@ def get_best_solver_algo(self, algorithm): else: from gillespy2 import ODESolver return ODESolver + + elif algorithm == 'Tau-Hybrid': + if can_use_cpp and chybrid_check: + from gillespy2 import TauHybridCSolver + return TauHybridCSolver + else: + from gillespy2 import TauHybridSolver + return TauHybridSolver + else: raise ModelError("Invalid value for the argument 'algorithm' entered. " "Please enter 'SSA', 'ODE', or 'Tau-leaping'.") From fb4519e567cb34ef60a45b53802c11b526620ccf Mon Sep 17 00:00:00 2001 From: seanebum Date: Mon, 31 Jan 2022 15:09:27 -0500 Subject: [PATCH 27/54] updating error-catching in solver-selection --- gillespy2/core/model.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 6be30e0f6..c5f1c0afa 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -874,16 +874,15 @@ def get_best_solver(self): from gillespy2.solvers.cpp.build.build_engine import BuildEngine can_use_cpp = not len(BuildEngine.get_missing_dependencies()) + if not can_use_cpp and not can_use_numpy: + raise ModelError('Dependency Error, cannot run model.') + if can_use_cpp and hybrid_check and chybrid_check: from gillespy2 import TauHybridCSolver return TauHybridCSolver elif can_use_numpy and hybrid_check: from gillespy2 import TauHybridSolver return TauHybridSolver - elif not can_use_numpy and hybrid_check: - raise ModelError('TauHybridSolver is the only solver currently that supports ' - 'AssignmentRules, RateRules, FunctionDefinitions, or Events. ' - 'Please install Numpy.') if can_use_cpp is False and can_use_numpy and not hybrid_check: from gillespy2 import NumPySSASolver @@ -931,7 +930,7 @@ def get_best_solver_algo(self, algorithm): from gillespy2 import ODESolver return ODESolver - elif algorithm == 'Tau-Hybrid': + elif algorithm == 'Hybrid': if can_use_cpp and chybrid_check: from gillespy2 import TauHybridCSolver return TauHybridCSolver From 0780448eff4d3d814b99a6c14a5ba9082c2c8730 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 2 Feb 2022 14:58:42 -0500 Subject: [PATCH 28/54] Set BuildEngine.clean() to ignore errors. Prepended the tmpdir with 'gillespy2_build_'. Prepended the obj and template directories with 'gillespy2_'. Prepended the executable file with 'GillesPy2_'. --- gillespy2/solvers/cpp/build/build_engine.py | 10 ++++++---- gillespy2/solvers/cpp/build/make.py | 6 +++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/gillespy2/solvers/cpp/build/build_engine.py b/gillespy2/solvers/cpp/build/build_engine.py index 5a9423a59..f92d428a5 100644 --- a/gillespy2/solvers/cpp/build/build_engine.py +++ b/gillespy2/solvers/cpp/build/build_engine.py @@ -91,13 +91,15 @@ def prepare(self, model: "Union[Model, template_gen.SanitizedModel]", variable=F # If the output directory is None, create and set it to a temporary directory. if self.output_dir is None: - self.output_dir = Path(tempfile.mkdtemp()) + self.output_dir = Path(tempfile.mkdtemp( + prefix='gillespy2_build_', dir=os.environ.get('GILLESPY2_TMPDIR')) + ) # If object files haven't been compiled yet, go ahead and compile them with make. # Precompilation only happens if the cache is enabled but hasn't been built yet. # Make target for individual simulation will still succeed if prebuild() isn't called. - self.obj_dir = self.output_dir.joinpath("obj") - self.template_dir = self.output_dir.joinpath("template") + self.obj_dir = self.output_dir.joinpath("gillespy2_obj") + self.template_dir = self.output_dir.joinpath("gillespy2_template") if gillespy2.cache_enabled: self.obj_dir = self.__get_cache_dir() @@ -186,4 +188,4 @@ def clean(self): return if self.output_dir.exists(): - shutil.rmtree(self.output_dir) \ No newline at end of file + shutil.rmtree(self.output_dir, ignore_errors=True) \ No newline at end of file diff --git a/gillespy2/solvers/cpp/build/make.py b/gillespy2/solvers/cpp/build/make.py index 80113cafe..e6d7fb7d6 100644 --- a/gillespy2/solvers/cpp/build/make.py +++ b/gillespy2/solvers/cpp/build/make.py @@ -36,7 +36,7 @@ def __init__(self, makefile: str, output_dir: str, obj_dir: str = None): # If not supplied, it should be assumed ethereal and cleaned up # with the rest of the tmp directory. if obj_dir is None: - self.obj_dir = self.output_dir.joinpath("obj").resolve() + self.obj_dir = self.output_dir.joinpath("gillespy2_obj").resolve() else: self.obj_dir = Path(obj_dir).resolve() @@ -50,9 +50,9 @@ def __init__(self, makefile: str, output_dir: str, obj_dir: str = None): if not path.is_dir(): path.mkdir() - self.output_file = "Simulation.out" + self.output_file = "GillesPy2_Simulation.out" if os.name == "nt": - self.output_file = "Simulation.exe" + self.output_file = "GillesPy2_Simulation.exe" self.output_file = Path(self.output_dir, self.output_file) From e2c92beb00778ab9f9d0148d83e3297c82151f42 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 2 Feb 2022 15:05:54 -0500 Subject: [PATCH 29/54] Added script to manually cleanup temp files. --- gillespy2/core/__init__.py | 1 + gillespy2/core/cleanup.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 gillespy2/core/cleanup.py diff --git a/gillespy2/core/__init__.py b/gillespy2/core/__init__.py index fb0458621..2959a97ad 100644 --- a/gillespy2/core/__init__.py +++ b/gillespy2/core/__init__.py @@ -18,6 +18,7 @@ import logging from .assignmentrule import * +from .cleanup import * from .events import * from .functiondefinition import * from .gillespyError import * diff --git a/gillespy2/core/cleanup.py b/gillespy2/core/cleanup.py new file mode 100644 index 000000000..5da96f534 --- /dev/null +++ b/gillespy2/core/cleanup.py @@ -0,0 +1,30 @@ +""" +GillesPy2 is a modeling toolkit for biochemical simulation. +Copyright (C) 2019-2021 GillesPy2 developers. + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +""" + +import os +import shutil +import tempfile + +def cleanup_tempfiles(): + ''' + Cleanup all tempfiles in gillespy2 build. + ''' + tempdir = tempfile.gettempdir() + for file_obj in os.listdir(tempdir): + if file_obj.startswith("gillespy2_build"): + cleanup_build_files(build_dir=os.path.join(tempdir, file_obj)) From bcdaabcbd46802b871203d1f0802d128d900316f Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 2 Feb 2022 15:12:02 -0500 Subject: [PATCH 30/54] Fixed broken test. --- test/test_build_engine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_build_engine.py b/test/test_build_engine.py index c33db36fd..595fcd6c3 100644 --- a/test/test_build_engine.py +++ b/test/test_build_engine.py @@ -104,7 +104,7 @@ def test_build_output(self): Should result in an expected definitions file and an executable simulation. """ simulation_file = self.build_engine.output_dir.joinpath( - "Simulation.exe" if os.name == "nt" else "Simulation.out") + "GillesPy2_Simulation.exe" if os.name == "nt" else "GillesPy2_Simulation.out") for solver_name in self.solver_names: with self.subTest(solver=solver_name): From e16b181a308877feda59d6ee4b9242b1e4276894 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 4 Feb 2022 12:34:21 -0500 Subject: [PATCH 31/54] Add unit test for continuous-valued species --- test/test_hybrid_solver.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index a3c391433..206720098 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -204,5 +204,36 @@ def test_run_example__with_tspan_and_increment(self): results = TauHybridSolver.run(model=model, increment=0.2) +class TestAllHybridSolvers(unittest.TestCase): + from gillespy2.solvers import TauHybridCSolver + + class TruncatedStateModel(gillespy2.Model): + def __init__(self): + gillespy2.Model.__init__(self, name="TruncatedStateModel") + S1 = gillespy2.Species(name="S1", initial_value=0, mode="discrete") + rate = gillespy2.Species(name="rate", initial_value=0.9999, mode="continuous") + self.add_species([S1, rate]) + self.add_rate_rule(gillespy2.RateRule(variable="rate", formula="-1/((t+0.9999)**2)")) + self.add_reaction( + # Because S1 is a "discrete" species, our reaction will be marked "stochastic." + gillespy2.Reaction(products={S1: 1}, propensity_function="10.0*rate") + ) + solvers = [ + TauHybridSolver, + TauHybridCSolver, + ] + + def test_continuous_state_values(self): + """ + Continuous values should be evaluate appropriately, without being truncated/casted to an integer. + """ + model = TestAllHybridSolvers.TruncatedStateModel() + for solver in self.solvers: + with self.subTest(solver=solver.name): + result = model.run(solver=solver, seed=1) + self.assertGreater(result["S1"][-1], 0.0, + "Reaction never fired; indicates that continuous species is being truncated") + + if __name__ == '__main__': unittest.main() From 67c528e4a260d42ff441b308a05f43d6eb6f604a Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 4 Feb 2022 12:54:19 -0500 Subject: [PATCH 32/54] Modify Hybrid Tau select to accept `double` as state - Does not change `TauLeapingCSolver`; uses different template instantiation --- gillespy2/solvers/cpp/c_base/Tau/tau.cpp | 8 ++++---- gillespy2/solvers/cpp/c_base/Tau/tau.h | 8 ++++---- .../cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp | 4 ++-- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/Tau/tau.cpp b/gillespy2/solvers/cpp/c_base/Tau/tau.cpp index 5cbfe3627..b49e4ea4a 100644 --- a/gillespy2/solvers/cpp/c_base/Tau/tau.cpp +++ b/gillespy2/solvers/cpp/c_base/Tau/tau.cpp @@ -109,7 +109,7 @@ namespace Gillespy return tau_args; } - template + template double select( Gillespy::Model &model, TauArgs &tau_args, @@ -117,7 +117,7 @@ namespace Gillespy const double ¤t_time, const double &save_time, const std::vector &propensity_values, - const std::vector ¤t_state) + const std::vector ¤t_state) { bool critical = false; // system-wide flag, true when any reaction is critical @@ -261,14 +261,14 @@ namespace Gillespy // Explicitly instantiate initialize/select functions for DISCRETE simulations template TauArgs initialize(Model &model, double tau_tol); - template double select( + template double select( Model &model, TauArgs &tau_args, const double &tau_tol, const double ¤t_time, const double &save_time, const std::vector &propensity_values, - const std::vector ¤t_state); + const std::vector ¤t_state); // Explicitly instantiate initialize/select functions for HYBRID simulations template TauArgs initialize(Model &model, double tau_tol); diff --git a/gillespy2/solvers/cpp/c_base/Tau/tau.h b/gillespy2/solvers/cpp/c_base/Tau/tau.h index de3f5a54f..6a8cb22b8 100644 --- a/gillespy2/solvers/cpp/c_base/Tau/tau.h +++ b/gillespy2/solvers/cpp/c_base/Tau/tau.h @@ -47,7 +47,7 @@ namespace Gillespy template TauArgs initialize(Model &model, double tau_tol); - template + template double select( Model &model, TauArgs &tau_args, @@ -55,19 +55,19 @@ namespace Gillespy const double ¤t_time, const double &save_time, const std::vector &propensity_values, - const std::vector ¤t_state); + const std::vector ¤t_state); // Continuous tau args is used by hybrid solver template struct TauArgs; extern template TauArgs initialize(Model &model, double tau_tol); - extern template double select( + extern template double select( Model &model, TauArgs &tau_args, const double &tau_tol, const double ¤t_time, const double &save_time, const std::vector &propensity_values, - const std::vector ¤t_state); + const std::vector ¤t_state); // Discrete tau args is used by tau leaping solver template struct TauArgs; diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index e28b9f320..198e42d20 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -158,14 +158,14 @@ namespace Gillespy } // Expected tau step is determined. - tau_step = select( + tau_step = select( model, tau_args, tau_tol, simulation->current_time, save_time, sol.data.propensities, - current_populations + current_state ); partition_species( simulation->reaction_state, From c9188cbbc30681c73bf06ffe1634b5a544da799e Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 4 Feb 2022 13:51:09 -0500 Subject: [PATCH 33/54] Add explicit `ode` and `ssa` propensity methods - Resolved by each reaction based on which "mode" that reaction is in - Removed all instances/references to `populations` vector --- .../tau_hybrid_cpp_solver/HybridModel.cpp | 27 ++++--- .../tau_hybrid_cpp_solver/HybridModel.h | 54 ++++++++++++- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 17 +--- .../tau_hybrid_cpp_solver/integrator.cpp | 17 ++-- .../c_base/tau_hybrid_cpp_solver/integrator.h | 3 - .../solvers/cpp/c_base/template/template.cpp | 78 +++++++++++++++++++ .../solvers/cpp/c_base/template/template.h | 7 ++ 7 files changed, 161 insertions(+), 42 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp index 81ff83799..5c16f48e2 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp @@ -203,11 +203,19 @@ namespace Gillespy HybridReaction::HybridReaction() : mode(SimulationState::DISCRETE), - base_reaction(nullptr) + m_base_reaction(nullptr), + m_id(-1) { // Empty constructor body } + HybridReaction::HybridReaction(Reaction *base_reaction) + : HybridReaction() + { + base_reaction = base_reaction; + m_id = base_reaction->id; + } + HybridSpecies::HybridSpecies() : user_mode(SimulationState::DYNAMIC), partition_mode(SimulationState::DISCRETE), @@ -235,7 +243,7 @@ namespace Gillespy for (int rxn_i = 0; rxn_i < model.number_reactions; ++rxn_i) { - reaction_state[rxn_i].base_reaction = &model.reactions[rxn_i]; + reaction_state[rxn_i].set_base_reaction(&model.reactions[rxn_i]); } } @@ -272,9 +280,8 @@ namespace Gillespy spec.diff_equation.formulas.clear(); } - for (int rxn_i = 0; rxn_i < reactions.size(); ++rxn_i) + for (HybridReaction &rxn : reactions) { - HybridReaction rxn = reactions[rxn_i]; if (rxn.mode == SimulationState::DISCRETE) { continue; @@ -283,17 +290,17 @@ namespace Gillespy for (int spec_i = 0; spec_i < species.size(); ++spec_i) { // A species change of 0 indicates that this species is not a dependency for this reaction. - if (rxn.base_reaction->species_change[spec_i] == 0) + if (rxn.get_base_reaction()->species_change[spec_i] == 0) { continue; } HybridSpecies &spec = species[spec_i]; auto &formula_set = spec.diff_equation.formulas; - int spec_diff = rxn.base_reaction->species_change[spec_i]; + int spec_diff = rxn.get_base_reaction()->species_change[spec_i]; - formula_set.push_back([rxn_i, spec_diff](double *state) { - return spec_diff * Reaction::propensity(rxn_i, state); + formula_set.emplace_back([&rxn, spec_diff](double *state) { + return spec_diff * rxn.propensity(state); }); } } @@ -321,7 +328,7 @@ namespace Gillespy { // 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.base_reaction->species_change[spec_i] == 0) + if (rxn.get_base_reaction()->species_change[spec_i] == 0) { continue; } @@ -389,7 +396,7 @@ namespace Gillespy // 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.base_reaction->species_change[spec_i]; + int spec_dx = rxn.get_base_reaction()->species_change[spec_i]; if (spec_dx < 0) { // Selected species is a reactant. diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h index da680d44e..5bf24c5f9 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h @@ -227,10 +227,62 @@ namespace Gillespy struct HybridReaction { - Reaction *base_reaction; SimulationState mode; + inline double propensity(double *state) const + { + return propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get()); + } + + inline double propensity(double *state, double *parameters, const double *constants) const + { + switch (mode) + { + case SimulationState::CONTINUOUS: + return ode_propensity(state, parameters, constants); + case SimulationState::DISCRETE: + default: + return ssa_propensity(state, parameters, constants); + } + } + + inline double ode_propensity(double *state) const + { + return ode_propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get()); + } + + inline double ode_propensity(double *state, double *parameters, const double *constants) const + { + return map_ode_propensity(m_id, state, parameters, constants); + } + + inline double ssa_propensity(double *state) const + { + return ssa_propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get()); + } + + inline double ssa_propensity(double *state, double *parameters, const double *constants) const + { + return map_ssa_propensity(m_id, state, parameters, constants); + } + + inline void set_base_reaction(Reaction *reaction) + { + m_base_reaction = reaction; + m_id = reaction->id; + } + + inline Reaction *get_base_reaction() const + { + return m_base_reaction; + } + HybridReaction(); + explicit HybridReaction(Reaction *base_reaction); + + private: + Reaction *m_base_reaction; + ReactionId m_id; }; struct HybridSimulation : Simulation diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index 198e42d20..5847e0fce 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -88,13 +88,11 @@ namespace Gillespy // TODO: change back double -> hybrid_state, once we figure out how that works EventList event_list; std::vector current_state(num_species); - std::vector current_populations(num_species); // Initialize the species population for the trajectory. for (int spec_i = 0; spec_i < num_species; ++spec_i) { current_state[spec_i] = species[spec_i].initial_population; - current_populations[spec_i] = species[spec_i].initial_population; } // Check for initial event triggers at t=0 (based on initial_value of trigger) @@ -142,19 +140,7 @@ namespace Gillespy for (unsigned int rxn_i = 0; rxn_i < num_reactions; ++rxn_i) { HybridReaction &rxn = simulation->reaction_state[rxn_i]; - double propensity = 0.0; - switch (rxn.mode) - { - case SimulationState::CONTINUOUS: - propensity = Reaction::propensity(rxn_i, current_state.data()); - break; - case SimulationState::DISCRETE: - propensity = Reaction::propensity(rxn_i, current_populations.data()); - break; - default: - break; - } - sol.data.propensities[rxn_i] = propensity; + sol.data.propensities[rxn_i] = rxn.propensity(current_state.data()); } // Expected tau step is determined. @@ -301,7 +287,6 @@ namespace Gillespy current_state[p_i] += population_changes[p_i]; result.concentrations[p_i] = current_state[p_i]; } - current_populations[p_i] = (int) current_state[p_i]; } } } while (invalid_state); diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index a69974429..a80c8da91 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -27,8 +27,6 @@ IntegratorData::IntegratorData( int num_species, int num_reactions) : simulation(simulation), - concentrations(std::vector(num_species)), - populations(std::vector(num_species)), propensities(std::vector(num_reactions)), species_state(&simulation->species_state), reaction_state(&simulation->reaction_state) {} @@ -55,13 +53,6 @@ Integrator::Integrator(HybridSimulation *simulation, N_Vector y0, double reltol, NV_Ith_S(y, mem_i) = NV_Ith_S(this->y0, mem_i); } - for (int spec_i = 0; spec_i < num_species; ++spec_i) - { - data.populations[spec_i] - = data.concentrations[spec_i] - = simulation->model->species[spec_i].initial_population; - } - for (int rxn_i = 0; rxn_i < num_reactions; ++rxn_i) { data.propensities[rxn_i] = 0; @@ -197,7 +188,7 @@ void Integrator::use_reactions(const std::vector &reactions) // Reaction root-finder should only be used on discrete-valued reactions. // The required IDs are placed into a reference vector and are mapped back out // when the caller of integrate() retrieves them. - data.active_reaction_ids.push_back(reaction.base_reaction->id); + data.active_reaction_ids.push_back(reaction.get_base_reaction()->id); } } } @@ -320,10 +311,12 @@ int Gillespy::TauHybrid::rhs(realtype t, N_Vector y, N_Vector ydot, void *user_d // These updates get written directly to the integrator's concentration state for (unsigned int rxn_i = 0; rxn_i < num_reactions; ++rxn_i) { - switch ((*reactions)[rxn_i].mode) { + HybridReaction rxn = (*reactions)[rxn_i]; + + switch (rxn.mode) { case SimulationState::DISCRETE: // Process stochastic reaction state by updating the root offset for each reaction. - propensity = Reaction::propensity(rxn_i, Y); + propensity = rxn.ssa_propensity(Y); dydt_offsets[rxn_i] = propensity; propensities[rxn_i] = propensity; break; diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h index f0a3a0f88..6d8cd9d7f 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h @@ -58,9 +58,6 @@ namespace Gillespy // In `rootfn`, this means that gout[i] is the "output" of reaction active_reaction_ids[i]. // This is used to map the internal reaction number to the actual reaction id. std::vector active_reaction_ids; - - std::vector concentrations; - std::vector populations; std::vector propensities; IntegratorData(HybridSimulation *simulation); diff --git a/gillespy2/solvers/cpp/c_base/template/template.cpp b/gillespy2/solvers/cpp/c_base/template/template.cpp index 19418c342..b581a40b5 100644 --- a/gillespy2/solvers/cpp/c_base/template/template.cpp +++ b/gillespy2/solvers/cpp/c_base/template/template.cpp @@ -88,6 +88,84 @@ namespace Gillespy return constants; } + double map_ssa_propensity(unsigned int reaction_id, int *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + + double map_ssa_propensity(unsigned int reaction_id, unsigned int *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + + double map_ssa_propensity(unsigned int reaction_id, double *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + + double map_ode_propensity(unsigned int reaction_id, int *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_ODE_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + + double map_ode_propensity(unsigned int reaction_id, unsigned int *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_ODE_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + + double map_ode_propensity(unsigned int reaction_id, double *S, double *P, const double *C) + { + switch (reaction_id) + { + #define PROPENSITY(id, func) case(id): return(func); + GPY_ODE_PROPENSITIES + #undef PROPENSITY + + default: + return -1.0; + } + } + double map_propensity(unsigned int reaction_id, unsigned int *S, double *P, const double *C) { switch (reaction_id) diff --git a/gillespy2/solvers/cpp/c_base/template/template.h b/gillespy2/solvers/cpp/c_base/template/template.h index d64436515..9c4e7a01a 100644 --- a/gillespy2/solvers/cpp/c_base/template/template.h +++ b/gillespy2/solvers/cpp/c_base/template/template.h @@ -38,6 +38,13 @@ namespace Gillespy double map_propensity(unsigned int reaction_id, unsigned int *S, double *parameters, const double *constants); double map_propensity(unsigned int reaction_id, double *S, double *parameters, const double *constants); + double map_ssa_propensity(unsigned int reaction_id, int *state, double *parameters, const double *constants); + double map_ssa_propensity(unsigned int reaction_id, unsigned int *state, double *parameters, const double *constants); + double map_ssa_propensity(unsigned int reaction_id, double *state, double *parameters, const double *constants); + double map_ode_propensity(unsigned int reaction_id, int *state, double *parameters, const double *constants); + double map_ode_propensity(unsigned int reaction_id, unsigned int *state, double *parameters, const double *constants); + double map_ode_propensity(unsigned int reaction_id, double *state, double *parameters, const double *constants); + template void add_reactions(Model &model); From 902c10a3695077e133a5621c9e258b1778310110 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 10 Feb 2022 13:48:36 -0500 Subject: [PATCH 34/54] Fixed incorrect exception raising. --- gillespy2/core/model.py | 48 +++++++------------ gillespy2/solvers/numpy/CLE_solver.py | 4 +- gillespy2/solvers/numpy/ode_solver.py | 4 +- gillespy2/solvers/numpy/ssa_solver.py | 4 +- gillespy2/solvers/numpy/tau_hybrid_solver.py | 4 +- gillespy2/solvers/numpy/tau_leaping_solver.py | 4 +- 6 files changed, 32 insertions(+), 36 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 0c881d77e..40055de4d 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -281,28 +281,28 @@ def sanitized_species_names(self): def problem_with_name(self, name): if name in Model.reserved_names: - return ModelError( + raise ModelError( 'Name "{}" is unavailable. It is reserved for internal GillesPy use. Reserved Names: ({}).'.format(name, Model.reserved_names)) if name in self.listOfSpecies: - return ModelError('Name "{}" is unavailable. A species with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. A species with that name exists.'.format(name)) if name in self.listOfParameters: - return ModelError('Name "{}" is unavailable. A parameter with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. A parameter with that name exists.'.format(name)) if name in self.listOfReactions: - return ModelError('Name "{}" is unavailable. A reaction with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. A reaction with that name exists.'.format(name)) if name in self.listOfEvents: - return ModelError('Name "{}" is unavailable. An event with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. An event with that name exists.'.format(name)) if name in self.listOfRateRules: - return ModelError('Name "{}" is unavailable. A rate rule with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. A rate rule with that name exists.'.format(name)) if name in self.listOfAssignmentRules: - return ModelError('Name "{}" is unavailable. An assignment rule with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. An assignment rule with that name exists.'.format(name)) if name in self.listOfFunctionDefinitions: - return ModelError('Name "{}" is unavailable. A function definition with that name exists.'.format(name)) + raise ModelError('Name "{}" is unavailable. A function definition with that name exists.'.format(name)) if name.isdigit(): - return ModelError('Name "{}" is unavailable. Names must not be numeric strings.'.format(name)) + raise ModelError('Name "{}" is unavailable. Names must not be numeric strings.'.format(name)) for special_character in Model.special_characters: if special_character in name: - return ModelError( + raise ModelError( 'Name "{}" is unavailable. Names must not contain special characters: {}.'.format(name, Model.special_characters)) @@ -335,9 +335,7 @@ def add_species(self, obj): self.add_species(S) else: try: - problem = self.problem_with_name(obj.name) - if problem is not None: - raise problem + self.problem_with_name(obj.name) self.listOfSpecies[obj.name] = obj self._listOfSpecies[obj.name] = 'S{}'.format(len(self._listOfSpecies)) except Exception as e: @@ -418,9 +416,7 @@ def add_parameter(self, params): self.add_parameter(p) else: if isinstance(params, Parameter) or type(params).__name__ == 'Parameter': - problem = self.problem_with_name(params.name) - if problem is not None: - raise problem + self.problem_with_name(params.name) self.update_namespace() params._evaluate(self.namespace) self.listOfParameters[params.name] = params @@ -501,9 +497,7 @@ def add_reaction(self, reactions): self.add_reaction(r) else: try: - problem = self.problem_with_name(reactions.name) - if problem is not None: - raise problem + self.problem_with_name(reactions.name) reactions.verify() self.validate_reactants_and_products(reactions) if reactions.name is None or reactions.name == '': @@ -543,9 +537,7 @@ def add_rate_rule(self, rate_rules): self.add_rate_rule(rr) else: try: - problem = self.problem_with_name(rate_rules.name) - if problem is not None: - raise problem + self.problem_with_name(rate_rules.name) if len(self.listOfAssignmentRules) != 0: for i in self.listOfAssignmentRules.values(): if rate_rules.variable == i.variable: @@ -588,9 +580,7 @@ def add_event(self, event): self.add_event(e) else: try: - problem = self.problem_with_name(event.name) - if problem is not None: - raise problem + self.problem_with_name(event.name) if event.trigger is None or not hasattr(event.trigger, 'expression'): raise ModelError( 'An Event must contain a valid trigger.') @@ -615,9 +605,7 @@ def add_function_definition(self, function_definitions): self.add_function_definition(fd) else: try: - problem = self.problem_with_name(function_definitions.name) - if problem is not None: - raise problem + self.problem_with_name(function_definitions.name) self.listOfFunctionDefinitions[function_definitions.name] = function_definitions except Exception as e: raise ParameterError( @@ -635,9 +623,7 @@ def add_assignment_rule(self, assignment_rules): self.add_assignment_rule(ar) else: try: - problem = self.problem_with_name(assignment_rules.name) - if problem is not None: - raise problem + self.problem_with_name(assignment_rules.name) if len(self.listOfRateRules) != 0: for i in self.listOfRateRules.values(): if assignment_rules.variable == i.variable: diff --git a/gillespy2/solvers/numpy/CLE_solver.py b/gillespy2/solvers/numpy/CLE_solver.py index 511cfbbf5..7e25e53e6 100644 --- a/gillespy2/solvers/numpy/CLE_solver.py +++ b/gillespy2/solvers/numpy/CLE_solver.py @@ -250,7 +250,9 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N while self.result is None: pass if hasattr(self, 'has_raised_exception'): - raise self.has_raised_exception + raise SimulationError( + f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" + ) from self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) diff --git a/gillespy2/solvers/numpy/ode_solver.py b/gillespy2/solvers/numpy/ode_solver.py index 957346903..8b8419dda 100644 --- a/gillespy2/solvers/numpy/ode_solver.py +++ b/gillespy2/solvers/numpy/ode_solver.py @@ -205,7 +205,9 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, show_l while self.result is None: pass if hasattr(self, 'has_raised_exception'): - raise self.has_raised_exception + raise SimulationError( + f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" + ) from self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) diff --git a/gillespy2/solvers/numpy/ssa_solver.py b/gillespy2/solvers/numpy/ssa_solver.py index 46c1984cf..61fd5f6e9 100644 --- a/gillespy2/solvers/numpy/ssa_solver.py +++ b/gillespy2/solvers/numpy/ssa_solver.py @@ -169,7 +169,9 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N while self.result is None: pass if hasattr(self, 'has_raised_exception'): - raise self.has_raised_exception + raise SimulationError( + f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" + ) from self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 921e9becc..13d0599aa 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -938,7 +938,9 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N except: pass if hasattr(self, 'has_raised_exception'): - raise self.has_raised_exception + raise SimulationError( + f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" + ) from self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) diff --git a/gillespy2/solvers/numpy/tau_leaping_solver.py b/gillespy2/solvers/numpy/tau_leaping_solver.py index a54766a47..cbc0d8c7b 100644 --- a/gillespy2/solvers/numpy/tau_leaping_solver.py +++ b/gillespy2/solvers/numpy/tau_leaping_solver.py @@ -237,7 +237,9 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N while self.result is None: pass if hasattr(self, 'has_raised_exception'): - raise self.has_raised_exception + raise SimulationError( + f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" + ) from self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) From af0d25e12029ca13a2572878dd5725ed596c0f31 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 10 Feb 2022 12:37:55 -0800 Subject: [PATCH 35/54] better names --- .../run_version_comparison.py} | 0 .../tyson_oscillator.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename test/{version_performance_tests/run_version_test.py => version_performance_comparison/run_version_comparison.py} (100%) rename test/{version_performance_tests => version_performance_comparison}/tyson_oscillator.py (100%) diff --git a/test/version_performance_tests/run_version_test.py b/test/version_performance_comparison/run_version_comparison.py similarity index 100% rename from test/version_performance_tests/run_version_test.py rename to test/version_performance_comparison/run_version_comparison.py diff --git a/test/version_performance_tests/tyson_oscillator.py b/test/version_performance_comparison/tyson_oscillator.py similarity index 100% rename from test/version_performance_tests/tyson_oscillator.py rename to test/version_performance_comparison/tyson_oscillator.py From af7802d88b27474f002737f29013f065668bbc81 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 10 Feb 2022 12:45:48 -0800 Subject: [PATCH 36/54] code cleanup --- .../run_version_comparison.py | 60 +++++++++---------- 1 file changed, 27 insertions(+), 33 deletions(-) diff --git a/test/version_performance_comparison/run_version_comparison.py b/test/version_performance_comparison/run_version_comparison.py index 5022bdd0a..aa2f22380 100755 --- a/test/version_performance_comparison/run_version_comparison.py +++ b/test/version_performance_comparison/run_version_comparison.py @@ -7,41 +7,35 @@ versions = [ ] -cmd3 = f"git clone https://github.com/StochSS/GillesPy2" -print(cmd3) -output3 = subprocess.run(cmd3, shell=True, capture_output=True) - -cmd4 = f"cd GillesPy2; git tag" -print(cmd4) -output4 = subprocess.run(cmd4, shell=True, capture_output=True) -versions = output4.stdout.decode("utf-8").split("\n") - - def SortFn(ver): x = re.findall("\d+",ver) return int(x[0])*10000+int(x[1])*100+int(x[2]) -while("" in versions) : - versions.remove("") - -versions.sort(key=SortFn) - -#print(f"versions={versions}") - - -for ver in versions: - print(f"{ver}: ",end='') - cmd2 = f"cd GillesPy2; git checkout {ver}" - output2 = subprocess.run(cmd2, shell=True, capture_output=True) - cmd = f'export PYTHONPATH="{os.getcwd()}/GillesPy2"; python3 ./tyson_oscillator.py' - #print(cmd) - - output = subprocess.run(cmd, shell=True, capture_output=True) - out = output.stdout.decode("utf-8").strip() - try: - print(float(out)) - except: - print("error"); - #print(output) - +if __name__=="__main__": + cmd3 = f"git clone https://github.com/StochSS/GillesPy2" + print(cmd3) + output3 = subprocess.run(cmd3, shell=True, capture_output=True) + + cmd4 = f"cd GillesPy2; git tag" + print(cmd4) + output4 = subprocess.run(cmd4, shell=True, capture_output=True) + versions = output4.stdout.decode("utf-8").split("\n") + + while("" in versions) : + versions.remove("") + + versions.sort(key=SortFn) + + for ver in versions: + print(f"{ver}: ",end='') + cmd2 = f"cd GillesPy2; git checkout {ver}" + output2 = subprocess.run(cmd2, shell=True, capture_output=True) + cmd = f'export PYTHONPATH="{os.getcwd()}/GillesPy2"; python3 ./tyson_oscillator.py' + + output = subprocess.run(cmd, shell=True, capture_output=True) + out = output.stdout.decode("utf-8").strip() + try: + print(float(out)) + except: + print("error"); From 65c1aeb9297d999730a50ec32bb53e502d376da1 Mon Sep 17 00:00:00 2001 From: seanebum Date: Thu, 10 Feb 2022 15:59:30 -0500 Subject: [PATCH 37/54] set algorithm back to Tau-Hybrid and updated error string --- gillespy2/core/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index c5f1c0afa..df09c95bf 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -930,7 +930,7 @@ def get_best_solver_algo(self, algorithm): from gillespy2 import ODESolver return ODESolver - elif algorithm == 'Hybrid': + elif algorithm == 'Tau-Hybrid': if can_use_cpp and chybrid_check: from gillespy2 import TauHybridCSolver return TauHybridCSolver @@ -940,7 +940,7 @@ def get_best_solver_algo(self, algorithm): else: raise ModelError("Invalid value for the argument 'algorithm' entered. " - "Please enter 'SSA', 'ODE', or 'Tau-leaping'.") + "Please enter 'SSA', 'ODE', 'Tau-leaping', or 'Tau-Hybrid'.") def run(self, solver=None, timeout=0, t=None, increment=None, show_labels=True, cpp_support=False, algorithm=None, **solver_args): From e02831a5127b03b0ac8c2c61cddae0504509cb66 Mon Sep 17 00:00:00 2001 From: seanebum Date: Thu, 10 Feb 2022 16:08:17 -0500 Subject: [PATCH 38/54] updated unit tests to pass for either hybrid solver in test_hybrid_solver.py --- test/test_hybrid_solver.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index e583ecde7..f9d88954f 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -41,7 +41,8 @@ def test_add_rate_rule(self): model.add_species([species]) model.add_rate_rule([rule]) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridCSolver') + valid_solvers = ('TauHybridSolver', 'TauHybridCSolver') + self.assertIn(results[0].solver_name, valid_solvers) def test_add_rate_rule_dict(self): model = Example() @@ -88,7 +89,8 @@ def test_add_continuous_species_dependent_event(self): event1.add_assignment([ea1, ea2]) model.add_event(event1) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridCSolver') + valid_solvers = ('TauHybridSolver', 'TauHybridCSolver') + self.assertIn(results[0].solver_name, valid_solvers) self.assertEqual(results['Sp'][-1], 1000) def test_add_stochastic_species_dependent_event(self): @@ -174,14 +176,16 @@ def test_ensure_hybrid_dynamic_species(self): species1 = gillespy2.Species('test_species1', initial_value=1,mode='dynamic') model.add_species(species1) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridCSolver') + valid_solvers = ('TauHybridSolver', 'TauHybridCSolver') + self.assertIn(results[0].solver_name, valid_solvers) def test_ensure_hybrid_continuous_species(self): model = Example() species1 = gillespy2.Species('test_species1', initial_value=1,mode='continuous') model.add_species(species1) results = model.run() - self.assertEqual(results[0].solver_name, 'TauHybridCSolver') + valid_solvers = ('TauHybridSolver', 'TauHybridCSolver') + self.assertIn(results[0].solver_name, valid_solvers) def test_ensure_continuous_dynamic_timeout_warning(self): model = Example() From 60ab5f90e05b158a6a048ca2bae3e8799cb8094f Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Thu, 10 Feb 2022 20:33:26 -0500 Subject: [PATCH 39/54] Eval trigger state *after* events have been evaluated --- .../cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp index 5c16f48e2..e36691c4e 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp @@ -505,14 +505,11 @@ namespace Gillespy // Step 1: Identify any fired event triggers. for (auto &event : m_events) { - bool trigger = event.trigger(t, event_state); - if (m_trigger_state.at(event.get_event_id()) != trigger) + if (m_trigger_state.at(event.get_event_id()) != event.trigger(t, event_state)) { double delay = event.delay(t, event_state); EventExecution execution = event.get_execution(t + delay, event_state, output_size); - // Update trigger state to prevent repeated firings. - m_trigger_state.find(event.get_event_id())->second = trigger; if (delay <= 0) { // Immediately put EventExecution on "triggered" pile @@ -590,6 +587,13 @@ namespace Gillespy m_trigger_pool.erase(event.get_event_id()); } + // Step 5: Update trigger states based on the new simulation state. + // This is to account for events that re-assign values that the event triggers depend on. + for (auto &event : m_events) + { + m_trigger_state.at(event.get_event_id()) = event.trigger(t, event_state); + } + return has_active_events(); } } From 70813aa03f394e4b1f5c311b67ed1d690ef8c882 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 11 Feb 2022 18:46:10 -0500 Subject: [PATCH 40/54] Fix missing "invalid state" check for C++ events --- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index c50e195f4..367ef5741 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -251,10 +251,6 @@ namespace Gillespy { population_changes[spec_i] += model.reactions[rxn_i].species_change[spec_i]; - if (current_state[spec_i] + population_changes[spec_i] < 0) - { - invalid_state = true; - } } rxn_state += log(urn.next()); @@ -269,6 +265,16 @@ namespace Gillespy } } + // Explicitly check for invalid population state, now that changes have been tallied. + for (int spec_i = 0; spec_i < num_species; ++spec_i) + { + if (current_state[spec_i] + population_changes[spec_i] < 0) + { + invalid_state = true; + break; + } + } + // Positive reaction state means a negative population was detected. // Only update state with the given population changes if valid. if (invalid_state) From 0ae31ee050cd3990c4ec8b6e9d6c0f24617fc7c4 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 11 Feb 2022 18:57:38 -0500 Subject: [PATCH 41/54] Fix initialization of hybrid data structures - Removed all default constructors; should not be able to create a hybrid species/reaction/simulation without a base --- .../tau_hybrid_cpp_solver/HybridModel.cpp | 26 +++++-------------- .../tau_hybrid_cpp_solver/HybridModel.h | 10 ++++--- .../solvers/cpp/c_base/template/template.cpp | 1 + 3 files changed, 14 insertions(+), 23 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp index 5c16f48e2..139b7ae75 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp @@ -200,27 +200,17 @@ namespace Gillespy return m_execution_time > rhs.m_execution_time; } - - HybridReaction::HybridReaction() - : mode(SimulationState::DISCRETE), - m_base_reaction(nullptr), - m_id(-1) - { - // Empty constructor body - } - HybridReaction::HybridReaction(Reaction *base_reaction) - : HybridReaction() + : mode(SimulationState::DISCRETE), + m_base_reaction(base_reaction) { - base_reaction = base_reaction; m_id = base_reaction->id; } - HybridSpecies::HybridSpecies() + HybridSpecies::HybridSpecies(Species *base_species) : user_mode(SimulationState::DYNAMIC), partition_mode(SimulationState::DISCRETE), - switch_tol(0.03), - switch_min(0) + m_base_species(base_species) { // Empty constructor body } @@ -232,18 +222,16 @@ namespace Gillespy } HybridSimulation::HybridSimulation(const Model &model) - : Simulation(), - species_state(model.number_species), - reaction_state(model.number_reactions) + : Simulation() { for (int spec_i = 0; spec_i < model.number_species; ++spec_i) { - species_state[spec_i].base_species = &model.species[spec_i]; + species_state.emplace_back(HybridSpecies(model.species.get() + spec_i)); } for (int rxn_i = 0; rxn_i < model.number_reactions; ++rxn_i) { - reaction_state[rxn_i].set_base_reaction(&model.reactions[rxn_i]); + reaction_state.emplace_back(HybridReaction(model.reactions.get() + rxn_i)); } } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h index 5bf24c5f9..a14428988 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.h @@ -192,8 +192,6 @@ namespace Gillespy struct HybridSpecies { - Species *base_species; - // allows the user to specify if a species' population should definitely be modeled continuously or // discretely // CONTINUOUS or DISCRETE @@ -222,7 +220,11 @@ namespace Gillespy // If `boundary_condition` is true, then reactants are not consumed, and products are not produced. bool boundary_condition = false; - HybridSpecies(); + HybridSpecies() = delete; + explicit HybridSpecies(Species *species); + + private: + Species *m_base_species; }; struct HybridReaction @@ -277,7 +279,7 @@ namespace Gillespy return m_base_reaction; } - HybridReaction(); + HybridReaction() = delete; explicit HybridReaction(Reaction *base_reaction); private: diff --git a/gillespy2/solvers/cpp/c_base/template/template.cpp b/gillespy2/solvers/cpp/c_base/template/template.cpp index b581a40b5..1b73a9a40 100644 --- a/gillespy2/solvers/cpp/c_base/template/template.cpp +++ b/gillespy2/solvers/cpp/c_base/template/template.cpp @@ -237,6 +237,7 @@ namespace Gillespy for (rxn_i = 0; rxn_i < GPY_NUM_REACTIONS; ++rxn_i) { + model.reactions[rxn_i].id = rxn_i; for (spec_i = 0; spec_i < GPY_NUM_SPECIES; ++spec_i) { model.reactions[rxn_i].species_change[spec_i] = reactions[rxn_i][spec_i]; From 29751aaf7245bb567775e8b4ca2b4f75ca3fb503 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Fri, 11 Feb 2022 19:22:02 -0500 Subject: [PATCH 42/54] Silent CVODE errors in hybrid C++ solver - Can be re-enabled with `debug` flag in Python solver object --- gillespy2/solvers/cpp/c_base/arg_parser.cpp | 9 +++++++++ gillespy2/solvers/cpp/c_base/arg_parser.h | 2 ++ .../tau_hybrid_cpp_solver/TauHybridSimulation.cpp | 2 +- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 15 ++++++++++++++- .../tau_hybrid_cpp_solver/TauHybridSolver.h | 2 +- .../c_base/tau_hybrid_cpp_solver/integrator.cpp | 5 +++++ .../cpp/c_base/tau_hybrid_cpp_solver/integrator.h | 2 ++ gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 2 ++ 8 files changed, 36 insertions(+), 3 deletions(-) diff --git a/gillespy2/solvers/cpp/c_base/arg_parser.cpp b/gillespy2/solvers/cpp/c_base/arg_parser.cpp index ba4002239..14b993aac 100644 --- a/gillespy2/solvers/cpp/c_base/arg_parser.cpp +++ b/gillespy2/solvers/cpp/c_base/arg_parser.cpp @@ -25,6 +25,11 @@ char ArgParser::match_arg(std::string &token) { + if (!token.compare("--verbose")) + { + return 'v'; + } + if (!token.compare("--timesteps")) { return 't'; @@ -154,6 +159,10 @@ ArgParser::ArgParser(int argc, char *argv[]) arg_stream >> output_interval; break; + case 'v': + verbose = true; + break; + default: std::cerr << usage << std::endl; exit(0); diff --git a/gillespy2/solvers/cpp/c_base/arg_parser.h b/gillespy2/solvers/cpp/c_base/arg_parser.h index a4bd8748f..2a8ded5bd 100644 --- a/gillespy2/solvers/cpp/c_base/arg_parser.h +++ b/gillespy2/solvers/cpp/c_base/arg_parser.h @@ -47,6 +47,8 @@ class ArgParser double switch_tol = 0.0; double tau_tol = 0.03; + bool verbose = false; + ArgParser(int argc, char *argv[]); ~ArgParser(); }; \ No newline at end of file diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp index fde42c463..829a5bc87 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp @@ -70,7 +70,7 @@ int main(int argc, char* argv[]) std::vector events; Gillespy::TauHybrid::Event::use_events(events); - TauHybrid::TauHybridCSolver(&simulation, events, tau_tol); + TauHybrid::TauHybridCSolver(&simulation, events, tau_tol, parser.verbose); simulation.output_buffer_final(std::cout); return 0; } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp index 367ef5741..f74eae6c9 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp @@ -32,6 +32,9 @@ #include "integrator.h" #include "tau.h" +static void silent_error_handler(int error_code, const char *module, const char *function_name, + char *message, void *eh_data); + namespace Gillespy { static volatile bool interrupted = false; @@ -42,7 +45,7 @@ namespace Gillespy namespace TauHybrid { - void TauHybridCSolver(HybridSimulation *simulation, std::vector &events, const double tau_tol) + void TauHybridCSolver(HybridSimulation *simulation, std::vector &events, const double tau_tol, bool verbose) { if (simulation == NULL) { @@ -71,6 +74,10 @@ namespace Gillespy y = y0; } Integrator sol(simulation, y, GPY_HYBRID_RELTOL, GPY_HYBRID_ABSTOL); + if (!verbose) + { + sol.set_error_handler(silent_error_handler); + } // Tau selector initialization. Used to select a valid tau step. TauArgs tau_args = initialize(model, tau_tol); @@ -362,3 +369,9 @@ namespace Gillespy } } } + +void silent_error_handler(int error_code, const char *module, const char *function_name, + char *message, void *eh_data) +{ + // Do nothing +} diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.h index 319ceb59a..ad4b3f251 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.h @@ -25,6 +25,6 @@ namespace Gillespy { namespace TauHybrid { - void TauHybridCSolver(HybridSimulation* simulation, std::vector &events, double tau_tol); + void TauHybridCSolver(HybridSimulation* simulation, std::vector &events, double tau_tol, bool verbose); } } diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index a80c8da91..ad8d70e2b 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -212,6 +212,11 @@ bool Integrator::disable_root_finder() return validate(this, CVodeRootInit(cvode_mem, 0, NULL)); } +void Integrator::set_error_handler(CVErrHandlerFn error_handler) +{ + validate(this, CVodeSetErrHandlerFn(cvode_mem, error_handler, nullptr)); +} + URNGenerator::URNGenerator(unsigned long long seed) : uniform(0, 1), rng(seed) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h index 6d8cd9d7f..b570ae596 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.h @@ -151,6 +151,8 @@ namespace Gillespy /// and the underlying SBML event data and reaction data are removed. bool disable_root_finder(); + void set_error_handler(CVErrHandlerFn error_handler); + IntegrationResults integrate(double *t); IntegrationResults integrate(double *t, std::set &event_roots, std::set &reaction_roots); IntegratorData data; diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index 04d6372d4..825ba2751 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -219,6 +219,8 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int display_args = None args = self._make_args(args) + if debug: + args.append("--verbose") decoder = IterativeSimDecoder.create_default(number_of_trajectories, number_timesteps, len(self.model.listOfSpecies)) sim_exec = self._build(self.model, self.target, self.variable, False) From 534b32aef6f0b3d0b844a3a37bb4957464a97460 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Sat, 12 Feb 2022 14:10:28 -0500 Subject: [PATCH 43/54] Add unit test for SBML feature validation --- test/test_all_solvers.py | 57 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 3 deletions(-) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index b3b26339b..37be0d2b2 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -45,6 +45,39 @@ class TestAllSolvers(unittest.TestCase): TauHybridCSolver, ] + sbml_features = { + "AssignmentRule": lambda model, variable: + model.add_assignment_rule(gillespy2.AssignmentRule(variable=variable, formula="1/(t+1)")), + "RateRule": lambda model, variable: + model.add_rate_rule(gillespy2.RateRule(variable=variable, formula="2*t")), + "Event": lambda model, variable: + model.add_event(gillespy2.Event( + trigger=gillespy2.EventTrigger(expression="t>1"), + assignments=[gillespy2.EventAssignment(variable=variable, expression="100")] + )), + "FunctionDefinition": lambda model, variable: + model.add_function_definition( + gillespy2.FunctionDefinition(name="fn", function="variable", args=["variable"])), + } + + # List of `sbml_features` that each solver does NOT support. + sbml_support_check = { + # list(sbml_features.keys()) basically means "supports no SBML features, + # so ensure that all of them raise an error." + NumPySSASolver: list(sbml_features.keys()), + TauLeapingSolver: list(sbml_features.keys()), + ODESolver: list(sbml_features.keys()), + TauHybridSolver: [], + + SSACSolver: list(sbml_features.keys()), + ODECSolver: list(sbml_features.keys()), + TauLeapingCSolver: list(sbml_features.keys()), + TauHybridCSolver: [ + "AssignmentRule", + "FunctionDefinition", + ], + } + model = Example() for sp in model.listOfSpecies.values(): sp.mode = 'discrete' @@ -52,9 +85,9 @@ class TestAllSolvers(unittest.TestCase): labeled_results = {} labeled_results_more_trajectories = {} - for solver in solvers: - labeled_results[solver] = model.run(solver=solver, number_of_trajectories=1,seed=1) - labeled_results_more_trajectories[solver] = model.run(solver=solver, number_of_trajectories=2) + # for solver in solvers: + # labeled_results[solver] = model.run(solver=solver, number_of_trajectories=1,seed=1) + # labeled_results_more_trajectories[solver] = model.run(solver=solver, number_of_trajectories=2) def test_instantiated(self): for solver in self.solvers: @@ -129,6 +162,24 @@ def test_basic_solver_import(self): results3 = model.run(solver=BasicTauHybridSolver) self.assertTrue(results3[0].solver_name == 'TauHybridSolver') + @unittest.expectedFailure + def test_sbml_feature_validation(self): + class TestModel(gillespy2.Model): + def __init__(self): + gillespy2.Model.__init__(self, name="TestModel") + self.add_species(gillespy2.Species(name="S", initial_value=0)) + self.timespan(np.linspace(0, 10, 11)) + + for solver in self.solvers: + with self.subTest(solver=solver.name): + for sbml_feature_name in self.sbml_support_check.get(solver): + model = TestModel() + with self.subTest(sbml_feature=sbml_feature_name): + add_sbml_feature = self.sbml_features.get(sbml_feature_name) + add_sbml_feature(model, "S") + with self.assertRaises(Exception): + model.run(solver=solver) + if __name__ == '__main__': unittest.main() From 4b2122fa82c6dc1ea67f7ddff905cb415ea82940 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Sun, 13 Feb 2022 10:59:02 -0500 Subject: [PATCH 44/54] Fix unintentionally commented-out code in solver tests --- test/test_all_solvers.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 37be0d2b2..7d0c2e6b0 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -85,9 +85,9 @@ class TestAllSolvers(unittest.TestCase): labeled_results = {} labeled_results_more_trajectories = {} - # for solver in solvers: - # labeled_results[solver] = model.run(solver=solver, number_of_trajectories=1,seed=1) - # labeled_results_more_trajectories[solver] = model.run(solver=solver, number_of_trajectories=2) + for solver in solvers: + labeled_results[solver] = model.run(solver=solver, number_of_trajectories=1,seed=1) + labeled_results_more_trajectories[solver] = model.run(solver=solver, number_of_trajectories=2) def test_instantiated(self): for solver in self.solvers: From f0ae0155ba8e4f77f4682b772bf104a959d1906f Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Sun, 13 Feb 2022 12:50:48 -0500 Subject: [PATCH 45/54] Update solver validation test structure - Uses "supported" feature list, instead of "unsupported" list - Added check for supported feature validation --- test/test_all_solvers.py | 43 ++++++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 7d0c2e6b0..04aa97bdb 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -60,22 +60,26 @@ class TestAllSolvers(unittest.TestCase): gillespy2.FunctionDefinition(name="fn", function="variable", args=["variable"])), } - # List of `sbml_features` that each solver does NOT support. - sbml_support_check = { - # list(sbml_features.keys()) basically means "supports no SBML features, - # so ensure that all of them raise an error." - NumPySSASolver: list(sbml_features.keys()), - TauLeapingSolver: list(sbml_features.keys()), - ODESolver: list(sbml_features.keys()), - TauHybridSolver: [], - - SSACSolver: list(sbml_features.keys()), - ODECSolver: list(sbml_features.keys()), - TauLeapingCSolver: list(sbml_features.keys()), - TauHybridCSolver: [ + # List of supported SBML features for each solver. + # When a feature is implemented for a particular solver, add the feature to its list. + solver_supported_sbml_features = { + NumPySSASolver: [], + TauLeapingSolver: [], + ODESolver: [], + TauHybridSolver: [ "AssignmentRule", + "RateRule", + "Event", "FunctionDefinition", ], + + SSACSolver: [], + ODECSolver: [], + TauLeapingCSolver: [], + TauHybridCSolver: [ + "RateRule", + "Event", + ], } model = Example() @@ -170,16 +174,25 @@ def __init__(self): self.add_species(gillespy2.Species(name="S", initial_value=0)) self.timespan(np.linspace(0, 10, 11)) + all_features = set(self.sbml_features.keys()) for solver in self.solvers: + unsupported_features = all_features.difference(self.solver_supported_sbml_features.get(solver)) with self.subTest(solver=solver.name): - for sbml_feature_name in self.sbml_support_check.get(solver): + for sbml_feature_name in unsupported_features: model = TestModel() - with self.subTest(sbml_feature=sbml_feature_name): + with self.subTest("Unsupported model features raise an error", sbml_feature=sbml_feature_name): add_sbml_feature = self.sbml_features.get(sbml_feature_name) add_sbml_feature(model, "S") with self.assertRaises(Exception): model.run(solver=solver) + for sbml_feature_name in self.solver_supported_sbml_features.get(solver): + model = TestModel() + with self.subTest("Supported model features validate successfully", sbml_feature=sbml_feature_name): + add_sbml_feature = self.sbml_features.get(sbml_feature_name) + add_sbml_feature(model, "S") + model.run(solver=solver) + if __name__ == '__main__': unittest.main() From eac55854d266fe17c81d5605f0e71bdaa39d0088 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Sun, 13 Feb 2022 16:20:16 -0500 Subject: [PATCH 46/54] Add validation/support methods to solver/model superclass --- gillespy2/core/gillespySolver.py | 18 ++++++++++++++++-- gillespy2/core/model.py | 21 ++++++++++++++++++++- gillespy2/solvers/cpp/c_solver.py | 11 ----------- test/test_all_solvers.py | 6 +++--- 4 files changed, 39 insertions(+), 17 deletions(-) diff --git a/gillespy2/core/gillespySolver.py b/gillespy2/core/gillespySolver.py index 5926d341f..6dc19d21f 100644 --- a/gillespy2/core/gillespySolver.py +++ b/gillespy2/core/gillespySolver.py @@ -16,7 +16,8 @@ along with this program. If not, see . """ -from .gillespyError import SimulationError +from .gillespyError import SimulationError, ModelError +from typing import Set, Type class GillesPySolver: @@ -79,4 +80,17 @@ def get_increment(self, increment): `increment` argument from this `solver.run()` call. """ ) - return increment \ No newline at end of file + return increment + + @classmethod + def get_supported_features(cls) -> "Set[Type]": + return set() + + @classmethod + def validate_sbml_features(cls, model): + unsupported_features = model.get_model_features() - cls.get_supported_features() + if unsupported_features: + unsupported_features = [feature.__name__ for feature in unsupported_features] + raise ModelError(f"Could not run Model, " + f"SBML Features not supported by {cls.name}: " + + ", ".join(unsupported_features)) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 640010e49..427e7280e 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -15,7 +15,7 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . """ - +import gillespy2 from gillespy2.core.jsonify import TranslationTable from gillespy2.core.reaction import * from gillespy2.core.raterule import RateRule @@ -27,6 +27,7 @@ from collections import OrderedDict from gillespy2.core.gillespyError import * from .gillespyError import SimulationError +from typing import Set, Type try: import lxml.etree as eTree @@ -928,6 +929,24 @@ def get_best_solver_algo(self, algorithm): raise ModelError("Invalid value for the argument 'algorithm' entered. " "Please enter 'SSA', 'ODE', 'Tau-leaping', or 'Tau-Hybrid'.") + def get_model_features(self) -> "Set[Type]": + """ + Determine what solver-specific model features are present on the model. + Used to validate that the model is compatible with the given solver. + + :returns: Set containing the classes of every solver-specific feature present on the model. + """ + features = set() + if len(self.listOfEvents): + features.add(gillespy2.Event) + if len(self.listOfRateRules): + features.add(gillespy2.RateRule) + if len(self.listOfAssignmentRules): + features.add(gillespy2.AssignmentRule) + if len(self.listOfFunctionDefinitions): + features.add(gillespy2.FunctionDefinition) + return features + def run(self, solver=None, timeout=0, t=None, increment=None, show_labels=True, cpp_support=False, algorithm=None, **solver_args): """ diff --git a/gillespy2/solvers/cpp/c_solver.py b/gillespy2/solvers/cpp/c_solver.py index 0a403de27..bd33cd507 100644 --- a/gillespy2/solvers/cpp/c_solver.py +++ b/gillespy2/solvers/cpp/c_solver.py @@ -365,17 +365,6 @@ def _validate_kwargs(self, **kwargs): for key, val in kwargs.items(): log.warning(f"Unsupported keyword argument for solver {self.name}: {key}") - def _validate_sbml_features(self, unsupported_features: "dict[str, str]"): - detected = [ ] - for feature_name, count in unsupported_features.items(): - if count: - detected.append(feature_name) - - if len(detected): - raise gillespyError.ModelError(f"Could not run Model, " - f"SBML Features not supported by {self.name}: " - + ", ".join(detected)) - def _validate_seed(self, seed: int): if seed is None: return None diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 04aa97bdb..8adef5eb9 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -183,15 +183,15 @@ def __init__(self): with self.subTest("Unsupported model features raise an error", sbml_feature=sbml_feature_name): add_sbml_feature = self.sbml_features.get(sbml_feature_name) add_sbml_feature(model, "S") - with self.assertRaises(Exception): - model.run(solver=solver) + with self.assertRaises(gillespy2.ModelError): + solver.run(model=model) for sbml_feature_name in self.solver_supported_sbml_features.get(solver): model = TestModel() with self.subTest("Supported model features validate successfully", sbml_feature=sbml_feature_name): add_sbml_feature = self.sbml_features.get(sbml_feature_name) add_sbml_feature(model, "S") - model.run(solver=solver) + solver.run(model=model) if __name__ == '__main__': From 0e992eeba9a2effb15719c25c4829afc51453d84 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Sun, 13 Feb 2022 16:21:54 -0500 Subject: [PATCH 47/54] Add validation support to all solvers, update unit test --- gillespy2/solvers/cpp/ode_c_solver.py | 7 +------ gillespy2/solvers/cpp/ssa_c_solver.py | 7 +------ gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 16 ++++++++++------ gillespy2/solvers/cpp/tau_leaping_c_solver.py | 7 +------ gillespy2/solvers/numpy/CLE_solver.py | 1 + gillespy2/solvers/numpy/ode_solver.py | 1 + gillespy2/solvers/numpy/ssa_solver.py | 1 + gillespy2/solvers/numpy/tau_hybrid_solver.py | 12 +++++++++++- gillespy2/solvers/numpy/tau_leaping_solver.py | 1 + test/test_all_solvers.py | 1 - 10 files changed, 28 insertions(+), 26 deletions(-) diff --git a/gillespy2/solvers/cpp/ode_c_solver.py b/gillespy2/solvers/cpp/ode_c_solver.py index 4bde7118b..65d9a4532 100644 --- a/gillespy2/solvers/cpp/ode_c_solver.py +++ b/gillespy2/solvers/cpp/ode_c_solver.py @@ -48,6 +48,7 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal ODECSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) @@ -56,12 +57,6 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int self._validate_resume(t, resume) self._validate_kwargs(**kwargs) - self._validate_sbml_features({ - "Rate Rules": len(self.model.listOfRateRules), - "Assignment Rules": len(self.model.listOfAssignmentRules), - "Events": len(self.model.listOfEvents), - "Function Definitions": len(self.model.listOfFunctionDefinitions) - }) if resume is not None: t = abs(t - int(resume["time"][-1])) diff --git a/gillespy2/solvers/cpp/ssa_c_solver.py b/gillespy2/solvers/cpp/ssa_c_solver.py index 7aa81d2e7..9d647ee52 100644 --- a/gillespy2/solvers/cpp/ssa_c_solver.py +++ b/gillespy2/solvers/cpp/ssa_c_solver.py @@ -49,6 +49,7 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal SSACSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) @@ -57,12 +58,6 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int self._validate_resume(t, resume) self._validate_kwargs(**kwargs) - self._validate_sbml_features({ - "Rate Rules": len(self.model.listOfRateRules), - "Assignment Rules": len(self.model.listOfAssignmentRules), - "Events": len(self.model.listOfEvents), - "Function Definitions": len(self.model.listOfFunctionDefinitions) - }) if resume is not None: t = abs(t - int(resume["time"][-1])) diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index 04d6372d4..54c4723df 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -1,7 +1,7 @@ import gillespy2 from gillespy2.solvers.cpp.c_decoder import IterativeSimDecoder from gillespy2.solvers.utilities import solverutils as cutils -from gillespy2.core import GillesPySolver, Model +from gillespy2.core import GillesPySolver, Model, Event, RateRule from gillespy2.core.gillespyError import * from typing import Union from gillespy2.core import Results @@ -141,6 +141,13 @@ def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel" sanitized_model.options["GPY_HYBRID_NUM_EVENT_ASSIGNMENTS"] = str(len(event_assignment_list)) return sanitized_model + @classmethod + def get_supported_features(cls): + return { + Event, + RateRule, + } + 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 @@ -156,7 +163,7 @@ def get_solver_settings(self): return ('model', 't', 'number_of_trajectories', 'timeout', 'increment', 'seed', 'debug', 'profile') def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int = 1, timeout: int = 0, - increment: int = None, seed: int = None, debug: bool = False, profile: bool = False, variables={}, + increment: int = None, seed: int = None, debug: bool = False, profile: bool = False, variables={}, resume=None, live_output: str = None, live_output_options: dict = {}, tau_step: int = .03, tau_tol=0.03, **kwargs): if self is None: @@ -168,6 +175,7 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal TauHybridCSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) @@ -176,10 +184,6 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int self._validate_variables_in_set(variables, self.species + self.parameters) self._validate_resume(t, resume) self._validate_kwargs(**kwargs) - self._validate_sbml_features({ - "Assignment Rules": len(self.model.listOfAssignmentRules), - "Function Definitions": len(self.model.listOfFunctionDefinitions) - }) if resume is not None: t = abs(t - int(resume["time"][-1])) diff --git a/gillespy2/solvers/cpp/tau_leaping_c_solver.py b/gillespy2/solvers/cpp/tau_leaping_c_solver.py index 36d1d34aa..42110c2bb 100644 --- a/gillespy2/solvers/cpp/tau_leaping_c_solver.py +++ b/gillespy2/solvers/cpp/tau_leaping_c_solver.py @@ -48,6 +48,7 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal TauLeapingCSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) @@ -56,12 +57,6 @@ def run(self=None, model: Model = None, t: int = 20, number_of_trajectories: int self._validate_variables_in_set(variables, self.species + self.parameters) self._validate_resume(t, resume) self._validate_kwargs(**kwargs) - self._validate_sbml_features({ - "Rate Rules": len(self.model.listOfRateRules), - "Assignment Rules": len(self.model.listOfAssignmentRules), - "Events": len(self.model.listOfEvents), - "Function Definitions": len(self.model.listOfFunctionDefinitions) - }) if resume is not None: t = abs(t - int(resume["time"][-1])) diff --git a/gillespy2/solvers/numpy/CLE_solver.py b/gillespy2/solvers/numpy/CLE_solver.py index 7e25e53e6..96c9876dd 100644 --- a/gillespy2/solvers/numpy/CLE_solver.py +++ b/gillespy2/solvers/numpy/CLE_solver.py @@ -157,6 +157,7 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal CLESolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) diff --git a/gillespy2/solvers/numpy/ode_solver.py b/gillespy2/solvers/numpy/ode_solver.py index 8b8419dda..40fb7d495 100644 --- a/gillespy2/solvers/numpy/ode_solver.py +++ b/gillespy2/solvers/numpy/ode_solver.py @@ -116,6 +116,7 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, show_l if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal OSESolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) diff --git a/gillespy2/solvers/numpy/ssa_solver.py b/gillespy2/solvers/numpy/ssa_solver.py index 61fd5f6e9..a25bb6e10 100644 --- a/gillespy2/solvers/numpy/ssa_solver.py +++ b/gillespy2/solvers/numpy/ssa_solver.py @@ -81,6 +81,7 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal NumPySSASolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 13d0599aa..2559f392c 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -24,7 +24,7 @@ import threading import gillespy2 from gillespy2.solvers.utilities import Tau -from gillespy2.core import GillesPySolver, log +from gillespy2.core import GillesPySolver, log, Event, RateRule, AssignmentRule, FunctionDefinition from gillespy2.core.gillespyError import * from gillespy2.core.results import Results @@ -764,6 +764,15 @@ def get_solver_settings(self): return ('model', 't', 'number_of_trajectories', 'increment', 'seed', 'debug', 'profile', 'tau_tol', 'event_sensitivity', 'integrator', 'integrator_options', 'timeout') + @classmethod + def get_supported_features(cls): + return { + Event, + RateRule, + AssignmentRule, + FunctionDefinition, + } + @classmethod def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=None, debug=False, profile=False, tau_tol=0.03, event_sensitivity=100, integrator='LSODA', @@ -827,6 +836,7 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal TauHybridSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) diff --git a/gillespy2/solvers/numpy/tau_leaping_solver.py b/gillespy2/solvers/numpy/tau_leaping_solver.py index cbc0d8c7b..29d993df9 100644 --- a/gillespy2/solvers/numpy/tau_leaping_solver.py +++ b/gillespy2/solvers/numpy/tau_leaping_solver.py @@ -144,6 +144,7 @@ def run(self, model=None, t=20, number_of_trajectories=1, increment=None, seed=N if model is not None and model.get_json_hash() != self.model.get_json_hash(): raise SimulationError("Model must equal TauLeapingSolver.model.") self.model.resolve_parameters() + self.validate_sbml_features(model=model) increment = self.get_increment(increment=increment) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index 8adef5eb9..f5fb73959 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -166,7 +166,6 @@ def test_basic_solver_import(self): results3 = model.run(solver=BasicTauHybridSolver) self.assertTrue(results3[0].solver_name == 'TauHybridSolver') - @unittest.expectedFailure def test_sbml_feature_validation(self): class TestModel(gillespy2.Model): def __init__(self): From ace83f1f38ebaca8c57f53b42173cbbbea0d0ea8 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Mon, 14 Feb 2022 16:03:24 -0500 Subject: [PATCH 48/54] Update unit test to check feature validation w/o running --- test/test_all_solvers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_all_solvers.py b/test/test_all_solvers.py index f5fb73959..47804bb6b 100644 --- a/test/test_all_solvers.py +++ b/test/test_all_solvers.py @@ -183,14 +183,14 @@ def __init__(self): add_sbml_feature = self.sbml_features.get(sbml_feature_name) add_sbml_feature(model, "S") with self.assertRaises(gillespy2.ModelError): - solver.run(model=model) + solver.validate_sbml_features(model=model) for sbml_feature_name in self.solver_supported_sbml_features.get(solver): model = TestModel() with self.subTest("Supported model features validate successfully", sbml_feature=sbml_feature_name): add_sbml_feature = self.sbml_features.get(sbml_feature_name) add_sbml_feature(model, "S") - solver.run(model=model) + solver.validate_sbml_features(model=model) if __name__ == '__main__': From a94666b71c38aa3caf56b6b35ffff1200c9948d2 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Wed, 16 Feb 2022 16:54:04 -0500 Subject: [PATCH 49/54] Add root-finding directionality check - Reaction root-finder should only handle negative-to-positive roots --- .../solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp index ad8d70e2b..62a062a2c 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/integrator.cpp @@ -156,7 +156,7 @@ IntegrationResults Integrator::integrate(double *t, std::set &event_roots, for (; root_id < num_rxn_roots; ++root_id) { - if (root_results[root_id] != 0) + if (root_results[root_id] < 0) { reaction_roots.insert(data.active_reaction_ids[root_id]); } From 03a1f667377ce9e7fcfdd09e19db641d6938b4d5 Mon Sep 17 00:00:00 2001 From: Joshua Cooper Date: Thu, 17 Feb 2022 17:28:51 -0500 Subject: [PATCH 50/54] Add support for math functions in `Event` expressions --- gillespy2/solvers/cpp/build/template_gen.py | 1 + .../cpp/c_base/tau_hybrid_cpp_solver/hybrid_template.cpp | 3 +++ 2 files changed, 4 insertions(+) diff --git a/gillespy2/solvers/cpp/build/template_gen.py b/gillespy2/solvers/cpp/build/template_gen.py index ad71c379d..b75b697f5 100644 --- a/gillespy2/solvers/cpp/build/template_gen.py +++ b/gillespy2/solvers/cpp/build/template_gen.py @@ -43,6 +43,7 @@ class SanitizedModel: # as well as functions in Python that have a different name in C++. function_map = { "abs": "abs", + "round": "round", } def __init__(self, model: Model, variable=False): diff --git a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/hybrid_template.cpp b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/hybrid_template.cpp index 691a4acb4..cbcc38344 100644 --- a/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/hybrid_template.cpp +++ b/gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/hybrid_template.cpp @@ -18,11 +18,14 @@ #include "hybrid_template.h" #include "template_params.h" +#include // , cannot be overridden, so we can't use it as a delimiter // Use a separate macro to represent a delimiter #define AND , +using namespace std; + namespace Gillespy { namespace TauHybrid From fdffec490a92041ca75f59c7f03a70cf7932076f Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Fri, 18 Feb 2022 10:14:49 -0500 Subject: [PATCH 51/54] Fixed cleanup_tempfiles method. --- gillespy2/core/cleanup.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gillespy2/core/cleanup.py b/gillespy2/core/cleanup.py index 5da96f534..83e2d79e0 100644 --- a/gillespy2/core/cleanup.py +++ b/gillespy2/core/cleanup.py @@ -21,10 +21,10 @@ import tempfile def cleanup_tempfiles(): - ''' - Cleanup all tempfiles in gillespy2 build. - ''' - tempdir = tempfile.gettempdir() - for file_obj in os.listdir(tempdir): - if file_obj.startswith("gillespy2_build"): - cleanup_build_files(build_dir=os.path.join(tempdir, file_obj)) + ''' + Cleanup all tempfiles in gillespy2 build. + ''' + tempdir = tempfile.gettempdir() + for file_obj in os.listdir(tempdir): + if file_obj.startswith("gillespy2_build"): + shutil.rmtree(os.path.join(tempdir, file_obj)) From 6298eaf369c2761ca5f651747f31f5dd15262579 Mon Sep 17 00:00:00 2001 From: seanebum Date: Fri, 18 Feb 2022 10:31:22 -0500 Subject: [PATCH 52/54] added unit test for hybrid solvers based on multifiring event notebook --- test/example_models.py | 47 +++++++++++++++++++++++++++++++++++++- test/test_hybrid_solver.py | 15 ++++++++++-- 2 files changed, 59 insertions(+), 3 deletions(-) diff --git a/test/example_models.py b/test/example_models.py index 1ca2349a2..1711f0dd3 100644 --- a/test/example_models.py +++ b/test/example_models.py @@ -98,7 +98,7 @@ def __init__(self, parameter_values=None): class Dimerization(Model): def __init__(self, parameter_values=None): - # First call the gillespy2.Model initializer. + # First call the Model initializer. Model.__init__(self, name="Dimerization") # Define parameters for the rates of creation and dissociation. @@ -525,6 +525,51 @@ def __init__(self, parameter_values=None): # Timespan self.timespan(np.arange(0, 20, 0.05)) +class MultiFiringEvent(Model): + """ + This is a simple example for mass-action degradation of species S. We will add two events + to demonstrate the usage of events. The first event will assign the value '0' to our species + once time passes 20, and the second event will be triggered once time crosses 30, assigning + a value of "100" to our species and changing the value of our degradation rate parameter + "k1" from .01 to .1, causing the species to decay more quickly. + """ + + def __init__(self, parameter_values=None): + + # Initialize the model. + Model.__init__(self, name="Example") + + # Species + S = Species(name='Sp', initial_value=100, mode='discrete') +# ev_time1 = Species(name='ev_time1', initial_value=20) +# ev_time2 = Species(name='ev_time2', initial_value=30) + self.add_species([S]) + + # Parameters + k1 = Parameter(name='k1', expression=0.01) + ev_time1 = Parameter(name='ev_time1', expression=20) + ev_time2 = Parameter(name='ev_time2', expression=30) + self.add_parameter([k1, ev_time1, ev_time2]) + + # Events + et = EventTrigger(expression='t>ev_time1') + ea = EventAssignment(variable=S, expression='0') + ea1 = EventAssignment(variable=ev_time1, expression='ev_time1+15') + e = Event(name='event1', trigger=et, assignments=[ea, ea1]) + + et2 = EventTrigger(expression='t>ev_time2', persistent=True) + ea2 = EventAssignment(variable=S, expression='100') + ea3 = EventAssignment(variable=k1, expression='.1') + ea4 = EventAssignment(variable=ev_time2, expression='ev_time2+15') + e2 = Event(name='event2', trigger=et2, assignments=[ea2, ea3, ea4]) + self.add_event([e, e2]) + + #Reactions + r = Reaction(name='R', reactants={S:1}, products={}, rate=k1) #Multiple reactions + self.add_reaction([r]) + + self.timespan(np.linspace(0, 60, 175)) + __all__ = ['Trichloroethylene', 'LacOperon', 'Schlogl', 'MichaelisMenten', 'ToggleSwitch', 'Example', 'Tyson2StateOscillator', 'Oregonator', 'VilarOscillator', 'Dimerization', 'RobustModel'] diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index 2248cd3b1..75c7c143c 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -15,12 +15,11 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . """ - import unittest import numpy as np import gillespy2 from gillespy2.core.gillespyError import * -from example_models import Example, ExampleNoTspan +from example_models import Example, ExampleNoTspan, MultiFiringEvent from gillespy2 import TauHybridSolver @@ -238,6 +237,18 @@ def test_continuous_state_values(self): self.assertGreater(result["S1"][-1], 0.0, "Reaction never fired; indicates that continuous species is being truncated") + def test_multi_firing_event(self): + model = MultiFiringEvent() + for solver in self.solvers: + with self.subTest(solver=solver.name): + res = model.run(solver=solver, seed=1) + self.assertNotEqual(res['Sp'][15], 0) + self.assertEqual(res['Sp'][25], 0) + self.assertNotEqual(res['Sp'][32], 0) + self.assertEqual(res['Sp'][40], 0) + self.assertNotEqual(res['Sp'][48], 0) + self.assertEqual(res['Sp'][55], 0) + self.assertEqual(res['Sp'][60],100) if __name__ == '__main__': unittest.main() From 2df511c441dfdb76f1e203304b57b92180c6b88e Mon Sep 17 00:00:00 2001 From: seanebum Date: Fri, 18 Feb 2022 12:07:14 -0500 Subject: [PATCH 53/54] realigned asserts to fall within subtests --- test/test_hybrid_solver.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index 75c7c143c..4308d3db1 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -242,13 +242,13 @@ def test_multi_firing_event(self): for solver in self.solvers: with self.subTest(solver=solver.name): res = model.run(solver=solver, seed=1) - self.assertNotEqual(res['Sp'][15], 0) - self.assertEqual(res['Sp'][25], 0) - self.assertNotEqual(res['Sp'][32], 0) - self.assertEqual(res['Sp'][40], 0) - self.assertNotEqual(res['Sp'][48], 0) - self.assertEqual(res['Sp'][55], 0) - self.assertEqual(res['Sp'][60],100) + self.assertNotEqual(res['Sp'][15], 0) + self.assertEqual(res['Sp'][25], 0) + self.assertNotEqual(res['Sp'][32], 0) + self.assertEqual(res['Sp'][40], 0) + self.assertNotEqual(res['Sp'][48], 0) + self.assertEqual(res['Sp'][55], 0) + self.assertEqual(res['Sp'][60],100) if __name__ == '__main__': unittest.main() From 19985051b07108719c6cbb195b4c8ae32ec13aa8 Mon Sep 17 00:00:00 2001 From: seanebum Date: Fri, 18 Feb 2022 12:21:48 -0500 Subject: [PATCH 54/54] fixed time indexing of unit test to match event points --- test/example_models.py | 2 +- test/test_hybrid_solver.py | 13 ++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/test/example_models.py b/test/example_models.py index 1711f0dd3..bbcec04f5 100644 --- a/test/example_models.py +++ b/test/example_models.py @@ -568,7 +568,7 @@ def __init__(self, parameter_values=None): r = Reaction(name='R', reactants={S:1}, products={}, rate=k1) #Multiple reactions self.add_reaction([r]) - self.timespan(np.linspace(0, 60, 175)) + self.timespan(np.linspace(0, 60, 181)) __all__ = ['Trichloroethylene', 'LacOperon', 'Schlogl', 'MichaelisMenten', 'ToggleSwitch', 'Example', 'Tyson2StateOscillator', 'Oregonator', diff --git a/test/test_hybrid_solver.py b/test/test_hybrid_solver.py index 4308d3db1..c0d0c5fc2 100644 --- a/test/test_hybrid_solver.py +++ b/test/test_hybrid_solver.py @@ -242,13 +242,12 @@ def test_multi_firing_event(self): for solver in self.solvers: with self.subTest(solver=solver.name): res = model.run(solver=solver, seed=1) - self.assertNotEqual(res['Sp'][15], 0) - self.assertEqual(res['Sp'][25], 0) - self.assertNotEqual(res['Sp'][32], 0) - self.assertEqual(res['Sp'][40], 0) - self.assertNotEqual(res['Sp'][48], 0) - self.assertEqual(res['Sp'][55], 0) - self.assertEqual(res['Sp'][60],100) + self.assertNotEqual(res['Sp'][45], 0) + self.assertEqual(res['Sp'][75], 0) + self.assertNotEqual(res['Sp'][96], 0) + self.assertEqual(res['Sp'][120], 0) + self.assertNotEqual(res['Sp'][144], 0) + self.assertEqual(res['Sp'][165], 0) if __name__ == '__main__': unittest.main()