Integrator configuration

Integrator configuration¶

You can explicitly control which integrator is used for the simulations.
As all integrators have different settings, changing them is done on the integrator object itself, which can be found at Simulator().integrator.

By default, the first available of these integrators is used:

1. Assimulo (CVODE, installable from conda-forge)
2. Sundials (CVODE via the official sundials4py bindings, pip install mxlpy[sundials])
3. Scipy (LSODA)

Assimulo and Sundials run the same solver with the same default settings.
However, that package currenlty requires installation via conda-forge.
If mxlpy was installed from pypi, this package is not available and mxlpy will fall back to the Scipy package.
Thus, when setting integrator settings you should always check which integrator you are actually working with, as this otherwise might lead bugs on other systems.

In [1]:

from functools import partial

from example_models import get_linear_chain_2v
import mxlpy as mb
from mxlpy import Simulator
from mxlpy.integrators import Scipy

model = get_linear_chain_2v()

from functools import partial from example_models import get_linear_chain_2v import mxlpy as mb from mxlpy import Simulator from mxlpy.integrators import Scipy model = get_linear_chain_2v()

Scipy¶

In [2]:

# You can explicitly set an integrator for the simulation
s = Simulator(model, integrator=Scipy)

# You can also change the default settings of each integrator
s = Simulator(model, integrator=partial(Scipy, atol=1e-6, rtol=1e-6))
s.simulate(5)

# You can change integration settings mid simulation
# As not all integrators have the same settings, we recommend explicitly checking
# which integrator is currently set
if isinstance(s.integrator, Scipy):
    s.integrator.atol = 1e-8  # set absolute tolerance
    s.integrator.rtol = 1e-8  # set relative tolerance

_ = (s.simulate(10).get_result().unwrap_or_err()).variables.plot()

# You can explicitly set an integrator for the simulation s = Simulator(model, integrator=Scipy) # You can also change the default settings of each integrator s = Simulator(model, integrator=partial(Scipy, atol=1e-6, rtol=1e-6)) s.simulate(5) # You can change integration settings mid simulation # As not all integrators have the same settings, we recommend explicitly checking # which integrator is currently set if isinstance(s.integrator, Scipy): s.integrator.atol = 1e-8 # set absolute tolerance s.integrator.rtol = 1e-8 # set relative tolerance _ = (s.simulate(10).get_result().unwrap_or_err()).variables.plot()

Sundials¶

Sundials wraps CVODE from the SUNDIALS suite. It needs the optional sundials extra (pip install mxlpy[sundials]), which has wheels for Linux x86_64, macOS and Windows x64. Use method="BDF" (default) for stiff and method="Adams" for non-stiff problems.

In [3]:

from mxlpy import Sundials

s = Simulator(
    model, integrator=partial(Sundials, atol=1e-6, rtol=1e-6, method="BDF")
)
s.simulate(5)

_ = s.get_result().unwrap_or_err().variables.plot()

from mxlpy import Sundials s = Simulator( model, integrator=partial(Sundials, atol=1e-6, rtol=1e-6, method="BDF") ) s.simulate(5) _ = s.get_result().unwrap_or_err().variables.plot()

Scans, fits and other routines¶

Routines like scan.*, fit.*, mc.* and mca.* don't take a Simulator, but a runner, which describes how each model is simulated. A runner holds the integrator as well as the settings of the steady-state search, so you configure them in one place for every routine.

There is one runner per model builder: KineticRunner for the default Model, OdeRunner for an OdeModelBuilder and EvaluationRunner for a SteadyStateModelBuilder. If no runner is passed, the default runner of the model builder is used.

In [4]:

import numpy as np
import pandas as pd

from mxlpy import KineticRunner, scan

runner = KineticRunner(
    integrator=partial(Scipy, atol=1e-8, rtol=1e-8),
    tolerance=1e-8,  # steady-state tolerance
)

res = scan.steady_state(
    model,
    to_scan=pd.DataFrame({"k1": np.linspace(1, 2, 5)}),
    runner=runner,
)
res.variables

import numpy as np import pandas as pd from mxlpy import KineticRunner, scan runner = KineticRunner( integrator=partial(Scipy, atol=1e-8, rtol=1e-8), tolerance=1e-8, # steady-state tolerance ) res = scan.steady_state( model, to_scan=pd.DataFrame({"k1": np.linspace(1, 2, 5)}), runner=runner, ) res.variables

  0%|          | 0/5 [00:00<?, ?it/s]

 20%|██        | 1/5 [00:01<00:06,  1.57s/it]

100%|██████████| 5/5 [00:01<00:00,  3.12it/s]

Out[4]:

	x	y
k1
1.00	0.500	1.00
1.25	0.625	1.25
1.50	0.750	1.50
1.75	0.875	1.75
2.00	1.000	2.00

For behaviour the settings don't cover, subclass a runner and override the method you need. The custom runner is then used by every routine you pass it to.

In [5]:

from mxlpy.simulation import Simulation
from mxlpy.types import Result


class LoggingRunner(KineticRunner):
    def steady_state(self, model: mb.KineticModelBuilder) -> Result[Simulation]:
        res = KineticRunner.steady_state(self, model)
        if isinstance(res.value, Exception):
            print(f"Steady state failed: {res.value}")
        return res


res = scan.steady_state(
    model,
    to_scan=pd.DataFrame({"k1": np.linspace(1, 2, 5)}),
    runner=LoggingRunner(),
)

from mxlpy.simulation import Simulation from mxlpy.types import Result class LoggingRunner(KineticRunner): def steady_state(self, model: mb.KineticModelBuilder) -> Result[Simulation]: res = KineticRunner.steady_state(self, model) if isinstance(res.value, Exception): print(f"Steady state failed: {res.value}") return res res = scan.steady_state( model, to_scan=pd.DataFrame({"k1": np.linspace(1, 2, 5)}), runner=LoggingRunner(), )

  0%|          | 0/5 [00:00<?, ?it/s]

 20%|██        | 1/5 [00:01<00:07,  1.88s/it]

100%|██████████| 5/5 [00:01<00:00,  2.56it/s]

Diffrax¶

In [6]:

from diffrax import Tsit5

from mxlpy import Diffrax

model = get_linear_chain_2v()

# Use diffrax integrator suite
s = (
    Simulator(model, integrator=Diffrax)  # pyright: ignore[reportArgumentType]
    .simulate(5)
    .get_result()
    .unwrap_or_err()
).variables.plot()


# Explicit solver choice
s = (
    Simulator(model, integrator=partial(Diffrax, solver=Tsit5()))  # pyright: ignore[reportArgumentType]
    .simulate(5)
    .get_result()
    .unwrap_or_err()
).variables.plot()

from diffrax import Tsit5 from mxlpy import Diffrax model = get_linear_chain_2v() # Use diffrax integrator suite s = ( Simulator(model, integrator=Diffrax) # pyright: ignore[reportArgumentType] .simulate(5) .get_result() .unwrap_or_err() ).variables.plot() # Explicit solver choice s = ( Simulator(model, integrator=partial(Diffrax, solver=Tsit5())) # pyright: ignore[reportArgumentType] .simulate(5) .get_result() .unwrap_or_err() ).variables.plot()

INFO:jax._src.xla_bridge:Unable to initialize backend 'tpu': INTERNAL: Failed to open libtpu.so: libtpu.so: cannot open shared object file: No such file or directory

INFO:jax._src.xla_bridge:Unable to initialize backend 'proxy': The 'proxy' backend requires the address of an IFRT Proxy server. Set it via the JAX_BACKEND_TARGET environment variable or the jax_backend_target config option, e.g. JAX_BACKEND_TARGET=grpc://localhost:29000.

WARNING:jax._src.xla_bridge:An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.

In [ ]:

Native integrator¶

Native compiles a model and runs a C solver outside Python, which makes simulate typically 50-200x faster than with Scipy. Derived values, reactions and readouts are compiled too, on first read of the results, with the Python path as a fallback for any that can't be. It only works for models that can be represented symbolically, and raises NativeLoweringError for any other. It isn't installed by default, and there are two backends:

• pip install "mxlpy[native]": llvmlite and a compiled solver. The fastest, needing a wheel for your platform.
• pip install "mxlpy[native-wasm]": wasmtime only, more portable, about 1.2-3x slower than native.

from functools import partial
from mxlpy import Native, Simulator

Simulator(model, integrator=Native).simulate(100)
Simulator(model, integrator=partial(Native, rtol=1e-10, atol=1e-10)).simulate(100)

The first simulation of a model compiles it, usually in well under a second. Changing parameters doesn't recompile. Events aren't supported.

The default method is Radau, for stiff models. partial(Native, method="LSODA") uses LSODA instead, which switches between non-stiff and stiff methods by itself. It's the same C code Scipy(method="LSODA") runs, and on non-stiff models it can be faster than Radau. method="CVODE_BDF" (stiff) and method="CVODE_Adams" (non-stiff) use SUNDIALS' CVODE, the solver Assimulo runs, so they give the same results without needing conda.

partial(Native, reuse_jacobian=True) (Radau only) reuses Jacobians for longer. That can be about twice as fast on larger stiff models, but error can rise noticeably on some small models, and it's too rough for fits that use finite differences, so keep it off for fitting.