Merge branch 'release_1_1_2' of https://github.com/GillesPy2/GillesPy2 …

…into release_1_1_2

gillespy2/core/__init__.py

@@ -11,6 +11,6 @@
 log = logging.getLogger()
 log.setLevel(logging.WARN)
 log.addHandler(_handler)
-version = '1.1.1'
+version = '1.1.2'
 
 __all__ = [s for s in dir() if not s.startswith('_')]

gillespy2/core/gillespy2.py

@@ -421,9 +421,7 @@ def delete_reaction(self, obj):
     def delete_all_reactions(self):
         self.listOfReactions.clear()
 
-    def run(self, number_of_trajectories=1, seed=None, solver=None, show_labels=True, 
-            switch_tol=0.03, tau_tol=0.03, integrator='lsoda', integrator_options={},
-            stoch_kit_home=None):
+    def run(self, solver=None, **solver_args):
         """
         Function calling simulation of the model. There are a number of
         parameters to be set here.
@@ -439,11 +437,11 @@ def run(self, number_of_trajectories=1, seed=None, solver=None, show_labels=True
         solver : gillespy.GillesPySolver
             The solver by which to simulate the model. This solver object may
             be initialized separately to specify an algorithm. Optional, 
-            defaults to StochKitSolver SSA.
+            defaults to ssa solver.
         show_labels: bool (True)
             If true, simulation returns a list of trajectories, where each list entry is a dictionary containing key value pairs of species : trajectory.  If false, returns a numpy array with shape [traj_no, time, species]
         switch_tol: float
-            Relative error tolerance value for deterministic/stochastic switching condition between 0.0 and 1.0
+            Tolerance for Continuous/Stochastic representation of species, based on coefficient of variance for each step.
         tau_tol: float
             Relative error tolerance value for calculating tau step between 0.0 and 1.0
         integrator: String
@@ -455,19 +453,14 @@ def run(self, number_of_trajectories=1, seed=None, solver=None, show_labels=True
         if solver is not None:
             if ((isinstance(solver, type)
                     and issubclass(solver, GillesPySolver))) or issubclass(type(solver), GillesPySolver):
-                return solver.run(model=self, t=self.tspan[-1], increment=self.tspan[-1] - self.tspan[-2],
-                                  seed=seed, number_of_trajectories=number_of_trajectories, 
-                                  show_labels=show_labels, switch_tol=0.03, tau_tol=0.03, integrator='lsoda', 
-                                  integrator_options={}, stoch_kit_home=None)
+                return solver.run(model=self, t=self.tspan[-1], increment=self.tspan[-1] - self.tspan[-2], **solver_args)
             else:
                 raise SimulationError(
                     "argument 'solver' to run() must be a subclass of GillesPySolver")
         else:
             from gillespy2.solvers.auto import SSASolver
             return SSASolver.run(model=self, t=self.tspan[-1],
-                                      increment=self.tspan[-1] - self.tspan[-2], seed=seed,
-                                      number_of_trajectories=number_of_trajectories,
-                                      show_labels=show_labels)
+                                      increment=self.tspan[-1] - self.tspan[-2], **solver_args)
 
 
 

gillespy2/solvers/cpp/ssa_c_solver.py

@@ -1,5 +1,5 @@
 import gillespy2
-from gillespy2.core import Model, Reaction, gillespyError, GillesPySolver
+from gillespy2.core import Model, Reaction, gillespyError, GillesPySolver, log
 import os #for getting directories for C++ files
 import shutil #for deleting/copying files
 import subprocess #For calling make and executing c solver
@@ -163,6 +163,9 @@ def run(self=None, model=None, t=20, number_of_trajectories=1,
             increment=0.05, seed=None, debug=False, profile=False, show_labels=True, **kwargs):
         if self is None:
             self = SSACSolver(model)
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
         if self.compiled:
             self.simulation_data = None
             number_timesteps = int(t//increment + 1)                    

gillespy2/solvers/cython/cython_ssa_solver.pyx

@@ -1,5 +1,5 @@
 # encoding: utf-8
-from gillespy2.core import GillesPySolver, gillespyError
+from gillespy2.core import GillesPySolver, gillespyError, log
 import numpy as np
 import random
 cimport numpy as np
@@ -137,10 +137,20 @@ def convert_infix_prefix(equation):
 
 class CythonSSASolver(GillesPySolver):
     name = "CythonSSASolver"
+    def __init__(self):
+        name = "CythonSSASolver"
     #@cython.boundscheck(False)
     @classmethod
     def run(self, model, t=20, number_of_trajectories=1,
             increment=0.05, seed=None, debug=False, profile=False, show_labels=True, **kwargs):
+
+        if not isinstance(self, CythonSSASolver):
+            self = CythonSSASolver()
+
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
+
         self.simulation_data = []
         #convert dictionary of species to species array
         species = list(model.listOfSpecies.keys())

gillespy2/solvers/numpy/Tau.py

@@ -96,9 +96,9 @@ def select(*tau_args):
                 if reactant not in sigma_i:
                     sigma_i[str(reactant)] = 0
                 mu_i[str(reactant)] += model.listOfReactions[r].reactants[reactant] * propensities[
-                    r]  # Cao, Gillespie, Petzold 29a
+                    r]  # Cao, Gillespie, Petzold 32a
                 sigma_i[str(reactant)] += model.listOfReactions[r].reactants[reactant] ** 2 * propensities[
-                    r]  # Cao, Gillespie, Petzold 29b
+                    r]  # Cao, Gillespie, Petzold 32b
     for r in reactants:
         calculated_max = epsilon_i[str(r)] * curr_state[r.name]
         #print('calculated max: ', calculated_max)

gillespy2/solvers/numpy/basic_ode_solver.py

@@ -4,15 +4,17 @@
 from scipy.integrate import odeint
 from collections import OrderedDict
 import numpy as np
-from gillespy2.core import GillesPySolver
-from gillespy2.core import log
+from gillespy2.core import GillesPySolver, log
 
 
 class BasicODESolver(GillesPySolver):
     """
     This Solver produces the deterministic continuous solution via ODE.
     """
     name = "BasicODESolver"
+
+    def __init__(self):
+        name = "BasicODESolver"
 
     @staticmethod
     def __f(t, y, curr_state, model, c_prop):
@@ -41,7 +43,7 @@ def __f(t, y, curr_state, model, c_prop):
         return state_change
 
     @classmethod
-    def run(cls, model, t=20, number_of_trajectories=1, increment=0.05, 
+    def run(self, model, t=20, number_of_trajectories=1, increment=0.05, 
             show_labels=True, integrator='lsoda', integrator_options={}, **kwargs):
         """
 
@@ -57,6 +59,11 @@ def run(cls, model, t=20, number_of_trajectories=1, increment=0.05,
         :param kwargs:
         :return:
         """
+        if not isinstance(self, BasicODESolver):
+            self = BasicODESolver()
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
         if number_of_trajectories > 1:
             log.warning("Generating duplicate trajectories for model with ODE Solver. Consider running with only 1 trajectory.")
 

gillespy2/solvers/numpy/basic_tau_hybrid_solver.py

@@ -4,7 +4,7 @@
 import numpy as np
 import gillespy2
 from gillespy2.solvers.numpy import Tau
-from gillespy2.core import GillesPySolver
+from gillespy2.core import GillesPySolver, log
 from gillespy2.core.gillespyError import *
 
 eval_globals = math.__dict__
@@ -19,10 +19,10 @@ class BasicTauHybridSolver(GillesPySolver):
     is bounded by bounding the relative change in the state of the system, resulting in increased
     run-time performance with little accuracy trade-off.
     """
-    name = "Basic Tau Hybrid Solver"
+    name = "BasicTauHybridSolver"
 
-    def __init__(self, debug=False):
-        self.debug = debug
+    def __init__(self):
+        name = 'BasicTauHybridSolver'
 
 
     def toggle_reactions(self, model, all_compiled, deterministic_reactions, dependencies, curr_state, rxn_offset, det_spec):
@@ -32,6 +32,15 @@ def toggle_reactions(self, model, all_compiled, deterministic_reactions, depende
         rate_rules = all_compiled['rules']
         rxns = all_compiled['rxns']
 
+        #print(inactive_reactions.keys())
+        #print(rxns.keys())
+        '''
+        print('ALL evals for rate rule combos')
+        for comb in rate_rules:
+            print('current deterministic set: ', comb)
+            for species, rr in rate_rules[comb].items():
+                print('Species: ', species, 'evaluated: ', eval(rr, eval_globals, curr_state))
+        '''
         #If the set has changed, reactivate non-determinsitic reactions
         reactivate = []
         for r in inactive_reactions:
@@ -73,7 +82,6 @@ def create_diff_eqs(self, comb, model, dependencies):
         # loop through each det reaction and concatenate it's diff eq for each species
         for reaction in comb:
             factor = OrderedDict()
-            pure_continuous = True
             for dep in dependencies[reaction]:
                 if model.listOfSpecies[dep].mode != 'continuous':
                     pure_continuous = False
@@ -158,7 +166,6 @@ def __f(t, y, curr_state, reactions, rate_rules, propensities, compiled_reaction
         """
         curr_state['t'] = t
         state_change = []
-
         for i, rr in enumerate(compiled_rate_rules):
             curr_state[rr] = y[i]
         for i, rr in enumerate(compiled_rate_rules):
@@ -224,7 +231,6 @@ def __get_reactions(self, integrator, integrator_options, step, curr_state, y0,
                 urn = (math.log(random.uniform(0, 1)))
                 current[i+len(compiled_rate_rules)] += urn
                 rxn_offset[r] += urn
-
 
 
         if debug:
@@ -275,10 +281,13 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
             Example use: {max_step : 0, rtol : .01}
         """
 
-        if not sys.warnoptions:
-            warnings.simplefilter("ignore")
         if not isinstance(self, BasicTauHybridSolver):
             self = BasicTauHybridSolver()
+
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
+
         if debug:
             print("t = ", t)
             print("increment = ", increment)

gillespy2/solvers/numpy/basic_tau_leaping_solver.py

@@ -4,7 +4,7 @@
 import random, math, sys, warnings
 import numpy as np
 from gillespy2.solvers.numpy import Tau
-from gillespy2.core import GillesPySolver
+from gillespy2.core import GillesPySolver, log
 
 
 class BasicTauLeapingSolver(GillesPySolver):
@@ -14,12 +14,12 @@ class BasicTauLeapingSolver(GillesPySolver):
     over this step are bounded by bounding the relative change in state, yielding greatly improved
     run-time performance with very little trade-off in accuracy.
     """
-    name = "Basic Tau Leaping Solver"
+    name = "BasicTauLeapingSolver"
 
     def __init__(self, debug=False, profile=False):
+        name = "BasicTauLeapingSolver"
         self.debug = debug
         self.profile = profile
-        self.epsilon = 0.03
 
     def get_reactions(self, seed, step, curr_state, curr_time, save_time, propensities, reactions):
         """
@@ -84,10 +84,13 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
                 show_labels : bool (True)
                     Use names of species as index of result object rather than position numbers.
                 """
-        if not sys.warnoptions:
-            warnings.simplefilter("ignore")
+
         if not isinstance(self, BasicTauLeapingSolver):
-            self = BasicTauLeapingSolver()
+            self = BasicTauLeapingSolver(debug=debug, profile=profile)
+
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
         if debug:
             print("t = ", t)
             print("increment = ", increment)

gillespy2/solvers/numpy/ssa_solver.py

@@ -1,4 +1,4 @@
-from gillespy2.core import GillesPySolver, Model, Reaction
+from gillespy2.core import GillesPySolver, Model, Reaction, log
 import random
 import math
 import numpy as np
@@ -7,8 +7,11 @@
 class NumPySSASolver(GillesPySolver):
     name = "NumPySSASolver"
 
-    @staticmethod
-    def run(model, t=20, number_of_trajectories=1, increment=0.05, seed=None, debug=False, show_labels=True, **kwargs):
+    def __init__(self):
+        name = 'NumPySSASolver'
+
+    @classmethod
+    def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None, debug=False, show_labels=True, **kwargs):
         """
         Run the SSA algorithm using a NumPy for storing the data in arrays and generating the timeline.
         :param model: The model on which the solver will operate.
@@ -22,6 +25,13 @@ def run(model, t=20, number_of_trajectories=1, increment=0.05, seed=None, debug=
         :param show_labels: Use names of species as index of result object rather than position numbers.
         :return: a list of each trajectory simulated.
         """
+
+        if not isinstance(self, NumPySSASolver):
+            self = NumPySSASolver()
+
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
         random.seed(seed)
         # create mapping of species dictionary to array indices
         species_mappings = model.sanitized_species_names()

gillespy2/solvers/stochkit/stochkit_solvers.py

@@ -5,7 +5,7 @@
 import uuid
 import subprocess
 import shutil
-from gillespy2.core import GillesPySolver, Model
+from gillespy2.core import GillesPySolver, Model, log
 from gillespy2.core.gillespyError import SimulationError, InvalidModelError
 
 
@@ -54,6 +54,9 @@ def run(cls, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
         """
         Call out and run the solver. Collect the results.
         """
+        if len(kwargs) > 0:
+            for key in kwargs:
+                log.warning('Unsupported keyword argument to solver: {0}'.format(key))
 
         if algorithm is None:
             raise SimulationError("No algorithm selected")

setup.py

@@ -54,7 +54,7 @@ class easy_install_new(easy_install):
     full_description = fh.read()
 
 setup(name="gillespy2",
-      version="1.1.1",
+      version="1.1.2",
       packages=find_packages('.'),
       include_package_data=True,
       description='Python interface for Gillespie style biochemical simulations',

test/test_all_solvers.py

@@ -40,6 +40,13 @@ def test_random_seed(self):
             diff_results = self.model.run(solver=solver, show_labels=False, seed=2)
             if solver.name != 'BasicODESolver':
                 self.assertFalse(np.array_equal(diff_results, self.results[solver]))
+
+    def test_extraneous_args(self):
+        for solver in self.solvers:
+            print(solver.name)
+            with self.assertWarns(Warning):
+                model = Example()
+                results = model.run(solver=solver, nonsense='ABC')
 
 if __name__ == '__main__':
     unittest.main()