model

model#

model

class pygom.model.CompileCanary(states=None)[source]#

Hold the need for (re-)compilation for various functions

A subclass of this should specify the states to watch

They may all be tripped to True using the trip() method An individual may be reset with the reset() method or with a direct assignment (they may not be tripped in this way).

reset(name)[source]#

Reset a canary

Parameters:

name (string) – the name of the canary to reset

Return type:

None

trip()[source]#

Trip all the canaries :rtype: None

pygom.model.DFE(ode, disease_state)[source]#

Returns the disease free equilibrium from an ode object

Parameters:

  • ode (BaseOdeModel) – a class object from pygom

  • diseaseState (array like) – name of the disease states

Returns:

e – disease free equilibrium

Return type:

array like

class pygom.model.DeterministicOde(state=None, param=None, derived_param=None, transition=None, event=None, birth_death=None, ode=None, backend='lambda')[source]#

This contains the interface and operation built above the already defined set of ode

Parameters:

  • state (list) – A list of states (string)

  • param (list) – A list of the parameters (string)

  • derived_param (list) – A list of the derived parameters (tuple of (string,string))

  • transition (list) – A list of transition (Transition)

  • birth_death (list) – A list of birth or death process (Transition)

  • ode (list) – A list of ode (Transition)

adjoint(state, t, state_param, func=None)[source]#

Compute the adjoint given the adjoint vector, time, state variable and the objective function. Note that this function is very restrictive in the sense that the (original) state variable changes through time but this assumes it is a constant, i.e. we assume that the original system is linear.

Parameters:

  • state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.

  • t (double) – The current time.

  • state_param (array like) – The state vector that is (or maybe) required to evaluate the jacobian of the original system

  • func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns:

output of the same length as the ode

Return type:

numpy.ndarray

Notes

The size of lambda should be the same as the state. The integral should be starting from T, the final time of the original system and is integrated backwards (for stability).

adjoint_T(t, state, state_param, func=None)[source]#

Same as adjoint() but with t as first parameter

adjoint_interpolate(state, t, interpolant, func=None)[source]#

Compute the adjoint given the adjoint vector, time, the functions which was used to interpolate the state variable

Parameters:

  • state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.

  • t (double) – The current time.

  • interpolant (list) – list of interpolating functions of the state

  • func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns:

output of the same length as the ode

Return type:

numpy.ndarray

adjoint_interpolate_T(t, state, interpolant, objInput=None)[source]#

Same as adjoint_interpolate() but with t as first parameter

adjoint_interpolate_jacobian(state, t, interpolant, func=None)[source]#

Compute the Jacobian of the adjoint given the adjoint vector, time, function of the interpolation on the state variables and the objective function. This is simply the same as the negative Jacobian of the ode transposed.

Parameters:

  • state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.

  • t (double) – The current time.

  • interpolant (list) – list of interpolating functions of the state

  • func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns:

output of is a two dimensional array of size [number of state x number of state]

Return type:

numpy.ndarray

Notes

Same as adjoint_jacobian() but takes a list of interpolating function instead of a single (vector) value

See also

adjoint_jacobian()

adjoint_interpolate_jacobian_T(t, state, interpolant, func=None)[source]#

Same as adjoint_interpolate_jacobian() but with t as first parameter

adjoint_jacobian(state, t, state_param, func=None)[source]#

Compute the jacobian of the adjoint given the adjoint vector, time, state variable and the objective function. This is simply the same as the negative jacobian of the ode transposed.

Parameters:

  • state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.

  • t (double) – The current time.

  • state_param (array like) – The state vector that is (or maybe) required to evaluate the jacobian of the original system

  • func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns:

output of is a two dimensional array of size [number of state x number of state]

Return type:

numpy.ndarray

Notes

It takes the same number of argument as the adjoint for simplicity when integrating.

See also

adjoint()

adjoint_jacobian_T(t, state, state_param, func=None)[source]#

Same as adjoint_jacobian_T() but with t being the first parameter

eval_forwardforward(FF, S, state, t)[source]#

Evaluate a single f(x) of the forward-forward sensitivities

Parameters:

  • FF (array like) – this is in fact a 3rd order Tensor, aka 3d array

  • S (array like) – sensitivities in matrix form

  • state (array like) – the current state

  • t (numeric) – time

Returns:

f(x) of size [number of state * (number of parameters * number of parameters)]

Return type:

numpy.ndarray

eval_sensitivity(S, t, state, by_state=False)[source]#

Evaluate the sensitivity given state and time

Parameters:

  • S (array like) – Which should be numpy.ndarray. The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.

  • t (double) – The current time

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

  • by_state (bool) – how we want the output to be arranged. Default is True so that we have a block diagonal structure

Return type:

numpy.ndarray

Notes

It is different to eval_ode() and eval_jacobian() in that the extra input argument is not a parameter

See also

sensitivity()

eval_sensitivityIV(S, IV, t, state)[source]#

Evaluate the sensitivity with initial values given state and time

Parameters:

  • S (array like) – Which should be numpy.ndarray. The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.

  • IV (array like) – sensitivities for the initial values

  • t (double) – The current time

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

Returns:

\(f(s(x,\theta))\) and \(f(s(x_{0}))\)

Return type:

numpy.ndarray

Notes

It is different to eval_ode() and eval_jacobian() in that the extra input argument is not a parameter.

See also

sensitivityIV()

forwardforward(ff, t, state, s)[source]#

Evaluate a single \(f(x)\) of the forward-forward sensitivities

Parameters:

  • ff (array like) – the forward-forward sensitivities in vector form

  • t (numeric) – time

  • state (array like) – the current state

  • s (array like) – forward sensitivities in vector form

Returns:

\(f(x)\) of size [number of state * (number of parameters * number of parameters)]

Return type:

numpy.ndarray

forwardforward_T(t, ff, s, state)[source]#

Same as forwardforward() but with t as the first parameter

get_transition_graph(file_name=None, show=True)[source]#

Returns the transition graph using graphviz

Parameters:

  • file_name (str, optional) – name of the output file, defaults to None

  • show (bool, optional) – If the graph should be plotted, defaults to True

Return type:

graphviz.Digraph

property initial_state#

Return the initial state values

property initial_time#

Return the initial time

property initial_values#

Returns the initial values, both time and state as a tuple (x0, t0)

integrate(t, full_output=False)[source]#

Integrate over a range of t when t is an array and a output at time t

Parameters:

  • t (array like) – the range of time points which we want to see the result of

  • full_output (bool) – if we want additional information

integrate2(t, full_output=False, method=None)[source]#

Integrate over a range of t when t is an array and a output at time t. Select a suitable method to integrate when method is None.

Parameters:

  • t (array like) – the range of time points which we want to see the result of

  • full_output (bool) – if we want additional information

  • method (str, optional) – the integration method. All those available in ode are allowed with ‘vode’ and ‘ivode’ representing the non-stiff and stiff version respectively. Defaults to None, which tries to choose the integration method via eigenvalue analysis (only one) using the initial conditions

is_stiff(state=None, t=None)[source]#

Test on the eigenvalues of the jacobian. We classify the problem as stiff if any of the eigenvalues are positive

Parameters:

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

  • t (double) – The current time

Returns:

eigenvalues of the system given input

Return type:

numpy.ndarray

jacobian_eigenvalue(state=None, t=None)[source]#

Find out the eigenvalues of the jacobian given state and time. If None is given, the initial values are used.

Parameters:

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

  • t (double) – The current time

Returns:

True if any eigenvalue is positive

Return type:

bool

linear_ode()[source]#

To check whether the input ode is linear

Returns:

True if it is linear, False otherwise

Return type:

bool

ode_and_forwardforward(state_param, t)[source]#

Evaluate a single f(x) of the ode and the forward-forward sensitivities

Parameters:

  • state_param (array like) – state and forward-forward sensitivities in vector form

  • t (numeric) – time

Returns:

same size as the state_param input

Return type:

numpy.ndarray

ode_and_forwardforward_T(t, state_param)[source]#

Same as odeAndForwardForward() but with time as the first input

ode_and_forwardforward_jacobian(state_param, t)[source]#

Return the jacobian after evaluation given the input of the state and the forward forward sensitivities

Parameters:

  • state_param (array like) – state and forward-forward sensitivities in vector form

  • t (numeric) – time

Returns:

size of (a,a) where a is the length of the state_param input

Return type:

numpy.ndarray

ode_and_forwardforward_jacobian_T(t, state_param)[source]#

Same as ode_and_forwardforward_jacobian() but with t being the first parameters

ode_and_sensitivity(state_param, t, by_state=False)[source]#

Evaluate the sensitivity given state and time

Parameters:

  • state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.

  • t (double) – The current time

  • by_state (bool) – Whether the output vector should be arranged by state or by parameters. If False, then it means that the vector of output is arranged according to looping i,j from Sensitivity_{i,j} with i being the state and j the param. This is the preferred way because it leds to a block diagonal Jacobian

Returns:

concatenation of 2 element. First contains the ode, second the sensitivity. Both are of type numpy.ndarray

Return type:

list

See also

sensitivity(), ode()

ode_and_sensitivityIV(state_param, t)[source]#

Evaluate the sensitivity given state and time

Parameters:

  • state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.

  • t (double) – The current time

Returns:

concatenation of 3 element. First contains the ode, second the sensitivity, then the sensitivity of the initial value. All of them are of type numpy.ndarray

Return type:

list

See also

sensitivity(), ode()

ode_and_sensitivityIV_T(t, state_param)[source]#

Same as ode_and_sensitivityIV() but with t as first parameter

ode_and_sensitivityIV_jacobian(state_param, t)[source]#

Evaluate the sensitivity given state and time. Output a block diagonal sparse matrix as default.

Parameters:

  • state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.

  • t (double) – The current time

  • byState (bool) – How the output is arranged, according to the vector of output. It can be in terms of state or parameters, where by state means that the jacobian is a block diagonal matrix.

Returns:

output of a square matrix of size: number of ode + 1 times number of parameters

Return type:

numpy.ndarray

See also

ode_and_sensitivity()

ode_and_sensitivityIV_jacobian_T(t, state_param)[source]#

Same as ode_and_sensitivityIV_jacobian() but with t as first parameter

ode_and_sensitivity_T(t, state_param, by_state=False)[source]#

Same as ode_and_sensitivity() but with t as first parameter

ode_and_sensitivity_jacobian(state_param, t, by_state=False)[source]#

Evaluate the sensitivity given state and time. Output a block diagonal sparse matrix as default.

Parameters:

  • state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.

  • t (double) – The current time

  • by_state (bool) – How the output is arranged, according to the vector of output. It can be in terms of state or parameters, where by state means that the jacobian is a block diagonal matrix.

Returns:

output of a square matrix of size: number of ode + 1 times number of parameters

Return type:

numpy.ndarray

See also

ode_and_sensitivity()

ode_and_sensitivity_jacobian_T(t, state_param, by_state=False)[source]#

Same as ode_and_sensitivity_jacobian() but with t as first parameter

plot()[source]#

Plot the results of the integration

Notes

If we have 3 states or more, it will always be arrange such that it has 3 columns. Uses the operation from odeutils

print_ode(latex_output=False)[source]#

Prints the ode in symbolic form onto the screen/console in actual symbols rather than the word of the symbol.

Parameters:

latex_output (bool, optional) – Defaults to false which prints the equation in terms of symbols, if set to yes then the formula in terms of latex equations will be printed onto the screen.

sens_jacobian_state(state_param, t)[source]#

Evaluate the jacobian of the sensitivity w.r.t. the state given state and time

Parameters:

  • state_param (array like) – The current numerical value for the states as well as the sensitivities, which can be numpy.ndarray or list

  • t (double) – The current time

Returns:

Matrix of dimension [number of state * number of parameters x number of state]

Return type:

numpy.ndarray

sensitivity(sens, t, state, by_state=False)[source]#

Evaluate the sensitivity given state and time. The default is to output the values by parameters, i.e. \(s_{i},\ldots,s_{i+n}\) are partial derivatives w.r.t. the states for \(i \in {1,1+p,1+2p,1+3p, \ldots, 1+(n-1)p}\). This is to take advantage of the fact that we have a block diagonal jacobian that was already evaluated

Parameters:

  • sens (array like) – The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.

  • t (double) – The current time

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

  • by_state (bool) – how we want the output to be arranged. Default is True so that we have a block diagonal structure

Return type:

numpy.ndarray

sensitivityIV(sensIV, t, state)[source]#

Evaluate the sensitivity which include the initial values as our parameters given state and time. The default is to output the values by parameters, i.e. \(s_{i},\ldots,s_{i+n}\) are partial derivatives w.r.t. the states for \(i \in {1,1+p,1+2p,1+3p, \ldots, 1+(n-1)p}\). This is to take advantage of the fact that we have a block diagonal Jacobian that was already evaluated.

Parameters:

  • sensIV (array like) – The starting sensitivity of size [number of state x number of parameters] + [number of state x number of state] for the initial condition. The latter is an identity matrix at time zero.

  • t (double) – The current time

  • state (array like) – The current numerical value for the states which can be numpy.ndarray or list

Returns:

output of the same length as the ode

Return type:

numpy.ndarray

sensitivityIV_T(t, sensIV, state)[source]#

Same as sensitivityIV() but with t as first parameter

sensitivity_T(t, sens, state, by_state=False)[source]#

Same as sensitivity() but with t as first parameter

class pygom.model.Event(transition_list, rate=None)[source]#

Class to contain transitions

pygom.model.FitzHugh(param=None)[source]#

The standard FitzHugh model without external input [FitzHugh1961]

\[\begin{split}\frac{dV}{dt} &= c ( V - \frac{V^{3}}{3} + R) \\ \frac{dR}{dt} &= -\frac{1}{c}(V - a + bR).\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.FitzHugh({'a':0.2, 'b':0.2, 'c':3.0})
>>> t = np.linspace(0, 20, 101)
>>> x0 = [1.0, -1.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
pygom.model.Influenza_SLIARD(param=None)[source]#

A simple influenza model from [Brauer2008], page 323.

\[\begin{split}\frac{dS}{dt} &= -S \beta (I + \delta A) \\ \frac{dL}{dt} &= S \beta (I + \delta A) - \kappa L \\ \frac{dI}{dt} &= p \kappa L - \alpha I \\ \frac{dA}{dt} &= (1 - p) \kappa L - \eta A \\ \frac{dR}{dt} &= f \alpha I + \eta A \\ \frac{dD}{dt} &= (1 - f) \alpha I\end{split}\]
exception pygom.model.InputError[source]#

As the name suggest… the type of input to the function or method is not of the expected types

pygom.model.Legrand_Ebola_SEIHFR(param=None)[source]#

The Legrand Ebola model [Legrand2007] with 6 compartments that includes the H = hospitalization and F = funeral state. Note that because this is an non-autonomous system, there are in fact a total of 7 states after conversion. The set of equations that describes the model are

\[\begin{split}\frac{dS}{dt} &= -(\beta_{I}SI + \beta_{H}SH + \beta_{F}SF) \\ \frac{dE}{dt} &= (\beta_{I}SI + \beta_{H}SH + \beta_{F}SF) - \alpha E \\ \frac{dI}{dt} &= \alpha E - (\gamma_{H} \theta_{1} + \gamma_{I}(1-\theta_{1})(1-\delta_{1}) + \gamma_{D}(1-\theta_{1})\delta_{1})I \\ \frac{dH}{dt} &= \gamma_{H}\theta_{1}I - (\gamma_{DH}\delta_{2} + \gamma_{IH}(1-\delta_{2}))H \\ \frac{dF}{dt} &= \gamma_{D}(1-\theta_{1})\delta_{1}I + \gamma_{DH}\delta_{2}H - \gamma_{F}F \\ \frac{dR}{dt} &= \gamma_{I}(1-\theta_{1})(1-\delta_{1})I + \gamma_{IH}(1-\delta_{2})H + \gamma_{F}F.\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> x0 = [1.0, 3.0/200000.0, 0.0, 0.0, 0.0, 0.0, 0.0]
>>> t = np.linspace(0, 25, 100)
>>> ode = common_models.Legrand_Ebola_SEIHFR([('beta_I',0.588),('beta_H',0.794),('beta_F',7.653),('omega_I',10.0/7.0),('omega_D',9.6/7.0),('omega_H',5.0/7.0),('omega_F',2.0/7.0),('alphaInv',7.0/7.0),('delta',0.81),('theta',0.80),('kappa',300.0),('interventionTime',7.0)])
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
pygom.model.Lorenz(param=None)[source]#

Lorenz attractor define by three parameters, \(\beta,\sigma,\rho\) as per [Lorenz1963].

\[\begin{split}\frac{dx}{dt} &= \sigma (y-x) \\ \frac{dy}{dt} &= x (\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z\end{split}\]

Examples

>>> import matplotlib.pyplot as plt
>>> import numpy
>>> from pygom import common_models
>>> t = numpy.linspace(0, 20, 101)
>>> params = {'beta':8.0/3.0, 'sigma':10.0, 'rho':28.0}
>>> ode = common_models.Lorenz(params)
>>> ode.initial_values = ([1., 1., 1.], t[0])
>>> solution = ode.integrate(t[1::])
>>> plt.plot(solution[:,0], solution[:,2])
>>> plt.show()
pygom.model.Lotka_Volterra(param=None)[source]#

Standard Lotka-Volterra model with two states and four parameters [Lotka1920]

\[\begin{split}\frac{dx}{dt} &= \alpha x - cxy \\ \frac{dy}{dt} &= -\delta y + \gamma xy\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> params = {'alpha':1, 'delta':3, 'c':2, 'gamma':6}
>>> ode = common_models.Lotka_Volterra(params)
>>> ode.initial_values = ([2.0, 6.0], 0)
>>> t = np.linspace(0.1, 100, 10000)
>>> ode.integrate(t)
>>> ode.plot()
class pygom.model.ODEVariable(ID: str, symbol: Union[None, Symbol] = None, value: Union[None, str] = None, units: Union[None, Quantity] = None, real: bool = False, limits: Union[None, tuple] = None)[source]#

A class that defines the variables in our ODE

Parameters:

  • ID (str) – identifier of the variable

  • name (str, optional) – name of the variable in human readable format. Defaults to None, which then takes the ID as the name

  • units (str, optional) – what unit the variable takes. Defaults to None.

  • real (bool, optional) – if the variable can only be a real number, defaults to True

class pygom.model.OrderedDict[source]#

Dictionary that remembers insertion order

clear() → None.  Remove all items from od.#
copy() → a shallow copy of od#
fromkeys(value=None)#

Create a new ordered dictionary with keys from iterable and values set to value.

items() → a set-like object providing a view on D's items#
keys() → a set-like object providing a view on D's keys#
move_to_end(key, last=True)#

Move an existing element to the end (or beginning if last is false).

Raise KeyError if the element does not exist.

pop(key[, default]) → v, remove specified key and return the corresponding value.#

If the key is not found, return the default if given; otherwise, raise a KeyError.

popitem(last=True)#

Remove and return a (key, value) pair from the dictionary.

Pairs are returned in LIFO order if last is true or FIFO order if false.

setdefault(key, default=None)#

Insert key with a value of default if key is not in the dictionary.

Return the value for key if key is in the dictionary, else default.

update([E, ]**F) → None.  Update D from dict/iterable E and F.#

If E is present and has a .keys() method, then does: for k in E: D[k] = E[k] If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v In either case, this is followed by: for k in F: D[k] = F[k]

values() → an object providing a view on D's values#
class pygom.model.ParameterStore[source]#

A class to store parameters of an ODE system

This is a specialised version of VariableStore which is able to handle values of a parameter that are draws from a stochatic distribution.

property has_stochastic_parameters: bool#

Simple check to see if there are any stochastic parameters in the store

new_realisation() → None[source]#

Generate a new realiasation of the parameters

set_value(variable, value)[source]#

Sets the value of a variable

property stochastic_parameters: dict[slice(<class 'str'>, <class 'scipy.stats._distn_infrastructure.rv_frozen'>, None)]#

Provides a dict of stochastic parameters (i.e. ones where the variable) has been defined as an instance of rv_frozen.

Return type:

Dict keyed on parameter name with value = the distribution

property values: list[float]#

Provides the values for the parameters

If there are stochastic parameters then a draw will be made and stored and returned on subsequent calls to this method.

To generate a new realisation call new_realisation.

Returns:

  • A list of numeric values.

  • For each element in the list if a parameter isstochastic then a new

  • value is drawn, if it is deterministic then the value is simply added.

class pygom.model.Quantity(name, abbrev=None, latex_repr=None, pretty_unicode_repr=None, pretty_ascii_repr=None, mathml_presentation_repr=None, is_prefixed=False, **assumptions)[source]#

Physical quantity: can be a unit of measure, a constant or a generic quantity.

property abbrev#

Symbol representing the unit name.

Prepend the abbreviation with the prefix symbol if it is defines.

convert_to(other, unit_system='SI')[source]#

Convert the quantity to another quantity of same dimensions.

Examples

>>> from sympy.physics.units import speed_of_light, meter, second
>>> speed_of_light
speed_of_light
>>> speed_of_light.convert_to(meter/second)
299792458*meter/second

>>> from sympy.physics.units import liter
>>> liter.convert_to(meter**3)
meter**3/1000
property free_symbols#

Return free symbols from quantity.

property is_prefixed#

Whether or not the quantity is prefixed. Eg. kilogram is prefixed, but gram is not.

property scale_factor#

Overall magnitude of the quantity as compared to the canonical units.

set_global_relative_scale_factor(scale_factor, reference_quantity)[source]#

Setting a scale factor that is valid across all unit system.

pygom.model.R0(ode, disease_state)[source]#

Returns the basic reproduction number, in symbolic form when the parameter values are not available

Parameters:

  • ode (BaseOdeModel) – a class object from pygom

  • diseaseStateIndex (array like) – name of the disease states

Returns:

e – R0

Return type:

array like

See also

getDiseaseProgressionMatrices(), getR0GivenMatrix()

pygom.model.R0_from_matrix(F, V, disease_state=None)[source]#

Returns the symbolic form of the basic reproduction number. This will include the states symbols which is different from getR0() where the states is replaced by the values of the disease-free equilibrium.

Parameters:

  • F (sympy.matrices.MatrixBase) – secondary infection rates

  • V (sympy.matrices.MatrixBase) – disease progression rates

  • disease_state (list like, optional) – list of the disease state as sympy.Symbol. Defaults to None which assumes that \(F,V\) had been differentiated

Returns:

e – the eigenvalues of \(FV^{-1}\) for the disease states

Return type:

sympy.matrices.MatrixBase

See also

getDiseaseProgressionMatrices(), getR0()

pygom.model.Robertson(param=None)[source]#

The so called Robertson problem [Robertson1966], which is a standard example used to test stiff integrator.

\[\begin{split}\frac{dy_{1}}{dt} &= -0.04 y_{1} + 1 \cdot 10^{4} y_{2} y_{3} \\ \frac{dy_{2}}{dt} &= 0.04 y_{1} - 1 \cdot 10^{4} y_{2} y_{3} - 3 \cdot 10^{7} y_{2}^{2}\\ \frac{dy_{3}}{dt} &= 3 \cdot 10^{7} y_{2}^{2}\end{split}\]

Examples

>>> from pygom import common_models
>>> import numpy
>>> t = numpy.append(0, 4*numpy.logspace(-6, 6, 1000))
>>> ode = common_models.Robertson()
>>> ode.initial_values = ([1.0,0.0,0.0], t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot() # note that this is not being plotted in the log scale
pygom.model.SEIR(param=None, init=None)[source]#

A standard SIR model [Brauer2008] with population N

pygom.model.SEIR_Birth_Death(param=None)[source]#

A standard SEIR model with birth and death [Aron1984], defined by the ode

\[\begin{split}\frac{dS}{dt} &= \mu - \beta SI - \mu S \\ \frac{dE}{dt} &= \beta SI - (\mu + \alpha) E \\ \frac{dI}{dt} &= \alpha E - (\mu + \gamma) I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

Uses the same set of parameters as the examples in SEIR() apart from \(\mu\) which is new.

>>> import numpy as np
>>> from pygom import common_models
>>> params = {'beta':1800, 'gamma':100, 'alpha':35.84, 'mu':0.02}
>>> ode = common_models.SEIR_Birth_Death(params)
>>> t = np.linspace(0, 50, 1001)
>>> x0 = [0.0658, 0.0007, 0.0002, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()

See also

SEIR()

pygom.model.SEIR_Birth_Death_Periodic(param=None)[source]#

A SEIR birth death model with periodic contact [Aron1984], defined by the ode

\[\begin{split}\frac{dS}{dt} &= \mu - \beta(t)SI - \mu S \\ \frac{dE}{dt} &= \beta(t)SI - (\mu + \alpha) E \\ \frac{dI}{dt} &= \alpha E - (\mu + \gamma) I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

where

\[\beta(t) = \beta_{0} (1 + \beta_{1} \cos(2 \pi t)).\]

An extension of an SEIR birth death model by varying the contact rate through time.

Examples

Uses the same set of parameters as the examples in SEIR_Birth_Death() but now we have two beta parameters instead of one.

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from pygom import common_models
>>> params = {'beta0':1800, 'beta1':0.2, 'gamma':100, 'alpha':35.84, 'mu':0.02}
>>> ode = common_models.SEIR_Birth_Death_Periodic(params)
>>> t = np.linspace(0, 50, 1001)
>>> x0 = [0.0658, 0.0007, 0.0002, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()
>>> plt.plot(np.log(solution[:,0]), np.log(solution[:,1]))
>>> plt.show()
>>> plt.plot(np.log(solution[:,0]), np.log(solution[:,2]))
>>> plt.show()

See also

SEIR(), SEIR_Birth_Death(), SIR_Periodic()

pygom.model.SEIR_Birth_Death_Periodic_Waning_Intro(param=None)[source]#

SEIR model with vital dynamics, periodic infectivity, immune waning and external introductions

pygom.model.SEIR_Multiple(n=2, param=None)[source]#

An SEIR model that describe spatial heterogeneity [Brauer2008], page 180. The model originated from [Lloyd1996] and notations used here follows [Brauer2008].

\[\begin{split}\frac{dS_{i}}{dt} &= dN_{i} - dS_{i} - \lambda_{i} S_{i} \\ \frac{dE_{i}}{dt} &= \lambda_{i}S_{i} - (d + \epsilon) E_{i} \\ \frac{dI_{i}}{dt} &= \epsilon E_{i} - (d + \gamma) I_{i} \\ \frac{dR_{i}}{dt} &= \gamma I_{i} - dR_{i}\end{split}\]

where

\[\lambda_{i} = \sum_{j=1}^{n} \beta_{i,j} I_{j} (1\{i \neq j\} p)\]

with \(n\) being the number of patch and \(p\) the coupled factor.

Examples

Use the initial conditions that were derived from the stationary condition specified in [Brauer2008].

>>> import numpy as np
>>> from pygom import common_models
>>> paramEval = {'beta_00':0.0010107, 'beta_01':0.0010107,
>>>              'beta_10':0.0010107, 'beta_11':0.0010107,
>>>              'd':0.02,'epsilon':45.6, 'gamma':73.0,
>>>              'N_0':10**6,'N_1':10**6,'p':0.01}
>>> x0 = [36139.3224081278, 422.560577637822,
>>>       263.883351688369, 963174.233662546]
>>> ode = common_models.SEIR_Multiple()
>>> t = np.linspace(0, 40, 100)
>>> x01 = []
>>> for s in x0:
>>>     x01 += [s]
>>>     x01 += [s]
>>> ode.parameters = paramEval
>>> ode.initial_values = (x01, t[0])
>>> solution, output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()
pygom.model.SIR(param=None)[source]#

A standard SIR model as per [Brauer2008]

\[\begin{split}\frac{dS}{dt} &= -\ \frac{beta SI}{N} \\ \frac{dI}{dt} &= \ \frac{beta SI}{N} - \gamma I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced top two graph in Figure 1.3 of the reference above. First, when everyone is susceptible and only one individual was infected.

>>> import numpy as np
>>> from pygom import common_models
>>> N=1e5
>>> ode = common_models.SIR({'beta':0.5, 'gamma':0.2, 'N':N})
>>> t = np.linspace(0, 730, 1001)
>>> i0=1
>>> x0 = [N-i0, i0, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
pygom.model.SIR_Birth_Death(param=None)[source]#

Extension of the standard SIR model [Brauer2008] to also include birth and death

\[\begin{split}\frac{dS}{dt} &= B -\beta SI - \mu S \\ \frac{dI}{dt} &= \beta SI - \gamma I - \mu I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced bottom graph in Figure 1.3 of the reference above.

>>> import numpy as np
>>> from pygom import common_models
>>> B = 126372.0/365.0
>>> N = 7781984.0
>>> params = {'beta':3.6, 'gamma':0.2, 'B':B/N, 'mu':B/N}
>>> ode = common_models.SIR_Birth_Death(params)
>>> t = np.linspace(0, 35*365, 10001)
>>> x0 = [0.065, 123.0*(5.0/30.0)/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()

See also

SIR()

pygom.model.SIR_norm(param=None)[source]#

A normalized SIR model:

\[\begin{split}\frac{dS}{dt} &= -\beta SI \\ \frac{dI}{dt} &= \beta SI - \gamma I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced top two graph in Figure 1.3 of the reference above. First, when everyone is susceptible and only one individual was infected.

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIR({'beta':3.6, 'gamma':0.2})
>>> t = np.linspace(0, 730, 1001)
>>> N = 7781984.0
>>> x0 = [1.0, 10/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

Second model with a more realistic scenario

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIR({'beta':3.6, 'gamma':0.2})
>>> t = np.linspace(0, 730, 1001)
>>> N = 7781984.0
>>> x0 = [0.065, 123*(5.0/30.0)/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
pygom.model.SIS(param=None)[source]#

Susceptible Infected Susceptible model

\[\begin{split}\frac{dS}{dt} &= -\beta SI + \gamma I \\ \frac{dI}{dt} &= \beta SI - \gamma I\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIS({'beta':0.5, 'gamma':0.2})
>>> t = np.linspace(0, 20, 100)
>>> x0 = [1.0, 0.1]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
pygom.model.SIS_Periodic(param=None)[source]#

A SIS model with periodic contact, defined by the ode as per [Hethcote1973]

\[\frac{dI}{dt} = (\beta(t)N - \alpha) I - \beta(t)I^{2}\]

where

\[\beta(t) = 2 - 1.8 \cos(5t).\]

As the name suggests, it achieves a (stable) periodic solution.

Examples

>>> from pygom import common_models
>>> import numpy as np
>>> ode = common_models.SIS_Periodic({'alpha':1.0})
>>> t = np.linspace(0, 10, 101)
>>> x0 = [0.1, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()
class pygom.model.SimulateOde(state=None, param=None, derived_param=None, transition=None, event=None, birth_death=None, ode=None, backend='lambda')[source]#

This builds on top of DeterministicOde which we simulate the outcome instead of solving it deterministically

Parameters:

  • state (list) – A list of states (string) or (string, (numeric, numeric)) if specifying limits

  • param (list) – A list of the parameters (string)

  • derived_param (list) – A list of the derived parameters (tuple of (string,string))

  • transition (list) – A list of transition (Transition) #TODO Now this might actually only be deterministic ODE objects. Check.

  • transition – A list of events (Event)

  • birth_death (list #TODO Now these are wrapped in events. Should try to make work for back compat.) – A list of birth or death process (Transition)

  • ode (list) – A list of ode (Transition)

get_bd_from_ode(A=None)[source]#

Returns a list of:class:Transition from this object by unrolling the odes. All the elements are of TransitionType.B or TransitionType.D

get_transitions_from_ode()[source]#

Returns a list of Transition from this object by unrolling the odes. All the elements are of TransitionType.T

get_unrolled_obj()[source]#

Returns a SimulateOde with the same state and parameters as the current object but with the equations defined by a set of transitions and birth death process instead of say, odes

plot(sim_X=None, sim_T=None)[source]#

Plot the results of a simulation

Takes the output of a function like solve_stochast

Parameters:

  • sim_X (list) – of length iteration each with (len(t),len(state)) if t is a vector, else it outputs unequal shape that was record of all the jumps

  • sim_T (list or numpy.ndarray) – if t is a single value, it outputs unequal shape that was record of all the jumps. if t is a vector, it outputs t so that it is a numpy.ndarray instead

Notes

If either sim_X or sim_T are None the this function will attempt to plot the deterministic ODE

If we have 3 states or more, it will always be arrange such that it has 3 columns. Uses the operation from odeutils

class pygom.model.StateStore[source]#

A class to store parameters of an ODE system

class pygom.model.Transition(origin=None, equation=None, transition_type=None, destination=None, magnitude='1', ID=None, name=None)[source]#

This class carries the information for transitions defined for an ode, which includes the ode itself, a birth death process where only one state is involved and also a transition between two states

Parameters:

  • origin (str) – Origin state.

  • equation (str) – Equation defining the transition

  • transition_type (enum or str, optional) – of type TransitionType or one of (‘ODE’, ‘T’, ‘B’, ‘D’) defaults to ‘ODE’

  • destination (str, optional) – Destination State. If the transition is not between state, such as a birth or death process, then this is is not required. If it is stated as a birth, death or an ode then it throws an error

property destination#

Return the destination state

Returns:

The destination state

Return type:

string

property equation#

Return the transition _equation

Returns:

The transition _equation

Return type:

string

is_between_state()[source]#

Return whether it is a transition between two state

Returns:

True if it is a transition between two state False if it is only related to the origin state

Return type:

bool

property origin#

Return the origin state

Returns:

The origin state

Return type:

string

property stochastic#

Return the secondary effects

Returns:

The destination state

Return type:

string

property transition_type#

Return the type of transition

Returns:

One of the four type available from transition_type

Return type:

transition_type

class pygom.model.TransitionType(value)[source]#

This is an Enum describing the four feasible type of transitions use to define the ode model BaseOdeModel

The following four types of transitions are available.

B = Birth process

D = Death process

T = Transition between states

ODE = ODE _equation

pygom.model.check_array_type(x, accept_booleans=False)[source]#

Check to see if the type of input is suitable. Only operate on one or two dimension arrays

Parameters:

  • x (array like) – which can be either a numpy.ndarray or list or tuple

  • accept_boolean (boolean) – If true boolean elements are accepted, else they are not.

Returns:

x – checked and converted array

Return type:

numpy.ndarray

pygom.model.check_dimension(x, y)[source]#

Compare the length of two array like objects. Converting both to a numpy array in the process if they are not already one.

Parameters:

  • x (array like) – first array

  • y (array like) – second array

Returns:

  • x (numpy.array) – checked and converted first array

  • y (numpy.array) – checked and converted second array

class pygom.model.compileCode(backend=None)[source]#

A class that compiles an algebraic expression in sympy to a faster numerical file using the appropriate backend.

compileExpr(inputSymb, inputExpr, backend=None, compileType=False)[source]#

Compiles the expression given the symbols. Determines the backend if required.

Parameters:

  • inputSymb (list) – the set of symbols for the input expression

  • inputExpr (expr) – expression in sympy

  • backend (optional) – the backend we want to use to compile

  • compileType (optional) – defaults to False. If True, return an extra output that informs the end user of the method used to compile the equation, can be one of (np, mpmath, sympy)

Return type:

Compiled function taking arguments of the input symbols

pygom.model.disease_progression_matrices(ode, disease_states, diff=True)[source]#

Returns (F,V), the secondary infection rates and disease progression rate respectively.

Parameters:

  • ode (BaseOdeModel) – an ode class in pygom

  • disease_states (array like) – the name of the disease states

  • diff (bool, optional) – if the first derivative of the matrices are return, defaults to true

Returns:

(F, V) – The progression matrices. If diff=False, then we return the \(F_{i}\) and \(V_{i}\) matrices as per [Brauer2008].

Return type:

tuple

pygom.model.is_list_like(x)[source]#

Test whether the input is a type that behaves like a list, such as (list,tuple,np.ndarray)

Parameters:

x – anything

Returns:

True if it belongs to one of the three expected type (list,tuple,np.ndarray)

Return type:

bool

pygom.model.plot_det(solution, t, stateList=None, y=None, yStateList=None)[source]#

Plot the results of the integration

Parameters:

  • solution (numpy.ndarray) – solution from the integration

  • t (array like) – the vector of time where the integration output correspond to

  • stateList (list) – name of the states, if available

  • Notes –

  • more (If we have 5 states or) –

  • such (it will always be arrange) –

  • columns. (that it has 3) –

pygom.model.plot_stoc(solution, t, stochastic_model)[source]#

Plot the results of a stocastic simulation

Parameters:

  • solution – results of the stochastic simulation

  • t (array like) – the vector of time where the integration output correspond to

  • stochastic_model – the model from which this simulation was generated

  • Notes –

  • more (If we have 5 states or) –

  • such (it will always be arrange) –

  • columns. (that it has 3) –

class pygom.model.shapeAdjust(numState, numParam, numTarget=None)[source]#

A class that change vector into matrix and vice versa for vectors used in DeterministicOde

Parameters:

  • numState (int) – number of states

  • numParam (int) – number of parameters

  • numTarget (int, optional) – number of targeted states, default assumes that this is the same as numState

kronParam(A, pre=False)[source]#

A sparse multiplication with an identity matrix of size equal to the number of parameters as initialized

Parameters:

  • A (array like) – a 2d array

  • pre (bool, optional) – If True, then returns \(I \otimes A\). If False then \(A \otimes I\), where \(A\) is the input matrix, \(I\) is the identity matrix and \(\otimes\) is the kron operator

kronState(A, pre=False)[source]#

A sparse multiplication with an identity matrix of size equal to the number of state as initialized

Parameters:

  • A (array like) – a 2d array

  • pre (bool, optional) – If True, then returns \(I \otimes A\). If False then \(A \otimes I\), where \(A\) is the input matrix, \(I\) is the identity matrix and \(\otimes\) is the kron operator

matToVecFF(FF)[source]#

Transforms the forward forward sensitivity matrix to a vector

matToVecSens(S)[source]#

Transforms the sensitivity matrix to a vector

vecToMatFF(ff)[source]#

Transforms the forward forward sensitivity vector to a matrix

vecToMatSens(s)[source]#

Transforms the sensitivity vector to a matrix

pygom.model.str_or_list(x)[source]#

Test to see whether input is a string or a list. If it is a string, then we convert it to a list.

Parameters:

x – str or list

Returns:

x in list form

Return type:

x

pygom.model.vanDerPol(param=None)[source]#

The van der Pol equation [vanderpol1926], a second order ode

\[y^{\prime\prime} - \mu (1-y^{2}) y^{\prime} + y = 0\]

where \(\mu > 0\). This can be converted to a first order ode by equating \(x = y^{\prime}\)

\[x^{\prime} - \mu (1 - y^{2}) x + y = 0\]

which result in a coupled ode

\[\begin{split}x^{\prime} &= \mu (1 - y^{2}) x - y \\ y^{\prime} &= x\end{split}\]

and this can be solved via standard method.

Examples

>>> from pygom import common_models
>>> import numpy
>>> t = numpy.linspace(0, 20, 1000)
>>> ode = common_models.vanDerPol({'mu':1.0})
>>> ode.initial_values = ([2.0,0.0], t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()