odetoolbox package

odetoolbox.analysis(indict, disable_stiffness_check: bool = False, disable_analytic_solver: bool = False, disable_singularity_detection: bool = False, preserve_expressions: bool | Iterable[str] = False, log_level: str | int = 30) → List[Dict]

The main entry point of the ODE-toolbox API.

Parameters:

  • indict – Input dictionary for the analysis. For details, see https://ode-toolbox.readthedocs.io/en/main/#input

  • disable_stiffness_check – Whether to perform stiffness checking.

  • disable_analytic_solver – Set to True to return numerical solver recommendations, and no propagators, even for ODEs that are analytically tractable.

  • preserve_expressions – Set to True, or a list of strings corresponding to individual variable names, to disable internal rewriting of expressions, and return same output as input expression where possible. Only applies to variables specified as first-order differential equations.

  • log_level – Sets the logging threshold. Logging messages which are less severe than log_level will be ignored. Log levels can be provided as an integer or string, for example “INFO” (more messages) or “WARN” (fewer messages). For a list of valid logging levels, see https://docs.python.org/3/library/logging.html#logging-levels

Returns:

The result of the analysis. For details, see https://ode-toolbox.readthedocs.io/en/latest/index.html#output

Submodules

odetoolbox.analytic_integrator module

class odetoolbox.analytic_integrator.AnalyticIntegrator(solver_dict, spike_times: Dict[str, List[float]] | None = None, enable_caching: bool = True)

Bases: Integrator

Integrate a dynamical system by means of the propagators returned by ODE-toolbox.

__init__(solver_dict, spike_times: Dict[str, List[float]] | None = None, enable_caching: bool = True)
Parameters:

  • solve_dict – The results dictionary returned by a call to odetoolbox.analysis().

  • spike_times – For each variable, used as a key, the list of times at which a spike occurs.

  • enable_caching – Allow caching of results between requested times.

disable_cache_update()

Disallow caching of results between requested times.

enable_cache_update()

Allow caching of results between requested times.

get_all_variable_symbols()
get_value(t)

Get numerical solution of the dynamical system at time t.

Parameters:

t – The time to compute the solution for.

reset()

Reset time to zero and state to initial values.

set_initial_values(vals)

Set initial values, i.e. the state of the system at \(t = 0\). This will additionally cause the system state to be reset to \(t = 0\) and the new initial conditions.

Parameters:

vals – New initial values.

odetoolbox.config module

class odetoolbox.config.Config

Bases: object

Static class to store global configuration options.

Options are stored in the static dict config. Access using either Config().key or Config()[key] or Config.config[key].

config = {'differential_order_symbol': '__d', 'expression_simplification_threshold': 1000, 'forbidden_names': ['oo', 'zoo', 'nan', 'NaN', '__h'], 'input_time_symbol': 't', 'integration_accuracy_abs': 1e-06, 'integration_accuracy_rel': 1e-06, 'max_step_size': 999.0, 'output_timestep_symbol': '__h', 'propagators_prefix': '__P', 'sim_time': 0.1, 'simplify_expression': 'sympy.simplify(expr)'}
keys()

odetoolbox.dependency_graph_plotter module

class odetoolbox.dependency_graph_plotter.DependencyGraphPlotter

Bases: object

Use graphviz to plot a dependency graph between state variables.

classmethod plot_graph(shapes, dependency_edges, node_is_lin, fn=None)

Plot graph and write to file.

Parameters:

  • shapes – List of Shape instances.

  • dependency_edges – List of edges returned from dependency analysis.

  • node_is_lin – List of Booleans returned from dependency analysis.

  • fn – Filename to write PNG image as.

odetoolbox.integrator module

class odetoolbox.integrator.Integrator

Bases: object

Integrate a dynamical system by means of the propagators returned by ODE-toolbox (base class).

all_variable_symbols: List[Symbol] = []
get_sorted_spike_times()

Returns a global, sorted list of spike times.

Return all_spike_times:

A sorted list of all spike times for all variables.

Return all_spike_times_sym:

For the spike at time all_spike_times[i], the variables to which that spike applies are listed in all_spike_times_sym[i].

get_spike_times()

Get spike times.

Return spike_times:

For each variable, used as a key, the list of spike times associated with it.

set_spike_times(spike_times: Dict[str, List[float]] | None)

Internally converts to a global, sorted list of spike times.

Parameters:

spike_times – For each variable, used as a key, the list of spike times associated with it.

odetoolbox.mixed_integrator module

class odetoolbox.mixed_integrator.MixedIntegrator(numeric_integrator, system_of_shapes, shapes, analytic_solver_dict=None, parameters=None, spike_times=None, random_seed=123, max_step_size=inf, integration_accuracy_abs=1e-06, integration_accuracy_rel=1e-06, sim_time=1.0, alias_spikes=False, debug_plot_dir: str | None = None)

Bases: Integrator

Mixed numeric+analytic integrator. Supply with a result from ODE-toolbox analysis; calculates numeric approximation of the solution.

__init__(numeric_integrator, system_of_shapes, shapes, analytic_solver_dict=None, parameters=None, spike_times=None, random_seed=123, max_step_size=inf, integration_accuracy_abs=1e-06, integration_accuracy_rel=1e-06, sim_time=1.0, alias_spikes=False, debug_plot_dir: str | None = None)
Parameters:

  • numeric_integrator – A method from the GSL library for evolving ODEs, e.g. odeiv.step_rk4

  • system_of_shapes – Dynamical system to solve.

  • shapes – List of shapes in the dynamical system.

  • analytic_solver_dict – Analytic solver dictionary from ODE-toolbox analysis result.

  • parameters – Dictionary mapping parameter name (as string) to value expression.

  • spike_times – For each variable, used as a key, the list of times at which a spike occurs.

  • random_seed – Random number generator seed.

  • max_step_size – The maximum step size taken by the integrator.

  • integration_accuracy_abs – Absolute integration accuracy.

  • integration_accuracy_rel – Relative integration accuracy.

  • sim_time – How long to simulate for.

  • alias_spikes – Whether to alias spike times to the numerical integration grid. False means that precise integration will be used for spike times whenever possible. True means that after taking a timestep \(dt\) and arriving at \(t\), spikes from \(\langle t - dt, t]\) will only be processed at time \(t\).

  • debug_plot_dir – If given, enable debug plotting to this directory. If enabled, matplotlib is imported and used for plotting.

integrate_ode(initial_values=None, h_min_lower_bound=5e-09, raise_errors=True, debug=False)

This function computes the average step size and the minimal step size that a given integration method from GSL uses to evolve a certain system of ODEs during a certain simulation time, integration method from GSL and spike train.

Parameters:

  • initial_values – A dictionary mapping variable names (as strings) to initial value expressions.

  • h_min_lower_bound – The minimum acceptable step size.

  • raise_errors – Stop and raise exception when error occurs, or try to continue.

  • debug – Return extra values useful for debugging.

Returns:

Average and minimal step size, and elapsed wall clock time.

integrator_debug_plot(t_log, h_log, y_log, dir)
numerical_jacobian(t, y, params)

Compute the numerical values of the Jacobian matrix at the current time t and state y.

If the dynamics of variables \(x_1, \ldots, x_N\) is defined as \(x_i' = f_i\), then row \(i\) of the Jacobian matrix \(\mathbf{J}_i = \left[\begin{matrix}\frac{\partial f_i}{\partial x_0} & \cdots & \frac{\partial f_i}{\partial x_N}\end{matrix}\right]\).

Parameters:

  • t – Current time.

  • y – Current state vector of the dynamical system.

  • params – GSL parameters (not used here).

Returns:

Jacobian matrix \(\mathbf{J}\).

step(t, y, params)

“Stepping function”: compute the (numerical) value of the derivative of y over time, at the current time t and state y.

Parameters:

  • t – Current time.

  • y – Current state vector of the dynamical system.

  • params – GSL parameters (not used here).

Returns:

Updated state vector

exception odetoolbox.mixed_integrator.ParametersIncompleteException

Bases: Exception

Thrown in case not all parameters are assigned a numerical value before integration was attempted.

odetoolbox.plot_helper module

odetoolbox.plot_helper.import_matplotlib()

Try to import and configure matplotlib. Returns the “mpl” and “plt” packages if the import was successful, or return (None, None) if unsuccessful.

odetoolbox.shapes module

exception odetoolbox.shapes.MalformedInputException

Bases: Exception

Thrown in case an error occurred while processing malformed input.

class odetoolbox.shapes.Shape(symbol, order, initial_values, derivative_factors, inhom_term=0.0, nonlin_term=0.0, lower_bound=None, upper_bound=None)

Bases: object

This class provides a canonical representation of a shape function independently of the way in which the user specified the shape. It assumes a differential equation of the general form (where bracketed superscript \({}^{(n)}\) indicates the \(n\)-th derivative with respect to time):

\[x^{(n)} = N + \sum_{i=0}^{n-1} c_i x^{(i)}\]

Any constant or nonlinear part is here contained in the term N.

In the input and output, derivatives are indicated by adding one prime (single quotation mark) for each derivative order. For example, in the expression

x''' = c0 + c1*x + c2*x' + c3*x'' + x*y + x**2

the symbol of the ODE would be x (i.e. without any qualifiers), order would be 3, derivative_factors would contain the linear part in the form of the list [c1, c2, c3], the inhomogeneous term would be c0 and the nonlinear part x*y + x**2 is stored in nonlin_term.

__init__(symbol, order, initial_values, derivative_factors, inhom_term=0.0, nonlin_term=0.0, lower_bound=None, upper_bound=None)

Perform type and consistency checks and assign arguments to member variables.

Parameters:

  • symbol – Symbolic name of the shape without additional qualifiers like prime symbols or similar.

  • order – Order of the ODE representing the shape.

  • initial_values – Initial values of the ODE representing the shape. The dict contains order many key-value pairs: one for each derivative that occurs in the ODE. The keys are strings created by concatenating the variable symbol with as many single quotation marks (’) as the derivation order. The values are SymPy expressions.

  • derivative_factors – The factors for the derivatives that occur in the ODE. This list has to contain order many values, i.e. one for each derivative that occurs in the ODE. The values have to be in ascending order, i.e. [c1, c2, c3] for the given example.

  • inhom_term – Inhomogeneous part of the ODE representing the shape, i.e. c0 for the given example.

  • nonlin_term – Nonlinear part of the ODE representing the shape, i.e. x*y + x**2 for the given example.

__str__()

Return str(self).

classmethod from_function(symbol: str, definition, max_t=100, max_order=4, all_variable_symbols=None, debug=False) → Shape

Create a Shape object given a function of time.

Only functions of time that have a homogeneous ODE equivalent are supported (inhomogeneous ODE functions are not supported).

For a complete description of the algorithm, see https://ode-toolbox.readthedocs.io/en/latest/index.html#converting-direct-functions-of-time Note that this function uses sympy.simplify() rather than the custom simplification expression, because the latter risks that we fail to recognize certain shapes.

Parameters:

  • symbol – The variable name of the shape (e.g. "alpha", "I")

  • definition – The definition of the shape (e.g. "(e/tau_syn_in) * t * exp(-t/tau_syn_in)")

  • all_variable_symbols – All known variable symbols. None or list of string.

Returns:

The canonical representation of the postsynaptic shape

Example:

>>> Shape.from_function("I_in", "(e/tau) * t * exp(-t/tau)")
classmethod from_json(indict, all_variable_symbols=None, parameters=None, _debug=False)

Create a Shape instance from an input dictionary.

Parameters:

  • indict – Input dictionary, i.e. one element of the "dynamics" list supplied in the ODE-toolbox input dictionary.

  • all_variable_symbols – All known variable symbols. None or list of string.

  • parameters – An optional dictionary of parameters to their defining expressions.

classmethod from_ode(symbol: str, definition: str, initial_values: dict, all_variable_symbols=None, lower_bound=None, upper_bound=None, parameters=None, debug=False, **kwargs) → Shape

Create a Shape object given an ODE and initial values.

Note that shapes are only aware of their own state variables: if an equation for \(x\) depends on another state variable \(y\) of another shape, then \(y\) will be assumed to be a parameter and the term will appear in the inhomogeneous component of \(x\).

Parameters:

  • symbol – The symbol (variable name) of the ODE.

  • definition – The definition of the ODE.

  • initial_values – A dictionary mapping initial values to expressions.

  • all_variable_symbols – All known variable symbols. None or list of string.

Example:

>>> Shape.from_ode("alpha",
                  "-1/tau**2 * shape_alpha -2/tau * shape_alpha'",
                  {"alpha" : "0", "alpha'" : "e/tau", "0"]})
get_all_variable_symbols(shapes=None, derivative_symbol="'") → List[Symbol]

Get all variable symbols for this shape and all other shapes in shapes, without duplicates, in no particular order.

Returns:

all_symbols

get_initial_value(sym: str)

Get the initial value corresponding to the variable symbol.

Parameters:

sym – String representation of a sympy symbol, e.g. "V_m'"

get_state_variables(derivative_symbol="'") → List[Symbol]

Get all variable symbols for this shape, ordered according to derivative order, up to the shape’s order \(N\): [sym, dsym/dt, d^2sym/dt^2, ..., d^(N-1)sym/dt^(N-1)]

is_homogeneous() → bool
Returns:

False if and only if the shape has a nonzero right-hand side.

is_lin_const_coeff() → bool
Returns:

True if and only if the shape is linear and constant coefficient.

is_lin_const_coeff_in(symbols, parameters=None)
Returns:

True if and only if the shape is linear and constant coefficient in those variables passed in symbols.

reconstitute_expr() → Expr

Recreate right-hand side expression from internal representation (linear coefficients, inhomogeneous, and nonlinear parts).

static split_lin_inhom_nonlin(expr, x, parameters=None)

Split an expression into a linear, inhomogeneous and nonlinear part.

For example, in the expression

x''' = c0 + c1*x + c2*x' + c3*x'' + x*y + x**2

lin_factors would contain the linear part in the form of the list [c1, c2, c3], the inhomogeneous term would be c0 and the nonlinear part x*y + x**2 is returned as nonlin_term.

All parameters in parameters are assumed to be constants.

odetoolbox.singularity_detection module

class odetoolbox.singularity_detection.SingularityDetection

Bases: object

Singularity detection for generated propagator matrix.

Some ordinary differential equations (ODEs) can be solved analytically: an expression for the solution can be readily derived by algebraic manipulation. This allows us to formulate an “exact integrator”, that yields the next state of the system given the current state and the timestep Δt, to floating point (machine) precision [1].

In some cases, an ODE is analytically tractable, but vulnerable to an edge case condition in the generated propagator matrices. Consider the following example: Let the system of ODEs be given by

\[y' = A \cdot y\]

Then the propagator matrix for a timestep \(\Delta t\) is

\[P = \exp(A \cdot \Delta t)\]

which we can use to advance the system

\[y(t + \Delta t) = P \cdot y(t)\]

If \(A\) is of the form:

\[\begin{split}\begin{bmatrix} -a & 0 & 0\\ 1 & -a & 0\\ 0 & 1 & -b \end{bmatrix}\end{split}\]

Then the generated propagator matrix contains denominators that include the factor \(a - b\). When the parameters are chosen such that \(a = b\), a singularity (division by zero fault) occurs. However, the singularity is readily avoided if we assume that \(a = b\) before generating the propagator, i.e. we start out with the matrix

\[\begin{split}\begin{bmatrix} -a & 0 & 0\\ 1 & -a & 0\\ 0 & 1 & -a \end{bmatrix}\end{split}\]

The resulting propagator contains no singularities.

This class detects the potential occurrence of such singularities (potential division by zero) in the generated propagator matrix, which occur under certain choices of parameter values. These choices are reported as “conditions” by the find_singularities() function.

References

static find_inhomogeneous_singularities(expr) → Set[SymmetricEq]

Find singularities in the inhomogeneous part of the update equations.

Returns:

a set with equations, where the left-hand side of each equation is the variable that is to be subsituted, and the right-hand side is the expression to put in its place

Return type:

conditions

static find_propagator_singularities(P: MutableDenseMatrix, A: MutableDenseMatrix) → Set[SymmetricEq]

Find singularities in the propagator matrix \(P\) given the system matrix \(A\).

Parameters:

  • P (sympy.Matrix) – propagator matrix to check for singularities

  • A (sympy.Matrix) – system matrix

Returns:

a set with equations, where the left-hand side of each equation is the variable that is to be subsituted, and the right-hand side is the expression to put in its place

Return type:

conditions

static find_singularity_conditions_in_expression_(expr, symbols) → Set[SymmetricEq]
exception odetoolbox.singularity_detection.SingularityDetectionException

Bases: Exception

Thrown in case an error occurs while detecting singularities.

odetoolbox.spike_generator module

class odetoolbox.spike_generator.SpikeGenerator

Bases: object

classmethod spike_times_from_json(stimuli, sim_time) → Mapping[str, List[float]]

Read or generate spike times according to a JSON specification.

Returns:

spike_times: Each variable symbol is a key in the dictionary, and the list of spike times is the corresponding value. Symbol names use derivative_symbol to indicate differential order.

odetoolbox.stiffness module

class odetoolbox.stiffness.StiffnessTester(system_of_shapes, shapes, analytic_solver_dict=None, parameters=None, stimuli=None, random_seed=123, max_step_size=inf, integration_accuracy_abs=1e-06, integration_accuracy_rel=1e-06, sim_time=100.0, alias_spikes=False)

Bases: object

__init__(system_of_shapes, shapes, analytic_solver_dict=None, parameters=None, stimuli=None, random_seed=123, max_step_size=inf, integration_accuracy_abs=1e-06, integration_accuracy_rel=1e-06, sim_time=100.0, alias_spikes=False)
Parameters:

  • system_of_shapes – Dynamical system to solve.

  • shapes – List of shapes in the dynamical system.

  • analytic_solver_dict – Analytic solver dictionary from ODE-toolbox analysis result.

  • parameters – Dictionary mapping parameter name (as string) to value expression.

  • stimuli – Dictionary containing spiking stimuli.

  • random_seed – Random number generator seed.

  • max_step_size – The maximum step size taken by the integrator.

  • integration_accuracy_abs – Absolute integration accuracy.

  • integration_accuracy_rel – Relative integration accuracy.

  • sim_time – How long to simulate for.

  • alias_spikes – Whether to alias spike times to the numerical integration grid. False means that precise integration will be used for spike times whenever possible. True means that after taking a timestep \(dt\) and arriving at \(t\), spikes from \(\langle t - dt, t]\) will only be processed at time \(t\).

check_stiffness(raise_errors=False)

Perform stiffness testing: use implicit and explicit solvers to simulate the dynamical system, then decide which is the better solver to use.

For details, see https://ode-toolbox.readthedocs.io/en/latest/index.html#numeric-solver-selection-criteria

Returns:

Either "implicit", "explicit" or "warning".

Return type:

str

property random_seed

odetoolbox.sympy_helpers module

exception odetoolbox.sympy_helpers.NumericalIssueException

Bases: Exception

Thrown in case of numerical issues, like division by zero.

class odetoolbox.sympy_helpers.SymmetricEq(lhs, rhs, **options)

Bases: Equality

The sympy Eq class is by default not symmetric. This subclasses makes Eqs equivalent under symmetry. This helps prevent duplications when they are added to a set.

default_assumptions: ClassVar[StdFactKB] = {}
class odetoolbox.sympy_helpers.SympyPrinter(settings=None)

Bases: StrPrinter

odetoolbox.sympy_helpers.symbol_in_expression(list_of_symbols, expr) → bool

odetoolbox.system_of_shapes module

exception odetoolbox.system_of_shapes.GetBlockDiagonalException

Bases: Exception

Thrown in case an error occurs while block diagonalising a matrix.

exception odetoolbox.system_of_shapes.PropagatorGenerationException

Bases: Exception

Thrown in case an error occurs while generating propagators.

class odetoolbox.system_of_shapes.SystemOfShapes(x: MutableDenseMatrix, A: MutableDenseMatrix, b: MutableDenseMatrix, c: MutableDenseMatrix, shapes: List[Shape])

Bases: object

Represent a dynamical system in the canonical form \(\mathbf{x}' = \mathbf{Ax} + \mathbf{b} + \mathbf{c}\).

__init__(x: MutableDenseMatrix, A: MutableDenseMatrix, b: MutableDenseMatrix, c: MutableDenseMatrix, shapes: List[Shape])

Initialize a dynamical system in the canonical form \(\mathbf{x}' = \mathbf{Ax} + \mathbf{b} + \mathbf{c}\).

Parameters:

  • x – Vector containing variable symbols.

  • A – Matrix containing linear part.

  • b – Vector containing inhomogeneous part (constant term).

  • c – Vector containing nonlinear part.

classmethod from_shapes(shapes: List[Shape], parameters=None)

Construct the global system matrix \(\mathbf{A}\) and inhomogeneous part (constant term) \(\mathbf{b}\) and nonlinear part \(\mathbf{c}\) on the basis of all shapes in shapes, and such that

\[\mathbf{x}' = \mathbf{Ax} + \mathbf{b} + \mathbf{c}\]
generate_numeric_solver(state_variables=None)

Generate the symbolic expressions for numeric integration state updates; return as JSON.

generate_propagator_solver(disable_singularity_detection: bool = False)

Generate the propagator matrix and symbolic expressions for propagator-based updates; return as JSON.

get_connected_symbols(idx: int) → List[Symbol]

Extract all symbols belonging to a shape with symbol self.x_[idx] from the system matrix.

For example, if symbol i is x, and symbol j is y, and the system is:

\[\begin{split}\frac{dx}{dt} &= y\\ \frac{dy}{dt} &= y' = -\frac{1}{\tau^2} x - \frac{2}{\tau} y\end{split}\]

Then get_connected_symbols() for symbol x would return [x, y], and get_connected_symbols() for y would return the same.

get_dependency_edges()
get_initial_value(sym)
get_jacobian_matrix()

Get the Jacobian matrix as symbolic expressions. Entries in the matrix are sympy expressions.

If the dynamics of variables \(x_1, \ldots, x_N\) is defined as \(x_i' = f_i\), then row \(i\) of the Jacobian matrix \(\mathbf{J}_i = \left[\begin{matrix}\frac{\partial f_i}{\partial x_0} & \cdots & \frac{\partial f_i}{\partial x_N}\end{matrix}\right]\).

get_lin_cc_symbols(E, parameters=None)

Retrieve the variable symbols of those shapes that are linear and constant coefficient. In the case of a higher-order shape, will return all the variable symbols with "__d" suffixes up to the order of the shape.

get_shape_by_symbol(sym: str) → Shape | None
get_sub_system(symbols)

Return a new SystemOfShapes instance which discards all symbols and equations except for those in symbols. This is probably only sensible when the elements in symbols do not dependend on any of the other symbols that will be thrown away.

propagate_lin_cc_judgements(node_is_lin, E)

Propagate: if a node depends on a node that is not linear and constant coefficient, it cannot be linear and constant coefficient.

Parameters:

  • node_is_lin – Initial assumption about whether node is linear and constant coefficient.

  • E – List of edges returned from dependency analysis.

reconstitute_expr(state_variables=None)

Reconstitute a sympy expression from a system of shapes (which is internally encoded in the form \(\mathbf{x}' = \mathbf{Ax} + \mathbf{b} + \mathbf{c}\)).

Before returning, the expression is simplified using a custom series of steps, passed via the simplify_expression argument (see the ODE-toolbox documentation for more details).

shape_order_from_system_matrix(idx: int) → int

Determine shape differential order from system matrix of symbol self.x_[idx]

odetoolbox.system_of_shapes.get_block_diagonal_blocks(A)