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