Skip to content
Browse files

Small modifications to support BNG's 'stop_if' method (#407)

* Some small modifications to support the 'stop_if' method in BNG, which
allows one to define a logical function that if True will terminate a
simulation. First problem is that BNG creates a new internal function
called '_stop_if' that gets output to the .gdat file. This was causing
problems in PySB since the number of defined Expressions wasn't matching
the number of functions in the .gdat file. Fixed this by explicitly
using the length of 'model.expressions' when reading in the .gdat file.
Second problem is that if 'stop_if' is satisfied at the outset BNG will
immediately quit, only outputting the initial state to the .cdat and
.gdat files. This meant that PySB was reading this in as a single numpy
array, rather than an array of arrays, which was causing a downstream
error. Fixed this by requiring the number of dimensions 'ndmin=2' in the
'numpy.loadtxt' calls. Also added a small test model to the validation
test suite called 'test_stop_if'.

* Fix BNG gdat parsing when constant expressions in model

* Fix test so works consistently across platforms
  • Loading branch information...
lh64 authored and alubbock committed Jan 25, 2019
1 parent c267f1e commit 470945e9c9f27a5c7b8c24cc4204d7e0533f3266
Showing with 25 additions and 4 deletions.
  1. +2 −2 pysb/
  2. +4 −1 pysb/simulator/
  3. +19 −1 pysb/tests/
@@ -225,7 +225,7 @@ def read_simulation_results_multi(base_filenames):

# Read concentrations data
cdat_arr = numpy.loadtxt(base_filename + '.cdat', skiprows=1)
cdat_arr = numpy.loadtxt(base_filename + '.cdat', skiprows=1, ndmin=2)
# -1 for time column
names += ['__s%d' % i for i in range(cdat_arr.shape[1] - 1)]
except IOError:
@@ -237,7 +237,7 @@ def read_simulation_results_multi(base_filenames):
# Exclude \# and time column
names += f.readline().split()[2:]
# Exclude first column (time)
gdat_arr = numpy.loadtxt(f)
gdat_arr = numpy.loadtxt(f, ndmin=2)
if cdat_arr is None:
cdat_arr = numpy.ndarray((len(gdat_arr), 0))
@@ -228,7 +228,10 @@ def run(self, tspan=None, initials=None, param_values=None, n_runs=1,
1:(len(self.model.species) + 1)])
if len(self.model.observables) or len(self.model.expressions):
(len(self.model.species) + 1):])
(len(self.model.species) + 1):
(len(self.model.species) + 1) +
len(self.model.observables) +

return SimulationResult(self, tout=tout, trajectories=species_out,
@@ -1,6 +1,6 @@
from pysb.testing import *
import numpy as np
from pysb import Monomer, Parameter, Initial, Observable, Rule
from pysb import Monomer, Parameter, Initial, Observable, Rule, Expression
from pysb.simulator.bng import BngSimulator, PopulationMap
from pysb.bng import generate_equations
from pysb.examples import robertson, expression_observables, earm_1_0
@@ -122,3 +122,21 @@ def test_hpp():
observables = np.array(x.observables)
assert len(observables) == 50

def test_stop_if():
Rule('A_synth', None >> A(), Parameter('k', 1))
Observable('Atot', A())
Expression('exp_const', k + 1)
Expression('exp_dyn', Atot + 1)
sim = BngSimulator(model, verbose=5)
tspan = np.linspace(0, 100, 101)
x =, stop_if='Atot>9', seed=_BNG_SEED)
# All except the last Atot value should be <=9
assert all(x.observables['Atot'][:-1] <= 9)
assert x.observables['Atot'][-1] > 9
# Starting with Atot > 9 should terminate simulation immediately
y =, initials=x.species[-1], stop_if='Atot>9')
assert len(y.observables) == 1

0 comments on commit 470945e

Please sign in to comment.
You can’t perform that action at this time.