From 557856cd5ed7f6d75f5eee9617e04cdec4c8ac1d Mon Sep 17 00:00:00 2001 From: BryanRumsey <44621966+BryanRumsey@users.noreply.github.com> Date: Wed, 5 Oct 2022 16:42:09 -0400 Subject: [PATCH 01/47] Updated version to v1.7.2 --- gillespy2/__version__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/__version__.py b/gillespy2/__version__.py index 629322003..dfa98cb65 100644 --- a/gillespy2/__version__.py +++ b/gillespy2/__version__.py @@ -22,7 +22,7 @@ # ============================================================================= -__version__ = '1.7.1' +__version__ = '1.7.2' __title__ = 'GillesPy2' __description__ = 'Python interface for Gillespie-style biochemical simulations' From 3573e3ac14e935eddf8052475b8d176711bba352 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 6 Oct 2022 11:29:14 -0400 Subject: [PATCH 02/47] Refactored the reaction constructed, add_reactant, and add_product functions to handle duplicate stoich species. --- gillespy2/core/reaction.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/gillespy2/core/reaction.py b/gillespy2/core/reaction.py index fb76e21e0..c0c733d47 100644 --- a/gillespy2/core/reaction.py +++ b/gillespy2/core/reaction.py @@ -106,18 +106,20 @@ def __init__(self, name=None, reactants=None, products=None, propensity_function if reactants is not None: for r in reactants: rtype = type(r).__name__ - if rtype == 'Species': - self.reactants[r.name] = reactants[r] + name = r.name if rtype == 'Species' else r + if name in self.reactants: + self.reactants[name] += reactants[r] else: - self.reactants[r] = reactants[r] + self.reactants[name] = reactants[r] if products is not None: for p in products: - rtype = type(p).__name__ - if rtype == 'Species': - self.products[p.name] = products[p] + ptype = type(p).__name__ + name = p.name if ptype == 'Species' else p + if name in self.products: + self.products[name] += products[p] else: - self.products[p] = products[p] + self.products[name] = products[p] if self.marate is not None: rtype = type(self.marate).__name__ @@ -365,7 +367,10 @@ def add_product(self, species, stoichiometry): except TypeError as err: raise ReactionError(f"Failed to validate product. Reason given: {err}") from err - self.products[name] = stoichiometry + if name in self.products: + self.products[name] += stoichiometry + else: + self.products[name] = stoichiometry def addReactant(self, *args, **kwargs): """ @@ -404,7 +409,10 @@ def add_reactant(self, species, stoichiometry): except TypeError as err: raise ReactionError(f"Failed to validate reactant. Reason given: {err}") from err - self.reactants[name] = stoichiometry + if name in self.reactants: + self.reactants[name] += stoichiometry + else: + self.reactants[name] = stoichiometry if self.massaction and self.type == "mass-action": self._create_mass_action() From 6f345a9b657dc97896c7b102a89760f2cc7d8dd9 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 11 Oct 2022 09:55:55 -0400 Subject: [PATCH 03/47] Removed unused imports. --- gillespy2/core/gillespySolver.py | 1 - gillespy2/core/model.py | 1 - 2 files changed, 2 deletions(-) diff --git a/gillespy2/core/gillespySolver.py b/gillespy2/core/gillespySolver.py index 79b9ddef1..183c29fa1 100644 --- a/gillespy2/core/gillespySolver.py +++ b/gillespy2/core/gillespySolver.py @@ -14,7 +14,6 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . import copy -import numpy from .timespan import TimeSpan from .gillespyError import SimulationError, ModelError diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 498cb6c65..c987f5f13 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -14,7 +14,6 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . import platform -import numpy as np from typing import Set, Type from collections import OrderedDict From 4fb3e9bbc98231fba8005e3dd23a9f59fa2552b5 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 11 Oct 2022 09:56:37 -0400 Subject: [PATCH 04/47] Added numpy to the requirements file. --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index 105027372..bb2793fde 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,3 +2,4 @@ matplotlib>=3.0.3 plotly>=4.1.0 scipy>=1.3.1 setuptools>=41.0.1 +numpy>=1.18.5 From e45672949d7c4ba7993c8bc83d01a5cc87700512 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 11 Oct 2022 10:09:06 -0400 Subject: [PATCH 05/47] Removed the try catch block from the numpy module init. --- gillespy2/solvers/numpy/__init__.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/gillespy2/solvers/numpy/__init__.py b/gillespy2/solvers/numpy/__init__.py index aa395d8e9..069fd20c4 100644 --- a/gillespy2/solvers/numpy/__init__.py +++ b/gillespy2/solvers/numpy/__init__.py @@ -22,19 +22,12 @@ log.setLevel(logging.WARN) log.addHandler(_handler) -try: - import numpy as np - can_use_numpy = True - from .ssa_solver import NumPySSASolver - from .ode_solver import ODESolver - from .tau_leaping_solver import TauLeapingSolver - from .tau_hybrid_solver import TauHybridSolver - from .CLE_solver import CLESolver - log.debug("Successful Import of NumPy solvers.") - -except Exception as e: - log.warn(" Unable to use NumPy: {0}. The performance of this package can be significantly increased if you install NumPy.".format(e)) - can_use_numpy = False - +import numpy as np +from .ssa_solver import NumPySSASolver +from .ode_solver import ODESolver +from .tau_leaping_solver import TauLeapingSolver +from .tau_hybrid_solver import TauHybridSolver +from .CLE_solver import CLESolver +log.debug("Successful Import of NumPy solvers.") __all__ = ['NumPySSASolver', 'ODESolver', 'TauLeapingSolver', 'TauHybridSolver', 'CLESolver'] if can_use_numpy else [] From 75edbe46f2082a758abf463270f3e4ebe051195f Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 11 Oct 2022 10:55:53 -0400 Subject: [PATCH 06/47] Removed uses of 'can_use_numpy' and refactored the get_best_solver and get_best_solver_algo. --- gillespy2/core/model.py | 76 +++++++++++++---------------- gillespy2/solvers/numpy/__init__.py | 2 +- 2 files changed, 34 insertions(+), 44 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index c987f5f13..a86ce1ac6 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -1110,15 +1110,10 @@ def get_best_solver(self): :returns: gillespy2.gillespySolver """ - from gillespy2.solvers.numpy import can_use_numpy hybrid_check = False - 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(): tempMode = self.get_species(i).mode @@ -1126,84 +1121,79 @@ def get_best_solver(self): hybrid_check = True break + chybrid_check = hybrid_check + if len(self.get_all_assignment_rules()) or len(self.get_all_function_definitions()): + hybrid_check = True + chybrid_check = False + from gillespy2.solvers.cpp.build.build_engine import BuildEngine missing_deps = BuildEngine.get_missing_dependencies() windows_space = platform.system() == "Windows" and " " in gillespy2.__file__ can_use_cpp = len(missing_deps) == 0 and not windows_space - if not can_use_cpp and not can_use_numpy: - raise ModelError('Dependency Error, cannot run model.') + if hybrid_check: + if can_use_cpp and chybrid_check: + from gillespy2 import TauHybridCSolver + return TauHybridCSolver - 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 - if can_use_cpp is False and can_use_numpy and not hybrid_check: - from gillespy2 import NumPySSASolver - return NumPySSASolver - - else: + if can_use_cpp: from gillespy2 import SSACSolver return SSACSolver + from gillespy2 import NumPySSASolver + return NumPySSASolver + def get_best_solver_algo(self, algorithm): """ If user has specified a particular algorithm, we return either the Python or C++ version of that algorithm """ - from gillespy2.solvers.numpy import can_use_numpy from gillespy2.solvers.cpp.build.build_engine import BuildEngine missing_deps = BuildEngine.get_missing_dependencies() windows_space = platform.system() == "Windows" and " " in gillespy2.__file__ can_use_cpp = len(missing_deps) == 0 and not windows_space - 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.") + chybrid_check = not (len(self.get_all_assignment_rules()) or len(self.get_all_function_definitions())) if algorithm == 'Tau-Leaping': if can_use_cpp: from gillespy2 import TauLeapingCSolver return TauLeapingCSolver - else: - from gillespy2 import TauLeapingSolver - return TauLeapingSolver - elif algorithm == 'SSA': + from gillespy2 import TauLeapingSolver + return TauLeapingSolver + + if algorithm == 'SSA': if can_use_cpp: from gillespy2 import SSACSolver return SSACSolver - else: - from gillespy2 import NumPySSASolver - return NumPySSASolver - elif algorithm == 'ODE': + from gillespy2 import NumPySSASolver + return NumPySSASolver + + if algorithm == 'ODE': if can_use_cpp: from gillespy2 import ODECSolver return ODECSolver - else: - from gillespy2 import ODESolver - return ODESolver - elif algorithm == 'Tau-Hybrid': + from gillespy2 import ODESolver + return ODESolver + + if algorithm == 'Tau-Hybrid': if can_use_cpp and chybrid_check: from gillespy2 import TauHybridCSolver return TauHybridCSolver - else: - from gillespy2 import TauHybridSolver - return TauHybridSolver - elif algorithm == 'CLE': + from gillespy2 import TauHybridSolver + return TauHybridSolver + + if algorithm == 'CLE': from gillespy2 import CLESolver return CLESolver - else: - raise ModelError("Invalid value for the argument 'algorithm' entered. " - "Please enter 'SSA', 'ODE', 'CLE', 'Tau-Leaping', or 'Tau-Hybrid'.") + raise ModelError("Invalid value for the argument 'algorithm' entered. " + "Please enter 'SSA', 'ODE', 'CLE', 'Tau-Leaping', or 'Tau-Hybrid'.") def get_model_features(self) -> "Set[Type]": """ diff --git a/gillespy2/solvers/numpy/__init__.py b/gillespy2/solvers/numpy/__init__.py index 069fd20c4..b409e2c11 100644 --- a/gillespy2/solvers/numpy/__init__.py +++ b/gillespy2/solvers/numpy/__init__.py @@ -30,4 +30,4 @@ from .CLE_solver import CLESolver log.debug("Successful Import of NumPy solvers.") -__all__ = ['NumPySSASolver', 'ODESolver', 'TauLeapingSolver', 'TauHybridSolver', 'CLESolver'] if can_use_numpy else [] +__all__ = ['NumPySSASolver', 'ODESolver', 'TauLeapingSolver', 'TauHybridSolver', 'CLESolver'] From 5315bfe4df567c0e731db7052d9acdf4fb08ce99 Mon Sep 17 00:00:00 2001 From: Matthew Dippel Date: Wed, 12 Oct 2022 14:33:10 -0400 Subject: [PATCH 07/47] didn't realize __ changes __dict__ --- gillespy2/core/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 498cb6c65..0d3ee4746 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -1064,7 +1064,7 @@ def make_translation_table(self): events = self.listOfEvents.values() for event in events: for assignment in event.__dict__['assignments']: - del assignment.__dict__['__name_deprecated'] + del assignment.__dict__['_EventAssignment__name_deprecated'] functions = self.listOfFunctionDefinitions.values() # A translation table is used to anonymize user-defined variable names and formulas into generic counterparts. From e402cdb9379013a8c735ccd80cbb1acb7e88de71 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 25 Oct 2022 09:23:11 -0400 Subject: [PATCH 08/47] Fixed event assignment instantiation when building sbml models. --- gillespy2/sbml/SBMLimport.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/sbml/SBMLimport.py b/gillespy2/sbml/SBMLimport.py index e59c71e3b..43631e708 100644 --- a/gillespy2/sbml/SBMLimport.py +++ b/gillespy2/sbml/SBMLimport.py @@ -360,8 +360,8 @@ def __get_events(sbml_model, gillespy_model): gillespy_model.delete_parameter(a.getVariable()) gillespy_model.add_species([gillespy_species]) - gillespy_assignment = gillespy2.EventAssignment(a.getVariable(), - __get_math(a.getMath())) + gillespy_assignment = gillespy2.EventAssignment(variable=a.getVariable(), + expression=__get_math(a.getMath())) gillespy_assignments.append(gillespy_assignment) gillespy_event = gillespy2.Event( name=event.getId(), trigger=gillespy_trigger, From b24e26acd666733e36f0dc45210b43465ae963a4 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Fri, 18 Nov 2022 15:35:22 -0500 Subject: [PATCH 09/47] Added make rules to support update docs yaml. --- docs/Makefile | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/Makefile b/docs/Makefile index e73a85a33..d061a67f8 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -84,6 +84,17 @@ autobuild: $(FORMAT) clean: rm -rf "$(BUILDDIR)" +build-docs: + git pull + make clean + git stash + git checkout main + make html + +commit-docs: + -git add -f "$(BUILDDIR)" + -git commit -m "Latest docs build." "$(BUILDDIR)" + publish: make clean git stash From 2516dbfd5bca580cdbea1f25a2721744f04ac921 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Fri, 18 Nov 2022 15:36:21 -0500 Subject: [PATCH 10/47] Added a github action to automatically update the docs. --- .github/workflows/update-docs.yaml | 46 ++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 .github/workflows/update-docs.yaml diff --git a/.github/workflows/update-docs.yaml b/.github/workflows/update-docs.yaml new file mode 100644 index 000000000..ec16de74b --- /dev/null +++ b/.github/workflows/update-docs.yaml @@ -0,0 +1,46 @@ +name: Publish GillesPy2 Docs + +on: + push: + branches: [main] + +jobs: + publish: + runs-on: ubuntu-latest + + steps: + - name: Initialize environment + uses: actions/checkout@v2 + with: + ref: main + fetch-depth: 0 + + - name: Install Python + uses: actions/setup-python@v2 + with: + python-version: '3.7' + + - name: Install Sphinx Dependency + run: | + python3 -m pip install --upgrade pip + python3 -m pip install -U sphinx + python3 -m pip install numpy + python3 -m pip install -r requirements.txt + + - name: Update the Docs + working-directory: docs + run: | + make build-docs + + - name: Commit Changes + working-directory: docs + run: | + git config --local user.email "41898282+github-actions[bot]@users.noreply.github.com" + git config --local user.name "github-actions[bot]" + make commit-docs + + - name: Push Changes + uses: ad-m/github-push-action@master + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + branch: main From 67ea771df8043acd11977ca6da83f107cab1a76c Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 1 Dec 2022 15:01:18 -0800 Subject: [PATCH 11/47] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a8ce93853..c138daf48 100644 --- a/README.md +++ b/README.md @@ -87,7 +87,7 @@ In GillesPy2, a model is expressed as an object. Components, such as the reacti ```python def create_dimerization(parameter_values=None): # First call the gillespy2.Model initializer. - model = gillespy2.model(name='Dimerization') + model = gillespy2.Model(name='Dimerization') # Define parameters for the rates of creation and dissociation. k_c = gillespy2.Parameter(name='k_c', expression=0.005) From 4f7ed4957e74688a5a2fbe3d845dca20c343007f Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 5 Dec 2022 15:20:13 -0800 Subject: [PATCH 12/47] Fixes #874 for TauHybridSolver --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 5a4722514..2ca7439fe 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -87,6 +87,12 @@ def __init__(self, model=None, profile_reactions=False): self.profile_data['time'] = [] for k in list(self.model.listOfSpecies)+list(self.model.listOfReactions): self.profile_data[k] = [] + self.non_negative_species = set() + for reaction, _ in model.listOfReactions.items(): + for key, value in model.listOfReactions[reaction].reactants.items(): + self.non_negative_species.add(key.name) + for key, value in model.listOfReactions[reaction].products.items(): + self.non_negative_species.add(key.name) def __toggle_reactions(self, all_compiled, deterministic_reactions, dependencies, curr_state, det_spec, rr_sets): @@ -593,6 +599,12 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, return sol, curr_time + def __simulate_negative_state_check(self, curr_state): + # check each species to see if they are negative + for s in self.non_negative_species: + if curr_state[s] < 0: + raise SimulationError(f"Negative State detected at begining of step.") + def __simulate_invalid_state_check(self, species_modified, curr_state, compiled_reactions): invalid_state = False err_message="" @@ -627,6 +639,9 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, data. """ + # first check if we have a valid state: + self.__simulate_negative_state_check(curr_state) + event_queue = [] prev_y0 = copy.deepcopy(y0) prev_curr_state = copy.deepcopy(curr_state) @@ -652,7 +667,6 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, if curr_state[r] >= 0 and propensities[r] == 0: curr_state[r] = math.log(random.uniform(0, 1)) - sol, curr_time = self.__integrate(integrator_options, curr_state, y0, curr_time, propensities, y_map, compiled_reactions, @@ -709,7 +723,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, if min_tau is None or min_tau > rxn_times[rname]: min_tau = rxn_times[rname] rxn_selected = rname - if rxn_selected is None: raise Exception(f"Negative State detected in step, and no reaction found to fire.\n\n error_message={invalid_err_message}\n curr_time={curr_time}\n tau_step={tau_step}\n curr_state={curr_state}\n\nstarting_curr_state={starting_curr_state}\n\n starting_tau_step={starting_tau_step}\nspecies_modified={species_modified}\nrxn_count={rxn_count}\n propensities={propensities}\nrxn_times={rxn_times}\ncompiled_reactions={compiled_reactions}\ncurr_state_after={curr_state_after}\n propensities_after={propensities_after}\nstarting_propensities={starting_propensities}\nfloored_curr_state={floored_curr_state}\nfloored_propensities={floored_propensities}\n ") + if rxn_selected is None: raise SimulationError(f"Negative State detected in step, and no reaction found to fire.\n\n error_message={invalid_err_message}\n curr_time={curr_time}\n tau_step={tau_step}\n curr_state={curr_state}\n\nstarting_curr_state={starting_curr_state}\n\n starting_tau_step={starting_tau_step}\nspecies_modified={species_modified}\nrxn_count={rxn_count}\n propensities={propensities}\nrxn_times={rxn_times}\ncompiled_reactions={compiled_reactions}\ncurr_state_after={curr_state_after}\n propensities_after={propensities_after}\nstarting_propensities={starting_propensities}\nfloored_curr_state={floored_curr_state}\nfloored_propensities={floored_propensities}\n ") tau_step = min_tau #estimated time to the first stochatic reaction From 6dbfa6a8de3c7f9f70b819474aff633160919213 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 6 Dec 2022 12:49:11 -0800 Subject: [PATCH 13/47] Fixes #874 for TauHybridCSolver --- .../tau_hybrid_cpp_solver/HybridModel.h | 1 + .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 30 +++++++++++++++++++ gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 4 ++- 3 files changed, 34 insertions(+), 1 deletion(-) 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 c3de19763..4c2b33845 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 @@ -300,6 +300,7 @@ namespace Gillespy INTEGRATOR_FAILED = 3, INVALID_AFTER_SSA = 4, NEGATIVE_STATE_NO_SSA_REACTION = 5, + NEGATIVE_STATE_AT_BEGINING_OF_STEP = 6, }; }; 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 735e66e92..148b0063f 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 @@ -133,6 +133,16 @@ namespace Gillespy + bool IsStateValidNonNegativeSpecies(int num_species, std::vector current_state, std::set> non_negative_species){ + // Explicitly check for invalid population state, now that changes have been tallied. + // Note: this should only check species that are reactants or products + for (const auto &r : non_negative_species) { + if (current_state[r.id] < 0) { + return false; + } + } + return true; + } bool IsStateNegativeCheck(int num_species, int*population_changes, std::vector current_state, std::set>tau_args_reactants){ // Explicitly check for invalid population state, now that changes have been tallied. @@ -169,6 +179,19 @@ namespace Gillespy bool in_event_handling = false; unsigned int neg_state_loop_cnt = 0; + std::set> non_negative_species; + + for (int r = 0; r < model.number_reactions; r++) { + for (int spec = 0; spec < model.number_species; spec++) { + if (model.reactions[r].products_change[spec] > 0) { + } else if (model.reactions[r].reactants_change[spec] > 0) { + non_negative_species.insert(model.species[spec]); + } + } + } + + + generator = std::mt19937_64(simulation->random_seed); URNGenerator urn(simulation->random_seed); Integrator sol(simulation, model, urn, config.rel_tol, config.abs_tol); @@ -281,6 +304,13 @@ namespace Gillespy break; } + // check if the state is valid + if( ! IsStateValidNonNegativeSpecies(num_species, non_negative_species) ) { + // throw an error + simulation->set_status(HybridSimulation::NEGATIVE_STATE_AT_BEGINING_OF_STEP); + return; + } + // Expected tau step is determined. tau_step = select( model, diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index e2f794e9d..a1401ecd0 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -28,6 +28,7 @@ class ErrorStatus(IntEnum): INTEGRATOR_FAILED = 3 INVALID_AFTER_SSA = 4 NEGATIVE_STATE_NO_SSA_REACTION = 5 + NEGATIVE_STATE_AT_BEGINING_OF_STEP = 6 @classmethod def __create_options(cls, sanitized_model: "SanitizedModel") -> "SanitizedModel": @@ -187,7 +188,8 @@ def _handle_return_code(self, return_code: "int") -> "int": raise ExecutionError("Invalid state after single SSA step") if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_NO_SSA_REACTION: raise ExecutionError("Negative State detected in step, and no reaction found to fire.") - + if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_AT_BEGINING_OF_STEP: + raise ExecutionError("Negative State detected at beginning of step.") return super()._handle_return_code(return_code) @classmethod From ea4edd297cc6b537a08dca55731c96fa25cc0fff Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 6 Dec 2022 14:49:20 -0800 Subject: [PATCH 14/47] Fixing bugs --- .../tau_hybrid_cpp_solver/TauHybridSolver.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 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 148b0063f..6503dd867 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 @@ -133,11 +133,11 @@ namespace Gillespy - bool IsStateValidNonNegativeSpecies(int num_species, std::vector current_state, std::set> non_negative_species){ + bool IsStateValidNonNegativeSpecies(int num_species, std::vector current_state, std::vector non_negative_species){ // Explicitly check for invalid population state, now that changes have been tallied. // Note: this should only check species that are reactants or products for (const auto &r : non_negative_species) { - if (current_state[r.id] < 0) { + if (current_state[r] < 0) { return false; } } @@ -179,13 +179,14 @@ namespace Gillespy bool in_event_handling = false; unsigned int neg_state_loop_cnt = 0; - std::set> non_negative_species; + std::vector non_negative_species; - for (int r = 0; r < model.number_reactions; r++) { - for (int spec = 0; spec < model.number_species; spec++) { - if (model.reactions[r].products_change[spec] > 0) { - } else if (model.reactions[r].reactants_change[spec] > 0) { - non_negative_species.insert(model.species[spec]); + for (int spec = 0; spec < model.number_species; spec++) { + for (int r = 0; r < model.number_reactions; r++) { + if (model.reactions[r].products_change[spec] > 0 || + model.reactions[r].reactants_change[spec] > 0) { + non_negative_species.push_back(model.species[spec].id); + break;// once we flagged it, skip to the next species } } } @@ -305,7 +306,7 @@ namespace Gillespy } // check if the state is valid - if( ! IsStateValidNonNegativeSpecies(num_species, non_negative_species) ) { + if( ! IsStateValidNonNegativeSpecies(num_species, current_state, non_negative_species) ) { // throw an error simulation->set_status(HybridSimulation::NEGATIVE_STATE_AT_BEGINING_OF_STEP); return; From 7b3d5240966c18cd42ac9796cf911e75cccfe74e Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 10:34:24 -0500 Subject: [PATCH 15/47] Fixed incorrect error raised in constructor. --- gillespy2/core/functiondefinition.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/gillespy2/core/functiondefinition.py b/gillespy2/core/functiondefinition.py index 85b54876e..6d5d445af 100644 --- a/gillespy2/core/functiondefinition.py +++ b/gillespy2/core/functiondefinition.py @@ -19,6 +19,8 @@ from gillespy2.core.sortableobject import SortableObject from gillespy2.core.jsonify import Jsonify +from gillespy2.core.gillespyError import FunctionDefinitionError + class FunctionDefinition(SortableObject, Jsonify): """ Object representation defining an evaluable function to be used during @@ -36,7 +38,7 @@ class FunctionDefinition(SortableObject, Jsonify): def __init__(self, name="", function=None, args=[]): if function is None: - raise TypeError("Function string provided for FunctionDefinition cannot be None") + raise FunctionDefinitionError("Function string provided for FunctionDefinition cannot be None") if name in (None, ""): self.name = f'fd{uuid.uuid4()}'.replace('-', '_') @@ -65,4 +67,4 @@ def sanitized_function(self, species_mappings, parameter_mappings): sanitized_function = self.function_string for id, name in enumerate(names): sanitized_function = sanitized_function.replace(name, "{" + str(id) + "}") - return sanitized_function.format(*replacements) \ No newline at end of file + return sanitized_function.format(*replacements) From 683cd4a8f4d065e79634ef1e30fa6311584e350d Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 10:44:04 -0500 Subject: [PATCH 16/47] Fixed incorrect error in Events. closes #790. --- gillespy2/core/events.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/gillespy2/core/events.py b/gillespy2/core/events.py index a03e201ae..e0260bc3f 100644 --- a/gillespy2/core/events.py +++ b/gillespy2/core/events.py @@ -19,6 +19,7 @@ from gillespy2.core.gillespyError import * from gillespy2.core.jsonify import Jsonify +from gillespy2.core.gillespyError import EventError class EventAssignment(Jsonify): """ @@ -257,17 +258,12 @@ def add_assignment(self, assignment): :param assignment: The event or list of events to be added to this event. :type assignment: EventAssignment or list[EventAssignment] """ - - if hasattr(assignment, 'variable'): - self.assignments.append(assignment) - elif isinstance(assignment, list): + if isinstance(assignment, list): for assign in assignment: - if hasattr(assign, 'variable'): - self.assignments.append(assign) - else: - raise EventError('add_assignment failed to add EventAssignment. Assignment to be added must be of ' - 'type EventAssignment or list of EventAssignment objects.') + self.add_assignment(assign) + elif hasattr(assignment, 'variable'): + self.assignments.append(assignment) else: - raise ModelError("Unexpected parameter for add_assignment. Parameter must be EventAssignment or list of " - "EventAssignments") + raise EventError("Unexpected parameter for add_assignment. Assignment must be an EventAssignment or list of " + "EventAssignments objects") return assignment From e8f5c91e781696e860594dfe2a47296e97571f16 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 11:02:52 -0500 Subject: [PATCH 17/47] Added pylintrc file. --- .pylintrc | 571 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 571 insertions(+) create mode 100644 .pylintrc diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 000000000..27e483833 --- /dev/null +++ b/.pylintrc @@ -0,0 +1,571 @@ +[MASTER] + +# A comma-separated list of package or module names from where C extensions may +# be loaded. Extensions are loading into the active Python interpreter and may +# run arbitrary code. +extension-pkg-allow-list= + +# A comma-separated list of package or module names from where C extensions may +# be loaded. Extensions are loading into the active Python interpreter and may +# run arbitrary code. (This is an alternative name to extension-pkg-allow-list +# for backward compatibility.) +extension-pkg-whitelist= + +# Return non-zero exit code if any of these messages/categories are detected, +# even if score is above --fail-under value. Syntax same as enable. Messages +# specified are enabled, while categories only check already-enabled messages. +fail-on= + +# Specify a score threshold to be exceeded before program exits with error. +fail-under=10.0 + +# Files or directories to be skipped. They should be base names, not paths. +ignore=CVS + +# Add files or directories matching the regex patterns to the ignore-list. The +# regex matches against paths and can be in Posix or Windows format. +ignore-paths= + +# Files or directories matching the regex patterns are skipped. The regex +# matches against base names, not paths. +ignore-patterns= + +# Python code to execute, usually for sys.path manipulation such as +# pygtk.require(). +#init-hook= + +# Use multiple processes to speed up Pylint. Specifying 0 will auto-detect the +# number of processors available to use. +jobs=1 + +# Control the amount of potential inferred values when inferring a single +# object. This can help the performance when dealing with large functions or +# complex, nested conditions. +limit-inference-results=100 + +# List of plugins (as comma separated values of python module names) to load, +# usually to register additional checkers. +load-plugins= + +# Pickle collected data for later comparisons. +persistent=yes + +# Minimum Python version to use for version dependent checks. Will default to +# the version used to run pylint. +py-version=3.10 + +# When enabled, pylint would attempt to guess common misconfiguration and emit +# user-friendly hints instead of false-positive error messages. +suggestion-mode=yes + +# Allow loading of arbitrary C extensions. Extensions are imported into the +# active Python interpreter and may run arbitrary code. +unsafe-load-any-extension=no + + +[MESSAGES CONTROL] + +# Only show warnings with the listed confidence levels. Leave empty to show +# all. Valid levels: HIGH, INFERENCE, INFERENCE_FAILURE, UNDEFINED. +confidence= + +# Disable the message, report, category or checker with the given id(s). You +# can either give multiple identifiers separated by comma (,) or put this +# option multiple times (only on the command line, not in the configuration +# file where it should appear only once). You can also use "--disable=all" to +# disable everything first and then reenable specific checks. For example, if +# you want to run only the similarities checker, you can use "--disable=all +# --enable=similarities". If you want to run only the classes checker, but have +# no Warning level messages displayed, use "--disable=all --enable=classes +# --disable=W". +disable=raw-checker-failed, + bad-inline-option, + locally-disabled, + file-ignored, + suppressed-message, + useless-suppression, + deprecated-pragma, + use-symbolic-message-instead + +# Enable the message, report, category or checker with the given id(s). You can +# either give multiple identifier separated by comma (,) or put this option +# multiple time (only on the command line, not in the configuration file where +# it should appear only once). See also the "--disable" option for examples. +enable=c-extension-no-member + + +[REPORTS] + +# Python expression which should return a score less than or equal to 10. You +# have access to the variables 'error', 'warning', 'refactor', and 'convention' +# which contain the number of messages in each category, as well as 'statement' +# which is the total number of statements analyzed. This score is used by the +# global evaluation report (RP0004). +evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / statement) * 10) + +# Template used to display messages. This is a python new-style format string +# used to format the message information. See doc for all details. +#msg-template= + +# Set the output format. Available formats are text, parseable, colorized, json +# and msvs (visual studio). You can also give a reporter class, e.g. +# mypackage.mymodule.MyReporterClass. +output-format=text + +# Tells whether to display a full report or only the messages. +reports=no + +# Activate the evaluation score. +score=yes + + +[REFACTORING] + +# Maximum number of nested blocks for function / method body +max-nested-blocks=5 + +# Complete name of functions that never returns. When checking for +# inconsistent-return-statements if a never returning function is called then +# it will be considered as an explicit return statement and no message will be +# printed. +never-returning-functions=sys.exit,argparse.parse_error + + +[LOGGING] + +# The type of string formatting that logging methods do. `old` means using % +# formatting, `new` is for `{}` formatting. +logging-format-style=old + +# Logging modules to check that the string format arguments are in logging +# function parameter format. +logging-modules=logging + + +[SPELLING] + +# Limits count of emitted suggestions for spelling mistakes. +max-spelling-suggestions=4 + +# Spelling dictionary name. Available dictionaries: none. To make it work, +# install the 'python-enchant' package. +spelling-dict= + +# List of comma separated words that should be considered directives if they +# appear and the beginning of a comment and should not be checked. +spelling-ignore-comment-directives=fmt: on,fmt: off,noqa:,noqa,nosec,isort:skip,mypy: + +# List of comma separated words that should not be checked. +spelling-ignore-words= + +# A path to a file that contains the private dictionary; one word per line. +spelling-private-dict-file= + +# Tells whether to store unknown words to the private dictionary (see the +# --spelling-private-dict-file option) instead of raising a message. +spelling-store-unknown-words=no + + +[MISCELLANEOUS] + +# List of note tags to take in consideration, separated by a comma. +notes=FIXME, + XXX, + TODO + +# Regular expression of note tags to take in consideration. +#notes-rgx= + + +[TYPECHECK] + +# List of decorators that produce context managers, such as +# contextlib.contextmanager. Add to this list to register other decorators that +# produce valid context managers. +contextmanager-decorators=contextlib.contextmanager + +# List of members which are set dynamically and missed by pylint inference +# system, and so shouldn't trigger E1101 when accessed. Python regular +# expressions are accepted. +generated-members= + +# Tells whether missing members accessed in mixin class should be ignored. A +# class is considered mixin if its name matches the mixin-class-rgx option. +ignore-mixin-members=yes + +# Tells whether to warn about missing members when the owner of the attribute +# is inferred to be None. +ignore-none=yes + +# This flag controls whether pylint should warn about no-member and similar +# checks whenever an opaque object is returned when inferring. The inference +# can return multiple potential results while evaluating a Python object, but +# some branches might not be evaluated, which results in partial inference. In +# that case, it might be useful to still emit no-member and other checks for +# the rest of the inferred objects. +ignore-on-opaque-inference=yes + +# List of class names for which member attributes should not be checked (useful +# for classes with dynamically set attributes). This supports the use of +# qualified names. +ignored-classes=optparse.Values,thread._local,_thread._local + +# List of module names for which member attributes should not be checked +# (useful for modules/projects where namespaces are manipulated during runtime +# and thus existing member attributes cannot be deduced by static analysis). It +# supports qualified module names, as well as Unix pattern matching. +ignored-modules= + +# Show a hint with possible names when a member name was not found. The aspect +# of finding the hint is based on edit distance. +missing-member-hint=yes + +# The minimum edit distance a name should have in order to be considered a +# similar match for a missing member name. +missing-member-hint-distance=1 + +# The total number of similar names that should be taken in consideration when +# showing a hint for a missing member. +missing-member-max-choices=1 + +# Regex pattern to define which classes are considered mixins ignore-mixin- +# members is set to 'yes' +mixin-class-rgx=.*[Mm]ixin + +# List of decorators that change the signature of a decorated function. +signature-mutators= + + +[VARIABLES] + +# List of additional names supposed to be defined in builtins. Remember that +# you should avoid defining new builtins when possible. +additional-builtins= + +# Tells whether unused global variables should be treated as a violation. +allow-global-unused-variables=yes + +# List of names allowed to shadow builtins +allowed-redefined-builtins= + +# List of strings which can identify a callback function by name. A callback +# name must start or end with one of those strings. +callbacks=cb_, + _cb + +# A regular expression matching the name of dummy variables (i.e. expected to +# not be used). +dummy-variables-rgx=_+$|(_[a-zA-Z0-9_]*[a-zA-Z0-9]+?$)|dummy|^ignored_|^unused_ + +# Argument names that match this expression will be ignored. Default to name +# with leading underscore. +ignored-argument-names=_.*|^ignored_|^unused_ + +# Tells whether we should check for unused import in __init__ files. +init-import=no + +# List of qualified module names which can have objects that can redefine +# builtins. +redefining-builtins-modules=six.moves,past.builtins,future.builtins,builtins,io + + +[FORMAT] + +# Expected format of line ending, e.g. empty (any line ending), LF or CRLF. +expected-line-ending-format= + +# Regexp for a line that is allowed to be longer than the limit. +ignore-long-lines=^\s*(# )??$ + +# Number of spaces of indent required inside a hanging or continued line. +indent-after-paren=4 + +# String used as indentation unit. This is usually " " (4 spaces) or "\t" (1 +# tab). +indent-string=' ' + +# Maximum number of characters on a single line. +max-line-length=120 + +# Maximum number of lines in a module. +max-module-lines=1000 + +# Allow the body of a class to be on the same line as the declaration if body +# contains single statement. +single-line-class-stmt=no + +# Allow the body of an if to be on the same line as the test if there is no +# else. +single-line-if-stmt=no + + +[SIMILARITIES] + +# Comments are removed from the similarity computation +ignore-comments=yes + +# Docstrings are removed from the similarity computation +ignore-docstrings=yes + +# Imports are removed from the similarity computation +ignore-imports=no + +# Signatures are removed from the similarity computation +ignore-signatures=no + +# Minimum lines number of a similarity. +min-similarity-lines=4 + + +[BASIC] + +# Naming style matching correct argument names. +argument-naming-style=snake_case + +# Regular expression matching correct argument names. Overrides argument- +# naming-style. +#argument-rgx= + +# Naming style matching correct attribute names. +attr-naming-style=snake_case + +# Regular expression matching correct attribute names. Overrides attr-naming- +# style. +#attr-rgx= + +# Bad variable names which should always be refused, separated by a comma. +bad-names=foo, + bar, + baz, + toto, + tutu, + tata + +# Bad variable names regexes, separated by a comma. If names match any regex, +# they will always be refused +bad-names-rgxs= + +# Naming style matching correct class attribute names. +class-attribute-naming-style=any + +# Regular expression matching correct class attribute names. Overrides class- +# attribute-naming-style. +#class-attribute-rgx= + +# Naming style matching correct class constant names. +class-const-naming-style=UPPER_CASE + +# Regular expression matching correct class constant names. Overrides class- +# const-naming-style. +#class-const-rgx= + +# Naming style matching correct class names. +class-naming-style=PascalCase + +# Regular expression matching correct class names. Overrides class-naming- +# style. +#class-rgx= + +# Naming style matching correct constant names. +const-naming-style=UPPER_CASE + +# Regular expression matching correct constant names. Overrides const-naming- +# style. +#const-rgx= + +# Minimum line length for functions/classes that require docstrings, shorter +# ones are exempt. +docstring-min-length=-1 + +# Naming style matching correct function names. +function-naming-style=snake_case + +# Regular expression matching correct function names. Overrides function- +# naming-style. +#function-rgx= + +# Good variable names which should always be accepted, separated by a comma. +good-names=i, + j, + k, + ex, + Run, + _ + +# Good variable names regexes, separated by a comma. If names match any regex, +# they will always be accepted +good-names-rgxs=listOfParameters, listOfSpecies, listOfReaction, listOfAssignmentRules, + listOfRateRules, listOfEvents, listOfFunctionDefinitions, + _listOfParameters, _listOfSpecies, _listOfReaction, _listOfAssignmentRules, + _listOfRateRules, _listOfEvents, _listOfFunctionDefinitions + +# Include a hint for the correct naming format with invalid-name. +include-naming-hint=no + +# Naming style matching correct inline iteration names. +inlinevar-naming-style=any + +# Regular expression matching correct inline iteration names. Overrides +# inlinevar-naming-style. +#inlinevar-rgx= + +# Naming style matching correct method names. +method-naming-style=snake_case + +# Regular expression matching correct method names. Overrides method-naming- +# style. +#method-rgx= + +# Naming style matching correct module names. +module-naming-style=snake_case + +# Regular expression matching correct module names. Overrides module-naming- +# style. +#module-rgx= + +# Colon-delimited sets of names that determine each other's naming style when +# the name regexes allow several styles. +name-group= + +# Regular expression which should only match function or class names that do +# not require a docstring. +no-docstring-rgx=^_ + +# List of decorators that produce properties, such as abc.abstractproperty. Add +# to this list to register other decorators that produce valid properties. +# These decorators are taken in consideration only for invalid-name. +property-classes=abc.abstractproperty + +# Naming style matching correct variable names. +variable-naming-style=snake_case + +# Regular expression matching correct variable names. Overrides variable- +# naming-style. +#variable-rgx= + + +[STRING] + +# This flag controls whether inconsistent-quotes generates a warning when the +# character used as a quote delimiter is used inconsistently within a module. +check-quote-consistency=no + +# This flag controls whether the implicit-str-concat should generate a warning +# on implicit string concatenation in sequences defined over several lines. +check-str-concat-over-line-jumps=no + + +[IMPORTS] + +# List of modules that can be imported at any level, not just the top level +# one. +allow-any-import-level= + +# Allow wildcard imports from modules that define __all__. +allow-wildcard-with-all=no + +# Analyse import fallback blocks. This can be used to support both Python 2 and +# 3 compatible code, which means that the block might have code that exists +# only in one or another interpreter, leading to false positives when analysed. +analyse-fallback-blocks=no + +# Deprecated modules which should not be used, separated by a comma. +deprecated-modules= + +# Output a graph (.gv or any supported image format) of external dependencies +# to the given file (report RP0402 must not be disabled). +ext-import-graph= + +# Output a graph (.gv or any supported image format) of all (i.e. internal and +# external) dependencies to the given file (report RP0402 must not be +# disabled). +import-graph= + +# Output a graph (.gv or any supported image format) of internal dependencies +# to the given file (report RP0402 must not be disabled). +int-import-graph= + +# Force import order to recognize a module as part of the standard +# compatibility libraries. +known-standard-library= + +# Force import order to recognize a module as part of a third party library. +known-third-party=enchant + +# Couples of modules and preferred modules, separated by a comma. +preferred-modules= + + +[CLASSES] + +# Warn about protected attribute access inside special methods +check-protected-access-in-special-methods=no + +# List of method names used to declare (i.e. assign) instance attributes. +defining-attr-methods=__init__, + __new__, + setUp, + __post_init__ + +# List of member names, which should be excluded from the protected access +# warning. +exclude-protected=_asdict, + _fields, + _replace, + _source, + _make + +# List of valid names for the first argument in a class method. +valid-classmethod-first-arg=cls + +# List of valid names for the first argument in a metaclass class method. +valid-metaclass-classmethod-first-arg=cls + + +[DESIGN] + +# List of regular expressions of class ancestor names to ignore when counting +# public methods (see R0903) +exclude-too-few-public-methods= + +# List of qualified class names to ignore when counting class parents (see +# R0901) +ignored-parents= + +# Maximum number of arguments for function / method. +max-args=20 + +# Maximum number of attributes for a class (see R0902). +max-attributes=25 + +# Maximum number of boolean expressions in an if statement (see R0916). +max-bool-expr=5 + +# Maximum number of branch for function / method body. +max-branches=20 + +# Maximum number of locals for function / method body. +max-locals=25 + +# Maximum number of parents for a class (see R0901). +max-parents=7 + +# Maximum number of public methods for a class (see R0904). +max-public-methods=30 + +# Maximum number of return / yield for function / method body. +max-returns=6 + +# Maximum number of statements in function / method body. +max-statements=60 + +# Minimum number of public methods for a class (see R0903). +min-public-methods=2 + + +[EXCEPTIONS] + +# Exceptions that will emit a warning when being caught. Defaults to +# "BaseException, Exception". +overgeneral-exceptions=BaseException, + Exception From 3cd3e7efd255b79a0453b3586e4dd1009ec11bac Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 16:15:55 -0500 Subject: [PATCH 18/47] Pylint changes to core. New pylint score == 9.23 > --- gillespy2/core/__init__.py | 3 +- gillespy2/core/assignmentrule.py | 19 +- gillespy2/core/events.py | 51 ++-- gillespy2/core/functiondefinition.py | 16 +- gillespy2/core/gillespyError.py | 4 - gillespy2/core/gillespySolver.py | 69 +++-- gillespy2/core/jsonify.py | 10 +- gillespy2/core/liveGraphing.py | 91 ++++-- gillespy2/core/model.py | 418 ++++++++++++++------------- gillespy2/core/raterule.py | 16 +- gillespy2/core/reaction.py | 221 +++++++++----- gillespy2/core/results.py | 256 ++++++++++------ gillespy2/core/sortableobject.py | 2 +- gillespy2/core/species.py | 15 +- gillespy2/core/timespan.py | 9 +- 15 files changed, 731 insertions(+), 469 deletions(-) diff --git a/gillespy2/core/__init__.py b/gillespy2/core/__init__.py index d5e70b30c..8182306ac 100644 --- a/gillespy2/core/__init__.py +++ b/gillespy2/core/__init__.py @@ -15,6 +15,8 @@ # along with this program. If not, see . import logging + +from gillespy2.__version__ import __version__ from .assignmentrule import * from .cleanup import * from .events import * @@ -29,7 +31,6 @@ from .sortableobject import * from .species import * from .timespan import TimeSpan -from gillespy2.__version__ import __version__ version = __version__ diff --git a/gillespy2/core/assignmentrule.py b/gillespy2/core/assignmentrule.py index d12e77724..281d9c01c 100644 --- a/gillespy2/core/assignmentrule.py +++ b/gillespy2/core/assignmentrule.py @@ -19,7 +19,6 @@ from gillespy2.core.sortableobject import SortableObject from gillespy2.core.jsonify import Jsonify - class AssignmentRule(SortableObject, Jsonify): """ An AssignmentRule is used to express equations that set the values of @@ -48,11 +47,23 @@ def __str__(self): return f"{self.name}: Var: {var_name}: {self.formula}" def sanitized_formula(self, species_mappings, parameter_mappings): + ''' + Sanitize the assignment rule formula. + + :param species_mappings: Mapping of species names to sanitized species names. + :type species_mappings: dict + + :param parameter_mappings: Mapping of parameter names to sanitized parameter names. + :type parameter_mappings: dict + + :returns: The sanitized formula. + :rtype: str + ''' names = sorted(list(species_mappings.keys()) + list(parameter_mappings.keys()), key=lambda x: len(x), reverse=True) replacements = [parameter_mappings[name] if name in parameter_mappings else species_mappings[name] for name in names] sanitized_formula = self.formula - for id, name in enumerate(names): - sanitized_formula = sanitized_formula.replace(name, "{" + str(id) + "}") - return sanitized_formula.format(*replacements) \ No newline at end of file + for i, name in enumerate(names): + sanitized_formula = sanitized_formula.replace(name, "{" + str(i) + "}") + return sanitized_formula.format(*replacements) diff --git a/gillespy2/core/events.py b/gillespy2/core/events.py index e0260bc3f..b3fcacb03 100644 --- a/gillespy2/core/events.py +++ b/gillespy2/core/events.py @@ -16,7 +16,8 @@ import uuid -from gillespy2.core.gillespyError import * +from gillespy2.core.parameter import Parameter +from gillespy2.core.species import Species from gillespy2.core.jsonify import Jsonify from gillespy2.core.gillespyError import EventError @@ -39,14 +40,12 @@ class EventAssignment(Jsonify): :type expression: str """ - - def __init__(self, name=None, variable=None, expression=None): if name in (None, ""): name = f'evn{uuid.uuid4()}'.replace('-', '_') else: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning("EventAssignment.name has been deprecated.") self.__name_deprecated = name @@ -56,8 +55,6 @@ def __init__(self, name=None, variable=None, expression=None): if expression is not None: self.expression = str(expression) - from gillespy2.core.parameter import Parameter - from gillespy2.core.species import Species #TODO: ADD Compartment to valid variable types once implemented valid_variable_types = [Species, Parameter, str] @@ -76,11 +73,10 @@ def __str__(self): def __getattr__(self, key): if key == 'name': - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning('EventAssignment.name has been deprecated.') return self.__name_deprecated - else: - raise AttributeError + raise AttributeError class EventTrigger(Jsonify): """ @@ -121,16 +117,27 @@ def __str__(self): return self.expression def sanitized_expression(self, species_mappings, parameter_mappings): + ''' + Sanitize the event trigger expression. + + :param species_mappings: Mapping of species names to sanitized species names. + :type species_mappings: dict + + :param parameter_mappings: Mapping of parameter names to sanitized parameter names. + :type parameter_mappings: dict + + :returns: The sanitized expression. + :rtype: str + ''' names = sorted(list(species_mappings.keys()) + list(parameter_mappings.keys()), key=lambda x: len(x), reverse=True) replacements = [parameter_mappings[name] if name in parameter_mappings else species_mappings[name] for name in names] sanitized_expression = self.expression - for id, name in enumerate(names): - sanitized_expression = sanitized_expression.replace(name, "{"+str(id)+"}") + for i, name in enumerate(names): + sanitized_expression = sanitized_expression.replace(name, "{"+str(i)+"}") return sanitized_expression.format(*replacements) - class Event(Jsonify): """ An Event can be given as an assignment_expression (function) or directly @@ -201,14 +208,14 @@ def __init__(self, name="", delay=None, assignments=[], priority="0", trigger=No else: raise EventError( 'name must be a valid string') - + # Trigger if hasattr(trigger, 'expression'): self.trigger = trigger else: raise EventError( 'trigger must be set to a valid EventTrigger') - + # Delay if delay is None or isinstance(delay, str): self.delay = delay @@ -242,13 +249,13 @@ def __init__(self, name="", delay=None, assignments=[], priority="0", trigger=No def __str__(self): print_string = self.name print_string += '\n\tTrigger: ' + str(self.trigger) - if len(self.assignments): + if len(self.assignments) > 0: print_string += '\n\tAssignments:' - for a in self.assignments: - if isinstance(a.variable, str): - print_string += '\n\t\t' + a.variable + ': ' + a.expression + for assign in self.assignments: + if isinstance(assign.variable, str): + print_string += '\n\t\t' + assign.variable + ': ' + assign.expression else: - print_string += '\n\t\t' + a.variable.name + ': ' + a.expression + print_string += '\n\t\t' + assign.variable.name + ': ' + assign.expression return print_string def add_assignment(self, assignment): @@ -264,6 +271,8 @@ def add_assignment(self, assignment): elif hasattr(assignment, 'variable'): self.assignments.append(assignment) else: - raise EventError("Unexpected parameter for add_assignment. Assignment must be an EventAssignment or list of " - "EventAssignments objects") + raise EventError( + "Unexpected parameter for add_assignment. Assignment must be an EventAssignment or list of " + "EventAssignments objects" + ) return assignment diff --git a/gillespy2/core/functiondefinition.py b/gillespy2/core/functiondefinition.py index 6d5d445af..153cb13ff 100644 --- a/gillespy2/core/functiondefinition.py +++ b/gillespy2/core/functiondefinition.py @@ -60,11 +60,23 @@ def get_arg_string(self) -> str: return ','.join(self.args) def sanitized_function(self, species_mappings, parameter_mappings): + ''' + Sanitize the function definition function. + + :param species_mappings: Mapping of species names to sanitized species names. + :type species_mappings: dict + + :param parameter_mappings: Mapping of parameter names to sanitized parameter names. + :type parameter_mappings: dict + + :returns: The sanitized function. + :rtype: str + ''' names = sorted(list(species_mappings.keys()) + list(parameter_mappings.keys()), key=lambda x: len(x), reverse=True) replacements = [parameter_mappings[name] if name in parameter_mappings else species_mappings[name] for name in names] sanitized_function = self.function_string - for id, name in enumerate(names): - sanitized_function = sanitized_function.replace(name, "{" + str(id) + "}") + for i, name in enumerate(names): + sanitized_function = sanitized_function.replace(name, "{" + str(i) + "}") return sanitized_function.format(*replacements) diff --git a/gillespy2/core/gillespyError.py b/gillespy2/core/gillespyError.py index 463f601f1..91d04d3b3 100644 --- a/gillespy2/core/gillespyError.py +++ b/gillespy2/core/gillespyError.py @@ -88,10 +88,6 @@ class SimulationTimeoutError(SimulationError): pass -class EventError(ModelError): - pass - - # Results errors class ResultsError(Exception): pass diff --git a/gillespy2/core/gillespySolver.py b/gillespy2/core/gillespySolver.py index 183c29fa1..ab856bfb4 100644 --- a/gillespy2/core/gillespySolver.py +++ b/gillespy2/core/gillespySolver.py @@ -17,11 +17,9 @@ from .timespan import TimeSpan from .gillespyError import SimulationError, ModelError -from typing import Set, Type - class GillesPySolver: - """ + """ Abstract class for a solver. """ @@ -29,34 +27,36 @@ class GillesPySolver: def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None, debug=False, profile=False, **kwargs): - """ + """ Call out and run the solver. Collect the results. :param model: The model on which the solver will operate. :type model: gillespy.Model - + :param t: The end time of the solver :type t: int - + :param number_of_trajectories: The number of times to sample the chemical master equation. Each trajectory will be returned at the end of the simulation. :type number_of_trajectories: int :param increment: The time step of the solution :type increment: float - - :param seed: The random seed for the simulation. Defaults to None. + + :param seed: The random seed for the simulation. Defaults to None. :type seed: int - + :param debug: Set to True to provide additional debug information about the simulation. :type debug: bool - + :returns: Simulation trajectories. """ raise SimulationError("This abstract solver class cannot be used directly.") def get_solver_settings(self): - + ''' + Abstract function for getting a solvers settings. + ''' raise SimulationError("This abstract solver class cannot be used directly") def validate_tspan(self, increment, t): @@ -76,8 +76,8 @@ def validate_tspan(self, increment, t): """ Failed while preparing to run the model. Neither increment or timespan are set. - To continue either add a `timespan` definition to your Model or add the - `increment` and `t` arguments to this `solver.run()` call. + To continue either add a `timespan` definition to your Model or add the + `increment` and `t` arguments to this `solver.run()` call. """ ) @@ -86,8 +86,8 @@ def validate_tspan(self, increment, t): """ Failed while preparing to run the model. Both increment and timespan are set. - To continue either remove your `timespan` definition from your Model or remove the - `increment` argument from this `solver.run()` call. + To continue either remove your `timespan` definition from your Model or remove the + `increment` argument from this `solver.run()` call. """ ) @@ -100,14 +100,31 @@ def validate_tspan(self, increment, t): @classmethod def get_supported_features(cls) -> "Set[Type]": + ''' + Get the set of supported features. + ''' return set() @classmethod def get_supported_integrator_options(cls) -> "Set[str]": + ''' + Get the supported integrator options + ''' return set() @classmethod def validate_model(cls, sol_model, model): + ''' + Validate that the solvers model is the same as the model passed to solver.run. + + :param sol_model: Solver model object. + :type sol_model: gillespy2.Model + + :param model: Model object passed to solver.run. + :type model: gillespy2.Model + + :raises SimulationError: If the models are not equal. + ''' if model is not None: model.compile_prep() if model.tspan is None: @@ -118,6 +135,14 @@ def validate_model(cls, sol_model, model): @classmethod def validate_sbml_features(cls, model): + ''' + Validate the models SBML features. + + :param model: The model object to be checked. + :type model: gillespy2.Model + + :raises ModelError: If the model contains unsupported SBML features. + ''' unsupported_features = model.get_model_features() - cls.get_supported_features() if unsupported_features: unsupported_features = [feature.__name__ for feature in unsupported_features] @@ -127,5 +152,15 @@ def validate_sbml_features(cls, model): @classmethod def validate_integrator_options(cls, options: "dict[str, float]") -> "dict[str, float]": - return { option: value for option, value in options.items() if option in cls.get_supported_integrator_options() } - + ''' + Validate the integrator options. + + :param options: Integrator options to validate. + :type options: dict + + :returns: Dict of supported integrator options. + :rtype: dict + ''' + return { + option: value for option, value in options.items() if option in cls.get_supported_integrator_options() + } diff --git a/gillespy2/core/jsonify.py b/gillespy2/core/jsonify.py index da9461a04..b5a1a2021 100644 --- a/gillespy2/core/jsonify.py +++ b/gillespy2/core/jsonify.py @@ -145,7 +145,7 @@ def get_json_hash(self, ignore_whitespace=True, hash_private_vals=False) -> str: prior to being hashed. :type ignore_whitespace: bool - :param hash_private_vals: If set to True all private and non-private variables will + :param hash_private_vals: If set to True all private and non-private variables will be included in the hash. :type hash_private_vals: bool @@ -161,17 +161,17 @@ def get_json_hash(self, ignore_whitespace=True, hash_private_vals=False) -> str: return hashlib.md5(str.encode(model_json)).hexdigest() - def __eq__(self, o: "Jsonify"): + def __eq__(self, obj: "Jsonify"): """ Overload to compare the json of two objects that derive from Jsonify. This method will not do any additional translation. - :param o: The Jsonify object to compare against. - :type o: Jsonify + :param obj: The Jsonify object to compare against. + :type obj: Jsonify :returns: True if equal, False if not. """ - return self.get_json_hash() == o.get_json_hash() + return self.get_json_hash() == obj.get_json_hash() class ComplexJsonCoder(JSONEncoder): """ diff --git a/gillespy2/core/liveGraphing.py b/gillespy2/core/liveGraphing.py index 07b59e450..3622fae35 100644 --- a/gillespy2/core/liveGraphing.py +++ b/gillespy2/core/liveGraphing.py @@ -16,10 +16,14 @@ import sys import json -import time import threading -from gillespy2.core import log +from math import floor +from IPython.display import clear_output + +from gillespy2.core.results import common_rgb_values +from gillespy2.core.gillespyError import SimulationError +from gillespy2.core import log class CRepeatTimer(threading.Timer): """ @@ -30,12 +34,11 @@ class CRepeatTimer(threading.Timer): pause = False def run(self): - from IPython.display import clear_output - type = str.join('', [*self.args[1]]) + _ = str.join('', [*self.args[1]]) while not self.finished.wait(self.interval): args = self.args[0].get() self.function(*args, **self.kwargs) - + if not self.pause: args = self.args[0].get() self.kwargs['finished'] = True @@ -49,8 +52,7 @@ class RepeatTimer(threading.Timer): pause = False def run(self): - from IPython.display import clear_output - type = str.join('', [*self.args[3]]) + _ = str.join('', [*self.args[3]]) self.args = self.args[:3] while not self.finished.wait(self.interval): self.function(*self.args, **self.kwargs) @@ -62,28 +64,41 @@ def run(self): def display_types(): + ''' + Get the list of supported display types. + + :returns: Supported display types. + :rtype: list + ''' return ["graph", "text", "progress"] def valid_graph_params(live_output_options): + ''' + Validated the live output options. + + :param live_output_options: Options to be validated. + :type live_output_options: dict + + :raises SimulationError: If the display type is invalid. + ''' if live_output_options['type'] not in ['progress', 'graph', 'text']: - from gillespy2.core.gillespyError import SimulationError raise SimulationError("Invalid input to 'live_output', please check spelling and ensure input is" " lower case.") if 'interval' not in live_output_options: live_output_options['interval'] = 1 elif live_output_options['interval'] < 0: - - log.warning("In LiveGraphing live_output_options, got \"interval\" = \"{0}\". setting interval = 1" - .format(live_output_options['interval'])) + message = f"In LiveGraphing live_output_options, got 'interval' = '{live_output_options['interval']}'." + message += " setting interval = 1" + log.warning(message) live_output_options['interval'] = 1 if live_output_options['type'] == "graph" and live_output_options['interval'] < 1: - log.warning("In LiveGraphing live_output_options, got \"interval\" = \"{0}\". Consider using an interval >= 1 " - "when displaying graphs".format(live_output_options['interval'])) + message = f"In LiveGraphing live_output_options, got 'interval' = '{live_output_options['interval']}'." + message += "Consider using an interval >= 1 when displaying graphs." + log.warning(message) if 'clear_output' not in live_output_options: - if live_output_options['type'] == "graph" or live_output_options['type'] == "progress": live_output_options['clear_output'] = True else: @@ -119,14 +134,25 @@ def __init__(self, model=None, timeline=None, number_of_trajectories=1, live_out self.header_printed = False def trajectory_header(self): + ''' + Create the trajectory header for the output. + ''' return "Trajectory (" + str(self.current_trajectory) + "/" + str(self.number_of_trajectories) + ")" def increment_trajectory(self, trajectory_num): + ''' + Increment the trejectory counter. + ''' self.current_trajectory = trajectory_num + 1 self.header_printed = False def print_text_header(self, file_obj): + ''' + Print the header for text display type. + :param file_obj: File object to write text output. + :type file_obj: file object + ''' self.header_printed = True if self.number_of_trajectories > 1: print(self.trajectory_header(), file=file_obj) @@ -136,12 +162,22 @@ def print_text_header(self, file_obj): print(species[:10].ljust(10), end="|", file=file_obj) print("", file=file_obj) - ''' - curr_state and curr_time should be list of len 1 to get reference - ''' def display(self, curr_state, curr_time, trajectory_base, finished=False): - from IPython.display import clear_output - from math import floor + ''' + Display the output for the live grapher. + + :param curr_state: Current state of the simulation. Should be a list of len 1 to get reference. + :type curr_state: list + + :param curr_time: Current time of the simulation. Should be a list of len 1 to get reference. + :type curr_time: list + + :param trajectory_base: Current results of the simulation. + :type trajectory_base: list + + :param finished: Indicates whether or not the simulation has finished. + :type finished: bool + ''' curr_time = curr_time[0] curr_state = curr_state[0] @@ -159,7 +195,7 @@ def display(self, curr_state, curr_time, trajectory_base, finished=False): file_obj = sys.stdout else: mode = "w" if self.clear_output else "a" - file_obj = open(self.file_path, mode) + file_obj = open(self.file_path, mode, encoding="utf-8") if self.display_type == "text": @@ -181,15 +217,15 @@ def display(self, curr_state, curr_time, trajectory_base, finished=False): if self.resume is True: print(f"progress = {round(((curr_time-self.x_shift)/self.timeline_len)*100, 2)} %\n", file=file_obj) else: - print(f"progress = {round((curr_time / (self.timeline_len + self.x_shift)) * 100, 2) }%\n", file=file_obj) + print( + f"progress = {round((curr_time / (self.timeline_len + self.x_shift)) * 100, 2) }%\n", file=file_obj + ) elif self.display_type == "graph": - if finished: return - import matplotlib.pyplot as plt - from gillespy2.core.results import common_rgb_values + import matplotlib.pyplot as plt # pylint: disable=import-outside-toplevel entry_count = floor(curr_time) - self.x_shift @@ -210,12 +246,11 @@ def display(self, curr_state, curr_time, trajectory_base, finished=False): plt.show() elif self.display_type == "figure": - import plotly - import plotly.graph_objs as go - from gillespy2.core.results import common_rgb_values + import plotly # pylint: disable=import-outside-toplevel + import plotly.graph_objs as go # pylint: disable=import-outside-toplevel entry_count = floor(curr_time) - self.x_shift - + trace_list = [] for i, species in enumerate(self.species): line_dict = {"color": common_rgb_values()[(i) % len(common_rgb_values())]} diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 43e374ccb..0531b17a2 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -14,8 +14,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . import platform -from typing import Set, Type -from collections import OrderedDict +from collections import OrderedDict, ChainMap import gillespy2 from gillespy2.core.assignmentrule import AssignmentRule @@ -28,7 +27,6 @@ from gillespy2.core.timespan import TimeSpan from gillespy2.core.sortableobject import SortableObject from gillespy2.core.jsonify import Jsonify, TranslationTable -from gillespy2.core.results import Trajectory, Results from gillespy2.core.gillespyError import ( ModelError, SimulationError, @@ -41,7 +39,7 @@ no_pretty_print = False -except: +except ImportError: import xml.etree.ElementTree as eTree import xml.dom.minidom import re @@ -63,15 +61,15 @@ def import_SBML(filename, name=None, gillespy_model=None, report_silently_with_s :param gillespy_model: If desired, the SBML model may be added to an existing GillesPy model :type gillespy_model: gillespy.Model - :param report_silently_with_sbml_error: SBML import will fail silently and + :param report_silently_with_sbml_error: SBML import will fail silently and SBML errors will not output to the user. :type report_silently_with_sbml_error: bool """ try: - from gillespy2.sbml.SBMLimport import convert - except ImportError: - raise ImportError('SBML conversion not imported successfully') + from gillespy2.sbml.SBMLimport import convert # pylint: disable=import-outside-toplevel + except ImportError as err: + raise ImportError('SBML conversion not imported successfully') from err return convert( filename, model_name=name, gillespy_model=gillespy_model, @@ -90,9 +88,9 @@ def export_SBML(gillespy_model, filename=None): :type filename: str """ try: - from gillespy2.sbml.SBMLexport import export - except ImportError: - raise ImportError('SBML export conversion not imported successfully') + from gillespy2.sbml.SBMLexport import export # pylint: disable=import-outside-toplevel + except ImportError as err: + raise ImportError('SBML export conversion not imported successfully') from err return export(gillespy_model, path=filename) @@ -108,9 +106,9 @@ def export_StochSS(gillespy_model, filename=None, return_stochss_model=False): :type filename: str """ try: - from gillespy2.stochss.StochSSexport import export - except ImportError: - raise ImportError('StochSS export conversion not imported successfully') + from gillespy2.stochss.StochSSexport import export # pylint: disable=import-outside-toplevel + except ImportError as err: + raise ImportError('StochSS export conversion not imported successfully') from err return export(gillespy_model, path=filename, return_stochss_model=return_stochss_model) @@ -214,32 +212,32 @@ def decorate(header): print_string = self.name if len(self.listOfSpecies): print_string += decorate('Species') - for s in sorted(self.listOfSpecies.values()): - print_string += '\n' + str(s) + for species in sorted(self.listOfSpecies.values()): + print_string += '\n' + str(species) if len(self.listOfParameters): print_string += decorate('Parameters') - for p in sorted(self.listOfParameters.values()): - print_string += '\n' + str(p) + for parameter in sorted(self.listOfParameters.values()): + print_string += '\n' + str(parameter) if len(self.listOfReactions): print_string += decorate('Reactions') - for r in sorted(self.listOfReactions.values()): - print_string += '\n' + str(r) + for reaction in sorted(self.listOfReactions.values()): + print_string += '\n' + str(reaction) if len(self.listOfEvents): print_string += decorate('Events') - for e in sorted(self.listOfEvents.values()): - print_string += '\n' + str(e) + for event in sorted(self.listOfEvents.values()): + print_string += '\n' + str(event) if len(self.listOfAssignmentRules): print_string += decorate('Assignment Rules') - for ar in sorted(self.listOfAssignmentRules.values()): - print_string += '\n' + str(ar) + for assign_rule in sorted(self.listOfAssignmentRules.values()): + print_string += '\n' + str(assign_rule) if len(self.listOfRateRules): print_string += decorate('Rate Rules') - for rr in sorted(self.listOfRateRules.values()): - print_string += '\n' + str(rr) + for rate_rule in sorted(self.listOfRateRules.values()): + print_string += '\n' + str(rate_rule) if len(self.listOfFunctionDefinitions): print_string += decorate('Function Definitions') - for fd in sorted(self.listOfFunctionDefinitions.values()): - print_string += '\n' + str(fd) + for func_def in sorted(self.listOfFunctionDefinitions.values()): + print_string += '\n' + str(func_def) if self.tspan is not None: print_string += decorate('Timespan') print_string += str(self.tspan) @@ -283,7 +281,7 @@ def _problem_with_name(self, name): def _resolve_event(self, event): def validate(event): - from gillespy2.core.gillespyError import EventError + from gillespy2.core.gillespyError import EventError # pylint: disable=import-outside-toplevel if event.trigger is None or not hasattr(event.trigger, 'expression'): raise EventError('An Event must contain a valid trigger.') try: @@ -318,10 +316,10 @@ def _resolve_all_parameters(self): def _resolve_rule(self, rule): def validate(rule): - from gillespy2.core.gillespyError import RateRuleError, AssignmentRuleError + from gillespy2.core.gillespyError import RateRuleError, AssignmentRuleError # pylint: disable=import-outside-toplevel errors = {"RateRule": RateRuleError, "AssignmentRule": AssignmentRuleError} error_class = errors[type(rule).__name__] - if rule.variable == None: + if rule.variable is None: raise error_class('A GillesPy2 Rate/Assignment Rule must be associated with a valid variable') if rule.formula == '': raise error_class('Invalid Rate/Assignment Rule. Expression must be a non-empty string value') @@ -551,7 +549,7 @@ def sanitized_species_names(self): """ species_name_mapping = OrderedDict([]) for i, name in enumerate(self.listOfSpecies.keys()): - species_name_mapping[name] = 'S[{}]'.format(i) + species_name_mapping[name] = f'S[{i}]' return species_name_mapping def add_parameter(self, parameters): @@ -623,7 +621,7 @@ def get_parameter(self, name): def get_all_parameters(self): """ Get all of the parameters in the model object. - + :returns: A dict of all parameters in the model, in the form: {name : parameter object} :rtype: OrderedDict """ @@ -640,7 +638,7 @@ def sanitized_parameter_names(self): parameter_name_mapping['vol'] = 'V' for i, name in enumerate(self.listOfParameters.keys()): if name not in parameter_name_mapping: - parameter_name_mapping[name] = 'P{}'.format(i) + parameter_name_mapping[name] = f'P{i}' return parameter_name_mapping def add_reaction(self, reactions): @@ -716,7 +714,7 @@ def get_reaction(self, name): def get_all_reactions(self): """ Get all of the reactions in the model object. - + :returns: A dict of all reactions in the model, in the form: {name : reaction object}. :rtype: OrderedDict """ @@ -745,7 +743,7 @@ def add_rate_rule(self, rate_rule): raise ModelError( f"Duplicate variable in rate_rules AND assignment_rules: {rate_rule.variable}." ) - elif rate_rule.variable in rr_vars: + if rate_rule.variable in rr_vars: raise ModelError(f"Duplicate variable in rate_rules: {rate_rule.variable}.") self._resolve_rule(rate_rule) self.listOfRateRules[rate_rule.name] = rate_rule @@ -833,7 +831,7 @@ def add_assignment_rule(self, assignment_rule): raise ModelError( f"Duplicate variable in rate_rules AND assignment_rules: {assignment_rule.variable}." ) - elif assignment_rule.variable in ar_vars: + if assignment_rule.variable in ar_vars: raise ModelError(f"Duplicate variable in assignments_rules: {assignment_rule.variable}.") self._resolve_rule(assignment_rule) self.listOfAssignmentRules[assignment_rule.name] = assignment_rule @@ -844,7 +842,8 @@ def add_assignment_rule(self, assignment_rule): ) self._listOfAssignmentRules[assignment_rule.name] = sanitized_assignment_rule else: - errmsg = f"assignment_rule must be of type AssignmentRule or list of AssignmentRules not {type(assignment_rule)}" + errmsg = "assignment_rule must be of type AssignmentRule or " + errmsg += f"list of AssignmentRules not {type(assignment_rule)}." raise ModelError(errmsg) return assignment_rule @@ -989,7 +988,8 @@ def add_function_definition(self, function_definition): self._problem_with_name(function_definition.name) self.listOfFunctionDefinitions[function_definition.name] = function_definition else: - errmsg = f"function_definition must be of type FunctionDefinition or list of FunctionDefinitions not {type(function_definition)}" + errmsg = "function_definition must be of type FunctionDefinition or " + errmsg += f"list of FunctionDefinitions not {type(function_definition)}." raise ModelError(errmsg) def delete_function_definition(self, name): @@ -1006,7 +1006,7 @@ def delete_function_definition(self, name): except KeyError as err: raise ModelError( f"{self.name} does not contain a function definition named {name}." - ) + ) from err def delete_all_function_definitions(self): """ @@ -1041,20 +1041,19 @@ def get_all_function_definitions(self): def timespan(self, time_span): """ Set the time span of simulation. StochKit does not support non-uniform - timespans. + timespans. - :param time_span: Evenly-spaced list of times at which to sample the species populations during the simulation. - Best to use the form gillespy2.TimeSpan(np.linspace(, , )) + :param time_span: Evenly-spaced list of times at which to sample the species populations + during the simulation. Best to use the form gillespy2.TimeSpan(np.linspace( + , , )) :type time_span: gillespy2.TimeSpan | iterator - """ + """ if isinstance(time_span, TimeSpan) or type(time_span).__name__ == "TimeSpan": self.tspan = time_span else: self.tspan = TimeSpan(time_span) def make_translation_table(self): - from collections import ChainMap - species = self.listOfSpecies.values() reactions = self.listOfReactions.values() parameters = self.listOfParameters.values() @@ -1113,11 +1112,11 @@ def get_best_solver(self): hybrid_check = False if len(self.get_all_rate_rules()) or len(self.get_all_events()): hybrid_check = True - - if len(self.get_all_species()) and hybrid_check == False: + + if len(self.get_all_species()) and not hybrid_check: for i in self.get_all_species(): tempMode = self.get_species(i).mode - if tempMode == 'dynamic' or tempMode == 'continuous': + if tempMode in ('dynamic', 'continuous'): hybrid_check = True break @@ -1126,31 +1125,31 @@ def get_best_solver(self): hybrid_check = True chybrid_check = False - from gillespy2.solvers.cpp.build.build_engine import BuildEngine + from gillespy2.solvers.cpp.build.build_engine import BuildEngine # pylint: disable=import-outside-toplevel missing_deps = BuildEngine.get_missing_dependencies() windows_space = platform.system() == "Windows" and " " in gillespy2.__file__ can_use_cpp = len(missing_deps) == 0 and not windows_space if hybrid_check: if can_use_cpp and chybrid_check: - from gillespy2 import TauHybridCSolver + from gillespy2 import TauHybridCSolver # pylint: disable=import-outside-toplevel return TauHybridCSolver - from gillespy2 import TauHybridSolver + from gillespy2 import TauHybridSolver # pylint: disable=import-outside-toplevel return TauHybridSolver - + if can_use_cpp: - from gillespy2 import SSACSolver + from gillespy2 import SSACSolver # pylint: disable=import-outside-toplevel return SSACSolver - from gillespy2 import NumPySSASolver + from gillespy2 import NumPySSASolver # pylint: disable=import-outside-toplevel return NumPySSASolver def get_best_solver_algo(self, algorithm): """ If user has specified a particular algorithm, we return either the Python or C++ version of that algorithm """ - from gillespy2.solvers.cpp.build.build_engine import BuildEngine + from gillespy2.solvers.cpp.build.build_engine import BuildEngine # pylint: disable=import-outside-toplevel missing_deps = BuildEngine.get_missing_dependencies() windows_space = platform.system() == "Windows" and " " in gillespy2.__file__ can_use_cpp = len(missing_deps) == 0 and not windows_space @@ -1158,40 +1157,40 @@ def get_best_solver_algo(self, algorithm): if algorithm == 'Tau-Leaping': if can_use_cpp: - from gillespy2 import TauLeapingCSolver + from gillespy2 import TauLeapingCSolver # pylint: disable=import-outside-toplevel return TauLeapingCSolver - from gillespy2 import TauLeapingSolver + from gillespy2 import TauLeapingSolver # pylint: disable=import-outside-toplevel return TauLeapingSolver if algorithm == 'SSA': if can_use_cpp: - from gillespy2 import SSACSolver + from gillespy2 import SSACSolver # pylint: disable=import-outside-toplevel return SSACSolver - from gillespy2 import NumPySSASolver + from gillespy2 import NumPySSASolver # pylint: disable=import-outside-toplevel return NumPySSASolver if algorithm == 'ODE': if can_use_cpp: - from gillespy2 import ODECSolver + from gillespy2 import ODECSolver # pylint: disable=import-outside-toplevel return ODECSolver - from gillespy2 import ODESolver + from gillespy2 import ODESolver # pylint: disable=import-outside-toplevel return ODESolver if algorithm == 'Tau-Hybrid': if can_use_cpp and chybrid_check: - from gillespy2 import TauHybridCSolver + from gillespy2 import TauHybridCSolver # pylint: disable=import-outside-toplevel return TauHybridCSolver - from gillespy2 import TauHybridSolver + from gillespy2 import TauHybridSolver # pylint: disable=import-outside-toplevel return TauHybridSolver if algorithm == 'CLE': - from gillespy2 import CLESolver + from gillespy2 import CLESolver # pylint: disable=import-outside-toplevel return CLESolver - + raise ModelError("Invalid value for the argument 'algorithm' entered. " "Please enter 'SSA', 'ODE', 'CLE', 'Tau-Leaping', or 'Tau-Hybrid'.") @@ -1241,7 +1240,8 @@ def run(self, solver=None, timeout=0, t=None, increment=None, algorithm=None, ** be initialized separately to specify an algorithm. Optional, defaults to ssa solver. :type solver: gillespy.GillesPySolver - :param timeout: Allows a time_out value in seconds to be sent to a signal handler, restricting simulation run-time + :param timeout: Allows a time_out value in seconds to be sent to a signal handler, + restricting simulation run-time :type timeout: int :param t: End time of simulation @@ -1249,7 +1249,8 @@ def run(self, solver=None, timeout=0, t=None, increment=None, algorithm=None, ** :param solver_args: Solver-specific arguments to be passed to solver.run() - :param algorithm: Specify algorithm ('ODE', 'Tau-Leaping', or 'SSA') for GillesPy2 to automatically pick best solver using that algorithm. + :param algorithm: Specify algorithm ('ODE', 'Tau-Leaping', or 'SSA') for GillesPy2 to automatically + pick best solver using that algorithm. :type algorithm: str :returns: Returns a Results object that inherits UserList and contains one or more Trajectory objects that @@ -1258,8 +1259,9 @@ def run(self, solver=None, timeout=0, t=None, increment=None, algorithm=None, ** To pause a simulation and retrieve data before the simulation, keyboard interrupt the simulation by pressing control+c or pressing stop on a jupyter notebook. To resume a simulation, pass your previously ran results into the run method, and set t = to the time you wish the resuming simulation to end (run(resume=results, t=x)). - - **Pause/Resume is only supported for SINGLE TRAJECTORY simulations. T MUST BE SET OR UNEXPECTED BEHAVIOR MAY OCCUR.** + + **Pause/Resume is only supported for SINGLE TRAJECTORY simulations. + T MUST BE SET OR UNEXPECTED BEHAVIOR MAY OCCUR.** """ if solver is None: @@ -1280,16 +1282,16 @@ def run(self, solver=None, timeout=0, t=None, increment=None, algorithm=None, ** try: return solver.run(t=t, increment=increment, timeout=timeout, **solver_args) - except Exception as e: + except Exception as err: raise SimulationError( - "argument 'solver={}' to run() failed. Reason Given: {}".format(solver, e) - ) from e + f"argument 'solver={solver}' to run() failed. Reason Given: {err}" + ) from err def problem_with_name(self, *args): """ (deprecated) """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Model.problem_with_name has been deprecated. Future releases of GillesPy2 may @@ -1309,7 +1311,7 @@ def set_parameter(self, p_name, expression): parameter. May reference other parameters by name. (e.g. "k1*4") :type expression: str """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Model.set_parameter has been deprecated. Future releases of GillesPy2 may @@ -1332,7 +1334,7 @@ def validate_reactants_and_products(self, reactions): :raises ModelError: If the reaction can't be resolved. """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Model.validate_reactants_and_products has been deprecated. Future releases of @@ -1360,8 +1362,8 @@ def from_model(cls, model): Note, this method is intended to be used internally by the models 'serialization' function, which performs additional operations and - tests on the model prior to writing out the XML file. - + tests on the model prior to writing out the XML file. + You should NOT do: .. code-block:: python @@ -1378,54 +1380,54 @@ def from_model(cls, model): """ # Description - md = cls() + sml_model = cls() - d = eTree.Element('Description') + descript = eTree.Element('Description') # Prepare model for export model.compile_prep() if model.units.lower() == "concentration": - d.set('units', model.units.lower()) + descript.set('units', model.units.lower()) - d.text = model.annotation - md.document.append(d) + descript.text = model.annotation + sml_model.document.append(descript) # Number of Reactions - nr = eTree.Element('NumberOfReactions') - nr.text = str(len(model.listOfReactions)) - md.document.append(nr) + num_reacs = eTree.Element('NumberOfReactions') + num_reacs.text = str(len(model.listOfReactions)) + sml_model.document.append(num_reacs) # Number of Species - ns = eTree.Element('NumberOfSpecies') - ns.text = str(len(model.listOfSpecies)) - md.document.append(ns) + num_specs = eTree.Element('NumberOfSpecies') + num_specs.text = str(len(model.listOfSpecies)) + sml_model.document.append(num_specs) # Species spec = eTree.Element('SpeciesList') for sname in model.listOfSpecies: - spec.append(md.__species_to_element(model.listOfSpecies[sname])) - md.document.append(spec) + spec.append(sml_model.__species_to_element(model.listOfSpecies[sname])) + sml_model.document.append(spec) # Parameters params = eTree.Element('ParametersList') for pname in model.listOfParameters: - params.append(md.__parameter_to_element( + params.append(sml_model.__parameter_to_element( model.listOfParameters[pname])) vol = Parameter(name='vol', expression=model.volume) vol._evaluate() - params.append(md.__parameter_to_element(vol)) + params.append(sml_model.__parameter_to_element(vol)) - md.document.append(params) + sml_model.document.append(params) # Reactions reacs = eTree.Element('ReactionsList') for rname in model.listOfReactions: - reacs.append(md.__reaction_to_element(model.listOfReactions[rname], model.volume)) - md.document.append(reacs) + reacs.append(sml_model.__reaction_to_element(model.listOfReactions[rname], model.volume)) + sml_model.document.append(reacs) - return md + return sml_model @classmethod def from_file(cls, filepath): @@ -1433,9 +1435,9 @@ def from_file(cls, filepath): file read from disk. """ tree = eTree.parse(filepath) root = tree.getroot() - md = cls() - md.document = root - return md + sml_model = cls() + sml_model.document = root + return sml_model @classmethod def from_string(cls, string): @@ -1443,9 +1445,9 @@ def from_string(cls, string): file read from disk. """ root = eTree.fromString(string) - md = cls() - md.document = root - return md + sml_model = cls() + sml_model.document = root + return sml_model def to_model(self, name): """ Instantiates a Model object from a StochMLDocument. """ @@ -1459,8 +1461,7 @@ def to_model(self, name): name = root.find('Name') if name.text is None: raise NameError("The Name cannot be none") - else: - model.name = name.text + model.name = name.text # Set annotiation ann = root.find('Description') @@ -1493,28 +1494,28 @@ def to_model(self, name): model.units = "population" # Create parameters - for px in root.iter('Parameter'): - name = px.find('Id').text - expr = px.find('Expression').text + for sml_param in root.iter('Parameter'): + name = sml_param.find('Id').text + expr = sml_param.find('Expression').text if name.lower() == 'vol' or name.lower() == 'volume': model.volume = float(expr) else: - p = Parameter(name, expression=expr) + param = Parameter(name, expression=expr) # Try to evaluate the expression in the empty namespace # (if the expr is a scalar value) - p._evaluate() - model.add_parameter(p) + param._evaluate() + model.add_parameter(param) # Create species - for spec in root.iter('Species'): - name = spec.find('Id').text - val = spec.find('InitialPopulation').text + for sml_spec in root.iter('Species'): + name = sml_spec.find('Id').text + val = sml_spec.find('InitialPopulation').text if '.' in val: val = float(val) else: val = int(val) - s = Species(name, initial_value=val) - model.add_species([s]) + spec = Species(name, initial_value=val) + model.add_species([spec]) # The namespace_propensity for evaluating the propensity function # for reactions must contain all the species and parameters. @@ -1532,50 +1533,50 @@ def to_model(self, name): for reac in root.iter('Reaction'): try: name = reac.find('Id').text - except: - raise InvalidStochMLError("Reaction has no name.") + except Exception as err: + raise InvalidStochMLError("Reaction has no name.") from err # Type may be 'mass-action','customized' try: r_type = reac.find('Type').text - except: - raise InvalidStochMLError("No reaction type specified.") + except Exception as err: + raise InvalidStochMLError("No reaction type specified.") from err g_reactants = {} reactants = reac.find('Reactants') try: - for ss in reactants.iter('SpeciesReference'): - specname = ss.get('id') + for stoich_spec in reactants.iter('SpeciesReference'): + specname = stoich_spec.get('id') # The stochiometry should be an integer value, but some # exising StoxhKit models have them as floats. This is # why we need the slightly odd conversion below. - stoch = int(float(ss.get('stoichiometry'))) + stoch = int(float(stoich_spec.get('stoichiometry'))) # Select a reference to species with name specname sref = model.listOfSpecies[specname] try: # The sref list should only contain one element if # the XML file is valid. g_reactants[sref] = stoch - except Exception as e: - StochMLImportError(e) - except: + except Exception as err: + raise StochMLImportError(f"Reason given: {err}") from err + except Exception: # Yes, this is correct. 'reactants' can be None pass g_products = {} products = reac.find('Products') try: - for ss in products.iter('SpeciesReference'): - specname = ss.get('id') - stoch = int(float(ss.get('stoichiometry'))) + for stoich_spec in products.iter('SpeciesReference'): + specname = stoich_spec.get('id') + stoch = int(float(stoich_spec.get('stoichiometry'))) sref = model.listOfSpecies[specname] try: # The sref list should only contain one element if # the XML file is valid. g_products[sref] = stoch - except Exception as e: - raise StochMLImportError(e) - except: + except Exception as err: + raise StochMLImportError(f"Reason given: {err}") from err + except Exception: # Yes, this is correct. 'products' can be None pass @@ -1589,26 +1590,27 @@ def to_model(self, name): ratename = reac.find('Rate').text try: kwargs['rate'] = model.listOfParameters[ratename] - except KeyError as k: + except KeyError: # No paramter name is given. This is a valid use case # in StochKit. We generate a name for the paramter, # and create a new parameter instance. The parameter's # value should now be found in 'ratename'. generated_rate_name = "Reaction_" + name + \ "_rate_constant" - p = Parameter(name=generated_rate_name, + param = Parameter(name=generated_rate_name, expression=ratename) # Try to evaluate the parameter to set its value - model.add_parameter(p) + model.add_parameter(param) kwargs['rate'] = model.listOfParameters[generated_rate_name] - except Exception as e: - raise + except Exception as err: + raise StochMLImportError(f"Reason given: {err}") from err elif r_type == 'customized': try: propfunc = reac.find('PropensityFunction').text - except Exception as e: + except Exception as err: raise InvalidStochMLError( - "Found a customized propensity function, but no expression was given. {}".format(e)) + f"Found a customized propensity function, but no expression was given. {err}" + ) from err kwargs['propensity_function'] = propfunc else: raise InvalidStochMLError( @@ -1627,88 +1629,88 @@ def to_string(self): try: doc = eTree.tostring(self.document, pretty_print=True) return doc.decode("utf-8") - except: + except Exception: # Hack to print pretty xml without pretty-print # (requires the lxml module). doc = eTree.tostring(self.document) xmldoc = xml.dom.minidom.parseString(doc) - uglyXml = xmldoc.toprettyxml(indent=' ') + ugly_xml = xmldoc.toprettyxml(indent=' ') text_re = re.compile(">\n\s+([^<>\s].*?)\n\s+\g<1>\g<1>. import ast import uuid -from json.encoder import JSONEncoder import numpy as np @@ -70,13 +69,13 @@ class Reaction(SortableObject, Jsonify): rate represents the mass-action constant rate independent of volume. if a name is not provided for reactions, the name will be populated by the - model based on the order it was added. This could impact seeded simulation + model based on the order it was added. This could impact seeded simulation results if the order of addition is not preserved. """ def __init__(self, name=None, reactants=None, products=None, propensity_function=None, ode_propensity_function=None, rate=None, annotation=None, massaction=None): if massaction is not None: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ massaction has been deprecated. Future releases of GillesPy2 may not support this feature. @@ -92,7 +91,7 @@ def __init__(self, name=None, reactants=None, products=None, propensity_function ode_propensity_function = str(ode_propensity_function) if isinstance(rate, (int, float)): rate = str(rate) - + self.name = name self.reactants = {} self.products = {} @@ -100,27 +99,27 @@ def __init__(self, name=None, reactants=None, products=None, propensity_function self.annotation = annotation self.propensity_function = propensity_function self.ode_propensity_function = ode_propensity_function - + self.validate(reactants=reactants, products=products) - + if reactants is not None: - for r in reactants: - rtype = type(r).__name__ - name = r.name if rtype == 'Species' else r + for reactant in reactants: + rtype = type(reactant).__name__ + name = reactant.name if rtype == 'Species' else reactant if name in self.reactants: - self.reactants[name] += reactants[r] + self.reactants[name] += reactants[reactant] else: - self.reactants[name] = reactants[r] - + self.reactants[name] = reactants[reactant] + if products is not None: - for p in products: - ptype = type(p).__name__ - name = p.name if ptype == 'Species' else p + for product in products: + ptype = type(product).__name__ + name = product.name if ptype == 'Species' else product if name in self.products: - self.products[name] += products[p] + self.products[name] += products[product] else: - self.products[name] = products[p] - + self.products[name] = products[product] + if self.marate is not None: rtype = type(self.marate).__name__ if rtype == 'Parameter': @@ -146,29 +145,35 @@ def __str__(self): print_string = self.name if len(self.reactants): print_string += '\n\tReactants' - for r, stoich in self.reactants.items(): + for reactant, stoich in self.reactants.items(): try: - if isinstance(r, str): - print_string += '\n\t\t' + r + ': ' + str(stoich) + if isinstance(reactant, str): + print_string += '\n\t\t' + reactant + ': ' + str(stoich) else: - print_string += '\n\t\t' + r.name + ': ' + str(stoich) - except Exception as e: - print_string += '\n\t\t' + r + ': ' + 'INVALID - ' + str(e) + print_string += '\n\t\t' + reactant.name + ': ' + str(stoich) + except Exception as err: + print_string += '\n\t\t' + reactant + ': ' + 'INVALID - ' + str(err) if len(self.products): print_string += '\n\tProducts' - for p, stoich in self.products.items(): + for product, stoich in self.products.items(): try: - if isinstance(p, str): - print_string += '\n\t\t' + p + ': ' + str(stoich) + if isinstance(product, str): + print_string += '\n\t\t' + product + ': ' + str(stoich) else: - print_string += '\n\t\t' + p.name + ': ' + str(stoich) - except Exception as e: - print_string += '\n\t\t' + p + ': ' + 'INVALID - ' + str(e) + print_string += '\n\t\t' + product.name + ': ' + str(stoich) + except Exception as err: + print_string += '\n\t\t' + product + ': ' + 'INVALID - ' + str(err) print_string += '\n\tPropensity Function: ' + self.propensity_function return print_string - + class __ExpressionParser(ast.NodeTransformer): def visit_BinOp(self, node): + ''' + Visit binary operator node. + + :param node: AST node object. + :type node: ast.Node + ''' node.left = self.visit(node.left) node.right = self.visit(node.right) if isinstance(node.op, (ast.BitXor, ast.Pow)): @@ -185,16 +190,21 @@ def visit_BinOp(self, node): return call # No modification to node, classes extending NodeTransformer methods # Always return node or value - else: - return node + return node def visit_Name(self, node): + ''' + Visit name node. + + :param node: AST node object. + :type node: ast.Node + ''' # Visits Name nodes, if the name nodes "id" value is 'e', replace with numerical constant if node.id == 'e': - nameToConstant = ast.copy_location(ast.Num(float(np.e), ctx=node.ctx), node) - return nameToConstant + name_to_constant = ast.copy_location(ast.Num(float(np.e), ctx=node.ctx), node) + return name_to_constant return node - + class __ToString(ast.NodeVisitor): substitutions = { "ln": "log", @@ -207,6 +217,12 @@ def _string_changer(self, addition): self.string += addition def visit_BinOp(self, node): + ''' + Visit binary operator node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('(') self.visit(node.left) self.visit(node.op) @@ -214,14 +230,32 @@ def visit_BinOp(self, node): self._string_changer(')') def visit_Name(self, node): + ''' + Visit name name. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer(node.id) self.generic_visit(node) def visit_Num(self, node): + ''' + Visit numerical node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer(str(node.n)) self.generic_visit(node) def visit_Call(self, node): + ''' + Visit function call node. + + :param node: AST node object. + :type node: ast.Node + ''' func_name = self.substitutions.get(node.func.id) \ if node.func.id in self.substitutions \ else node.func.id @@ -235,30 +269,66 @@ def visit_Call(self, node): self._string_changer(')') def visit_Add(self, node): + ''' + Visit addition node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('+') self.generic_visit(node) def visit_Div(self, node): + ''' + Visit division node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('/') self.generic_visit(node) def visit_Mult(self, node): + ''' + Visit multiplication node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('*') self.generic_visit(node) def visit_UnaryOp(self, node): + ''' + Visit unary operator node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('(') self.visit_Usub(node) self._string_changer(')') def visit_Sub(self, node): + ''' + Visit subtraction node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('-') self.generic_visit(node) def visit_Usub(self, node): + ''' + Visit unary subtraction node. + + :param node: AST node object. + :type node: ast.Node + ''' self._string_changer('-') self.generic_visit(node) - + def _create_mass_action(self): """ Initializes the mass action propensity function given @@ -269,8 +339,8 @@ def _create_mass_action(self): # Users can still create such propensities if they really want to, # but should then use a custom propensity. total_stoch = 0 - for r in sorted(self.reactants): - total_stoch += self.reactants[r] + for reactant in sorted(self.reactants): + total_stoch += self.reactants[reactant] if total_stoch > 2: raise ReactionError( """ @@ -313,16 +383,16 @@ def _create_mass_action(self): self.propensity_function = self._create_custom_propensity(propensity_function=propensity_function) self.ode_propensity_function = self._create_custom_propensity(propensity_function=ode_propensity_function) - + def _create_custom_propensity(self, propensity_function): expr = propensity_function.replace('^', '**') expr = ast.parse(expr, mode='eval') expr = self.__ExpressionParser().visit(expr) - - newFunc = self.__ToString() - newFunc.visit(expr) - return newFunc.string - + + new_func = self.__ToString() + new_func.visit(expr) + return new_func.string + def _create_sanitized_reaction(self, n_ndx, species_mappings, parameter_mappings): name = f"R{n_ndx}" reactants = {species_mappings[species.name]: self.reactants[species] for species in self.reactants} @@ -340,7 +410,7 @@ def addProduct(self, *args, **kwargs): :param stoichiometry: Stoichiometry of this product. :type stoichiometry: int """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.addProduct has been deprecated. Future releases of GillesPy2 may @@ -371,7 +441,7 @@ def add_product(self, species, stoichiometry): self.products[name] += stoichiometry else: self.products[name] = stoichiometry - + def addReactant(self, *args, **kwargs): """ Add a reactant to this reaction (deprecated) @@ -382,7 +452,7 @@ def addReactant(self, *args, **kwargs): :param stoichiometry: Stoichiometry of this participant reactant :type stoichiometry: int """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.addReactant has been deprecated. Future releases of GillesPy2 may @@ -417,7 +487,7 @@ def add_reactant(self, species, stoichiometry): self._create_mass_action() self.validate(coverage="initialized") - + def Annotate(self, *args, **kwargs): """ Add an annotation to this reaction (deprecated). @@ -425,7 +495,7 @@ def Annotate(self, *args, **kwargs): :param annotation: Annotation note to be added to reaction :type annotation: str """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.Annotate has been deprecated. Future releases of GillesPy2 may @@ -434,13 +504,13 @@ def Annotate(self, *args, **kwargs): ) self.set_annotation(*args, **kwargs) - + def create_mass_action(self, *args, **kwargs): """ Initializes the mass action propensity function given self.reactants and a single parameter value. """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.create_mass_action has been deprecated. Future releases of GillesPy2 may @@ -465,15 +535,27 @@ def from_json(cls, json_object): new.validate(coverage="all") return new - + def sanitized_propensity_function(self, species_mappings, parameter_mappings, ode=False): + ''' + Sanitize the reaction propensity. + + :param species_mappings: Mapping of species names to sanitized species names. + :type species_mappings: dict + + :param parameter_mappings: Mapping of parameter names to sanitized parameter names. + :type parameter_mappings: dict + + :returns: The sanitized propensity. + :rtype: str + ''' names = sorted(list(species_mappings.keys()) + list(parameter_mappings.keys()), key=lambda x: len(x), reverse=True) replacements = [parameter_mappings[name] if name in parameter_mappings else species_mappings[name] for name in names] sanitized_propensity = self.ode_propensity_function if ode else self.propensity_function - for id, name in enumerate(names): - sanitized_propensity = sanitized_propensity.replace(name, "{" + str(id) + "}") + for i, name in enumerate(names): + sanitized_propensity = sanitized_propensity.replace(name, "{" + str(i) + "}") return sanitized_propensity.format(*replacements) def set_annotation(self, annotation): @@ -506,8 +588,11 @@ def set_propensities(self, propensity_function=None, ode_propensity_function=Non propensity_function = str(propensity_function) if isinstance(ode_propensity_function, (int, float)): ode_propensity_function = str(ode_propensity_function) - - self.validate(propensity_function=propensity_function, ode_propensity_function=ode_propensity_function, coverage="propensities") + + self.validate( + propensity_function=propensity_function, ode_propensity_function=ode_propensity_function, + coverage="propensities" + ) self.propensity_function = propensity_function self.ode_propensity_function = ode_propensity_function @@ -535,7 +620,7 @@ def set_rate(self, rate): """ if rate is None: raise ReactionError("rate can't be None type") - + if isinstance(rate, Parameter) or type(rate).__name__ == "Parameter": rate = rate.name elif isinstance(rate, (int, float)): @@ -557,7 +642,7 @@ def setType(self, rxntype): :param rxntype: Either "mass-action" or "customized" :type rxntype: str """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.setType has been deprecated. Future releases of GillesPy2 may not support this feature. @@ -569,13 +654,13 @@ def setType(self, rxntype): raise ReactionError("Invalid reaction type.") self.type = rxntype.lower() - self.massaction = False if self.type == 'customized' else True - + self.massaction = self.type != 'customized' + def to_dict(self): temp = vars(self).copy() # This to_dict overload is needed because, while Species is hashable (thanks to its inheretence), - # objects are not valid key values in the JSON spec. To fix this, we set the Species equal to some key 'key', + # objects are not valid key values in the JSON spec. To fix this, we set the Species equal to some key 'key', # and it's value equal to some key 'value'. temp["products"] = list({ "key": k, "value": v} for k, v in self.products.items() ) @@ -644,7 +729,7 @@ def validate(self, reactants=None, products=None, propensity_function=None, raise ReactionError("stoichiometry in products can't be None type.") if not isinstance(stoichiometry, int) or stoichiometry <= 0: raise ReactionError("stoichiometry in products must greater than 0 and of type int.") - + if coverage in ("all", "build", "propensities"): if propensity_function is None: propensity_function = self.propensity_function @@ -700,21 +785,23 @@ def validate(self, reactants=None, products=None, propensity_function=None, if self.massaction and self.type == "mass-action": raise ReactionError("You must specify either a mass-action rate or propensity functions") if self.massaction or self.type == "mass-action": - errmsg = "Invalid customized reaction. Customized reactions require massaction=False and type='customized'" + errmsg = "Invalid customized reaction. Customized reactions require " + errmsg += "massaction=False and type='customized'" raise ReactionError(errmsg) else: if not self.massaction and self.type == "customized": raise ReactionError("You cannot set the rate and simultaneously set propensity functions.") if not self.massaction or self.type == "customized": - errmsg = "Invalid mass-action reaction. Mass-action reactions require massaction=True and type='mass-action'" + errmsg = "Invalid mass-action reaction. Mass-action reactions require " + errmsg += "massaction=True and type='mass-action'" raise ReactionError(errmsg) - + def verify(self, *args, **kwargs): """ Check if the reaction is properly formatted. Does nothing on sucesss, raises and error on failure. """ - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning( """ Reaction.verify has been deprecated. Future releases of GillesPy2 may diff --git a/gillespy2/core/results.py b/gillespy2/core/results.py index fa4ddc7fa..07b8ed6d9 100644 --- a/gillespy2/core/results.py +++ b/gillespy2/core/results.py @@ -13,26 +13,33 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . - +import os +import csv from datetime import datetime -from gillespy2.core.gillespyError import * -from gillespy2.core.jsonify import Jsonify from collections import UserDict, UserList -# List of 50 hex color values used for plotting graphs -def common_rgb_values(): - return ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', - '#bcbd22', '#17becf', '#ff0000', '#00ff00', '#0000ff', '#ffff00', '#00ffff', '#ff00ff', - '#800000', '#808000', '#008000', '#800080', '#008080', '#000080', '#ff9999', '#ffcc99', - '#ccff99', '#cc99ff', '#ffccff', '#62666a', '#8896bb', '#77a096', '#9d5a6c', '#9d5a6c', - '#eabc75', '#ff9600', '#885300', '#9172ad', '#a1b9c4', '#18749b', '#dadecf', '#c5b8a8', - '#000117', '#13a8fe', '#cf0060', '#04354b', '#0297a0', '#037665', '#eed284', '#442244', - '#ffddee', '#702afb'] +import numpy as np +from gillespy2.core.jsonify import Jsonify +from gillespy2.core.gillespyError import ValidationError + +def common_rgb_values(): + ''' + List of 50 hex color values used for plotting graphs + ''' + return [ + '#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', + '#bcbd22', '#17becf', '#ff0000', '#00ff00', '#0000ff', '#ffff00', '#00ffff', '#ff00ff', + '#800000', '#808000', '#008000', '#800080', '#008080', '#000080', '#ff9999', '#ffcc99', + '#ccff99', '#cc99ff', '#ffccff', '#62666a', '#8896bb', '#77a096', '#9d5a6c', '#9d5a6c', + '#eabc75', '#ff9600', '#885300', '#9172ad', '#a1b9c4', '#18749b', '#dadecf', '#c5b8a8', + '#000117', '#13a8fe', '#cf0060', '#04354b', '#0297a0', '#037665', '#eed284', '#442244', + '#ffddee', '#702afb' + ] def _plot_iterate(self, show_labels=True, included_species_list=[]): - import matplotlib.pyplot as plt + import matplotlib.pyplot as plt # pylint: disable=import-outside-toplevel for i, species in enumerate(self.data): if species != 'time': @@ -57,7 +64,7 @@ def _plotplotly_iterate(trajectory, show_labels=True, trace_list=None, line_dict if trace_list is None: trace_list = [] - import plotly.graph_objs as go + import plotly.graph_objs as go # pylint: disable=import-outside-toplevel for i, species in enumerate(trajectory.data): if species != 'time': @@ -126,8 +133,8 @@ def __init__(self, data, model=None, solver_name="Undefined solver name", rc=0): self.status = status_list[rc] def __getitem__(self, key): - if type(key) is int: - from gillespy2.core import log + if isinstance(key, int): + from gillespy2.core import log # pylint: disable=import-outside-toplevel species = list(self.data.keys())[key] msg = "Trajectory is of type dictionary." msg += f"Use trajectory['[{species}]'] instead of trajectory[{key}]['{species}']" @@ -155,25 +162,22 @@ def __init__(self, data): def __getattribute__(self, key): if key in ('model', 'solver_name', 'rc', 'status'): if len(self.data) > 1: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel msg = f"Results is of type list. Use results[i]['{key}'] instead of results['{key}']" log.warning(msg) return getattr(Results.__getattribute__(self, key='data')[0], key) - else: - return UserList.__getattribute__(self, key) + return UserList.__getattribute__(self, key) def __getitem__(self, key): if key == 'data': return UserList.__getitem__(self, key) if isinstance(key, str): if len(self.data) > 1: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel msg = f"Results is of type list. Use results[i]['{key}'] instead of results['{key}']" log.warning(msg) return self.data[0][key] - else: - return(UserList.__getitem__(self,key)) - raise KeyError(key) + return UserList.__getitem__(self,key) def __add__(self, other): combined_data = Results(data=(self.data + other.data)) @@ -181,7 +185,7 @@ def __add__(self, other): consistent_model = combined_data._validate_model() if not consistent_solver: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning("Results objects contain Trajectory objects from multiple solvers.") if not consistent_model: @@ -192,8 +196,7 @@ def __add__(self, other): def __radd__(self, other): if other == 0: return self - else: - return self.__add__(other) + return self.__add__(other) def _validate_model(self, reference=None): is_valid = True @@ -233,16 +236,21 @@ def _validate_title(self, show_title): return title def to_array(self): - import numpy as np + ''' + Convert the results object into a numpy array. + + :returns: Array containing the result of the simulation. + :rtype: numpy.ndarray + ''' results = [] size1 = len(self.data[0]['time']) size2 = len(self.data[0]) - for trajectory in range(0, len(self.data)): - newArray = np.zeros((size1, size2)) - for i, key in enumerate(self.data[trajectory]): - newArray[:, i] = self.data[trajectory][key] - results.append(newArray) + for trajectory in self.data: + new_array = np.zeros((size1, size2)) + for j, key in enumerate(trajectory): + new_array[:, j] = trajectory[key] + results.append(new_array) return np.array(results) @classmethod @@ -259,22 +267,21 @@ def build_from_solver_results(cls, solver, live_output_options): :type live_output_options: dict """ if solver.rc == 33: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning('GillesPy2 simulation exceeded timeout.') if hasattr(solver.result[0], 'shape'): return solver.result if len(solver.result) > 0: results_list = [] - for i in range(0, len(solver.result)): - temp = Trajectory(data=solver.result[i], model=solver.model, solver_name=solver.name, rc=solver.rc) + for result in solver.result: + temp = Trajectory(data=result, model=solver.model, solver_name=solver.name, rc=solver.rc) results_list.append(temp) results = Results(results_list) if "type" in live_output_options.keys() and live_output_options['type'] == "graph": results.plot() return results - else: - raise ValueError("number_of_trajectories must be non-negative and non-zero") + raise ValueError("number_of_trajectories must be non-negative and non-zero") def to_csv(self, path=None, nametag=None, stamp=None): """ @@ -290,9 +297,6 @@ def to_csv(self, path=None, nametag=None, stamp=None): :param stamp: Allows the user to optionally "tag" the directory (not included files). Default is timestamp. :type stamp: str """ - import csv - import os - if stamp is None: now = datetime.now() stamp = datetime.timestamp(now) @@ -312,13 +316,13 @@ def to_csv(self, path=None, nametag=None, stamp=None): field_names = [] for species in trajectory: # build the header field_names.append(species) - with open(filename, 'w', newline='') as csv_file: + with open(filename, 'w', newline='', encoding="utf-8") as csv_file: csv_writer = csv.writer(csv_file) csv_writer.writerow(field_names) # write the header - for n,time in enumerate(trajectory['time']): # write all lines of the CSV file + for j, _ in enumerate(trajectory['time']): # write all lines of the CSV file this_line=[] for species in trajectory: # build one line of the CSV file - this_line.append(trajectory[species][n]) + this_line.append(trajectory[species][j]) csv_writer.writerow(this_line) # write one line of the CSV file def plot(self, index=None, xaxis_label="Time", xscale='linear', yscale='linear', yaxis_label="Value", @@ -354,14 +358,14 @@ def plot(self, index=None, xaxis_label="Time", xscale='linear', yscale='linear', :type figsize: tuple of ints (x,y) """ - import matplotlib.pyplot as plt - from collections.abc import Iterable + import matplotlib.pyplot as plt # pylint: disable=import-outside-toplevel + from collections.abc import Iterable # pylint: disable=import-outside-toplevel trajectory_list = [] if isinstance(index, Iterable): for i in index: trajectory_list.append(self.data[i]) elif isinstance(index, int): - trajectory_list.append(self.data[index]) + trajectory_list.append(self.data[index]) else: trajectory_list = self.data @@ -369,7 +373,7 @@ def plot(self, index=None, xaxis_label="Time", xscale='linear', yscale='linear', title = self._validate_title(show_title) if len(trajectory_list) < 2: - multiple_graphs = False + multiple_graphs = False if multiple_graphs: for i, trajectory in enumerate(trajectory_list): @@ -386,8 +390,8 @@ def plot(self, index=None, xaxis_label="Time", xscale='linear', yscale='linear', else: try: plt.style.use(style) - except: - from gillespy2.core import log + except Exception: + from gillespy2.core import log # pylint: disable=import-outside-toplevel msg = f"Invalid matplotlib style. Try using one of the following {plt.style.available}" log.warning(msg) plt.style.use("default") @@ -453,8 +457,8 @@ def plotplotly(self, index=None, xaxis_label="Time", yaxis_label="Value", title= :type **layout_args: dict """ - from plotly.offline import init_notebook_mode, iplot - import plotly.graph_objs as go + from plotly.offline import init_notebook_mode, iplot # pylint: disable=import-outside-toplevel + import plotly.graph_objs as go # pylint: disable=import-outside-toplevel # Backwards compatibility with xaxis_label argument (which duplicates plotly's xaxis_title argument) if layout_args.get('xaxis_title') is not None: @@ -466,7 +470,7 @@ def plotplotly(self, index=None, xaxis_label="Time", yaxis_label="Value", title= init_notebook_mode(connected=True) - from collections.abc import Iterable + from collections.abc import Iterable # pylint: disable=import-outside-toplevel trajectory_list = [] if isinstance(index, Iterable): for i in index: @@ -488,7 +492,7 @@ def plotplotly(self, index=None, xaxis_label="Time", yaxis_label="Value", title= if multiple_graphs: - from plotly import subplots + from plotly import subplots # pylint: disable=import-outside-toplevel fig = subplots.make_subplots(print_grid=False, rows=int(number_of_trajectories/2) + int(number_of_trajectories % 2), cols=2) @@ -501,14 +505,15 @@ def plotplotly(self, index=None, xaxis_label="Time", yaxis_label="Value", title= trace_list = _plotplotly_iterate(trajectory, trace_list=[], included_species_list= included_species_list) - for k in range(0, len(trace_list)): + for trace in trace_list: if i % 2 == 0: - fig.append_trace(trace_list[k], int(i/2) + 1, 1) + fig.append_trace(trace, int(i/2) + 1, 1) else: - fig.append_trace(trace_list[k], int(i/2) + 1, 2) + fig.append_trace(trace, int(i/2) + 1, 2) - fig['layout'].update(autosize=True, height=400*len(trajectory_list), showlegend=show_legend, title=title - ) + fig['layout'].update( + autosize=True, height=400*len(trajectory_list), showlegend=show_legend, title=title + ) else: trace_list = [] for i, trajectory in enumerate(trajectory_list): @@ -531,10 +536,10 @@ def plotplotly(self, index=None, xaxis_label="Time", yaxis_label="Value", title= fig['data'] = trace_list fig['layout'] = layout - if return_plotly_figure: - return fig - else: + if not return_plotly_figure: iplot(fig) + return None + return fig def average_ensemble(self): """ @@ -554,7 +559,8 @@ def average_ensemble(self): output_trajectory['time'] = trajectory_list[0]['time'] - for i in range(0, number_of_trajectories): # Add every value of every Trajectory Dict into one output Trajectory + # Add every value of every Trajectory Dict into one output Trajectory + for i in range(0, number_of_trajectories): trajectory_dict = trajectory_list[i] for species in trajectory_dict: if species == 'time': @@ -578,20 +584,20 @@ def stddev_ensemble(self, ddof=0): trajectories' outputs. :param ddof: Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents - the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard deviation - where ddof is 0. + the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard + deviation where ddof is 0. :type ddof: int :returns: the Results object """ - from math import sqrt + from math import sqrt # pylint: disable=import-outside-toplevel trajectory_list = self.data number_of_trajectories = len(trajectory_list) if ddof == number_of_trajectories: - from gillespy2.core import log + from gillespy2.core import log # pylint: disable=import-outside-toplevel log.warning("ddof must be less than the number of trajectories. Using ddof of 0") ddof = 0 @@ -655,8 +661,8 @@ def plotplotly_mean_stdev(self, xaxis_label="Time", yaxis_label="Value", title=N :type return_plotly_figure: bool :param ddof: Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents - the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard deviation - where ddof is 0. + the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard + deviation where ddof is 0. :type ddof: int :param **layout_args: Optional additional arguments to be passed to plotlys layout constructor. @@ -673,8 +679,8 @@ def plotplotly_mean_stdev(self, xaxis_label="Time", yaxis_label="Value", title=N average_trajectory = self.average_ensemble().data[0] stddev_trajectory = self.stddev_ensemble(ddof=ddof).data[0] - from plotly.offline import init_notebook_mode, iplot - import plotly.graph_objs as go + from plotly.offline import init_notebook_mode, iplot # pylint: disable=import-outside-toplevel + import plotly.graph_objs as go # pylint: disable=import-outside-toplevel init_notebook_mode(connected=True) @@ -747,11 +753,10 @@ def plotplotly_mean_stdev(self, xaxis_label="Time", yaxis_label="Value", title=N ) fig = dict(data=trace_list, layout=layout) - if return_plotly_figure: - return fig - else: + if not return_plotly_figure: iplot(fig) - + return None + return fig def plot_mean_stdev(self, xscale='linear', yscale='linear', xaxis_label="Time", yaxis_label="Value" , title=None, show_title=False, style="default", show_legend=True, included_species_list=[], @@ -782,8 +787,8 @@ def plot_mean_stdev(self, xscale='linear', yscale='linear', xaxis_label="Time", :type included_species_list: list :param ddof: Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents - the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard deviation - where ddof is 0. + the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard + deviation where ddof is 0. :type ddof: int :type save_png: bool or str @@ -795,12 +800,12 @@ def plot_mean_stdev(self, xscale='linear', yscale='linear', xaxis_label="Time", average_result = self.average_ensemble().data[0] stddev_trajectory = self.stddev_ensemble(ddof=ddof).data[0] - import matplotlib.pyplot as plt + import matplotlib.pyplot as plt # pylint: disable=import-outside-toplevel try: plt.style.use(style) - except: - from gillespy2.core import log + except Exception: + from gillespy2.core import log # pylint: disable=import-outside-toplevel msg = f"Invalid matplotlib style. Try using one of the following {plt.style.available}" log.warning(msg) plt.style.use("default") @@ -814,12 +819,12 @@ def plot_mean_stdev(self, xscale='linear', yscale='linear', xaxis_label="Time", if species not in included_species_list and included_species_list: continue - lowerBound = [a-b for a, b in zip(average_result[species], stddev_trajectory[species])] - upperBound = [a+b for a, b in zip(average_result[species], stddev_trajectory[species])] + lower_bound = [a-b for a, b in zip(average_result[species], stddev_trajectory[species])] + upper_bound = [a+b for a, b in zip(average_result[species], stddev_trajectory[species])] - plt.fill_between(average_result['time'], lowerBound, upperBound, color='whitesmoke') - plt.plot(average_result['time'], upperBound, color='grey', linestyle='dashed') - plt.plot(average_result['time'], lowerBound, color='grey', linestyle='dashed') + plt.fill_between(average_result['time'], lower_bound, upper_bound, color='whitesmoke') + plt.plot(average_result['time'], upper_bound, color='grey', linestyle='dashed') + plt.plot(average_result['time'], lower_bound, color='grey', linestyle='dashed') plt.plot(average_result['time'], average_result[species], label=species) if not show_title: @@ -846,12 +851,87 @@ def plot_mean_stdev(self, xscale='linear', yscale='linear', xaxis_label="Time", # for backwards compatability, we need to keep the old name around def plotplotly_std_dev_range(self, **kwargs): - from gillespy2.core import log - log.warning("The plotplotly_std_dev_range function has been deprecated. This function will be removed in a future release. Please use plotplotly_mean_stdev instead.") + """ + Plot a plotly graph depicting the mean and standard deviation of a results object + + :param xaxis_label: The label for the x-axis + :type xaxis_label: str + + :param yaxis_label: The label for the y-axis + :type yaxis_label: str + + :param title: The title of the graph + :type title: str + + :param show_title: If True, title will be shown on graph. + :type show_title: bool + + :param show_legend: Default True, if False, legend will not be shown on graph. + :type show_legend: bool + + :param included_species_list: A list of strings describing which species to include. By default displays all + species. + :type included_species_list: list + + :param return_plotly_figure: Whether or not to return a figure dicctionary of data(graph object traces) and + layout which may be edited by the user + :type return_plotly_figure: bool + + :param ddof: Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents + the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard + deviation where ddof is 0. + :type ddof: int + + :param **layout_args: Optional additional arguments to be passed to plotlys layout constructor. + :type **layout_args: dict + """ + from gillespy2.core import log # pylint: disable=import-outside-toplevel + log.warning( + "The plotplotly_std_dev_range function has been deprecated. This function will be removed in a " + "future release. Please use plotplotly_mean_stdev instead." + ) self.plotplotly_mean_stdev(**kwargs) # for backwards compatability, we need to keep the old name around def plot_std_dev_range(self, **kwargs): - from gillespy2.core import log - log.warning("The plot_std_dev_range function has been deprecated. This function will be removed in a future release. Please use plot_mean_stdev instead.") + """ + Plot a matplotlib graph depicting mean and standard deviation of a results object. + + :param xaxis_label: The label for the x-axis + :type xaxis_label: str + + :param yaxis_label: The label for the y-axis + :type yaxis_label: str + + :param title: The title of the graph + :type title: str + + :param show_title: Default False, if True, title will be displayed on the graph. + :type show_title: bool + + :param style: Matplotlib style to be displayed on graph. + :type style: str + + :param show_legend: Default to True, if False, legend will not be shown on graph. + :type show_legend: bool + + :param included_species_list: A list of strings describing which species to include. By default displays all + species. + :type included_species_list: list + + :param ddof: Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents + the number of trajectories. Sample standard deviation uses ddof of 1. Defaults to population standard + deviation where ddof is 0. + :type ddof: int + + :type save_png: bool or str + + :param figsize: The size of the graph. A tuple of the form (width,height). Is (18,10) by default. + :type figsize: tuple of ints (x,y) + """ + from gillespy2.core import log # pylint: disable=import-outside-toplevel + log.warning( + "The plot_std_dev_range function has been deprecated. This function will be removed in a " + "future release. Please use plot_mean_stdev instead." + ) self.plot_mean_stdev(**kwargs) diff --git a/gillespy2/core/sortableobject.py b/gillespy2/core/sortableobject.py index 1492c91db..a723ae58b 100644 --- a/gillespy2/core/sortableobject.py +++ b/gillespy2/core/sortableobject.py @@ -47,4 +47,4 @@ def __hash__(self): self._hash = hash(self.name) else: self._hash = hash(self) - return self._hash \ No newline at end of file + return self._hash diff --git a/gillespy2/core/species.py b/gillespy2/core/species.py index 37f12c6ae..8ffbb97ae 100644 --- a/gillespy2/core/species.py +++ b/gillespy2/core/species.py @@ -54,8 +54,8 @@ class Species(SortableObject, Jsonify): :type switch_tol: float :param switch_min: ***FOR USE WITH TauHybridSolver ONLY*** - Minimum population value at which species will be represented as continuous. If a value is given, switch_min will be - used instead of switch_tol + Minimum population value at which species will be represented as continuous. + If a value is given, switch_min will be used instead of switch_tol :type switch_min: float :raises SpeciesError: Arg is of invalid type. Required arg set to None. Arg value is outside of accepted bounds. @@ -86,13 +86,6 @@ def __init__(self, name=None, initial_value=0, constant=False, boundary_conditio def __str__(self): print_string = self.name print_string += ': ' + str(self.initial_value) - ''' - print_string += '\n\tInitial Value: ' + str(self.initial_value) - print_string += '\n\tConstant: ' + str(self.constant) - print_string += '\n\tBoundary Condition: ' + str(self.boundary_condition) - print_string += '\n\tMode: ' + self.mode - print_string += '\n\tAllow Negative Populations: ' + str(self.allow_negative_populations) - ''' return print_string def set_initial_value(self, initial_value): @@ -157,7 +150,9 @@ def validate(self, initial_value=None, coverage="all"): """ ) if not self.allow_negative_populations and initial_value < 0: - raise SpeciesError('A species initial value must be non-negative unless allow_negative_populations=True') + raise SpeciesError( + 'A species initial value must be non-negative unless allow_negative_populations=True' + ) # Check constant if coverage in ("all", "constant"): diff --git a/gillespy2/core/timespan.py b/gillespy2/core/timespan.py index 4797299d4..e1a4a3691 100644 --- a/gillespy2/core/timespan.py +++ b/gillespy2/core/timespan.py @@ -13,9 +13,10 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import numpy as np from collections.abc import Iterator +import numpy as np + from gillespy2.core.jsonify import Jsonify from .gillespyError import TimespanError @@ -24,10 +25,10 @@ class TimeSpan(Iterator, Jsonify): Model timespan that describes the duration to run the simulation and at which timepoint to sample the species populations during the simulation. - :param items: Evenly-spaced list of times at which to sample the species populations during the simulation. + :param items: Evenly-spaced list of times at which to sample the species populations during the simulation. Best to use the form np.linspace(, , ) :type items: list, tuple, range, or numpy.ndarray - + :raises TimespanError: items is an invalid type. """ def __init__(self, items): @@ -124,7 +125,7 @@ def validate(self): raise TimespanError("Timespans must contain values.") if self.items[0] < 0: raise TimespanError("Simulation must run from t=0 to end time (t must always be positive).") - + first_diff = self.items[1] - self.items[0] other_diff = self.items[2:] - self.items[1:-1] isuniform = np.isclose(other_diff, first_diff).all() From f55d67011fd75c75f1c00a2bf0e6b596d60bc3fe Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 16:16:39 -0500 Subject: [PATCH 19/47] Updates to the pylintrc file. --- .pylintrc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.pylintrc b/.pylintrc index 27e483833..da4869507 100644 --- a/.pylintrc +++ b/.pylintrc @@ -394,10 +394,12 @@ good-names=i, # Good variable names regexes, separated by a comma. If names match any regex, # they will always be accepted -good-names-rgxs=listOfParameters, listOfSpecies, listOfReaction, listOfAssignmentRules, +good-names-rgxs=t, rc, listOfParameters, listOfSpecies, listOfReaction, listOfAssignmentRules, listOfRateRules, listOfEvents, listOfFunctionDefinitions, _listOfParameters, _listOfSpecies, _listOfReaction, _listOfAssignmentRules, - _listOfRateRules, _listOfEvents, _listOfFunctionDefinitions + _listOfRateRules, _listOfEvents, _listOfFunctionDefinitions, + __ExpressionParser, __ToString, visit_BinOp, visit_Name, visit_Num, visit_Call, + visit_Add, visit_Div, visit_Mult, visit_UnaryOp, visit_Sub, visit_Usub # Include a hint for the correct naming format with invalid-name. include-naming-hint=no From f87d693b0dcec3d92920b68ff9ee4ed7fa9e09b4 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 16:58:13 -0500 Subject: [PATCH 20/47] Pylint changes to the SBML directory. New pylint score == 9.56. --- gillespy2/sbml/SBMLexport.py | 8 +- gillespy2/sbml/SBMLimport.py | 219 +++++++++++++++++------------------ 2 files changed, 109 insertions(+), 118 deletions(-) diff --git a/gillespy2/sbml/SBMLexport.py b/gillespy2/sbml/SBMLexport.py index 8370c15d7..a254ea607 100644 --- a/gillespy2/sbml/SBMLexport.py +++ b/gillespy2/sbml/SBMLexport.py @@ -17,8 +17,8 @@ import ast try: import libsbml -except ImportError: - raise ImportError('libsbml is required to convert GillesPy models to SBML files.') +except ImportError as err: + raise ImportError('libsbml is required to convert GillesPy models to SBML files.') from err def __get_math(math): @@ -82,7 +82,7 @@ def __add_events(event_list, model): evt = model.createEvent() evt.setId(name) evt.setUseValuesFromTriggerTime(event.use_values_from_trigger_time) - + if event.delay is not None: delay = __get_math(event.delay) dly = evt.createDelay() @@ -136,7 +136,7 @@ def __add_function_definitions(function_definition_list, model): def __write_to_file(document, path): writer = libsbml.SBMLWriter() - with open(path, "w") as sbml_file: + with open(path, "w", encoding="utf-8") as sbml_file: sbml_file.write(writer.writeSBMLToString(document)) def export(model, path=None): diff --git a/gillespy2/sbml/SBMLimport.py b/gillespy2/sbml/SBMLimport.py index 43631e708..ac47d84cc 100644 --- a/gillespy2/sbml/SBMLimport.py +++ b/gillespy2/sbml/SBMLimport.py @@ -13,20 +13,18 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . - -import os -import gillespy2 -import numpy as np -import math import re +import math try: import libsbml -except ImportError: - raise ImportError('libsbml is required to convert SBML files for GillesPy.') +except ImportError as err: + raise ImportError('libsbml is required to convert SBML files for GillesPy.') from err -from gillespy2.core import gillespyError +import numpy as np +import gillespy2 +from gillespy2.core.gillespyError import InvalidModelError, SBMLError init_state = {'INF': np.inf, 'NaN': np.nan} postponed_evals = {} @@ -34,9 +32,11 @@ def piecewise(*args): args = list(args) sol = None - if len(args) % 2: args.append(True) + if len(args) % 2: + args.append(True) for i, arg in enumerate(args): - if not i % 2: continue + if not i % 2: + continue if arg: sol = args[i-1] break @@ -54,17 +54,17 @@ def __read_sbml_model(filename): converter_code = 0 converter_code = -10 - errors.append(["SBML {0}, code {1}, line {2}: {3}".format(error.getSeverityAsString(), error.getErrorId(), - error.getLine(), error.getMessage()), - converter_code]) + errmsg = f"SBML {error.getSeverityAsString()}, code {error.getErrorId()}, " + errmsg += f"line {error.getLine()}: {error.getMessage()}" + errors.append([errmsg, converter_code]) if min([code for error, code in errors] + [0]) < 0: return None, errors sbml_model = document.getModel() return sbml_model, errors -def __get_math(math): - math_str = libsbml.formulaToL3String(math) +def __get_math(formula): + math_str = libsbml.formulaToL3String(formula) replacements = { r'\bln\b': 'log', r'\^': '**', @@ -98,17 +98,15 @@ def __get_species(sbml_model, gillespy_model, errors): if population == int_value: population = int_value else: # Else convert concentration to population - postponed_evals[name] = '{} * {}'.format(name, cid) + postponed_evals[name] = f'{name} * {cid}' value = concentration else: # Treat as concentration if concentration is not None: # If concentration is provided value = concentration else: # Else convert population to concentration - postponed_evals[name] = '{} / {}'.format(name, cid) + postponed_evals[name] = f'{name} / {cid}' value = population - - constant = species.getConstant() boundary_condition = species.getBoundaryCondition() is_negative = value < 0 @@ -117,7 +115,7 @@ def __get_species(sbml_model, gillespy_model, errors): constant=constant, boundary_condition=boundary_condition) gillespy_model.add_species([gillespy_species]) init_state[name] = value - + def __get_parameters(sbml_model, gillespy_model): for i in range(sbml_model.getNumParameters()): @@ -153,32 +151,21 @@ def __get_compartments(sbml_model, gillespy_model): init_state[name] = value gillespy_model.add_parameter([gillespy_parameter]) - ''' - for i in range(sbml_model.getNumCompartments()): - compartment = sbml_model.getCompartment(i) - vol = compartment.getSize() - gillespy_model.volume = vol - - errors.append([ - "Compartment '{0}' found on line '{1}' with volume '{2}' and dimension '{3}'. gillespy " - "assumes a single well-mixed, reaction volume".format( - compartment.getId(), compartment.getLine(), compartment.getVolume(), - compartment.getSpatialDimensions()), -5]) - ''' -def traverse_math(node, old_id, new_id): - if node is None: return +def __traverse_math(node, old_id, new_id): + if node is None: + return for i in range(node.getNumChildren()): if node.getChild(i).getName() == old_id: new_node = libsbml.ASTNode() new_node.setName(new_id) node.replaceChild(i, new_node) - traverse_math(node.getChild(i), old_id, new_id) + __traverse_math(node.getChild(i), old_id, new_id) def __get_kinetic_law(sbml_model, gillespy_model, reaction): kinetic_law = reaction.getKineticLaw() if kinetic_law is None: - raise gillespyError.InvalidModelError( + raise InvalidModelError( f"Failed to load SBML model: Reaction '{reaction}' is missing its propensity function." ) @@ -186,20 +173,20 @@ def __get_kinetic_law(sbml_model, gillespy_model, reaction): params = kinetic_law.getListOfParameters() local_params = kinetic_law.getListOfLocalParameters() for i in range(kinetic_law.getNumLocalParameters()): - lp = local_params.get(i) - old_id = lp.getId() - new_id = ('{}_{}'.format(reaction.getId(), lp.getId())) - traverse_math(tree, old_id, new_id) - lp.setId(new_id) - gillespy_parameter = gillespy2.Parameter(name=new_id, expression=lp.getValue()) + local_param = local_params.get(i) + old_id = local_param.getId() + new_id = (f'{reaction.getId()}_{local_param.getId()}') + __traverse_math(tree, old_id, new_id) + local_param.setId(new_id) + gillespy_parameter = gillespy2.Parameter(name=new_id, expression=local_param.getValue()) gillespy_model.add_parameter([gillespy_parameter]) for i in range(kinetic_law.getNumParameters()): - p = params.get(i) - if not p.getId() in gillespy_model.listOfParameters: - gillespy_parameter = gillespy2.Parameter(name=p.getId(), expression=p.getValue()) + param = params.get(i) + if not param.getId() in gillespy_model.listOfParameters: + gillespy_parameter = gillespy2.Parameter(name=param.getId(), expression=param.getValue()) gillespy_model.add_parameter([gillespy_parameter]) return tree - + def __get_reactions(sbml_model, gillespy_model, errors): # reactions for i in range(sbml_model.getNumReactions()): @@ -217,42 +204,42 @@ def __get_reactions(sbml_model, gillespy_model, errors): reactant = reaction.getReactant(j) species = reactant.getSpecies() - if species == "EmptySet": continue + if species == "EmptySet": + continue + stoichiometry = reactant.getStoichiometry() + if isinstance(stoichiometry, float): + if int(stoichiometry) != stoichiometry: + logmsg = f"Reaction {name} contains a float stoichiometry for reactant {species}. " + logmsg += "Please check your model as this may cause inaccuracies in the results." + gillespy2.log.warning(logmsg) + stoichiometry = int(stoichiometry) + + if species in r_set: + reactants[species] += stoichiometry else: - stoichiometry = reactant.getStoichiometry() - if isinstance(stoichiometry, float): - if int(stoichiometry) != stoichiometry: - logmsg = f"Reaction {name} contains a float stoichiometry for reactant {species}. " - logmsg += "Please check your model as this may cause inaccuracies in the results." - gillespy2.log.warning(logmsg) - stoichiometry = int(stoichiometry) - - if species in r_set: - reactants[species] += stoichiometry - else: - r_set.add(species) - reactants[species] = stoichiometry + r_set.add(species) + reactants[species] = stoichiometry # get products for j in range(reaction.getNumProducts()): product = reaction.getProduct(j) species = product.getSpecies() - if species == "EmptySet": continue + if species == "EmptySet": + continue + stoichiometry = product.getStoichiometry() + if isinstance(stoichiometry, float): + if int(stoichiometry) != stoichiometry: + logmsg = f"Reaction {name} contains a float stoichiometry for product {species}. " + logmsg += "Please check your model as this may cause inaccuracies in the results." + gillespy2.log.warning(logmsg) + stoichiometry = int(stoichiometry) + + if species in p_set: + products[species] += stoichiometry else: - stoichiometry = product.getStoichiometry() - if isinstance(stoichiometry, float): - if int(stoichiometry) != stoichiometry: - logmsg = f"Reaction {name} contains a float stoichiometry for product {species}. " - logmsg += "Please check your model as this may cause inaccuracies in the results." - gillespy2.log.warning(logmsg) - stoichiometry = int(stoichiometry) - - if species in p_set: - products[species] += stoichiometry - else: - p_set.add(species) - products[species] = stoichiometry + p_set.add(species) + products[species] = stoichiometry gillespy_reaction = gillespy2.Reaction(name=name, reactants=reactants, products=products, propensity_function=propensity) @@ -281,9 +268,9 @@ def __get_rules(sbml_model, gillespy_model, errors): gillespy_model.add_species([species]) t = [] - + if rule.isAssignment(): - assign_value = eval(rule_string, {**eval_globals, **init_state}) + _ = eval(rule_string, {**eval_globals, **init_state}) postponed_evals[rule_variable] = rule_string gillespy_rule = gillespy2.AssignmentRule(name=rule_name, variable=rule_variable, formula=rule_string) @@ -306,22 +293,21 @@ def __get_rules(sbml_model, gillespy_model, errors): else: msg = "Rule" - errors.append(["{0} '{1}' found on line '{2}' with equation '{3}'. gillespy does not support SBML Algebraic Rules".format( - msg, rule.getId(), rule.getLine(), libsbml.formulaToString(rule.getMath())), -5]) + errmsg = f"{msg} '{rule.getId()}' found on line '{rule.getLine()}' with equation " + errmsg += f"'{libsbml.formulaToString(rule.getMath())}'. gillespy does not support SBML Algebraic Rules." + errors.append([errmsg, -5]) -def __get_constraints(sbml_model, gillespy_model): +def __get_constraints(sbml_model, gillespy_model, errors): for i in range(sbml_model.getNumConstraints()): constraint = sbml_model.getConstraint(i) - errors.append([ - "Constraint '{0}' found on line '{1}' with equation '{2}'. gillespy does not support SBML " - "Constraints".format( - constraint.getId(), constraint.getLine(), libsbml.formulaToString(constraint.getMath())), - -5]) + errmsg = f"Constraint '{constraint.getId()}' found on line '{constraint.getLine()}' with equation " + errmsg += f"'{libsbml.formulaToString(constraint.getMath())}'. gillespy does not support SBML Constraints" + errors.append([errmsg, -5]) def __get_function_definitions(sbml_model, gillespy_model): # TODO: - # DOES NOT CURRENTLY SUPPORT ALL MATHML + # DOES NOT CURRENTLY SUPPORT ALL MATHML # ALSO DOES NOT SUPPORT NON-MATHML for i in range(sbml_model.getNumFunctionDefinitions()): function = sbml_model.getFunctionDefinition(i) @@ -330,15 +316,17 @@ def __get_function_definitions(sbml_model, gillespy_model): num_nodes = function_tree.getNumChildren() function_args = [function_tree.getChild(i).getName() for i in range(num_nodes-1)] function_string = __get_math(function_tree.getChild(num_nodes-1)) - gillespy_function = gillespy2.FunctionDefinition(name=function_name, function=function_string, args=function_args) + gillespy_function = gillespy2.FunctionDefinition( + name=function_name, function=function_string, args=function_args + ) gillespy_model.add_function_definition(gillespy_function) - init_state[gillespy_function.name] = gillespy_function.function + init_state[gillespy_function.name] = gillespy_function.function_string def __get_events(sbml_model, gillespy_model): for i in range(sbml_model.getNumEvents()): event = sbml_model.getEvent(i) gillespy_assignments = [] - + trigger = event.getTrigger() delay = event.getDelay() if delay is not None: @@ -348,20 +336,23 @@ def __get_events(sbml_model, gillespy_model): initial_value = trigger.getInitialValue() persistent = trigger.getPersistent() use_values_from_trigger_time = event.getUseValuesFromTriggerTime() - gillespy_trigger = gillespy2.EventTrigger(expression=expression, - initial_value=initial_value, persistent=persistent) + gillespy_trigger = gillespy2.EventTrigger( + expression=expression, initial_value=initial_value, persistent=persistent + ) assignments = event.getListOfEventAssignments() - for a in assignments: + for assign in assignments: # Convert Non-Constant Parameter to Species - if a.getVariable() in gillespy_model.listOfParameters: - gillespy_species = gillespy2.Species(name=a.getVariable(), - initial_value=gillespy_model.listOfParameters[a.getVariable()].expression, - mode='continuous', allow_negative_populations=True) - gillespy_model.delete_parameter(a.getVariable()) + if assign.getVariable() in gillespy_model.listOfParameters: + gillespy_species = gillespy2.Species( + name=assign.getVariable(), + initial_value=gillespy_model.listOfParameters[assign.getVariable()].expression, + mode='continuous', allow_negative_populations=True + ) + gillespy_model.delete_parameter(assign.getVariable()) gillespy_model.add_species([gillespy_species]) - gillespy_assignment = gillespy2.EventAssignment(variable=a.getVariable(), - expression=__get_math(a.getMath())) + gillespy_assignment = gillespy2.EventAssignment(variable=assign.getVariable(), + expression=__get_math(assign.getMath())) gillespy_assignments.append(gillespy_assignment) gillespy_event = gillespy2.Event( name=event.getId(), trigger=gillespy_trigger, @@ -372,9 +363,9 @@ def __get_events(sbml_model, gillespy_model): def __get_initial_assignments(sbml_model, gillespy_model): for i in range(sbml_model.getNumInitialAssignments()): - ia = sbml_model.getInitialAssignment(i) - variable = ia.getId() - expression = __get_math(ia.getMath()) + init_assign = sbml_model.getInitialAssignment(i) + variable = init_assign.getId() + expression = __get_math(init_assign.getMath()) assigned_value = eval(expression, {**init_state, **eval_globals}) init_state[variable] = assigned_value if assigned_value != assigned_value: @@ -391,8 +382,10 @@ def __resolve_evals(gillespy_model, init_state): successful = [] if len(postponed_evals): for var, expr in postponed_evals.items(): - try: assigned_value = eval(expr, {**eval_globals, **init_state}) - except: assigned_value = np.nan + try: + assigned_value = eval(expr, {**eval_globals, **init_state}) + except Exception: + assigned_value = np.nan if assigned_value == assigned_value: successful.append(var) init_state[var] = assigned_value @@ -400,13 +393,14 @@ def __resolve_evals(gillespy_model, init_state): gillespy_model.listOfSpecies[var].initial_value = assigned_value elif var in gillespy_model.listOfParameters: gillespy_model.listOfParameters[var].value = assigned_value - if not len(successful): break - for var in successful: del postponed_evals[var] + if len(successful) <= 0: + break + for var in successful: + del postponed_evals[var] def convert(filename, model_name=None, gillespy_model=None, report_silently_with_sbml_error=False): - sbml_model, errors = __read_sbml_model(filename) - + if sbml_model is None: if report_silently_with_sbml_error: return None, errors @@ -414,9 +408,10 @@ def convert(filename, model_name=None, gillespy_model=None, report_silently_with raise SBMLError(f"SBML model import failed. Reason Given: \n\t{errs}") if len(errors) > 0 and not report_silently_with_sbml_error: - from gillespy2 import log + from gillespy2 import log # pylint: disable=import-outside-toplevel errs = '\n\t'.join(errors) - log.warning(f"Error were detected in the SBML model. Error: \n\t{errs}") + errmsg = f"Error were detected in the SBML model. Error: \n\t{errs}" + log.warning(errmsg) if model_name is None: model_name = sbml_model.getName() @@ -430,13 +425,9 @@ def convert(filename, model_name=None, gillespy_model=None, report_silently_with __get_compartments(sbml_model, gillespy_model) __get_reactions(sbml_model, gillespy_model, errors) __get_rules(sbml_model, gillespy_model, errors) - __get_constraints(sbml_model, gillespy_model) + __get_constraints(sbml_model, gillespy_model, errors) __get_events(sbml_model, gillespy_model) __get_initial_assignments(sbml_model, gillespy_model) __resolve_evals(gillespy_model, init_state) - return gillespy_model, errors - - - From 5ab594164b6eb1e6ecfbe18a771562927e243749 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 8 Dec 2022 17:04:12 -0500 Subject: [PATCH 21/47] Fixed live graphing imports. --- gillespy2/core/liveGraphing.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gillespy2/core/liveGraphing.py b/gillespy2/core/liveGraphing.py index 3622fae35..1b6690bf4 100644 --- a/gillespy2/core/liveGraphing.py +++ b/gillespy2/core/liveGraphing.py @@ -19,7 +19,6 @@ import threading from math import floor -from IPython.display import clear_output from gillespy2.core.results import common_rgb_values from gillespy2.core.gillespyError import SimulationError @@ -178,6 +177,8 @@ def display(self, curr_state, curr_time, trajectory_base, finished=False): :param finished: Indicates whether or not the simulation has finished. :type finished: bool ''' + from IPython.display import clear_output # pylint: disable=import-outside-toplevel + curr_time = curr_time[0] curr_state = curr_state[0] From 510001180d20ae906a784f9c80a6a5657e726ed5 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 9 Dec 2022 12:42:59 -0800 Subject: [PATCH 22/47] Updating error message --- gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 2 +- gillespy2/solvers/numpy/tau_hybrid_solver.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index a1401ecd0..9cb5023d9 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -189,7 +189,7 @@ def _handle_return_code(self, return_code: "int") -> "int": if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_NO_SSA_REACTION: raise ExecutionError("Negative State detected in step, and no reaction found to fire.") if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_AT_BEGINING_OF_STEP: - raise ExecutionError("Negative State detected at beginning of step.") + raise ExecutionError("Negative State detected at beginning of step. Species involved in reactions can not be negative.") return super()._handle_return_code(return_code) @classmethod diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 2ca7439fe..be88a1204 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -603,7 +603,7 @@ def __simulate_negative_state_check(self, curr_state): # check each species to see if they are negative for s in self.non_negative_species: if curr_state[s] < 0: - raise SimulationError(f"Negative State detected at begining of step.") + raise SimulationError(f"Negative State detected at begining of step. Species involved in reactions can not be negative.") def __simulate_invalid_state_check(self, species_modified, curr_state, compiled_reactions): invalid_state = False From fe5e843fc758e9de24201c03d894bf2710130dff Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 9 Dec 2022 13:16:36 -0800 Subject: [PATCH 23/47] solver should not hide error --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index be88a1204..04bb634a1 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -1115,9 +1115,7 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None, except: pass if hasattr(self, 'has_raised_exception'): - raise SimulationError( - f"Error encountered while running simulation:\nReturn code: {int(self.rc)}.\n" - ) from self.has_raised_exception + raise self.has_raised_exception return Results.build_from_solver_results(self, live_output_options) From 9c61b1c3af4261ae9b9541d59fcc63810b343aa5 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 9 Dec 2022 16:22:13 -0800 Subject: [PATCH 24/47] Unit tests for hybrid solvers negative state behavior --- test/test_hybrid_negative_state.py | 72 ++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 test/test_hybrid_negative_state.py diff --git a/test/test_hybrid_negative_state.py b/test/test_hybrid_negative_state.py new file mode 100644 index 000000000..da9d4463e --- /dev/null +++ b/test/test_hybrid_negative_state.py @@ -0,0 +1,72 @@ +# GillesPy2 is a modeling toolkit for biochemical simulation. +# Copyright (C) 2019-2022 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 unittest +import numpy +import gillespy2 + +class TestNegativeState(unittest.TestCase): + def setUp(self): + + def create_hybrid_model(include_reaction=False): + model = gillespy2.Model(name="test_negative_state") + # Parameters + k1 = gillespy2.Parameter(name="k1", expression="1e-05") + model.add_parameter([k1]) + # Variables + A = gillespy2.Species(name="A", initial_value=10, mode="discrete") + model.add_species([A]) + # Reactions + if include_reaction: + model.add_reaction(gillespy2.Reaction(name="r1", reactants={}, + products={'A': 1}, rate="k1")) + # Rate Rule + model.add_rate_rule(gillespy2.RateRule(name="rr1", formula="-1", + variable="A")) + # Timespan + model.timespan(gillespy2.TimeSpan.arange(1,t=20)) + return model + + self.rr_only_model = create_hybrid_model(include_reaction=False) + self.rr_rxn_model = create_hybrid_model(include_reaction=True) + + def test_hybrid_python_rr_only(self): + sol = gillespy2.TauHybridSolver(model=self.rr_only_model) + result = sol.run() + self.assertLess(result['A'][-1], -9.9) + + def test_hybrid_python_rr_rxn(self): + sol = gillespy2.TauHybridSolver(model=self.rr_rxn_model) + with self.assertRaises(gillespy2.SimulationError) as ex: + result = sol.run() + self.assertIn('Negative State detected at begining of step. Species involved in reactions can not be negative.', str(ex.exception)) + + + def test_hybrid_C_rr_only(self): + sol = gillespy2.TauHybridCSolver(model=self.rr_only_model) + result = sol.run() + self.assertLess(result['A'][-1], -9.9) + + def test_hybrid_C_rr_rxn(self): + sol = gillespy2.TauHybridCSolver(model=self.rr_rxn_model) + with self.assertRaises(gillespy2.SimulationError) as ex: + result = sol.run() + self.assertIn('Negative State detected at begining of step. Species involved in reactions can not be negative.', str(ex.exception)) + + +if __name__ == '__main__': + unittest.main() + From e6f52f2ec68786940dcfd00f76a3bf79bb83fa71 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 9 Dec 2022 16:27:35 -0800 Subject: [PATCH 25/47] activating tests --- test/run_integration_tests.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/run_integration_tests.py b/test/run_integration_tests.py index 00b57fd36..a03db5046 100644 --- a/test/run_integration_tests.py +++ b/test/run_integration_tests.py @@ -58,6 +58,7 @@ import test_jsonify import test_notebooks import test_c_decode + import test_hybrid_negative_state modules = [ test_empty_model, @@ -85,6 +86,7 @@ test_jsonify, test_notebooks, test_c_decode, + test_hybrid_negative_state, ] for module in modules: From d5464724aa32d36e86e085c6bf1b1a1781a2f8ab Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 9 Dec 2022 17:29:51 -0800 Subject: [PATCH 26/47] Fixes #895 --- .../c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp | 7 +------ gillespy2/solvers/cpp/tau_hybrid_c_solver.py | 12 ++++++------ test/test_hybrid_negative_state.py | 2 +- 3 files changed, 8 insertions(+), 13 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 6503dd867..bc5ff07f2 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 @@ -391,12 +391,7 @@ namespace Gillespy { if (!simulation->species_state[p_i].boundary_condition) { - HybridSpecies *spec = &simulation->species_state[p_i]; - if( spec->partition_mode == SimulationState::CONTINUOUS ){ - result.concentrations[p_i] = result.concentrations[p_i] + population_changes[p_i]; - }else if( spec->partition_mode == SimulationState::DISCRETE ){ - result.concentrations[p_i] = current_state[p_i] + population_changes[p_i]; - } + result.concentrations[p_i] += population_changes[p_i]; } } // ===== ===== diff --git a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py index 9cb5023d9..bd06d1b29 100644 --- a/gillespy2/solvers/cpp/tau_hybrid_c_solver.py +++ b/gillespy2/solvers/cpp/tau_hybrid_c_solver.py @@ -179,17 +179,17 @@ def _build(self, model: "Union[Model, SanitizedModel]", simulation_name: str, va def _handle_return_code(self, return_code: "int") -> "int": if return_code == TauHybridCSolver.ErrorStatus.UNKNOWN: - raise ExecutionError("C++ solver failed (no error code given).") + raise SimulationError("C++ solver failed (no error code given).") if return_code == TauHybridCSolver.ErrorStatus.LOOP_OVER_INTEGRATE: - raise ExecutionError("Loop over integrate exceeded, problem space is too stiff") + raise SimulationError("Loop over integrate exceeded, problem space is too stiff") if return_code == TauHybridCSolver.ErrorStatus.INTEGRATOR_FAILED: - raise ExecutionError("Sundials ODE solver failed with 'BAD_STEP_SIZE'") + raise SimulationError("Sundials ODE solver failed with 'BAD_STEP_SIZE'") if return_code == TauHybridCSolver.ErrorStatus.INVALID_AFTER_SSA: - raise ExecutionError("Invalid state after single SSA step") + raise SimulationError("Invalid state after single SSA step") if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_NO_SSA_REACTION: - raise ExecutionError("Negative State detected in step, and no reaction found to fire.") + raise SimulationError("Negative State detected in step, and no reaction found to fire.") if return_code == TauHybridCSolver.ErrorStatus.NEGATIVE_STATE_AT_BEGINING_OF_STEP: - raise ExecutionError("Negative State detected at beginning of step. Species involved in reactions can not be negative.") + raise SimulationError("Negative State detected at beginning of step. Species involved in reactions can not be negative.") return super()._handle_return_code(return_code) @classmethod diff --git a/test/test_hybrid_negative_state.py b/test/test_hybrid_negative_state.py index da9d4463e..508515553 100644 --- a/test/test_hybrid_negative_state.py +++ b/test/test_hybrid_negative_state.py @@ -64,7 +64,7 @@ def test_hybrid_C_rr_rxn(self): sol = gillespy2.TauHybridCSolver(model=self.rr_rxn_model) with self.assertRaises(gillespy2.SimulationError) as ex: result = sol.run() - self.assertIn('Negative State detected at begining of step. Species involved in reactions can not be negative.', str(ex.exception)) + self.assertEqual('Negative State detected at beginning of step. Species involved in reactions can not be negative.', str(ex.exception)) if __name__ == '__main__': From b3a04818f366886e2fb0aae8af6eedea4ddbc0fc Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 10:00:59 -0800 Subject: [PATCH 27/47] Updating tests to ensure proper simulation behavior --- test/run_integration_tests.py | 2 ++ test/test_c_decode.py | 17 --------- test/test_run_output.py | 67 +++++++++++++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 17 deletions(-) create mode 100644 test/test_run_output.py diff --git a/test/run_integration_tests.py b/test/run_integration_tests.py index a03db5046..da981b10c 100644 --- a/test/run_integration_tests.py +++ b/test/run_integration_tests.py @@ -59,6 +59,7 @@ import test_notebooks import test_c_decode import test_hybrid_negative_state + import test_run_output modules = [ test_empty_model, @@ -86,6 +87,7 @@ test_jsonify, test_notebooks, test_c_decode, + test_run_output, test_hybrid_negative_state, ] diff --git a/test/test_c_decode.py b/test/test_c_decode.py index 76fce359c..46c0037aa 100644 --- a/test/test_c_decode.py +++ b/test/test_c_decode.py @@ -43,23 +43,6 @@ def setUp(self): (TauHybridCSolver(model=self.model), [1, 3]), ] - def test_run_output(self): - expected_tspan = self.model.tspan - expected_init = self.model.get_species("A").initial_value - for solver, trajectory_counts in self.solvers: - for number_of_trajectories in trajectory_counts: - with self.subTest("Processing simulation output for each solver with different trajectory sizes", - number_of_trajectories=number_of_trajectories, - solver=solver): - results = self.model.run(solver=solver, number_of_trajectories=number_of_trajectories, seed=1024) - for result in results: - tspan, first_value, last_value = result["time"], result["A"][0], result["A"][-1] - result_diff = np.concatenate([np.array([first_value]), result["A"][1:]]) - result_diff = result["A"] - result_diff - self.assertTrue(np.allclose(tspan, expected_tspan), msg="C++ Simulation output contains unexpected timeline values" - f"\n Received timeline: {tspan}\n Expected timeline: {expected_tspan}") - self.assertEqual(first_value, expected_init, msg=f"C++ Simulation output begins with unexpected value: {first_value}") - self.assertAlmostEqual(last_value, 0, places=3, msg=f"C++ Simulation output converges on an unexpectedly large value: {last_value}") def test_c_decoder(self): """ diff --git a/test/test_run_output.py b/test/test_run_output.py new file mode 100644 index 000000000..7a52ec46f --- /dev/null +++ b/test/test_run_output.py @@ -0,0 +1,67 @@ +# GillesPy2 is a modeling toolkit for biochemical simulation. +# Copyright (C) 2019-2022 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 sys + +sys.path.append(os.path.join(os.path.dirname(__file__), '..')) + +from datetime import time +import unittest +import numpy as np +import io +import gillespy2 +from gillespy2.solvers import ODESolver, ODECSolver, NumPySSASolver, SSACSolver, TauLeapingCSolver, TauLeapingSolver, TauHybridSolver, TauHybridCSolver +from example_models import create_degradation +from gillespy2.solvers.cpp.c_decoder import BasicSimDecoder, IterativeSimDecoder + +class TestRunOutput(unittest.TestCase): + decoders = [BasicSimDecoder, IterativeSimDecoder] + + def setUp(self): + self.model = create_degradation() + tspan = gillespy2.TimeSpan.linspace(t=300, num_points=301) + self.model.timespan(tspan) + self.solvers = [ + # Solver and List of trajectory numbers to try + (ODECSolver(model=self.model), [1]), + (SSACSolver(model=self.model), [1, 3]), + (TauLeapingCSolver(model=self.model), [1, 3]), + (TauHybridCSolver(model=self.model), [1, 3]), + (ODESolver(model=self.model), [1]), + (NumPySSASolver(model=self.model), [1, 3]), + (TauLeapingSolver(model=self.model), [1, 3]), + (TauHybridSolver(model=self.model), [1, 3]), + ] + + def test_run_output(self): + expected_tspan = self.model.tspan + expected_init = self.model.get_species("A").initial_value + for solver, trajectory_counts in self.solvers: + for number_of_trajectories in trajectory_counts: + with self.subTest("Processing simulation output for each solver with different trajectory sizes", + number_of_trajectories=number_of_trajectories, + solver=solver): + results = self.model.run(solver=solver, number_of_trajectories=number_of_trajectories, seed=1024) + for n,result in enumerate(results): + tspan, first_value, last_value = result["time"], result["A"][0], result["A"][-1] + result_diff = np.concatenate([np.array([first_value]), result["A"][1:]]) + result_diff = result["A"] - result_diff + self.assertTrue(np.allclose(tspan, expected_tspan), msg="C++ Simulation output contains unexpected timeline values" + f"\n Received timeline: {tspan}\n Expected timeline: {expected_tspan}, trajectory_number={n}") + self.assertEqual(first_value, expected_init, msg=f"C++ Simulation output begins with unexpected value: {first_value}, trajectory_number={n}") + self.assertAlmostEqual(last_value, 0, places=3, msg=f"C++ Simulation output converges on an unexpectedly large value: {last_value}, trajectory_number={n}") + From 749017894ecf298c7e132773c45be13dcf41ff90 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 13:34:20 -0800 Subject: [PATCH 28/47] The state must be updated differently depending if the mode of the species --- .../cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp | 7 ++++++- test/test_run_output.py | 8 +++----- 2 files changed, 9 insertions(+), 6 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 bc5ff07f2..6503dd867 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 @@ -391,7 +391,12 @@ namespace Gillespy { if (!simulation->species_state[p_i].boundary_condition) { - result.concentrations[p_i] += population_changes[p_i]; + HybridSpecies *spec = &simulation->species_state[p_i]; + if( spec->partition_mode == SimulationState::CONTINUOUS ){ + result.concentrations[p_i] = result.concentrations[p_i] + population_changes[p_i]; + }else if( spec->partition_mode == SimulationState::DISCRETE ){ + result.concentrations[p_i] = current_state[p_i] + population_changes[p_i]; + } } } // ===== ===== diff --git a/test/test_run_output.py b/test/test_run_output.py index 7a52ec46f..b1280e657 100644 --- a/test/test_run_output.py +++ b/test/test_run_output.py @@ -26,10 +26,8 @@ import gillespy2 from gillespy2.solvers import ODESolver, ODECSolver, NumPySSASolver, SSACSolver, TauLeapingCSolver, TauLeapingSolver, TauHybridSolver, TauHybridCSolver from example_models import create_degradation -from gillespy2.solvers.cpp.c_decoder import BasicSimDecoder, IterativeSimDecoder class TestRunOutput(unittest.TestCase): - decoders = [BasicSimDecoder, IterativeSimDecoder] def setUp(self): self.model = create_degradation() @@ -60,8 +58,8 @@ def test_run_output(self): tspan, first_value, last_value = result["time"], result["A"][0], result["A"][-1] result_diff = np.concatenate([np.array([first_value]), result["A"][1:]]) result_diff = result["A"] - result_diff - self.assertTrue(np.allclose(tspan, expected_tspan), msg="C++ Simulation output contains unexpected timeline values" + self.assertTrue(np.allclose(tspan, expected_tspan), msg="Simulation output contains unexpected timeline values" f"\n Received timeline: {tspan}\n Expected timeline: {expected_tspan}, trajectory_number={n}") - self.assertEqual(first_value, expected_init, msg=f"C++ Simulation output begins with unexpected value: {first_value}, trajectory_number={n}") - self.assertAlmostEqual(last_value, 0, places=3, msg=f"C++ Simulation output converges on an unexpectedly large value: {last_value}, trajectory_number={n}") + self.assertEqual(first_value, expected_init, msg=f"Simulation output begins with unexpected value: {first_value}, trajectory_number={n}") + self.assertAlmostEqual(last_value, 0, places=3, msg=f"Simulation output converges on an unexpectedly large value: {last_value}, trajectory_number={n}") From 193188b3404577fc33f220c7bf754d511eaa1867 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 15:32:32 -0800 Subject: [PATCH 29/47] Fixes #896 --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 04bb634a1..e67cd633b 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -145,8 +145,8 @@ def __create_diff_eqs(self, comb, dependencies, rr_sets, rate_rules): # Initialize sample dict rr_vars = {} for n, rr in self.model.listOfRateRules.items(): - rr_vars[rr.variable] = n - for spec in self.model.listOfSpecies: + rr_vars[rr.variable.name] = n + for spec in self.model.listOfSpecies.keys(): if spec in rr_vars.keys(): diff_eqs[self.model.listOfSpecies[spec]] = self.model.listOfRateRules[rr_vars[spec]].formula else: From a0c5285565da9114fe85072e9690ef797f3db5f6 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 15:36:38 -0800 Subject: [PATCH 30/47] Fixes #897 --- gillespy2/core/model.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 498cb6c65..49f120b7b 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -749,6 +749,14 @@ def add_rate_rule(self, rate_rule): elif rate_rule.variable in rr_vars: raise ModelError(f"Duplicate variable in rate_rules: {rate_rule.variable}.") self._resolve_rule(rate_rule) + # check if the rate_rule's target's mode is continious + if rate_rule.variable.mode == 'discrete': + raise ModelError("RateRules can not target discrete species") + if rate_rule.variable.mode is None or rate_rule.variable.mode == 'dynamic': + rate_rule.variable.mode = 'continuous' + from gillespy2.core import log + log.warning(f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule {rate_rule.name}") + self.listOfRateRules[rate_rule.name] = rate_rule # Build the sanitized rate rule sanitized_rate_rule = RateRule(name=f'RR{len(self._listOfRateRules)}') From 1047499921c5be6b949eb3fbc37021a1cd456a18 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 16:46:45 -0800 Subject: [PATCH 31/47] Fixes #898 --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 64 ++++++++++++-------- 1 file changed, 39 insertions(+), 25 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index e67cd633b..28e1be81c 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -94,6 +94,41 @@ def __init__(self, model=None, profile_reactions=False): for key, value in model.listOfReactions[reaction].products.items(): self.non_negative_species.add(key.name) + def __save_state_to_output(self, curr_time, save_index, curr_state, species, + trajectory, save_times): + """ + Helper function to save the curr_state to the trajectory output + """ + + # Now update the step and trajectories for this step of the simulation. + # Here we make our final assignments for this step, and begin + # populating our results trajectory. + + num_saves = 0 + for time in save_times: + if time > curr_time: + break + # if a solution is given for it + trajectory_index = save_index + assignment_state = copy.deepcopy(curr_state) + for s,sname in enumerate(species): + # Get ODE Solutions + trajectory[trajectory_index][s + 1] = curr_state[sname] + # Update Assignment Rules for all processed time points + if len(self.model.listOfAssignmentRules): + # Copy ODE state for assignments + assignment_state[sname] = curr_state[sname] + assignment_state['t'] = time + for ar in self.model.listOfAssignmentRules.values(): + assignment_value = eval(ar.formula, {**eval_globals, **assignment_state}) + assignment_state[ar.variable] = assignment_value + trajectory[trajectory_index][species.index(ar.variable.name) + 1] = assignment_value + num_saves += 1 + save_index += 1 + save_times = save_times[num_saves:] # remove completed save times + return save_times, save_index + + def __toggle_reactions(self, all_compiled, deterministic_reactions, dependencies, curr_state, det_spec, rr_sets): """ @@ -641,6 +676,9 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, # first check if we have a valid state: self.__simulate_negative_state_check(curr_state) + if curr_time == 0.0: + # save state at beginning of simulation + save_times, save_index = self.__save_state_to_output(curr_time, save_index, curr_state, species, trajectory, save_times) event_queue = [] prev_y0 = copy.deepcopy(y0) @@ -748,31 +786,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, raise Exception(f"Negative State detected in step, after single SSA step.\n\n error_message={invalid_err_message}\n curr_time={curr_time}\n tau_step={tau_step}\n curr_state={curr_state}\n\nstarting_curr_state={starting_curr_state}\n\n starting_tau_step={starting_tau_step}\nspecies_modified={species_modified}\nrxn_count={rxn_count}\n propensities={propensities}\nrxn_selected={rxn_selected}\ncompiled_reactions={compiled_reactions}\ncurr_state_after={curr_state_after} \n propensities_after={propensities_after}\nstarting_propensities={starting_propensities}\nfirst_rxn_count={first_rxn_count}\n first_err_message={first_err_message}\n 2nd_tau_step={tau_step}\nfloored_propensities={floored_propensities} ") - # Now update the step and trajectories for this step of the simulation. - # Here we make our final assignments for this step, and begin - # populating our results trajectory. - num_saves = 0 - for time in save_times: - if time > curr_time: - break - # if a solution is given for it - trajectory_index = save_index - assignment_state = copy.deepcopy(curr_state) - for s in range(len(species)): - # Get ODE Solutions - trajectory[trajectory_index][s + 1] = sol.y[s] - # Update Assignment Rules for all processed time points - if len(self.model.listOfAssignmentRules): - # Copy ODE state for assignments - assignment_state[species[s]] = sol.y[s] - assignment_state['t'] = time - for ar in self.model.listOfAssignmentRules.values(): - assignment_value = eval(ar.formula, {**eval_globals, **assignment_state}) - assignment_state[ar.variable] = assignment_value - trajectory[trajectory_index][species.index(ar.variable.name) + 1] = assignment_value - num_saves += 1 - save_index += 1 - save_times = save_times[num_saves:] # remove completed save times + save_times, save_index = self.__save_state_to_output(curr_time, save_index, curr_state, species, trajectory, save_times) events_processed = self.__process_queued_events(event_queue, trigger_states, curr_state) From 7e80677d9f970e95e2ed290867e1719d3987ba1d Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 17:06:39 -0800 Subject: [PATCH 32/47] RateRules can target Parameter, which don't have a .mode --- gillespy2/core/model.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 49f120b7b..6e53f95eb 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -749,13 +749,14 @@ def add_rate_rule(self, rate_rule): elif rate_rule.variable in rr_vars: raise ModelError(f"Duplicate variable in rate_rules: {rate_rule.variable}.") self._resolve_rule(rate_rule) - # check if the rate_rule's target's mode is continious - if rate_rule.variable.mode == 'discrete': - raise ModelError("RateRules can not target discrete species") - if rate_rule.variable.mode is None or rate_rule.variable.mode == 'dynamic': - rate_rule.variable.mode = 'continuous' - from gillespy2.core import log - log.warning(f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule {rate_rule.name}") + if rate_rule.variable.name in self.listOfSpecies.keys(): + # check if the rate_rule's target's mode is continious + if rate_rule.variable.mode == 'discrete': + raise ModelError("RateRules can not target discrete species") + if rate_rule.variable.mode is None or rate_rule.variable.mode == 'dynamic': + rate_rule.variable.mode = 'continuous' + from gillespy2.core import log + log.warning(f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule {rate_rule.name}") self.listOfRateRules[rate_rule.name] = rate_rule # Build the sanitized rate rule From c9ec6ff6d474d030f78e6eb8d9af79d64919b943 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Tue, 13 Dec 2022 18:31:26 -0800 Subject: [PATCH 33/47] correcting test --- test/test_hybrid_negative_state.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_hybrid_negative_state.py b/test/test_hybrid_negative_state.py index 508515553..e12642334 100644 --- a/test/test_hybrid_negative_state.py +++ b/test/test_hybrid_negative_state.py @@ -27,7 +27,7 @@ def create_hybrid_model(include_reaction=False): k1 = gillespy2.Parameter(name="k1", expression="1e-05") model.add_parameter([k1]) # Variables - A = gillespy2.Species(name="A", initial_value=10, mode="discrete") + A = gillespy2.Species(name="A", initial_value=10) model.add_species([A]) # Reactions if include_reaction: From 234b520e5465e336b878b569d2e05355569f129c Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 15 Dec 2022 09:04:01 -0500 Subject: [PATCH 34/47] Added github action to run pylint checks on pull requests. --- .github/workflows/pylint_on_pull_request.yml | 75 ++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 .github/workflows/pylint_on_pull_request.yml diff --git a/.github/workflows/pylint_on_pull_request.yml b/.github/workflows/pylint_on_pull_request.yml new file mode 100644 index 000000000..d881a1689 --- /dev/null +++ b/.github/workflows/pylint_on_pull_request.yml @@ -0,0 +1,75 @@ +name: PyLint On Pull Request +on: [pull_request] +jobs: + build: + runs-on: ubuntu-20.04 + steps: + - name: Set Up Python + uses: actions/setup-python@v2 + - name: Install PyLint + run: pip install --upgrade pylint + - name: Checkout + uses: actions/checkout@v2 + with: + fetch-depth: 0 + - name: Checkout Head + run: git checkout $HEAD_REF + env: + HEAD_REF: ${{ github.event.pull_request.head.ref }} + - name: Checkout Base + run: git checkout $BASE_REF + env: + BASE_REF: ${{ github.event.pull_request.base.ref }} + - name: Get Base Lint Score + run: | + echo BASE_LINT=$(git diff --name-only --diff-filter=M $HEAD_REF | grep -E "\.py" | xargs pylint | grep -E -o "at [0-9.-]+" | grep -E -o [0-9.-]+) >> $GITHUB_ENV + env: + HEAD_REF: ${{ github.event.pull_request.head.ref }} + if: always() + - name: Set Base Lint to 0 + run: echo BASE_LINT=0 >> $GITHUB_ENV + if: env.BASE_LINT == '' + - name: Checkout Head + run: git checkout $HEAD_REF + env: + HEAD_REF: ${{ github.event.pull_request.head.ref }} + - name: Get Head Lint Score + run: | + echo HEAD_LINT=$(git diff --name-only --diff-filter=M $BASE_REF | grep -E "\.py" | xargs pylint | grep -E -o "at [0-9.-]+" | grep -E -o [0-9.-]+) >> $GITHUB_ENV + env: + BASE_REF: ${{ github.event.pull_request.base.ref }} + if: always() + - name: Set Head Lint to 0 + run: echo HEAD_LINT=0 >> $GITHUB_ENV + if: env.HEAD_LINT == '' + - name: Get Added Files Lint Score + run: | + echo ADDED_LINT=$(git diff --name-only --diff-filter=A $BASE_REF | grep -E "\.py" | xargs pylint | grep -E -o "at [0-9.-]+" | grep -E -o [0-9.-]+) >> $GITHUB_ENV + env: + BASE_REF: ${{ github.event.pull_request.base.ref }} + if: always() + - name: Get Delta + run: | + import os + base = float(os.environ['BASE_LINT']) + head = float(os.environ['HEAD_LINT']) + delta = head - base + os.popen(f"echo DELTA={round(delta, 2)} >> $GITHUB_ENV") + shell: python + - name: Display Results + run: | + echo "Lint of modified files in base:" + echo ${{ env.BASE_LINT }} + echo "Lint of modified files in head:" + echo ${{ env.HEAD_LINT }} + echo "Delta (+/-):" + echo ${{ env.DELTA }} + echo "Lint of files added by head:" + echo ${{ env.ADDED_LINT }} + if: always() + - name: Fail If Negative Delta + run: | + import os + if float(os.environ['HEAD_LINT']) < 9 and float(os.environ['DELTA']) < 0: + raise Exception("Head lint score < 9 and negative delta.") + shell: python From f2c73af2d694e78ad5bb05ed27e03f2184244d3d Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 16 Dec 2022 12:46:02 -0800 Subject: [PATCH 35/47] Creating test to ensure proper default behavior of event assignments with discrete/continuous species --- test/run_integration_tests.py | 2 + test/test_hybrid_event_round.py | 85 +++++++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+) create mode 100644 test/test_hybrid_event_round.py diff --git a/test/run_integration_tests.py b/test/run_integration_tests.py index 00b57fd36..d5ebe539e 100644 --- a/test/run_integration_tests.py +++ b/test/run_integration_tests.py @@ -58,6 +58,7 @@ import test_jsonify import test_notebooks import test_c_decode + import test_hybrid_event_round modules = [ test_empty_model, @@ -85,6 +86,7 @@ test_jsonify, test_notebooks, test_c_decode, + test_hybrid_event_round, ] for module in modules: diff --git a/test/test_hybrid_event_round.py b/test/test_hybrid_event_round.py new file mode 100644 index 000000000..6332817da --- /dev/null +++ b/test/test_hybrid_event_round.py @@ -0,0 +1,85 @@ +# GillesPy2 is a modeling toolkit for biochemical simulation. +# Copyright (C) 2019-2022 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 sys + +sys.path.append(os.path.join(os.path.dirname(__file__), '..')) + +from datetime import time +import unittest +import numpy as np +import io +import gillespy2 +from gillespy2 import Model, Species, Parameter, Reaction, Event, \ + EventTrigger, EventAssignment, RateRule, \ + AssignmentRule, FunctionDefinition, TimeSpan +from gillespy2.solvers import TauHybridCSolver, TauHybridSolver + + +class TestHybridEventRound(unittest.TestCase): + + def create_model(self, discrete=True): + model = Model(name="test_issue710_round_event_assign") + + # Parameters + k1 = Parameter(name="k1", expression="1e-05") + model.add_parameter([k1]) + + # Variables + if discrete: + A = Species(name="A", initial_value=5, mode='discrete') + else: + A = Species(name="A", initial_value=5, mode="continuous") + model.add_species([A]) + + # Reactions + r1 = Reaction(name="r1", reactants={}, products={'A': 1}, rate="k1") + model.add_reaction([r1]) + + # Event Triggers + e1_trig = EventTrigger(expression="t>10", initial_value=False, persistent=False) + + # Event Assignments + e1_assign_1 = EventAssignment(variable="A", expression="A + 0.5") + + # Events + model.add_event(Event(name="e1", trigger=e1_trig, assignments=[e1_assign_1], delay=None, priority="0", use_values_from_trigger_time=False)) + + # Timespan + tspan = TimeSpan(np.arange(0, 20.05, 0.05)) + model.timespan(tspan) + return model + + def setUp(self): + self.d_model = self.create_model(discrete=True) + self.c_model = self.create_model(discrete=False) + + def test_event_rounding(self): + number_of_trajectories = 3 + for model, expected_last in {self.d_model: 6, self.c_model: 5.5}.items(): + expected_init = model.get_species("A").initial_value + + for sname, sclass in {'TauHybridCSolver':TauHybridCSolver,'TauHybridSolver':TauHybridSolver}.items(): + with self.subTest(f"Checking event assignment rounding for {model.listOfSpecies['A'].mode} species with {sname} solver."): + solver = sclass(model=model) + results = solver.run(number_of_trajectories=number_of_trajectories) + for result in results: + first_value = result["A"][0] + last_value = result["A"][-1], + self.assertEqual(first_value, expected_init, msg=f"Simulation output begins with {first_value}, should be {expected_init}") + self.assertAlmostEqual(result["A"][-1], expected_last, places=3, msg=f"Simulation output ends with {last_value}, should be {expected_last}") + From b9825e2f3c57dded210c154fa461c631f1ebb6ed Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 16 Dec 2022 15:13:38 -0800 Subject: [PATCH 36/47] adding test for a dynamic species --- test/test_hybrid_event_round.py | 64 +++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/test/test_hybrid_event_round.py b/test/test_hybrid_event_round.py index 6332817da..fa480af04 100644 --- a/test/test_hybrid_event_round.py +++ b/test/test_hybrid_event_round.py @@ -64,10 +64,74 @@ def create_model(self, discrete=True): model.timespan(tspan) return model + def create_dynamic_model(self): + model = Model(name="test_issue710_round_event_assign_switch") + + # Parameters + k1 = Parameter(name="k1", expression="1e-05") + model.add_parameter([k1]) + + # Variables + A = Species(name="A", initial_value=100, mode='dynamic', + switch_min=50) + model.add_species([A]) + + # Reactions + r1 = Reaction(name="r1", reactants={}, products={'A': 1}, rate="k1") + model.add_reaction([r1]) + + # Event Triggers + t5_trig = EventTrigger(expression="t>5", initial_value=False, persistent=False) + t10_trig = EventTrigger(expression="t>10", initial_value=False, persistent=False) + t15_trig = EventTrigger(expression="t>15", initial_value=False, persistent=False) + + # Event Assignments + eassign_1 = EventAssignment(variable="A", expression="A + 0.33") + eassign_2 = EventAssignment(variable="A", expression="A - 90") + + # Events + model.add_event(Event(name="e1", trigger=t5_trig, + assignments=[eassign_1], + delay=None, priority="0", + use_values_from_trigger_time=False)) + model.add_event(Event(name="e2", trigger=t10_trig, + assignments=[eassign_2], + delay=None, priority="0", + use_values_from_trigger_time=False)) + model.add_event(Event(name="e3", trigger=t15_trig, + assignments=[eassign_1], + delay=None, priority="0", + use_values_from_trigger_time=False)) + + # Timespan + tspan = TimeSpan(np.arange(0, 20.05, 0.05)) + model.timespan(tspan) + return model + def setUp(self): self.d_model = self.create_model(discrete=True) self.c_model = self.create_model(discrete=False) + self.dyn_model = self.create_dynamic_model() + + def test_dynamic_event_rounding(self): + expected_times_values = { + 99: 100.0, + 101: 100.33, + 201: 10.0, + 301: 10.0, + } + number_of_trajectories = 3 + + for sname, sclass in {'TauHybridCSolver':TauHybridCSolver,'TauHybridSolver':TauHybridSolver}.items(): + with self.subTest(f"Checking event assignment rounding for dynamic species with {sname} solver."): + solver = sclass(model=self.dyn_model) + results = solver.run(number_of_trajectories=number_of_trajectories) + for result in results: + for time_value, expected_value in expected_times_values.items(): + self.assertAlmostEqual(result['A'][time_value], expected_value, places=2, msg=f"Simulation output at {time_value} is {result['A'][time_value]}, should be {expected_value}") + + def test_event_rounding(self): number_of_trajectories = 3 for model, expected_last in {self.d_model: 6, self.c_model: 5.5}.items(): From cb58bbb042bfd240aff9b5b2c86211464e724383 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 16 Dec 2022 16:12:27 -0800 Subject: [PATCH 37/47] removing pure_ode and pure_stochastic optimization fixes issue #709 --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 93 ++++++++------------ test/test_hybrid_event_round.py | 8 +- 2 files changed, 42 insertions(+), 59 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 5a4722514..a75676720 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -382,7 +382,7 @@ def __get_next_step(self, event_times, reaction_times, delayed_events, return next_step[curr_time], curr_time def __process_queued_events(self, event_queue, trigger_states, - curr_state): + curr_state, det_spec): """ Helper method which processes the events queue. Method is primarily for evaluating assignments at the designated state (trigger time or event @@ -483,7 +483,7 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, propensities, y_map, compiled_reactions, active_rr, event_queue, delayed_events, trigger_states, - event_sensitivity, tau_step, pure_ode): + event_sensitivity, tau_step ): """ Helper function to perform the ODE integration of one step. This method uses scipy.integrate.LSODA to get simulation data, and @@ -508,11 +508,7 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, tau_step = max(integrator_options['min_step'], tau_step) else: tau_step = max(1e-6, tau_step) - if pure_ode: - next_tau = curr_time+max_step_size - dense_output = True - else: - next_tau = curr_time + tau_step + next_tau = curr_time + tau_step curr_state['t'] = curr_time curr_state['time'] = curr_time @@ -607,7 +603,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, propensities, species, parameters, compiled_reactions, active_rr, y_map, trajectory, save_times, save_index, delayed_events, trigger_states, event_sensitivity, - tau_step, pure_ode, debug): + tau_step, debug, det_spec): """ Function to process simulation until next step, which can be a stochastic reaction firing, an event trigger or assignment, or end of @@ -661,8 +657,8 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, delayed_events, trigger_states, event_sensitivity, - tau_step, - pure_ode) + tau_step + ) species_modified,rxn_count = self.__update_stochastic_rxn_states(compiled_reactions, curr_state) @@ -721,8 +717,8 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, delayed_events, trigger_states, event_sensitivity, - tau_step, - pure_ode) + tau_step + ) # only update the selected reaction first_rxn_count = copy.deepcopy(rxn_count) @@ -760,7 +756,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, save_index += 1 save_times = save_times[num_saves:] # remove completed save times - events_processed = self.__process_queued_events(event_queue, trigger_states, curr_state) + events_processed = self.__process_queued_events(event_queue, trigger_states, curr_state, det_spec) # Finally, perform a final check on events after all non-ODE assignment # changes have been carried out on model. @@ -777,7 +773,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, elif not triggered: curr_state[e.name] = False - events_processed = self.__process_queued_events(event_queue, trigger_states, curr_state) + events_processed = self.__process_queued_events(event_queue, trigger_states, curr_state, det_spec) return sol, curr_state, curr_time, save_times, save_index @@ -1137,30 +1133,27 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, t0_delayed_events, species_modified_by_events = self.__check_t0_events(initial_state) # Create deterministic tracking data structures - det_spec = {species: True for (species, value) in self.model.listOfSpecies.items() if value.mode == 'dynamic'} + #det_spec = {species: True for (species, value) in self.model.listOfSpecies.items() if value.mode == 'dynamic'} + det_spec = {species: False for species in self.model.listOfSpecies.keys()} + for (species, value) in self.model.listOfSpecies.items(): + print(f"species={species} value.mode={value.mode}") + if value.mode == 'dynamic' or value.mode == 'continuous': + det_spec[species] = True + print(f"det_spec={det_spec} l1144") det_rxn = {rxn: False for (rxn, value) in self.model.listOfReactions.items()} # Determine if entire simulation is ODE or Stochastic, in order to # avoid unnecessary calculations during simulation - pure_ode = True - pure_stochastic = True - for spec in self.model.listOfSpecies.values(): - if spec.mode != 'discrete': - pure_stochastic = False - if spec.mode != 'continuous': - pure_ode = False - simulation_data = [] dependencies = OrderedDict() # If considering deterministic changes, create dependency data # structure for creating diff eqs later - if not pure_stochastic: - for reaction in self.model.listOfReactions: - dependencies[reaction] = set() - [dependencies[reaction].add(reactant.name) for reactant in self.model.listOfReactions[reaction].reactants] - [dependencies[reaction].add(product.name) for product in self.model.listOfReactions[reaction].products] + for reaction in self.model.listOfReactions: + dependencies[reaction] = set() + [dependencies[reaction].add(reactant.name) for reactant in self.model.listOfReactions[reaction].reactants] + [dependencies[reaction].add(product.name) for product in self.model.listOfReactions[reaction].products] # Main trajectory loop for trajectory_num in range(number_of_trajectories): @@ -1188,8 +1181,7 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, save_index = 0 # Record Highest Order reactant for each reaction and set error tolerance - if not pure_ode: - HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, critical_threshold = Tau.initialize(self.model, tau_tol) + HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, critical_threshold = Tau.initialize(self.model, tau_tol) # One-time compilations to reduce time spent with eval compiled_reactions, compiled_rate_rules, compiled_inactive_reactions, compiled_propensities = \ @@ -1221,31 +1213,25 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, self.rc = 33 break # Get current propensities - if not pure_ode: - for i, r in enumerate(self.model.listOfReactions): - try: - propensities[r] = eval(compiled_propensities[r], eval_globals, curr_state[0]) - if curr_state[0][r] > 0 and propensities[r]==0: - # This is an edge case, that might happen after a single SSA step. - curr_state[0][r] = math.log(random.uniform(0, 1)) - except Exception as e: - raise SimulationError('Error calculation propensity for {0}.\nReason: {1}\ncurr_state={2}'.format(r, e, curr_state)) + for i, r in enumerate(self.model.listOfReactions): + try: + propensities[r] = eval(compiled_propensities[r], eval_globals, curr_state[0]) + if curr_state[0][r] > 0 and propensities[r]==0: + # This is an edge case, that might happen after a single SSA step. + curr_state[0][r] = math.log(random.uniform(0, 1)) + except Exception as e: + raise SimulationError('Error calculation propensity for {0}.\nReason: {1}\ncurr_state={2}'.format(r, e, curr_state)) # Calculate Tau statistics and select a good tau step - if not pure_ode: - tau_args = [HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, tau_tol, critical_threshold, - self.model, propensities, curr_state[0], curr_time[0], save_times[0]] - tau_step = save_times[-1] - curr_time[0] if pure_ode else Tau.select(*tau_args) + tau_args = [HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, tau_tol, critical_threshold, + self.model, propensities, curr_state[0], curr_time[0], save_times[0]] + tau_step = Tau.select(*tau_args) # Process switching if used - if not pure_stochastic and not pure_ode: - mn, sd, CV = self.__calculate_statistics(curr_time[0], propensities, curr_state[0], tau_step, det_spec) + mn, sd, CV = self.__calculate_statistics(curr_time[0], propensities, curr_state[0], tau_step, det_spec) # Calculate sd and CV for hybrid switching and flag deterministic reactions - if pure_stochastic: - deterministic_reactions = frozenset() # Empty if non-det - else: - deterministic_reactions = self.__flag_det_reactions(det_spec, det_rxn, dependencies) + deterministic_reactions = self.__flag_det_reactions(det_spec, det_rxn, dependencies) if debug: print('mean: {0}'.format(mn)) @@ -1256,10 +1242,7 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, # Set active reactions and rate rules for this integration step rr_sets = {frozenset() : compiled_rate_rules} # base rr set - if pure_stochastic: - active_rr = rr_sets[frozenset()] - else: - active_rr = self.__toggle_reactions(all_compiled, deterministic_reactions, + active_rr = self.__toggle_reactions(all_compiled, deterministic_reactions, dependencies, curr_state[0], det_spec, rr_sets) # Create integration initial state vector @@ -1275,8 +1258,8 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, trajectory, save_times, save_index, delayed_events, trigger_states, - event_sensitivity, tau_step, pure_ode, - debug) + event_sensitivity, tau_step, + debug, det_spec) # End of trajectory, format results data = {'time': timeline} diff --git a/test/test_hybrid_event_round.py b/test/test_hybrid_event_round.py index fa480af04..7c3fc06fa 100644 --- a/test/test_hybrid_event_round.py +++ b/test/test_hybrid_event_round.py @@ -117,9 +117,9 @@ def setUp(self): def test_dynamic_event_rounding(self): expected_times_values = { 99: 100.0, - 101: 100.33, - 201: 10.0, - 301: 10.0, + 102: 100.33, + 202: 10.0, + 302: 10.0, } number_of_trajectories = 3 @@ -144,6 +144,6 @@ def test_event_rounding(self): for result in results: first_value = result["A"][0] last_value = result["A"][-1], - self.assertEqual(first_value, expected_init, msg=f"Simulation output begins with {first_value}, should be {expected_init}") + self.assertAlmostEqual(first_value, expected_init, places=3, msg=f"Simulation output begins with {first_value}, should be {expected_init}") self.assertAlmostEqual(result["A"][-1], expected_last, places=3, msg=f"Simulation output ends with {last_value}, should be {expected_last}") From e310053867e33a0bc0e54a47ebe03430ac0225c2 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Fri, 16 Dec 2022 16:30:51 -0800 Subject: [PATCH 38/47] bug fix --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index a75676720..31e8e3d6c 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -1135,11 +1135,9 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, # Create deterministic tracking data structures #det_spec = {species: True for (species, value) in self.model.listOfSpecies.items() if value.mode == 'dynamic'} det_spec = {species: False for species in self.model.listOfSpecies.keys()} - for (species, value) in self.model.listOfSpecies.items(): - print(f"species={species} value.mode={value.mode}") + for (sname, value) in self.model.listOfSpecies.items(): if value.mode == 'dynamic' or value.mode == 'continuous': - det_spec[species] = True - print(f"det_spec={det_spec} l1144") + det_spec[sname] = True det_rxn = {rxn: False for (rxn, value) in self.model.listOfReactions.items()} # Determine if entire simulation is ODE or Stochastic, in order to From 5fbdff4e14e808fbd5702e34f1a9192ffa7712f7 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 19 Dec 2022 18:43:00 -0800 Subject: [PATCH 39/47] Fixes #902 --- gillespy2/solvers/utilities/Tau.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/gillespy2/solvers/utilities/Tau.py b/gillespy2/solvers/utilities/Tau.py index 34b327d91..880da3c66 100644 --- a/gillespy2/solvers/utilities/Tau.py +++ b/gillespy2/solvers/utilities/Tau.py @@ -114,14 +114,12 @@ def select(*tau_args): for r in model.listOfReactions: # Calculate abs mean and standard deviation for each reactant for reactant in model.listOfReactions[r].reactants: - if reactant not in mu_i: - mu_i[reactant.name] = 0 - if reactant not in sigma_i: - sigma_i[reactant.name] = 0 - mu_i[reactant.name] += model.listOfReactions[r].reactants[reactant] * propensities[ - r] # Cao, Gillespie, Petzold 32a - sigma_i[reactant.name] += model.listOfReactions[r].reactants[reactant] ** 2 * propensities[ - r] # Cao, Gillespie, Petzold 32b + net_change = model.listOfReactions[r].reactants[reactant] + if reactant in model.listOfReactions[r].products: + net_change -= model.listOfReactions[r].products[reactant] + net_change = abs(net_change) + mu_i[reactant.name] += net_change * propensities[r] # Cao, Gillespie, Petzold 32a + sigma_i[reactant.name] += net_change ** 2 * propensities[r] # Cao, Gillespie, Petzold 32b for r in reactants: calculated_max = epsilon_i[r.name] * curr_state[r.name] From 6da0d7729774a14b4822be370b93be7f9dca7bc1 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Mon, 19 Dec 2022 19:07:30 -0800 Subject: [PATCH 40/47] adding tests --- test/test_hybrid_solver_opioid_model.py | 109 ++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 test/test_hybrid_solver_opioid_model.py diff --git a/test/test_hybrid_solver_opioid_model.py b/test/test_hybrid_solver_opioid_model.py new file mode 100644 index 000000000..d97dfad50 --- /dev/null +++ b/test/test_hybrid_solver_opioid_model.py @@ -0,0 +1,109 @@ +# GillesPy2 is a modeling toolkit for biochemical simulation. +# Copyright (C) 2019-2022 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 sys + +sys.path.append(os.path.join(os.path.dirname(__file__), '..')) + +from datetime import time +import unittest +import numpy as np +import io +import gillespy2 +from gillespy2.solvers import TauHybridCSolver, TauHybridSolver + +class TestHbridSolverOpioidModel(unittest.TestCase): + + def create_opioid(self): + # Initialize Model + model = gillespy2.Model(name="Opioid") + + # Define Variables (GillesPy2.Species) + S = gillespy2.Species(name='Susceptibles', initial_value=1000) + P = gillespy2.Species(name='Prescribed_Users', initial_value=0) + A = gillespy2.Species(name='Addicted', initial_value=0) + R = gillespy2.Species(name='Rehab', initial_value=0) + Natural_Deaths = gillespy2.Species(name='Natural_Deaths', initial_value=0) + Addiction_Deaths = gillespy2.Species(name='Addiction_Deaths', initial_value=0) + + # Add Variables to Model + model.add_species([S,P,A,R,Natural_Deaths,Addiction_Deaths]) + + # Define Parameters + alpha = gillespy2.Parameter(name='alpha', expression= 0.15) + epsilon = gillespy2.Parameter(name='epsilon', expression= 0.8) + beta_p = gillespy2.Parameter(name='beta_p', expression= 0.00266) + beta_a = gillespy2.Parameter(name='beta_a', expression= 0.00094) + gamma = gillespy2.Parameter(name='gamma', expression= 0.00744) + zeta = gillespy2.Parameter(name='zeta', expression= 0.2) + delta = gillespy2.Parameter(name='delta', expression= 0.1) + sigma = gillespy2.Parameter(name='sigma', expression= 0.9) + mu = gillespy2.Parameter(name='mu', expression= 0.00729) + mu_prime = gillespy2.Parameter(name='mu_prime', expression= 0.01159) + + # Add Parameters to Model + model.add_parameter([alpha, epsilon, beta_p, beta_a, gamma, zeta, delta, sigma, mu, mu_prime]) + + # Define Reactions + SP = gillespy2.Reaction( + name="SP", reactants={'Susceptibles': 1}, products={'Prescribed_Users': 1}, rate='alpha' + ) + SA_a = gillespy2.Reaction(name="SA_a", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_a') + SA_p = gillespy2.Reaction(name="SA_p", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_p') + mu_S = gillespy2.Reaction( + name="mu_S", reactants={'Susceptibles': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + PA = gillespy2.Reaction(name="PA", reactants={'Prescribed_Users': 1}, products={'Addicted': 1}, rate='gamma') + PS = gillespy2.Reaction( + name="PS", reactants={'Prescribed_Users': 1}, products={'Susceptibles': 1}, rate='epsilon' + ) + AR = gillespy2.Reaction(name="AR", reactants={'Addicted': 1}, products={'Rehab': 1}, rate='zeta') + RA = gillespy2.Reaction(name="RA", reactants={'Rehab': 1}, products={'Addicted': 1}, rate='delta') + RS = gillespy2.Reaction(name="RS", reactants={'Rehab': 1}, products={'Susceptibles': 1}, rate='sigma') + mu_P = gillespy2.Reaction( + name="mu_P", reactants={'Prescribed_Users': 1}, + products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + mu_R = gillespy2.Reaction( + name="m_R", reactants={'Rehab': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + mu_prime_A = gillespy2.Reaction( + name="mu_prime_A", reactants={'Addicted': 1}, + products={'Susceptibles': 1, 'Addiction_Deaths': 1}, rate='mu_prime' + ) + + # Add Reactions to Model + model.add_reaction([SP, PS, SA_a, SA_p, PA, AR, RA, RS, mu_P, mu_R, mu_prime_A, mu_S]) + + # Define Timespan + tspan = gillespy2.TimeSpan.linspace(t=10, num_points=11) + + # Set Model Timespan + model.timespan(tspan) + return model + + def test_run_output(self): + model = self.create_opioid() + for solver in [TauHybridSolver, TauHybridCSolver]: + results = model.run(solver=solver, number_of_trajectories=3, seed=1024) + for result in results: + with self.subTest("Processing simulation output for solver {solver.name} for trajectory={n}"): + min_s = min(result['Susceptibles']) + self.assertGreater(min_s, 500) + max_p = max(result['Prescribed_Users']) + self.assertLess(max_p, 500) + From 890800568fec6c3c8e459323814da734c633f1ed Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 09:10:49 -0800 Subject: [PATCH 41/47] Moving Opioid model to example_models file --- test/example_models.py | 74 ++++++++++++++++++++++++- test/test_hybrid_solver_opioid_model.py | 71 +----------------------- 2 files changed, 74 insertions(+), 71 deletions(-) diff --git a/test/example_models.py b/test/example_models.py index 7020e5910..f0055eb3b 100644 --- a/test/example_models.py +++ b/test/example_models.py @@ -13,7 +13,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 import ( Event, Model, @@ -564,9 +564,79 @@ def create_multi_firing_event(parameter_values=None): model.timespan(np.linspace(0, 60, 181)) return model + +def create_opioid(): + # Initialize Model + model = gillespy2.Model(name="Opioid") + + # Define Variables (GillesPy2.Species) + S = gillespy2.Species(name='Susceptibles', initial_value=1000) + P = gillespy2.Species(name='Prescribed_Users', initial_value=0) + A = gillespy2.Species(name='Addicted', initial_value=0) + R = gillespy2.Species(name='Rehab', initial_value=0) + Natural_Deaths = gillespy2.Species(name='Natural_Deaths', initial_value=0) + Addiction_Deaths = gillespy2.Species(name='Addiction_Deaths', initial_value=0) + + # Add Variables to Model + model.add_species([S,P,A,R,Natural_Deaths,Addiction_Deaths]) + + # Define Parameters + alpha = gillespy2.Parameter(name='alpha', expression= 0.15) + epsilon = gillespy2.Parameter(name='epsilon', expression= 0.8) + beta_p = gillespy2.Parameter(name='beta_p', expression= 0.00266) + beta_a = gillespy2.Parameter(name='beta_a', expression= 0.00094) + gamma = gillespy2.Parameter(name='gamma', expression= 0.00744) + zeta = gillespy2.Parameter(name='zeta', expression= 0.2) + delta = gillespy2.Parameter(name='delta', expression= 0.1) + sigma = gillespy2.Parameter(name='sigma', expression= 0.9) + mu = gillespy2.Parameter(name='mu', expression= 0.00729) + mu_prime = gillespy2.Parameter(name='mu_prime', expression= 0.01159) + + # Add Parameters to Model + model.add_parameter([alpha, epsilon, beta_p, beta_a, gamma, zeta, delta, sigma, mu, mu_prime]) + + # Define Reactions + SP = gillespy2.Reaction( + name="SP", reactants={'Susceptibles': 1}, products={'Prescribed_Users': 1}, rate='alpha' + ) + SA_a = gillespy2.Reaction(name="SA_a", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_a') + SA_p = gillespy2.Reaction(name="SA_p", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_p') + mu_S = gillespy2.Reaction( + name="mu_S", reactants={'Susceptibles': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + PA = gillespy2.Reaction(name="PA", reactants={'Prescribed_Users': 1}, products={'Addicted': 1}, rate='gamma') + PS = gillespy2.Reaction( + name="PS", reactants={'Prescribed_Users': 1}, products={'Susceptibles': 1}, rate='epsilon' + ) + AR = gillespy2.Reaction(name="AR", reactants={'Addicted': 1}, products={'Rehab': 1}, rate='zeta') + RA = gillespy2.Reaction(name="RA", reactants={'Rehab': 1}, products={'Addicted': 1}, rate='delta') + RS = gillespy2.Reaction(name="RS", reactants={'Rehab': 1}, products={'Susceptibles': 1}, rate='sigma') + mu_P = gillespy2.Reaction( + name="mu_P", reactants={'Prescribed_Users': 1}, + products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + mu_R = gillespy2.Reaction( + name="m_R", reactants={'Rehab': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' + ) + mu_prime_A = gillespy2.Reaction( + name="mu_prime_A", reactants={'Addicted': 1}, + products={'Susceptibles': 1, 'Addiction_Deaths': 1}, rate='mu_prime' + ) + + # Add Reactions to Model + model.add_reaction([SP, PS, SA_a, SA_p, PA, AR, RA, RS, mu_P, mu_R, mu_prime_A, mu_S]) + + # Define Timespan + tspan = gillespy2.TimeSpan.linspace(t=10, num_points=11) + + # Set Model Timespan + model.timespan(tspan) + return model + + __all__ = [ 'create_trichloroethylene', 'create_lac_operon', 'create_schlogl', 'create_michaelis_menten', 'create_decay', 'create_decay_no_tspan', 'create_tyson_2_state_oscillator', 'create_oregonator', 'create_vilar_oscillator', 'create_dimerization', 'create_degradation', 'create_robust_model', 'create_multi_firing_event', 'create_toggle_switch' -] \ No newline at end of file +] diff --git a/test/test_hybrid_solver_opioid_model.py b/test/test_hybrid_solver_opioid_model.py index d97dfad50..261cbb13c 100644 --- a/test/test_hybrid_solver_opioid_model.py +++ b/test/test_hybrid_solver_opioid_model.py @@ -25,79 +25,12 @@ import io import gillespy2 from gillespy2.solvers import TauHybridCSolver, TauHybridSolver +from example_models import create_opioid class TestHbridSolverOpioidModel(unittest.TestCase): - def create_opioid(self): - # Initialize Model - model = gillespy2.Model(name="Opioid") - - # Define Variables (GillesPy2.Species) - S = gillespy2.Species(name='Susceptibles', initial_value=1000) - P = gillespy2.Species(name='Prescribed_Users', initial_value=0) - A = gillespy2.Species(name='Addicted', initial_value=0) - R = gillespy2.Species(name='Rehab', initial_value=0) - Natural_Deaths = gillespy2.Species(name='Natural_Deaths', initial_value=0) - Addiction_Deaths = gillespy2.Species(name='Addiction_Deaths', initial_value=0) - - # Add Variables to Model - model.add_species([S,P,A,R,Natural_Deaths,Addiction_Deaths]) - - # Define Parameters - alpha = gillespy2.Parameter(name='alpha', expression= 0.15) - epsilon = gillespy2.Parameter(name='epsilon', expression= 0.8) - beta_p = gillespy2.Parameter(name='beta_p', expression= 0.00266) - beta_a = gillespy2.Parameter(name='beta_a', expression= 0.00094) - gamma = gillespy2.Parameter(name='gamma', expression= 0.00744) - zeta = gillespy2.Parameter(name='zeta', expression= 0.2) - delta = gillespy2.Parameter(name='delta', expression= 0.1) - sigma = gillespy2.Parameter(name='sigma', expression= 0.9) - mu = gillespy2.Parameter(name='mu', expression= 0.00729) - mu_prime = gillespy2.Parameter(name='mu_prime', expression= 0.01159) - - # Add Parameters to Model - model.add_parameter([alpha, epsilon, beta_p, beta_a, gamma, zeta, delta, sigma, mu, mu_prime]) - - # Define Reactions - SP = gillespy2.Reaction( - name="SP", reactants={'Susceptibles': 1}, products={'Prescribed_Users': 1}, rate='alpha' - ) - SA_a = gillespy2.Reaction(name="SA_a", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_a') - SA_p = gillespy2.Reaction(name="SA_p", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_p') - mu_S = gillespy2.Reaction( - name="mu_S", reactants={'Susceptibles': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' - ) - PA = gillespy2.Reaction(name="PA", reactants={'Prescribed_Users': 1}, products={'Addicted': 1}, rate='gamma') - PS = gillespy2.Reaction( - name="PS", reactants={'Prescribed_Users': 1}, products={'Susceptibles': 1}, rate='epsilon' - ) - AR = gillespy2.Reaction(name="AR", reactants={'Addicted': 1}, products={'Rehab': 1}, rate='zeta') - RA = gillespy2.Reaction(name="RA", reactants={'Rehab': 1}, products={'Addicted': 1}, rate='delta') - RS = gillespy2.Reaction(name="RS", reactants={'Rehab': 1}, products={'Susceptibles': 1}, rate='sigma') - mu_P = gillespy2.Reaction( - name="mu_P", reactants={'Prescribed_Users': 1}, - products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' - ) - mu_R = gillespy2.Reaction( - name="m_R", reactants={'Rehab': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu' - ) - mu_prime_A = gillespy2.Reaction( - name="mu_prime_A", reactants={'Addicted': 1}, - products={'Susceptibles': 1, 'Addiction_Deaths': 1}, rate='mu_prime' - ) - - # Add Reactions to Model - model.add_reaction([SP, PS, SA_a, SA_p, PA, AR, RA, RS, mu_P, mu_R, mu_prime_A, mu_S]) - - # Define Timespan - tspan = gillespy2.TimeSpan.linspace(t=10, num_points=11) - - # Set Model Timespan - model.timespan(tspan) - return model - def test_run_output(self): - model = self.create_opioid() + model = create_opioid() for solver in [TauHybridSolver, TauHybridCSolver]: results = model.run(solver=solver, number_of_trajectories=3, seed=1024) for result in results: From b9ea2ed912e0819bc40428c33f649b7564522354 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 10:36:37 -0800 Subject: [PATCH 42/47] Addressing review comments --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 31e8e3d6c..95c7f83e8 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -1133,11 +1133,8 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, t0_delayed_events, species_modified_by_events = self.__check_t0_events(initial_state) # Create deterministic tracking data structures - #det_spec = {species: True for (species, value) in self.model.listOfSpecies.items() if value.mode == 'dynamic'} - det_spec = {species: False for species in self.model.listOfSpecies.keys()} - for (sname, value) in self.model.listOfSpecies.items(): - if value.mode == 'dynamic' or value.mode == 'continuous': - det_spec[sname] = True + det_spec = {species: value.mode != 'discrete' for (species, value) in self.model.listOfSpecies.items()} + det_rxn = {rxn: False for (rxn, value) in self.model.listOfReactions.items()} # Determine if entire simulation is ODE or Stochastic, in order to @@ -1221,9 +1218,7 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, raise SimulationError('Error calculation propensity for {0}.\nReason: {1}\ncurr_state={2}'.format(r, e, curr_state)) # Calculate Tau statistics and select a good tau step - tau_args = [HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, tau_tol, critical_threshold, - self.model, propensities, curr_state[0], curr_time[0], save_times[0]] - tau_step = Tau.select(*tau_args) + tau_step = Tau.select(HOR, reactants, mu_i, sigma_i, g_i, epsilon_i, tau_tol, critical_threshold, self.model, propensities, curr_state[0], curr_time[0], save_times[0]) # Process switching if used mn, sd, CV = self.__calculate_statistics(curr_time[0], propensities, curr_state[0], tau_step, det_spec) From 5d356489c9ff7daaa22d2b7618e386d0b2a832f6 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 10:48:44 -0800 Subject: [PATCH 43/47] Fixingexceptions --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 95c7f83e8..2d3ebf254 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -705,7 +705,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, if min_tau is None or min_tau > rxn_times[rname]: min_tau = rxn_times[rname] rxn_selected = rname - if rxn_selected is None: raise Exception(f"Negative State detected in step, and no reaction found to fire.\n\n error_message={invalid_err_message}\n curr_time={curr_time}\n tau_step={tau_step}\n curr_state={curr_state}\n\nstarting_curr_state={starting_curr_state}\n\n starting_tau_step={starting_tau_step}\nspecies_modified={species_modified}\nrxn_count={rxn_count}\n propensities={propensities}\nrxn_times={rxn_times}\ncompiled_reactions={compiled_reactions}\ncurr_state_after={curr_state_after}\n propensities_after={propensities_after}\nstarting_propensities={starting_propensities}\nfloored_curr_state={floored_curr_state}\nfloored_propensities={floored_propensities}\n ") + if rxn_selected is None: raise SimulationError(f"Negative State detected in step, and no reaction found to fire. error_message={invalid_err_message}") tau_step = min_tau #estimated time to the first stochatic reaction @@ -727,7 +727,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, (invalid_state, invalid_err_message) = self.__simulate_invalid_state_check(species_modified, curr_state, compiled_reactions) if invalid_state: - raise Exception(f"Negative State detected in step, after single SSA step.\n\n error_message={invalid_err_message}\n curr_time={curr_time}\n tau_step={tau_step}\n curr_state={curr_state}\n\nstarting_curr_state={starting_curr_state}\n\n starting_tau_step={starting_tau_step}\nspecies_modified={species_modified}\nrxn_count={rxn_count}\n propensities={propensities}\nrxn_selected={rxn_selected}\ncompiled_reactions={compiled_reactions}\ncurr_state_after={curr_state_after} \n propensities_after={propensities_after}\nstarting_propensities={starting_propensities}\nfirst_rxn_count={first_rxn_count}\n first_err_message={first_err_message}\n 2nd_tau_step={tau_step}\nfloored_propensities={floored_propensities} ") + raise SimulationError(f"Negative State detected in step, after single SSA step. error_message={invalid_err_message}") # Now update the step and trajectories for this step of the simulation. From bf3f15636eef3605b5a7789ca354e4d9dc2a614f Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 10:57:33 -0800 Subject: [PATCH 44/47] Pylint fixes --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 56 ++++++++++---------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 2d3ebf254..796d46192 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -88,7 +88,7 @@ def __init__(self, model=None, profile_reactions=False): for k in list(self.model.listOfSpecies)+list(self.model.listOfReactions): self.profile_data[k] = [] - def __toggle_reactions(self, all_compiled, deterministic_reactions, dependencies, + def __toggle_reactions(self, all_compiled, deterministic_reactions, dependencies, curr_state, det_spec, rr_sets): """ Helper method which is used to convert reaction channels into @@ -134,7 +134,7 @@ def __create_diff_eqs(self, comb, dependencies, rr_sets, rate_rules): differential equations, used dynamically throught the simulation. """ diff_eqs = OrderedDict() - output_rules = copy.deepcopy(rate_rules); + output_rules = copy.deepcopy(rate_rules) # Initialize sample dict rr_vars = {} @@ -161,10 +161,10 @@ def __create_diff_eqs(self, comb, dependencies, rr_sets, rate_rules): if factor[dep] != 0: if self.model.listOfSpecies[dep].mode == 'continuous': diff_eqs[self.model.listOfSpecies[dep]] += ' + {0}*({1})'.format(factor[dep], - self.model.listOfReactions[reaction].ode_propensity_function) + self.model.listOfReactions[reaction].ode_propensity_function) else: diff_eqs[self.model.listOfSpecies[dep]] += ' + {0}*({1})'.format(factor[dep], - self.model.listOfReactions[reaction].propensity_function) + self.model.listOfReactions[reaction].propensity_function) for spec in self.model.listOfSpecies: if diff_eqs[self.model.listOfSpecies[spec]] == '0': @@ -187,7 +187,7 @@ def __flag_det_reactions(self, det_spec, det_rxn, dependencies): det_rxn[rxn] = True # iterate through the dependent species of this reaction for species in dependencies[rxn]: - # if any of the dependencies are discrete or (dynamic AND the + # if any of the dependencies are discrete or (dynamic AND the # species itself has not been flagged as deterministic) # then allow it to be modelled discretely if self.model.listOfSpecies[species].mode == 'discrete': @@ -213,7 +213,7 @@ def __calculate_statistics(self, curr_time, propensities, curr_state, tau_step, NOTE: the argument cv_history should not be passed in, this is modified by the function to keep a persistent data set. """ - #TODO: move configuration to solver init + #move configuration to solver init (issue #905) history_length = 12 #cite: "Statistical rules of thumb" G Van Belle if curr_time==0.0: #re-set cv_history @@ -240,7 +240,7 @@ def __calculate_statistics(self, curr_time, propensities, curr_state, tau_step, for species,value in self.model.listOfSpecies.items(): if value.mode == 'dynamic': if mn[species] > 0 and sd[species] > 0: - CV[species] = math.sqrt(sd[species]) / mn[species] + CV[species] = math.sqrt(sd[species]) / mn[species] else: CV[species] = 1 # value chosen to guarantee species will be discrete @@ -293,7 +293,7 @@ def __f(t, y, curr_state, species, reactions, rate_rules, propensities, state_change[y_map[r]] += propensities[r] for event in events: triggered = eval(event.trigger.expression, {**eval_globals, **curr_state}) - if triggered: + if triggered: state_change[y_map[event]] = 1 return state_change @@ -484,18 +484,18 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, active_rr, event_queue, delayed_events, trigger_states, event_sensitivity, tau_step ): - """ + """ Helper function to perform the ODE integration of one step. This method uses scipy.integrate.LSODA to get simulation data, and determines the next stopping point of the simulation. The state is updated and returned to __simulate along with curr_time and the - solution object. + solution object. """ max_step_size = self.model.tspan[1] - self.model.tspan[0] / 100 from functools import partial events = self.model.listOfEvents.values() - dense_output = False + dense_output = False int_args = [curr_state, self.model.listOfSpecies, self.model.listOfReactions, self.model.listOfRateRules, propensities, y_map, @@ -525,7 +525,7 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, # ODE processes. This will update all species whose mode is set to # 'continuous', as well as 'dynamic' mode species which have been # flagged as deterministic. - + for spec_name, species in self.model.listOfSpecies.items(): if not species.constant: curr_state[spec_name] = sol.y[y_map[spec_name]] @@ -590,14 +590,14 @@ def __integrate(self, integrator_options, curr_state, y0, curr_time, def __simulate_invalid_state_check(self, species_modified, curr_state, compiled_reactions): - invalid_state = False - err_message="" - # check each species to see if they are negative - for s in species_modified.keys(): - if curr_state[s] < 0: - invalid_state = True - err_message += f"'{s}' has negative state '{curr_state[s]}'" - return (invalid_state, err_message) + invalid_state = False + err_message="" + # check each species to see if they are negative + for s in species_modified.keys(): + if curr_state[s] < 0: + invalid_state = True + err_message += f"'{s}' has negative state '{curr_state[s]}'" + return (invalid_state, err_message) def __simulate(self, integrator_options, curr_state, y0, curr_time, propensities, species, parameters, compiled_reactions, @@ -630,7 +630,7 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, loop_count = 0 invalid_state = False - + starting_curr_state = copy.deepcopy(curr_state) starting_propensities = copy.deepcopy(propensities) starting_tau_step=tau_step @@ -664,10 +664,10 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, # Occasionally, a tau step can result in an overly-aggressive # forward step and cause a species population to fall below 0, - # which would result in an erroneous simulation. - # We estimate the time to the first + # which would result in an erroneous simulation. + # We estimate the time to the first # stochatic reaction firing (assume constant propensities) and - # simulate the ODE system until that time, fire that reaction + # simulate the ODE system until that time, fire that reaction # and continue the simulation. (invalid_state, invalid_err_message) = self.__simulate_invalid_state_check(species_modified, curr_state, compiled_reactions) @@ -952,7 +952,7 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None, :param timeout: If set, if simulation takes longer than timeout, will exit. :type timeout: int - + :returns: A result object containing the results of the simulation. :rtype: gillespy2.Results """ @@ -1100,7 +1100,7 @@ def run(self=None, model=None, t=None, number_of_trajectories=1, increment=None, 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) def ___run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, live_grapher, t=20, @@ -1235,7 +1235,7 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, # Set active reactions and rate rules for this integration step rr_sets = {frozenset() : compiled_rate_rules} # base rr set - active_rr = self.__toggle_reactions(all_compiled, deterministic_reactions, + active_rr = self.__toggle_reactions(all_compiled, deterministic_reactions, dependencies, curr_state[0], det_spec, rr_sets) # Create integration initial state vector @@ -1251,7 +1251,7 @@ def __run(self, curr_state, curr_time, timeline, trajectory_base, initial_state, trajectory, save_times, save_index, delayed_events, trigger_states, - event_sensitivity, tau_step, + event_sensitivity, tau_step, debug, det_spec) # End of trajectory, format results From 7cda9161c8269958832e6d644ae9bb5f5d327fe2 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 11:13:17 -0800 Subject: [PATCH 45/47] Fixing pylint issues --- 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 2712e780b..9daadcf4a 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -746,13 +746,13 @@ def add_rate_rule(self, rate_rule): if rate_rule.variable in rr_vars: raise ModelError(f"Duplicate variable in rate_rules: {rate_rule.variable}.") self._resolve_rule(rate_rule) - if rate_rule.variable.name in self.listOfSpecies.keys(): + if rate_rule.variable.name in self.listOfSpecies: # check if the rate_rule's target's mode is continious if rate_rule.variable.mode == 'discrete': raise ModelError("RateRules can not target discrete species") if rate_rule.variable.mode is None or rate_rule.variable.mode == 'dynamic': rate_rule.variable.mode = 'continuous' - from gillespy2.core import log + from gillespy2.core import log # pylint:disable=import-outside-toplevel log.warning(f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule {rate_rule.name}") self.listOfRateRules[rate_rule.name] = rate_rule From 6b5955ff416195abc89023998c24e4023d24a8e5 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 11:18:49 -0800 Subject: [PATCH 46/47] more pylint --- gillespy2/core/model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gillespy2/core/model.py b/gillespy2/core/model.py index 9daadcf4a..9a188809a 100644 --- a/gillespy2/core/model.py +++ b/gillespy2/core/model.py @@ -753,7 +753,9 @@ def add_rate_rule(self, rate_rule): if rate_rule.variable.mode is None or rate_rule.variable.mode == 'dynamic': rate_rule.variable.mode = 'continuous' from gillespy2.core import log # pylint:disable=import-outside-toplevel - log.warning(f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule {rate_rule.name}") + errmsg = f"Changing {rate_rule.variable.name}.mode='continuous' as it is the target of RateRule" + errmsg += f" {rate_rule.name}" + log.warning(errmsg) self.listOfRateRules[rate_rule.name] = rate_rule # Build the sanitized rate rule From d35ec78a45d617cd523166a3086da62f17da4566 Mon Sep 17 00:00:00 2001 From: Brian Drawert Date: Thu, 22 Dec 2022 11:24:42 -0800 Subject: [PATCH 47/47] pylint fixes --- gillespy2/solvers/numpy/tau_hybrid_solver.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/gillespy2/solvers/numpy/tau_hybrid_solver.py b/gillespy2/solvers/numpy/tau_hybrid_solver.py index 28e1be81c..10b3101ec 100644 --- a/gillespy2/solvers/numpy/tau_hybrid_solver.py +++ b/gillespy2/solvers/numpy/tau_hybrid_solver.py @@ -181,7 +181,7 @@ def __create_diff_eqs(self, comb, dependencies, rr_sets, rate_rules): rr_vars = {} for n, rr in self.model.listOfRateRules.items(): rr_vars[rr.variable.name] = n - for spec in self.model.listOfSpecies.keys(): + for spec in self.model.listOfSpecies: if spec in rr_vars.keys(): diff_eqs[self.model.listOfSpecies[spec]] = self.model.listOfRateRules[rr_vars[spec]].formula else: @@ -638,7 +638,8 @@ def __simulate_negative_state_check(self, curr_state): # check each species to see if they are negative for s in self.non_negative_species: if curr_state[s] < 0: - raise SimulationError(f"Negative State detected at begining of step. Species involved in reactions can not be negative.") + raise SimulationError(f"Negative State detected at begining of step." \ + " Species involved in reactions can not be negative.") def __simulate_invalid_state_check(self, species_modified, curr_state, compiled_reactions): invalid_state = False @@ -678,7 +679,9 @@ def __simulate(self, integrator_options, curr_state, y0, curr_time, self.__simulate_negative_state_check(curr_state) if curr_time == 0.0: # save state at beginning of simulation - save_times, save_index = self.__save_state_to_output(curr_time, save_index, curr_state, species, trajectory, save_times) + save_times, save_index = self.__save_state_to_output( + curr_time, save_index, curr_state, species, trajectory, save_times + ) event_queue = [] prev_y0 = copy.deepcopy(y0)