Merge pull request #322 from GillesPy2/hotfix-eval-order

Hotfix eval order

gillespy2/solvers/numpy/basic_ode_solver.py

@@ -76,16 +76,31 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05,
                 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.")
-        sim_thread = Thread(target=self.__run, args=(model,), kwargs={'t':t,
+        sim_thread = Thread(target=self.___run, args=(model,), kwargs={'t':t,
                                         'number_of_trajectories':number_of_trajectories,
                                         'increment':increment, 'show_labels':show_labels, 
                                         'timeout':timeout})
-        sim_thread.start()
-        sim_thread.join(timeout=timeout)
-        self.stop_event.set()
-        while self.result is None: pass
+        try:
+            sim_thread.start()
+            sim_thread.join(timeout=timeout)
+            self.stop_event.set()
+            while self.result is None: pass
+        except:
+            pass
+        if hasattr(self, 'has_raised_exception'):
+            raise self.has_raised_exception
         return self.result, self.rc
 
+    def ___run(self, model, t=20, number_of_trajectories=1, increment=0.05, timeout=None,
+            show_labels=True, integrator='lsoda', integrator_options={}, **kwargs):
+        try:
+            self.__run(model, t, number_of_trajectories, increment, timeout,
+                        show_labels, integrator, integrator_options, **kwargs)
+        except Exception as e:
+            self.has_raised_exception = e
+            self.result = []
+            return [], -1
+
     def __run(self, model, t=20, number_of_trajectories=1, increment=0.05, timeout=None,
             show_labels=True, integrator='lsoda', integrator_options={}, **kwargs):
 

gillespy2/solvers/numpy/basic_tau_hybrid_solver.py

@@ -211,19 +211,19 @@ def __f(t, y, curr_state, species, reactions, rate_rules, propensities,
         curr_state['time'] = t
         for item, index in y_map.items():
             if item in assignment_rules:
-                curr_state[item] = eval(assignment_rules[item].formula, {**curr_state, **eval_globals})
+                curr_state[item] = eval(assignment_rules[item].formula, {**eval_globals, **curr_state})
             else:
                 curr_state[item] = y[index]
         for rr in compiled_rate_rules:
             try:
-                state_change[y_map[rr]] += eval(compiled_rate_rules[rr], {**curr_state, **eval_globals})
+                state_change[y_map[rr]] += eval(compiled_rate_rules[rr], {**eval_globals, **curr_state})
             except ValueError:
                 pass
         for i, r in enumerate(compiled_reactions):
-            propensities[r] = eval(compiled_reactions[r],{**curr_state, **eval_globals})
+            propensities[r] = eval(compiled_reactions[r],{**eval_globals, **curr_state})
             state_change[y_map[r]] += propensities[r]
         for event in events:
-            triggered = eval(event.trigger.expression, {**curr_state, **eval_globals})
+            triggered = eval(event.trigger.expression, {**eval_globals, **curr_state})
             if triggered: state_change[y_map[event]] = 1
 
 
@@ -367,7 +367,7 @@ def __handle_event(self, event, curr_state, curr_time, event_queue,
         else:
             curr_state['t'] = curr_time
             curr_state['time'] = curr_time
-            execution_time = curr_time + eval(event.delay,eval_globals, curr_state)
+            execution_time = curr_time + eval(event.delay, {**eval_globals,**curr_state})
             curr_state[event.name] = True
             heapq.heappush(delayed_events, (execution_time, event.name))
             if event.use_values_from_trigger_time:
@@ -385,14 +385,14 @@ def __check_t0_events(self, model, initial_state):
         t0_delayed_events = {}
         for e in model.listOfEvents.values():
             if not e.trigger.value:
-                t0_firing = eval(e.trigger.expression, eval_globals, initial_state)
+                t0_firing = eval(e.trigger.expression, {**eval_globals,**initial_state})
                 if t0_firing:
                     if e.delay is None:
                         for a in e.assignments:
-                            initial_state[a.variable.name] = eval(a.expression, eval_globals, initial_state)
+                            initial_state[a.variable.name] = eval(a.expression,{**eval_globals, **initial_state})
                             species_modified_by_events.append(a.variable.name)
                     else:
-                        execution_time = eval(e.delay,eval_globals, initial_state)
+                        execution_time = eval(e.delay,{**eval_globals,**initial_state})
                         t0_delayed_events[e.name] = execution_time
         return t0_delayed_events, species_modified_by_events
 
@@ -403,7 +403,6 @@ def __update_stochastic_rxn_states(self, model, compiled_reactions, curr_state):
         rxn_count = OrderedDict()
         species_modified = OrderedDict()
         # Update stochastic reactions
-
         for rxn in compiled_reactions:
             rxn_count[rxn] = 0
             while curr_state[rxn] > 0:
@@ -604,7 +603,7 @@ def __simulate(self, integrator, integrator_options, curr_state, y0, model, curr
                     assignment_state[species[s]] = sol.sol(time)[s]
             assignment_state['t'] = time
             for spec, ar in model.listOfAssignmentRules.items():
-                assignment_value = eval(ar.formula, eval_globals, assignment_state)
+                assignment_value = eval(ar.formula, {**eval_globals,**assignment_state})
                 assignment_state[spec] = assignment_value
                 trajectory[trajectory_index][species.index(spec)+1] = assignment_value
             num_saves += 1
@@ -618,7 +617,7 @@ def __simulate(self, integrator, integrator_options, curr_state, y0, model, curr
         while event_cycle:
             event_cycle = False
             for i, e in enumerate(model.listOfEvents.values()):
-                triggered = eval(e.trigger.expression, eval_globals, curr_state)
+                triggered = eval(e.trigger.expression, {**eval_globals,**curr_state})
                 if triggered and not curr_state[e.name]:
                     curr_state[e.name] = True
                     self.__handle_event(e, curr_state, curr_time, 
@@ -706,7 +705,7 @@ def __map_state(self, species, parameters, compiled_reactions, events, curr_stat
         y0 = [0] * (len(species) + len(parameters) + len(compiled_reactions) + len(events))
         for i, spec in enumerate(species):
             if isinstance(curr_state[spec], str):
-                y0[i] = eval(curr_state[spec], eval_globals, curr_state)
+                y0[i] = eval(curr_state[spec], {**eval_globals, **curr_state})
             else:
                 y0[i] = curr_state[spec]
             y_map[spec] = i
@@ -780,22 +779,40 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
 
         if timeout is not None and timeout <= 0: timeout = None
 
-        sim_thread = threading.Thread(target=self.__run, args=(model,), kwargs={'t':t,
+        sim_thread = threading.Thread(target=self.___run, args=(model,), kwargs={'t':t,
                                         'number_of_trajectories':number_of_trajectories,
                                         'increment':increment, 'seed':seed,
                                         'debug':debug, 'profile':profile,'show_labels':show_labels,
                                         'timeout':timeout, 'tau_tol':tau_tol,
                                         'event_sensitivity':event_sensitivity,
                                         'integrator':integrator,
                                         'integrator_options':integrator_options})
-        sim_thread.start()
-        sim_thread.join(timeout=timeout)
-        self.stop_event.set()
-        while self.result is None: pass
+        try:
+            sim_thread.start()
+            sim_thread.join(timeout=timeout)
+            self.stop_event.set()
+            while self.result is None: pass
+        except:
+            pass
+        if hasattr(self,'has_raised_exception'):
+            raise self.has_raised_exception
         return self.result, self.rc
 
 
 
+    def ___run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None, 
+            debug=False, profile=False, show_labels=True,
+            tau_tol=0.03, event_sensitivity=100, integrator='LSODA',
+            integrator_options={}, **kwargs):
+            try:
+                self.__run(model,t,number_of_trajectories, increment, seed, debug,
+                           profile,show_labels, tau_tol, event_sensitivity, integrator,
+                           integrator_options, **kwargs)
+            except Exception as e:
+                self.has_raised_exception = e
+                self.result = []
+                return [], -1
+
     def __run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None, 
             debug=False, profile=False, show_labels=True,
             tau_tol=0.03, event_sensitivity=100, integrator='LSODA',
@@ -927,7 +944,7 @@ def __run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None
                 if not pure_ode:
                     for i, r in enumerate(model.listOfReactions):
                         try:
-                            propensities[r] = eval(compiled_propensities[r], eval_globals, curr_state)
+                            propensities[r] = eval(compiled_propensities[r],{**eval_globals, **curr_state})
                         except Exception as e:
                             raise SimulationError('Error calculation propensity for {0}.\nReason: {1}'.format(r, e))
 

gillespy2/solvers/numpy/basic_tau_leaping_solver.py

@@ -102,17 +102,34 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
             for key in kwargs:
                 log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
 
-        sim_thread = Thread(target=self.__run, args=(model,), kwargs={'t':t,
+        sim_thread = Thread(target=self.___run, args=(model,), kwargs={'t':t,
                                         'number_of_trajectories':number_of_trajectories,
                                         'increment':increment, 'seed':seed,
                                         'debug':debug, 'show_labels':show_labels,
                                         'timeout':timeout, 'tau_tol':tau_tol})
-        sim_thread.start()
-        sim_thread.join(timeout=timeout)
-        self.stop_event.set()
-        while self.result is None: pass
+        try:
+            sim_thread.start()
+            sim_thread.join(timeout=timeout)
+            self.stop_event.set()
+            while self.result is None: pass
+        except:
+            pass
+        if hasattr(self, 'has_raised_exception'):
+            raise self.has_raised_exception
         return self.result, self.rc
 
+    def ___run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
+                    debug=False, profile=False, show_labels=True, 
+                    timeout=None, tau_tol=0.03, **kwargs):
+        try:
+            self.__run(model, t, number_of_trajectories, increment, seed,
+                        debug, profile, show_labels, timeout, tau_tol, **kwargs)
+        except Exception as e:
+            self.has_raised_exception = e
+            self.result = []
+            return [], -1
+
+
     def __run(self, model, t=20, number_of_trajectories=1, increment=0.05, seed=None,
                     debug=False, profile=False, show_labels=True, 
                     timeout=None, tau_tol=0.03, **kwargs):

gillespy2/solvers/numpy/ssa_solver.py

@@ -44,17 +44,32 @@ def run(self, model, t=20, number_of_trajectories=1, increment=0.05,
         if len(kwargs) > 0:
             for key in kwargs:
                 log.warning('Unsupported keyword argument to {0} solver: {1}'.format(self.name, key))
-        sim_thread = Thread(target=self.__run, args=(model,), kwargs={'t':t,
+        sim_thread = Thread(target=self.___run, args=(model,), kwargs={'t':t,
                                         'number_of_trajectories':number_of_trajectories,
                                         'increment':increment, 'seed':seed,
                                         'debug':debug, 'show_labels':show_labels,
                                         'timeout':timeout})
-        sim_thread.start()
-        sim_thread.join(timeout=timeout)
-        self.stop_event.set()
-        while self.result is None: pass
+        try:
+            sim_thread.start()
+            sim_thread.join(timeout=timeout)
+            self.stop_event.set()
+            while self.result is None: pass
+        except:
+            pass
+        if hasattr(self, 'has_raised_exception'):
+            raise self.has_raised_exception
         return self.result, self.rc
 
+    def ___run(self, model, t=20, number_of_trajectories=1, increment=0.05,
+                    seed=None, debug=False, show_labels=True, timeout=None):
+        try:
+            self.__run(model, t, number_of_trajectories, increment, seed,
+                            debug, show_labels, timeout)
+        except Exception as e:
+            self.has_raised_exception = e
+            self.result = []
+            return [], -1
+
     def __run(self, model, t=20, number_of_trajectories=1, increment=0.05,
                     seed=None, debug=False, show_labels=True, timeout=None):
 

test/sbml_test_runner/wrappers/hybrid_wrapper_for_sts.sh

@@ -119,8 +119,8 @@ EXPECTED_RESULTS_FILE = os.path.join(TEST_DIR, '{0}-results.csv'.format(CASE_NO)
 model, errors = gillespy2.import_SBML(TEST_FILE)
 
 # Create absolute and relative tolerance defaults
-atol = 1e-9
-rtol = 1e-6
+atol = 1e-12
+rtol = 1e-12
 
 # Retrieve simulation settings from file.
 start, duration, steps = [0]*3
@@ -129,8 +129,8 @@ with open(SETTINGS_FILE, 'r') as settings:
         if 'start' in line: start = float(line.split(': ')[1])
         elif 'duration' in line: duration = float(line.split(': ')[1])
         elif 'steps' in line: steps = int(line.split(': ')[1])
-        elif 'absolute' in line: atol = float(line.split(': ')[1])
-        elif 'relative' in line: rtol = float(line.split(': ')[1])
+        #elif 'absolute' in line: atol = float(line.split(': ')[1])
+        #elif 'relative' in line: rtol = float(line.split(': ')[1])
 
 #Force Continuous Species
 for species in model.listOfSpecies.values():

test/test_hybrid_solver.py

@@ -32,6 +32,15 @@ def test_add_bad_species_rate_rule_dict(self):
         with self.assertRaises(ModelError):
             self.model.add_rate_rule(rule)
 
+    def test_math_name_overlap(self):
+        gamma = gillespy2.Species('gamma',initial_value=2, mode='continuous')
+        self.model.add_species([gamma])
+        k2 = gillespy2.Parameter(name='k2', expression=1)
+        self.model.add_parameter([k2])
+        gamma_react = gillespy2.Reaction(name='gamma_react', reactants={'gamma':1}, products={}, rate=k2)
+        self.model.add_reaction([gamma_react])
+        self.model.run(solver=BasicTauHybridSolver)
+
     def test_add_bad_expression_rate_rule_dict(self):
         species2 = gillespy2.Species('test_species2', initial_value=2, mode='continuous')
         rule = gillespy2.RateRule(variable=species2, formula='')