# Integrators¶

Contains integrator methods that will be used to integrate ODEs or SDEs.

class epipack.integrators.IntegrationMixin[source]

Bases: object

A helper MixIn that enables the base class to set initial conditions and integrate numerical ODEs.

Expects the base class to have the following methods and attributes:

• get_numerical_dydt()

• get_compartment_id()

• compartments

• N_comp

integrate(time_points, return_compartments=None, *args, **kwargs)[source]

Returns values of the given compartments at the demanded time points (as a dictionary).

Parameters
• time_points (np.ndarray) -- An array of time points at which the compartment values should be evaluated and returned.

• return_compartments (list, default = None) -- A list of compartments for which the result should be returned. If return_compartments is None, all compartments will be returned.

• integrator (str, default = 'dopri5') -- Which method to use for integration. Currently supported are 'euler' and 'dopri5'. If 'euler' is chosen, $$\delta t$$ will be determined by the difference of consecutive entries in time_points.

• adopt_final_state (bool, default = False) -- Whether or not to adopt the final state of the integration as new initial conditions.

• diffusion_constants (numpy.ndarray) -- Scalar and constant diffusion coefficients as prefactors for each compartment's Wiener process (has to be of same length as y0)

Returns

result -- Dictionary mapping compartments to time series.

Return type

dict

integrate_and_return_by_index(time_points, return_compartments=None, integrator='dopri5', adopt_final_state=False, diffusion_constants=None)[source]

Returns values of the given compartments at the demanded time points (as a numpy.ndarray of shape (return_compartments), len(time_points).

Parameters
• time_points (np.ndarray) -- An array of time points at which the compartment values should be evaluated and returned.

• return_compartments (list, default = None) -- A list of compartments for which the result should be returned. If return_compartments is None, all compartments will be returned.

• integrator (str, default = 'dopri5') -- Which method to use for integration. Currently supported are 'euler' and 'dopri5'. If 'euler' is chosen, $$\delta t$$ will be determined by the difference of consecutive entries in time_points.

• adopt_final_state (bool, default = False) -- Whether or not to adopt the final state of the integration

• diffusion_constants (numpy.ndarray) -- Scalar and constant diffusion coefficients as prefactors for each compartment's Wiener process (has to be of same length as y0)

Returns

result -- First axis accesses system states, second axis refers to the system state corresponding to the time points in time_points.

Return type

numpy.ndarray

integrate_and_return_by_index_until(t0, stop_condition, return_compartments=None, adopt_final_state=False, solve_ivp_kwargs={})[source]

Integrate the system until a stop condition is reached, using scipy.integrate.solve_ivp. Returns values of the given compartments at the demanded time points (as a numpy.ndarray of shape (return_compartments), len(time_points).

Parameters
• t0 (float) -- Initial time

• stop_condition (function) --

A function stop_condition(t,y) that flips from positive to negative values when a certain condition is reached. For instance, when an element reaches a certain threshold value:

stop_condition = lambda t, y: 0.3 - y[2]


• return_compartments (list, default = None) -- A list of compartments for which the result should be returned. If return_compartments is None, all compartments will be returned.

• adopt_final_state (bool, default = False) -- Whether or not to adopt the final state of the integration

• solve_ivp_kwargs (dict) -- Keyword arguments that will be passed to scipy.integrate.solve_ivp.

Returns

• t (float) -- Times when the system state was recorded The final time corresponds to the time when the stop condition was reached.

• y (numpy.ndarray) -- First axis accesses system states, second axis refers to the system state corresponding to the time points in t.

integrate_until(t0, stop_condition, return_compartments=None, adopt_final_state=False, solve_ivp_kwargs={})[source]

Integrate the system until a stop condition is reached, using scipy.integrate.solve_ivp.

Parameters
• t0 (float) -- Initial time

• stop_condition (function) --

A function stop_condition(t,y) that flips from positive to negative values when a certain condition is reached. For instance, when an element reaches a certain threshold value:

stop_condition = lambda t, y: 0.3 - y[2]


• return_compartments (list, default = None) -- A list of compartments for which the result should be returned. If return_compartments is None, all compartments will be returned.

• adopt_final_state (bool, default = False) -- Whether or not to adopt the final state of the integration

• solve_ivp_kwargs (dict) -- Keyword arguments that will be passed to scipy.integrate.solve_ivp.

Returns

• t (float) -- Times when the system state was recorded The final time corresponds to the time when the stop condition was reached.

• result (dict) -- Dictionary mapping compartments to time series.

set_initial_conditions(initial_conditions, initial_time=0.0, allow_nonzero_column_sums=False)[source]

Set the initial conditions for integration

Parameters
• initial_conditions (dict) -- A dictionary that maps compartments to a fraction of the population. Compartments that are not set in this dictionary are assumed to have an initial condition of zero.

• allow_nonzero_column_sums (bool, default = False) -- If True, an error is raised when the initial conditions do not sum to the population size.

epipack.integrators.integrate_SDE(dydt, t, y0, diffusion_constants, *args)[source]

Integrate an SDE system with Euler's method.

Parameters
• dydt (function) -- A function returning the momenta of the deterministic ODE system.

• t (numpy.ndarray of float) -- Array of time points at which the functions should be evaluated. Time steps must be equally spaced.

• y0 (numpy.ndarray) -- Initial conditions

• diffusion_constants (numpy.ndarray) --

Scalar and constant diffusion coefficients as prefactors for each compartment's Wiener process (has to be of same length as y0)

corresponds to $$D_i$$ in

$dY_i = f_i(\mathbf Y,t) dt + D_i dW_i$

• *args (list) -- List of parameters that will be passed to the momenta function.

epipack.integrators.integrate_dopri5(dydt, t, y0, *args)[source]

Integrate an ODE system with the Runge-Kutte Dormand Prince method (with step-size control).

Parameters
• dydt (function) -- A function returning the momenta of the ODE system.

• t (numpy.ndarray of float) -- Array of time points at which the functions should be evaluated.

• y0 (numpy.ndarray) -- Initial conditions

• *args (list) -- List of parameters that will be passed to the momenta function.

epipack.integrators.integrate_euler(dydt, t, y0, *args)[source]

Integrate an ODE system with Euler's method.

Parameters
• dydt (function) -- A function returning the momenta of the ODE system.

• t (numpy.ndarray of float) -- Array of time points at which the functions should be evaluated.

• y0 (numpy.ndarray) -- Initial conditions

• *args (list) -- List of parameters that will be passed to the momenta function.

epipack.integrators.integrate_ivp_until(dydt, t0, y0, stop_condition, solve_ivp_kwargs={})[source]

Integrate an ODE system with the scipy.integrate.solve_ivp method.

Parameters
• dydt (function) -- A function returning the momenta of the ODE system.

• t0 (float) -- Initial time

• y0 (numpy.ndarray) -- Initial conditions

• stop_condition (function) --

A function stop_condition(t,y) that flips from positive to negative values when a certain condition is reached. For instance, when an element reaches a certain threshold value:

stop_condition = lambda t, y: 0.3 - y[2]


• solve_ivp_kwargs (dict) -- Keyword arguments that will be passed to scipy.integrate.solve_ivp.

Returns

• t (float) -- Times when the system state was recorded The final time corresponds to the time when the stop condition was reached.

• y (numpy.ndarray) -- First axis accesses system states, second axis refers to the system state corresponding to the time points in t.

epipack.integrators.time_leap_ivp(t0, y0, get_event_rates, rand=None)[source]

Compute a time leap for time-varying rate functions based on Gillespie's SSA using Newton's method on a quad-integrator.

Parameters
• t0 (float) -- The current time

• y0 (numpy.ndarray) -- Current state of the system

• get_event_rates (function) -- A function that takes time t and state y as input and returns an array of corresponding event rates

• rand (float, default = None) -- A random number from the unit interval.

Returns

new_t -- A new time.

Return type

float

epipack.integrators.time_leap_newton(t0, y0, get_event_rates, rand=None)[source]

Compute a time leap for time-varying rate functions based on Gillespie's SSA using Newton's method on a quad-integrator.

Parameters
• t0 (float) -- The current time

• y0 (numpy.ndarray) -- Current state of the system

• get_event_rates (function) -- A function that takes time t and state y as input and returns an array of corresponding event rates

• rand (float, default = None) -- A random number from the unit interval.

Returns

new_t -- A new time.

Return type

float