smfsb package

Submodules

smfsb.models module

smfsb.models.bd(th=[1, 1.1])

Create a birth-death model

Create and return a Spn object representing a discrete stochastic birth-death model.

Parameters

th: array: array of length 2 containing the birth and death rates

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> bd = smfsb.models.bd()
>>> step = bd.step_gillespie()
>>> smfsb.sim_time_series(bd.m, 0, 50, 0.1, step)

smfsb.models.dimer(th=[0.00166, 0.2])

Create a dimerisation kinetics model

Create and return a Spn object representing a discrete stochastic dimerisation kinetics model.

Parameters

th: array: array of length 2 containing the rates of the bind and unbind reactions

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> dimer = smfsb.models.dimer()
>>> step = dimer.step_gillespie()
>>> smfsb.sim_time_series(dimer.m, 0, 50, 0.1, step)

smfsb.models.id(th=[1, 0.1])

Create an immigration-death model

Create and return a Spn object representing a discrete stochastic immigration-death model.

Parameters

th: array: array of length 2 containing the immigration and death rates

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> id = smfsb.models.id()
>>> step = id.step_gillespie()
>>> smfsb.sim_time_series(id.m, 0, 50, 0.1, step)

smfsb.models.lv(th=[1, 0.005, 0.6])

Create a Lotka-Volterra model

Create and return a Spn object representing a discrete stochastic Lotka-Volterra model.

Parameters

th: array: array of length 3 containing the rates of the three governing reactions, prey reproduction, predator-prey interaction, and predator death

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> lv = smfsb.models.lv()
>>> step = lv.step_gillespie()
>>> smfsb.sim_time_series(lv.m, 0, 50, 0.1, step)

smfsb.models.mm(th=[0.00166, 0.0001, 0.1])

Create a Michaelis-Menten enzyme kinetic model

Create and return a Spn object representing a discrete stochastic Michaelis-Menten enzyme kinetic model.

Parameters

th: array: array of length 3 containing the binding, unbinding and production rates

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> mm = smfsb.models.mm()
>>> step = mm.step_gillespie()
>>> smfsb.sim_time_series(mm.m, 0, 50, 0.1, step)

smfsb.models.sir(th=[0.0015, 0.1])

Create a basic SIR compartmental epidemic model

Create and return a Spn object representing a discrete stochastic SIR model.

Parameters

th: array: array of length 2 containing the rates of the two governing transitions

Returns

Spn model object with rates th

Examples

>>> import smfsb
>>> sir = smfsb.models.sir()
>>> step = sir.step_gillespie()
>>> smfsb.sim_time_series(sir.m, 0, 50, 0.1, step)

smfsb.data module

smfsb.spn module

class smfsb.spn.Spn(n, t, pre, post, h, m)

Bases: object

Class for stochastic Petri net models.

gillespie(n)

Simulate a sample path from a stochastic kinetic model described by a stochastic Petri net

This function simulates a single realisation from a discrete stochastic kinetic model described by a stochastic Petri net (SPN).

Parameters

n: int An integer representing the number of events to simulate, excluding the initial state

Returns

A tuple consisting of a first component, a vector of length n containing the event times, and a second component, a matrix with n+1 rows containing the state of the system. The first row is the intial state prior to the first event.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> lv.gillespie(1000)

gillespied(tt, dt=1)

Simulate a sample path from a stochastic kinetic model described by a stochastic Petri net

This function simulates a single realisation from a discrete stochastic kinetic model described by a stochastic Petri net and discretises the output onto a regular time grid.

Parameters

tt: float The required length of simulation time. dt: float The grid size for the output. Note that this parameter simply determines the volume of output. It has no bearing on the correctness of the simulation algorithm. Defaults to one time unit.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> lv.gillespied(30, 0.1)

step_cle(dt=0.01)

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama integration method for the associated CLE

This method returns a function for advancing the state of an SPN model using a simple Euler-Maruyama integration method method for the chemical Langevin equation form of the model.The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SPN models.

Parameters

dtfloat: The time step for the time-stepping integration method. Defaults to 0.01.

Returns

A function which can be used to advance the state of the SPN model by using an Euler-Maruyama method with step size ‘dt’. The function closure has interface ‘function(x0, t0, deltat)’, where ‘x0’ and ‘t0’ represent the initial state and time, and ‘deltat’ represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_cle(0.001)
>>> stepLv([50, 100], 0, 1)

step_cle_1d(d, dt=0.01)

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama discretisation of the CLE on a 1D regular grid

This method creates a function for advancing the state of an SPN model using a simple Euler-Maruyama discretisation of the CLE on a 1D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_1d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore twice this value (as it can leave to the left or the right).
dtfloat: Time step for the Euler-Maruyama discretisation.

Returns

A function which can be used to advance the state of the SPN model by using a simple Euler-Maruyama algorithm. The function closure has parameters x0, t0, deltat, where x0 is a matrix with rows corresponding to species and columns corresponding to voxels, representing the initial condition, t0 represents the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv1d = lv.step_cle_1d(np.array([0.6,0.6]))
>>> N = 20
>>> x0 = np.zeros((2,N))
>>> x0[:,int(N/2)] = lv.m
>>> stepLv1d(x0, 0, 1)

step_cle_2d(d, dt=0.01)

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama discretisation of the CLE on a 2D regular grid

This method creates a function for advancing the state of an SPN model using a simple Euler-Maruyama discretisation of the CLE on a 2D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_2d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore four times this value (as it can leave in one of 4 directions).
dtfloat: Time step for the Euler-Maruyama discretisation.

Returns

A function which can be used to advance the state of the SPN model by using a simple Euler-Maruyama algorithm. The function closure has parameters x0, t0, deltat, where x0 is a 3d array with indices species, then rows and columns corresponding to voxels, representing the initial condition, t0 represents the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv2d = lv.step_cle_2d(np.array([0.6,0.6]))
>>> M = 15
>>> N = 20
>>> x0 = np.zeros((2,M,N))
>>> x0[:,int(M/2),int(N/2)] = lv.m
>>> stepLv2d(x0, 0, 1)

step_euler(dt=0.01)

Create a function for advancing the state of an SPN by using a simple continuous deterministic Euler integration method

This method returns a function for advancing the state of an SPN model using a simple continuous deterministic Euler integration method. The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SPN models.

Parameters

dtfloat: The time step for the time-stepping integration method. Defaults to 0.01.

Returns

A function which can be used to advance the state of the SPN model by using an Euler method with step size ‘dt’. The function closure has interface ‘function(x0, t0, deltat)’, where ‘x0’ and ‘t0’ represent the initial state and time, and ‘deltat’ represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_euler(0.001)
>>> stepLv([50, 100], 0, 1)

step_euler_1d(d, dt=0.01)

Create a function for advancing the state of an SPN by using a simple forward Euler discretisation of the reaction-diffusion on a 1D regular grid

This method creates a function for advancing the state of an SPN model using a simple forward Euler discretisation of the reaction-diffusion on a 1D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_1d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore twice this value (as it can leave to the left or the right).
dtfloat: Time step for the Euler-Maruyama discretisation.

Returns

A function which can be used to advance the state of the SPN model by using a simple forward Euler algorithm. The function closure has parameters x0, t0, deltat, where x0 is a matrix with rows corresponding to species and columns corresponding to voxels, representing the initial condition, t0 represents the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv1d = lv.step_euler_1d(np.array([0.6,0.6]))
>>> N = 20
>>> x0 = np.zeros((2,N))
>>> x0[:,int(N/2)] = lv.m
>>> stepLv1d(x0, 0, 1)

step_euler_2d(d, dt=0.01)

Create a function for advancing the state of an SPN by using a simple forward Euler discretisation of the reaction-diffusion on a 2D regular grid

This method creates a function for advancing the state of an SPN model using a simple forward Euler discretisation of the reaction-diffusion on a 2D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_2d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore four times this value (as it can leave in one of 4 directions).
dtfloat: Time step for the Euler discretisation.

Returns

A function which can be used to advance the state of the SPN model by using a simple forward Euler algorithm. The function closure has parameters x0, t0, deltat, where x0 is a 3d array with indices species, then rows and columns corresponding to voxels, representing the initial condition, t0 represents the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv2d = lv.step_euler_2d(np.array([0.6,0.6]))
>>> M = 15
>>> N = 20
>>> x0 = np.zeros((2,M,N))
>>> x0[:,int(M/2),int(N/2)] = lv.m
>>> stepLv2d(x0, 0, 1)

step_first()

Create a function for advancing the state of an SPN by using Gillespie’s first reaction method

This function creates a function for advancing the state of an SPN model using Gillespie’s first reaction method. The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SPN models.

Returns

A function which can be used to advance the state of the SPN model by using Gillespie’s first reaction method. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_first()
>>> stepLv([50, 100], 0, 1)

step_gillespie(min_haz=1e-10, max_haz=10000000.0)

Create a function for advancing the state of a SPN by using the Gillespie algorithm

This method returns a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SPN models.

Parameters

min_hazfloat: Minimum hazard to consider before assuming 0. Defaults to 1e-10.
max_hazfloat: Maximum hazard to consider before assuming an explosion and bailing out. Defaults to 1e07.

Returns

A function which can be used to advance the state of the SPN model by using the Gillespie algorithm. The function closure has interface ‘function(x0, t0, deltat)’, where ‘x0’ and ‘t0’ represent the initial state and time, and ‘deltat’ represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_gillespie()
>>> stepLv([50, 100], 0, 1)

step_gillespie_1d(d, min_haz=1e-10, max_haz=10000000.0)

Create a function for advancing the state of an SPN by using the Gillespie algorithm on a 1D regular grid

This method creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_1d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore twice this value (as it can leave to the left or the right).
min_hazfloat: Minimum hazard to consider before assuming 0. Defaults to 1e-10.
max_hazfloat: Maximum hazard to consider before assuming an explosion and bailing out. Defaults to 1e07.

Returns

A function which can be used to advance the state of the SPN model by using the Gillespie algorithm. The function closure has arguments x0, t0, deltat, where x0 is a matrix with rows corresponding to species and columns corresponding to voxels, representing the initial condition, t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv1d = lv.step_gillespie_1d(np.array([0.6,0.6]))
>>> N = 20
>>> x0 = np.zeros((2,N))
>>> x0[:,int(N/2)] = lv.m
>>> stepLv1d(x0, 0, 1)

step_gillespie_2d(d, min_haz=1e-10, max_haz=10000000.0)

Create a function for advancing the state of an SPN by using the Gillespie algorithm on a 2D regular grid

This method creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as sim_time_series_2d) for simulating realisations of SPN models in space and time.

Parameters

darray: A vector of diffusion coefficients - one coefficient for each reacting species, in order. The coefficient is the reaction rate for a reaction for a molecule moving into an adjacent compartment. The hazard for a given molecule leaving the compartment is therefore four times this value (as it can leave in one of 4 directions).
min_hazfloat: Minimum hazard to consider before assuming 0. Defaults to 1e-10.
max_hazfloat: Maximum hazard to consider before assuming an explosion and bailing out. Defaults to 1e07.

Returns

A function which can be used to advance the state of the SPN model by using the Gillespie algorithm. The function closure has arguments x0, t0, deltat, where x0 is a 3d array with dimensions corresponding to species then two spatial dimensions, representing the initial condition, t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns an array representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv2d = lv.step_gillespie_2d(np.array([0.6, 0.6]))
>>> N = 20
>>> x0 = np.zeros((2, N, N))
>>> x0[:, int(N/2), int(N/2)] = lv.m
>>> stepLv2d(x0, 0, 1)

step_poisson(dt=0.01)

Create a function for advancing the state of an SPN by using a simple approximate Poisson time stepping method

This method returns a function for advancing the state of an SPN model using a simple approximate Poisson time stepping method. The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SPN models.

Parameters

dtfloat: The time step for the time-stepping integration method. Defaults to 0.01.

Returns

A function which can be used to advance the state of the SPN model by using a Poisson time stepping method with step size ‘dt’. The function closure has interface ‘function(x0, t0, deltat)’, where ‘x0’ and ‘t0’ represent the initial state and time, and ‘deltat’ represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_poisson(0.001)
>>> stepLv([50, 100], 0, 1)

smfsb.sim module

smfsb.sim.discretise(times, states, dt=1, start=0)

Discretise output from a discrete event simulation algorithm

This function discretises output from a discrete event simulation algorithm such as ‘gillespie’ onto a regular time grid, and returns the results as a matrix.

Parameters

times: array: A vector of event times.
states: array: A matrix of states. There should be one more row than the length of times.
dt: float: The time step required for the output of the discretisation process. Defaults to one time unit.
start: float: The start time for the output. Defaults to zero.

Returns

A matrix with rows corresponding to the state of the system on a regular grid.

Examples

>>> import smfsb
>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> times, states = lv.gillespie(1000)
>>> smfsb.discretise(times, states, 0.1)

smfsb.sim.imdeath(n=20, x0=0, lamb=1, mu=0.1)

Simulate a sample path from the homogeneous immigration-death process

This function simulates a single realisation from a time-homogeneous immigration-death process.

Parameters

n: int: The number of states to be sampled from the process, not including the initial state, ‘x0’
x0: int: The initial state of the process, which defaults to zero.
lamb: float: The rate at which new individual immigrate into the population. Defaults to 1.
mu: float: The rate at which individuals within the population die, independently of all other individuals. Defaults to 0.1.

Returns

A tuple, (tvec, xvec), where tvec is a vector of event times of length n and xvec is a vector of states, of length n+1.

Examples

>>> import smfsb
>>> smfsb.imdeath(100)

smfsb.sim.rcfmc(n, q_mat, pi0)

Simulate a continuous time finite state space Markov chain

This function simulates a single realisation from a continuous time Markov chain having a finite state space based on a given transition rate matrix.

Parameters

n: int: The number of states to be sampled from the Markov chain, including the initial state, which will be sampled using ‘pi0’.
q_mat: matrix: The transition rate matrix of the Markov chain, where each off-diagonal element ‘q_mat[i,j]’ represents the rate of transition from state ‘i’ to state ‘j’. This matrix is assumed to be square, having rows summing to zero.
pi0: array: A vector representing the probability distribution of the initial state of the Markov chain. If this vector is of length ‘r’, then the transition matrix ‘P’ is assumed to be ‘r x r’. The elements of this vector are assumed to be non-negative and sum to one.

Returns

A tuple, (tvec, xvec), where tvec is a vector of event times of length n and xvec is a vector of states, of length n+1.

Examples

>>> import smfsb
>>> import numpy as np
>>> smfsb.rcfmc(200, np.array([[-0.5,0.5],[1,-1]]), np.array([1,0]))

smfsb.sim.rdiff(a_fun, b_fun, x0=0, t=50, dt=0.01)

Simulate a sample path from a univariate diffusion process

This function simulates a single realisation from a time-homogeneous univariate diffusion process.

Parameters

a_fun: function: A scalar-valued function representing the infinitesimal mean (drift) of the diffusion process. The argument is the current state of the process.
b_fun: function: A scalar-valued function representing the infinitesimal standard deviation of the process. The argument is the current state of the process.
x0: float: The initial state of the diffusion process.
t: float: The length of the time interval over which the diffusion process is to be simulated.
dt: float: The step size to be used _both_ for the time step of the Euler integration method _and_ the recording interval for the output. It would probably be better to have separate parameters for these two things. Defaults to 0.01 time units.

Returns

A vector of states of the diffusion process on the required time grid.

Examples

>>> import smfsb
>>> import numpy as np
>>> smfsb.rdiff(lambda x: 1 - 0.1*x, lambda x: np.sqrt(1 + 0.1*x))

smfsb.sim.rfmc(n, p_mat, pi0)

Simulate a finite state space Markov chain

This function simulates a single realisation from a discrete time Markov chain having a finite state space based on a given transition matrix.

Parameters

n: int: The number of states to be sampled from the Markov chain, including the initial state, which will be sampled using ‘pi0’.
p_mat: matrix: The transition matrix of the Markov chain. This is assumed to be a stochastic matrix, having non-negative elements and rows summing to one.
pi0: array: A vector representing the probability distribution of the initial state of the Markov chain. If this vector is of length ‘r’, then the transition matrix ‘p_mat’ is assumed to be ‘r x r’. The elements of this vector are assumed to be non-negative and sum to one.

Returns

A numpy array containing the sampled values from the Markov chain.

Examples

>>> import smfsb
>>> import numpy as np
>>> p_mat = np.array([[0.9,0.1],[0.2,0.8]])
>>> pi0 = np.array([0.5,0.5])
>>> smfsb.rfmc(200, p_mat, pi0)

smfsb.sim.show_source(fun)

Print to console the source code of a function or method.

Called for the side-effect of printing the function source to standard output.

Parameters

fun: function or method: The function of interest

Returns

None

Examples

>>> import smfsb
>>> smfsb.show_source(smfsb.Spn.step_gillespie)

smfsb.sim.sim_sample(n, x0, t0, deltat, step_fun)

Simulate a many realisations of a model at a given fixed time in the future given an initial time and state, using a function (closure) for advancing the state of the model

This function simulates many realisations of a model at a given fixed time in the future given an initial time and state, using a function (closure) for advancing the state of the model , such as created by step_gillespie or step_euler.

Parameters

n: int: The number of samples required.
x0: array of numbers: The intial state of the system at time t0.
t0: float: The intial time to be associated with the initial state.
deltat: float: The amount of time in the future of t0 at which samples of the system state are required.
step_fun: function: A function (closure) for advancing the state of the process, such as produced by step_gillespie or step_euler.

Returns

A matrix with rows representing simulated states at time t0+deltat.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_gillespie()
>>> smfsb.sim_sample(10, [50, 100], 0, 30, stepLv)

smfsb.sim.sim_time_series(x0, t0, tt, dt, step_fun)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

This function simulates single realisation of a model on a regular grid of times using a function (closure) for advancing the state of the model, such as created by ‘step_gillespie’ or ‘step_euler’.

Parameters

x0: array of numbers: The intial state of the system at time t0
t0: float: This intial time to be associated with the intial state.
tt: float: The terminal time of the simulation.
dt: float: The time step of the output. Note that this time step relates only to the recorded output, and has no bearing on the accuracy of the simulation process.
step_fun: function: A function (closure) for advancing the state of the process, such as produced by ‘step_gillespie’ or ‘step_euler’.

Returns

A matrix with rows representing the state of the system at successive times.

Examples

>>> import smfsb.models
>>> lv = smfsb.models.lv()
>>> stepLv = lv.step_gillespie()
>>> smfsb.sim_time_series([50, 100], 0, 100, 0.1, stepLv)

smfsb.sim.simple_euler(rhs, ic, t=50, dt=0.001)

Simulate a sample path from an ODE model

This function integrates an Ordinary Differential Equation (ODE) model using a simple first order Euler method. The function is pedagogic and not intended for serious use. See scipy.integrate.solve_ivp for better, more robust ODE solvers.

Parameters

rhs: function: A vector-valued function representing the right hand side of the ODE model. The first argument is a vector representing the current state of the model, ‘x’. The second argument of ‘rhs’ is the current simulation time, ‘t’. In the case of a homogeneous ODE model, this argument will be unused within the function. The output of ‘rhs’ should be a vector of the same dimension as ‘x’.
ic: array: The initial conditions for the ODE model. This should be a vector of the same dimensions as the output from ‘rhs’, and the first argument of ‘rhs’.
t: float: The length of the time interval over which the ODE model is to be integrated. Defaults to 50 time units.
dt: float: The step size to be used both for the time step of the Euler integration method and the recording interval for the output. It would probably be better to have separate parameters for these two things. Defaults to 0.001 time units.

Returns

A matrix with rows representing the states at each time step.

Examples

>>> import smfsb
>>> import numpy as np
>>> smfsb.simple_euler(lambda x,t: 1-0.1*x[0], np.array([0]))

smfsb.sim.step_sde(drift, diffusion, dt=0.01)

Create a function for advancing the state of an SDE model by using a simple Euler-Maruyama integration method

This function creates a function for advancing the state of an SDE model using a simple Euler-Maruyama integration method. The resulting function (closure) can be used in conjunction with other functions (such as ‘sim_time_series’) for simulating realisations of SDE models.

Parameters

drift: function: A function representing the drift vector of the SDE model (corresponding roughly to the RHS of ante ODE model). ‘drift’ should have arguments x and t, with ‘x’ representing current system state and ‘t’ representing current system time. The value of the function should be a vector of the same dimension as ‘x’, representing the infinitesimal mean of the Ito SDE.
diffusion: function: A function representing the diffusion matrix of the SDE model (the square root of the infinitesimal variance matrix). ‘diffusion’ should have arguments x and t, with ‘x’ representing current system state and ‘t’ representing current system time. The value of the function should be a square matrix with both dimensions the same as the length of ‘x’.
dt: float: Time step to be used by the simple Euler-Maruyama integration method. Defaults to 0.01.

Returns

A function which can be used to advance the state of the SDE model with given drift vector and diffusion matrix, by using an Euler-Maruyama method with step size ‘dt’. The function closure returns a vector representing the simulated state of the system at the new time.

Examples

>>> import smfsb
>>> import numpy as np
>>> lamb = 2; alpha = 1; mu = 0.1; sig = 0.2
>>> def myDrift(x, t):
>>>     return np.array([lamb - x[0]*x[1],
>>>                      alpha*(mu - x[1])])
>>>
>>> def myDiff(x, t):
>>>     return np.array([[np.sqrt(lamb + x[0]*x[1]), 0],
>>>                      [0 ,sig*np.sqrt(x[1])]])
>>>
>>> stepProc = smfsb.step_sde(myDrift, myDiff, dt=0.001)
>>> smfsb.sim_time_series(np.array([1, 0.1]), 0, 30, 0.01, stepProc)

smfsb.spatial module

smfsb.spatial.sim_time_series_1d(x0, t0, tt, dt, step_fun, verb=False)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

This function simulates single realisation of a model on a 1D regular spatial grid and regular grid of times using a function (closure) for advancing the state of the model, such as created by step_gillespie_1d.

Parameters

x0array: The initial state of the process at time t0, a matrix with rows corresponding to reacting species and columns corresponding to spatial location.
t0float: The initial time to be associated with the initial state x0.
ttfloat: The terminal time of the simulation.
dtfloat: The time step of the output. Note that this time step relates only to the recorded output, and has no bearing on the accuracy of the simulation process.
step_funfunction: A function (closure) for advancing the state of the process, such as produced by step_gillespie_1d.
verbboolean: Output progress to the console (this function can be very slow).

Returns

A 3d array representing the simulated process. The dimensions are species, space, and time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv1d = lv.step_gillespie_1d(np.array([0.6,0.6]))
>>> N = 10
>>> T = 5
>>> x0 = np.zeros((2,N))
>>> x0[:,int(N/2)] = lv.m
>>> smfsb.sim_time_series_1d(x0, 0, T, 1, stepLv1d, True)

smfsb.spatial.sim_time_series_2d(x0, t0, tt, dt, step_fun, verb=False)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

This function simulates single realisation of a model on a 2D regular spatial grid and regular grid of times using a function (closure) for advancing the state of the model, such as created by step_gillespie_2d.

Parameters

x0array: The initial state of the process at time t0, a 3d array with dimensions corresponding to reacting species and then two corresponding to spatial location.
t0float: The initial time to be associated with the initial state x0.
ttfloat: The terminal time of the simulation.
dtfloat: The time step of the output. Note that this time step relates only to the recorded output, and has no bearing on the accuracy of the simulation process.
step_funfunction: A function (closure) for advancing the state of the process, such as produced by step_gillespie_2d.
verbboolean: Output progress to the console (this function can be very slow).

Returns

A 4d array representing the simulated process. The dimensions are species, two space, and time.

Examples

>>> import smfsb.models
>>> import numpy as np
>>> lv = smfsb.models.lv()
>>> stepLv2d = lv.step_gillespie_2d(np.array([0.6,0.6]))
>>> M = 10
>>> N = 15
>>> T = 5
>>> x0 = np.zeros((2,M,N))
>>> x0[:,int(M/2),int(N/2)] = lv.m
>>> smfsb.sim_time_series_2d(x0, 0, T, 1, stepLv2d, True)

smfsb.inference module

smfsb.inference.abc_run(n, rprior, rdist, verb=False)

Run a set of simulations initialised with parameters sampled from a given prior distribution, and compute statistics required for an ABC analaysis

Run a set of simulations initialised with parameters sampled from a given prior distribution, and compute statistics required for an ABC analaysis. Typically used to calculate “distances” of simulated synthetic data from observed data.

Parameters

nint: An integer representing the number of simulations to run.
rpriorfunction: A function without arguments generating a single parameter (vector) from prior distribution.
rdistfunction: A function taking a parameter (vector) as argument and returning the required statistic of interest. This will typically be computed by first using the parameter to run a forward model, then computing required summary statistics, then computing a distance. See the example for details.
verbboolean: Print progress information to console?

Returns

A tuple with first component a list of parameters and second component a list of corresponding distances.

Examples

>>> import smfsb
>>> import numpy as np
>>> import scipy as sp
>>> data = np.random.normal(5, 2, 250)
>>> def rpr():
>>>   return np.exp(np.random.uniform(-3, 3, 2))
>>>
>>> def rmod(th):
>>>   return np.random.normal(th[0], th[1], 250)
>>>
>>> def sumStats(dat):
>>>   return np.array([np.mean(dat), np.std(dat)])
>>>
>>> ssd = sumStats(data)
>>> def dist(ss):
>>>   diff = ss - ssd
>>>   return np.sqrt(np.sum(diff*diff))
>>>
>>> def rdis(th):
>>>   return dist(sumStats(rmod(th)))
>>>
>>> smfsb.abc_run(100, rpr, rdis)

smfsb.inference.abc_smc(n, rprior, dprior, rdist, rperturb, dperturb, factor=10, steps=15, verb=False, debug=False)

Run an ABC-SMC algorithm for infering the parameters of a forward model

Run an ABC-SMC algorithm for infering the parameters of a forward model. This sequential Monte Carlo algorithm often performs better than simple rejection-ABC in practice.

Parameters

nint: An integer representing the number of simulations to pass on at each stage of the SMC algorithm. Note that the TOTAL number of forward simulations required by the algorithm will be (roughly) ‘n*steps*factor’.
rpriorfunction: A function without arguments generating single parameter (vector) from the prior.
dpriorfunction: A function taking a parameter vector as argumnent and returning the log of the prior density.
rdistfunction: A function taking a parameter (vector) as argument and returning a scalar “distance” representing a measure of how good the chosen parameter is. This will typically be computed by first using the parameter to run a forward model, then computing required summary statistics, then computing a distance. See the example for details.
rperturbfunction: A function which takes a parameter as its argument and returns a perturbed parameter from an appropriate kernel.
dperturbfunction: A function which takes a pair of parameters as its first two arguments (new first and old second), and returns the log of the density associated with this perturbation kernel.
factorint: At each step of the algorithm, ‘n*factor’ proposals are generated and the best ‘n’ of these are weighted and passed on to the next stage. Note that the effective sample size of the parameters passed on to the next step may be (much) smaller than ‘n’, since some of the particles may be assigned small (or zero) weight.
stepsint: The number of steps of the ABC-SMC algorithm. Typically, somewhere between 5 and 100 steps seems to be used in practice.
verbboolean: Boolean indicating whether some progress should be printed to the console (the number of steps remaining).

Returns

A matrix with rows representing samples from the approximate posterior distribution.

Examples

>>> import smfsb
>>> import numpy as np
>>> import scipy as sp
>>> data = np.random.normal(5, 2, 250)
>>> def rpr():
>>>   return np.exp(np.random.uniform(-3, 3, 2))
>>>
>>> def rmod(th):
>>>   return np.random.normal(np.exp(th[0]), np.exp(th[1]), 250)
>>>
>>> def sumStats(dat):
>>>   return np.array([np.mean(dat), np.std(dat)])
>>>
>>> ssd = sumStats(data)
>>> def dist(ss):
>>>   diff = ss - ssd
>>>   return np.sqrt(np.sum(diff*diff))
>>>
>>> def rdis(th):
>>>   return dist(sumStats(rmod(th)))
>>>
>>> smfsb.abc_smc(100, rpr, lambda x: np.log(np.sum(((x<3)&(x>-3))/6)),
>>>                           rdis, lambda x: np.random.normal(x, 0.1),
>>>                           lambda x,y: np.sum(sp.stats.norm.logpdf(y, x, 0.1)))

smfsb.inference.abc_smc_step(dprior, prior_sample, prior_lw, rdist, rperturb, dperturb, factor)

Carry out one step of an ABC-SMC algorithm

Not meant to be directly called by users. See abc_smc.

smfsb.inference.metrop(n, alpha)

Run a simple Metropolis sampler with standard normal target and uniform innovations

This function runs a simple Metropolis sampler with standard normal target distribution and uniform innovations.

Parameters

nint: The number of iterations of the Metropolis sampler.
alpha: float: The tuning parameter of the sampler. The innovations of the sampelr are of the form U(-alpha, alpha).

Returns

A vector containing the output of the sampler.

Examples

>>> import smfsb
>>> smfsb.metrop(100, 1)

smfsb.inference.metropolis_hastings(init, log_lik, rprop, ldprop=<function <lambda>>, ldprior=<function <lambda>>, iters=10000, thin=10, verb=True, debug=False)

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution. Note that the algorithm carries over the old likelihood from the previous iteration, making it suitable for problems with expensive likelihoods, and also for “exact approximate” pseudo-marginal or particle marginal MH algorithms.

Parameters

initvector: A parameter vector with which to initialise the MCMC algorithm.
log_likfunction: A function which takes a parameter (such as init) as its only required argument and returns the log-likelihood of the data. Note that it is fine for this to return the log of an unbiased estimate of the likelihood, in which case the algorithm will be an “exact approximate” pseudo-marginal MH algorithm.
rpropstochastic function: A function which takes a parameter as its only required argument and returns a single sample from a proposal distribution.
ldpropfunction: A function which takes a new and old parameter as its first two required arguments and returns the log density of the new value conditional on the old. Defaults to a flat function which causes this term to drop out of the acceptance probability. It is fine to use the default for _any_ _symmetric_ proposal, since the term will also drop out for any symmetric proposal.
ldpriorfunction: A function which take a parameter as its only required argument and returns the log density of the parameter value under the prior. Defaults to a flat function which causes this term to drop out of the acceptance probability. People often use a flat prior when they are trying to be “uninformative” or “objective”, but this is slightly naive. In particular, what is “flat” is clearly dependent on the parametrisation of the model.
itersint: The number of MCMC iterations required (_after_ thinning).
thinint: The required thinning factor. eg. only store every thin iterations.
verbboolean: Boolean indicating whether some progress information should be printed to the console. Defaults to True.
debugboolean: Boolean indicating whether debugging information is required. Prints information about each iteration to console, to, eg., debug a crashing sampler.

Returns

A matrix with rows representing samples from the posterior distribution.

Examples

>>> import smfsb
>>> import numpy as np
>>> import scipy as sp
>>> data = np.random.normal(5, 2, 250)
>>> llik = lambda x: np.sum(sp.stats.norm.logpdf(data, x[0], x[1]))
>>> prop = lambda x: np.random.normal(x, 0.1, 2)
>>> smfsb.metropolis_hastings([1,1], llik, prop)

smfsb.inference.normal_gibbs(iters, n, a, b, c, d, xbar, ssquared)

A simple Gibbs sampler for Bayesian inference for the mean and precision of a normal random sample

This function runs a simple Gibbs sampler for the Bayesian posterior distribution of the mean and precision given a normal random sample.

Parameters

itersint: The number of iterations of the Gibbs sampler
nint: The sample size of the normal random sample
afloat: The shape parameter of the gamma prior on the sample precision.
bfloat: The scale parameter of the gamma prior on the sample precision.
cfloat: Th mean of the normal prior on the sample mean.
dfloat: The precision of the normal prior on the sample mean.
xbarfloat: The sample mean of the data.
ssquaredfloat: The sample variance of the data.

Returns

A matrix containing the samples of the Gibbs sampler in rows.

Examples

>>> import smfsb
>>> postmat = smfsb.normal_gibbs(iters=1100, n=15, a=3, b=11, c=10, d=1/100,
>>>   xbar=25, ssquared=20)
>>> postmat = postmat[range(100,1100),:]

smfsb.inference.pf_marginal_ll(n, sim_x0, t0, step_fun, data_ll, data, debug=False)

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set using a simple bootstrap particle filter.

Parameters

nint: An integer representing the number of particles to use in the particle filter.
sim_x0function: A function with arguments t0 and th, where ‘t0’ is a time at which to simulate from an initial distribution for the state of the particle filter and th is a vector of parameters. The return value should be a state vector randomly sampled from the prior distribution. The function therefore represents a prior distribution on the initial state of the Markov process.
t0float: The time corresponding to the starting point of the Markov process. Can be no bigger than the smallest observation time.
step_funfunction: A function for advancing the state of the Markov process, with arguments x, t0, deltat and th, with th representing a vector of parameters.
data_llfunction: A function with arguments x, t, y, th, where x and t represent the true state and time of the process, y is the observed data, and th is a parameter vector. The return value should be the log of the likelihood of the observation. The function therefore represents the observation model.
datamatrix: A matrix with first column an increasing set of times. The remaining columns represent the observed values of y at those times.

Returns

A function with single argument th, representing a parameter vector, which evaluates to the log of the particle filters unbiased estimate of the marginal likelihood of the data (for parameter th).

Examples

>>> import numpy as np
>>> import scipy as sp
>>> import smfsb
>>> def obsll(x, t, y, th):
>>>     return np.sum(sp.stats.norm.logpdf(y-x, scale=10)
>>>
>>> def simX(t0, th):
>>>     return np.array([np.random.poisson(50), np.random.poisson(100)])
>>>
>>> def step(x, t, dt, th):
>>>     sf = smfsb.models.lv(th).step_gillespie()
>>>     return sf(x, t, dt)
>>>
>>> mll = smfsb.pf_marginal_ll(80, simX, 0, step, obsll, smfsb.data.lv_noise_10)
>>> mll(np.array([1, 0.005, 0.6]))
>>> mll(np.array([2, 0.005, 0.6]))

smfsb.smfsbsbml module

smfsb.smfsbsbml.file_to_spn(filename, verb=False)

Convert an SBML model into a Spn object

Read a file containing a model in SBML and convert into an Spn object for simulation and analysis.

Parameters

filename: string: String name of file containing the model
verb: boolean: Output some debugging info

Returns

An Spn object

Examples

>>> import smfsb
>>> myMod = smfsb.file_to_spn("myModel.xml")
>>> step = myMod.step_gillespie()

smfsb.smfsbsbml.mod_to_spn(filename, verb=False)

Convert an SBML-shorthand model into a Spn object

Read a file containing a model in SBML-shorthand and convert into an Spn object for simulation and analysis.

Parameters

filename: string: String name of file containing the model
verb: boolean: Output some debugging info

Returns

An Spn object

Examples

>>> import smfsb
>>> myMod = smfsb.mod_to_spn("myModel.mod")
>>> step = myMod.step_gillespie()

smfsb.smfsbsbml.model_to_spn(m, verb=False)

Convert a libSBML model into a Spn object

Convert a libSBML model into a Spn object for simulation and analysis.

Parameters

m: model: A libsbml model (not document) object
verb: boolean: Output some debugging info

Returns

An Spn object

Examples

>>> import smfsb
>>> import libsbml
>>> d = libsbml.readSBML("myModel.xml")
>>> m = d.getModel()
>>> myMod = smfsb.model_to_spn(m)
>>> step = myMod.step_gillespie()

smfsb.smfsbsbml.shorthand_to_spn(sh_string, verb=False)

Convert an SBML-shorthand model string into a Spn object

Parse a string containing a model in SBML-shorthand and convert into an Spn object for simulation and analysis.

Parameters

sh_string: string: String containing the model
verb: boolean: Output some debugging info

Returns

An Spn object

Examples

>>> import smfsb
>>> file = open('myModel.mod', 'r')
>>> myModStr = file.read()
>>> file.close()
>>> myMod = smfsb.shorthand_to_spn(myModStr)
>>> step = myMod.step_gillespie()