API

CoupledSDEs <: ContinuousTimeDynamicalSystem
CoupledSDEs(f, u0 [, p]; kwargs...)

A stochastic continuous time dynamical system defined by a set of coupled stochastic differential equations (SDEs) as follows:

\[\text{d}\mathbf{u} = \mathbf{f}(\mathbf{u}, p, t) \text{d}t + \mathbf{g}(\mathbf{u}, p, t) \text{d}\mathcal{N}_t\]

where $\mathbf{u}(t)$ is the state vector at time $t$, $\mathbf{f}$ describes the deterministic dynamics, and the noise term $\mathbf{g}(\mathbf{u}, p, t) \text{d}\mathcal{N}_t$ describes the stochastic forcing in terms of a noise function (or diffusion function) $\mathbf{g}$ and a noise process $\mathcal{N}_t$. The parameters of the functions $\mathcal{f}$ and $\mathcal{g}$ are contained in the vector $p$.

There are multiple ways to construct a CoupledSDEs depending on the type of stochastic forcing.

The only required positional arguments are the deterministic dynamic rule f(u, p, t), the initial state u0, and optinally the parameter container p (by default p = nothing). For construction instructions regarding f, u0 see the DynamicalSystems.jl tutorial .

By default, the noise term is standard Brownian motion, i.e. additive Gaussian white noise with identity covariance matrix. To construct different noise structures, see below.

Noise term

The noise term can be specified via several keyword arguments. Based on these keyword arguments, the noise function g is constructed behind the scenes unless explicitly given.

• The noise strength (i.e. the magnitude of the stochastic forcing) can be scaled with noise_strength (defaults to 1.0). This factor is multiplied with the whole noise term.
• For non-diagonal and correlated noise, a covariance matrix can be provided via covariance (defaults to identity matrix of size length(u0).)
• For more complicated noise structures, including state- and time-dependent noise, the noise function g can be provided explicitly as a keyword argument (defaults to nothing). For construction instructions, continue reading.

The function g interfaces to the diffusion function specified in an SDEProblem of DynamicalSystems.jl. g must follow the same syntax as f, i.e., g(u, p, t) for out-of-place (oop) and g!(du, u, p, t) for in-place (iip).

Unless g is of vector form and describes diagonal noise, a prototype type instance for the output of g must be specified via the keyword argument noise_prototype. It can be of any type A that has the method LinearAlgebra.mul!(Y, A, B) -> Y defined. Commonly, this is a matrix or sparse matrix. If this is not given, it defaults to nothing, which means the g should be interpreted as being diagonal.

The noise process can be specified via noise_process. It defaults to a standard Wiener process (Gaussian white noise). For details on defining noise processes, see the docs of DiffEqNoiseProcess.jl . A complete list of the pre-defined processes can be found here. Note that DiffEqNoiseProcess.jl also has an interface for defining custom noise processes.

By combining g and noise_process, you can define different types of stochastic systems. Examples of different types of stochastic systems are listed on the StochasticDiffEq.jl tutorial page. A quick overview of common types of stochastic systems can also be found in the online docs for CoupledSDEs.

Keyword arguments

• g: noise function (default nothing)
• noise_strength: scaling factor for noise strength (default 1.0)
• covariance: noise covariance matrix (default nothing)
• noise_prototype: prototype instance for the output of g (default nothing)
• noise_process: stochastic process as provided by DiffEqNoiseProcess.jl (default nothing, i.e. standard Wiener process)
• t0: initial time (default 0.0)
• diffeq: DiffEq solver settings (see below)
• seed: random number seed (default rand(UInt64), so each CoupledSDEs produces a different noise realization unless a seed is given explicitly)

DifferentialEquations.jl interfacing

The CoupledSDEs is evolved using the solvers of DifferentialEquations.jl. To specify a solver via the diffeq keyword argument, use the flag alg, which can be accessed after loading StochasticDiffEq.jl (using StochasticDiffEq). The default diffeq is:

(alg = SOSRA(), abstol = 1.0e-2, reltol = 1.0e-2)

diffeq keywords can also include a callback for event handling .

Dev note: CoupledSDEs is a light wrapper of SDEIntegrator from StochasticDiffEq.jl. The integrator is available as the field integ, and the SDEProblem is integ.sol.prob. The convenience syntax SDEProblem(ds::CoupledSDEs, tspan = (t0, Inf)) is available to extract the problem.

Converting between CoupledSDEs and CoupledODEs

You can convert a CoupledSDEs system to CoupledODEs to analyze its deterministic part using the function CoupledODEs(ds::CoupledSDEs; diffeq, t0). Similarly, use CoupledSDEs(ds::CoupledODEs [, p]; kwargs...) to convert a CoupledODEs into a CoupledSDEs.

CoupledSDEs(ds::CoupledODEs, p = current_parameters(ds); kwargs...)

Converts a CoupledODEs into a CoupledSDEs. While p is optional, it is likely to change as the diffusion (noise) function g is added.

CoupledODEs(ds::CoupledSDEs; kwargs...)

Converts a CoupledSDEs into a CoupledODEs by extracting the deterministic part of ds.

FreidlinWentzellHamiltonian{IIP, D, ...}

Freidlin-Wentzell Hamiltonian for the small-noise SDE $\mathrm{d}X_t = b(X_t)\,\mathrm{d}t + \sigma\,\Sigma(X_t)\,\mathrm{d}W_t$, with diffusion tensor $a(x) = \Sigma(x)\Sigma(x)^\top$:

\[H(x, p) \;=\; \langle b(x),\, p\rangle \;+\; \tfrac{1}{2}\,\langle p,\, a(x)\,p\rangle.\]

p is the conjugate momentum (Legendre dual of $\dot x$); Hamilton's equations are $\dot x = b(x) + a(x)\,p$ and $\dot p = -\partial_x b^\top p - \tfrac{1}{2}\langle p, \partial_x a\,p\rangle$. Freidlin-Wentzell instantons (action minimizers between invariant sets) live on the zero-energy shell $H \equiv 0$; the simplified geometric MAM ((Grafke et al., 2017)) minimizes the action directly in (x, p) space.

The stored a(x) is trace-normalized (internal helper _trace_normalized_a) so that the action is invariant to the overall scale of noise_strength.

Fields

• H_x, H_p: callables (x, p) -> ∂_x H, (x, p) -> ∂_p H, returning D × N matrices.
• a: trace-normalized diffusion callable; a Base.Returns for constant noise.
• x_ref: reference state used by Base.show; nothing if not supplied.

Constructors

• FreidlinWentzellHamiltonian(ds::ContinuousTimeDynamicalSystem): builds H_x, H_p, a from ds. Central finite differences are used for $\partial_x a$ when a is state-dependent. Validation and diagonal-vs-coupled classification happen at cache build (internal helper build_sgmam_cache); rank-deficient a is rejected there.
• FreidlinWentzellHamiltonian{IIP, D}(H_x, H_p; a = Returns(Diagonal(ones(D))), x_ref = nothing): for hand-rolled Hamiltonians; you are responsible for matching the convention above.

See (Freidlin and Wentzell, 1998) for the underlying theory.

solver(ds::ContinuousTimeDynamicalSystem) -> Any

Returns the SDE solver specified in the diffeq settings of the CoupledSDEs.

drift(
    sys::Union{CoupledODEs{IIP}, CoupledSDEs{IIP}},
    x;
    t
) -> Any

Returns the deterministic drift $f(x)$ of the CoupledSDEs sys at state x. For time-dependent systems, the time can be specified as a keyword argument t (by default t=0).

div_drift(sys::ContinuousTimeDynamicalSystem, x) -> Any
div_drift(sys::ContinuousTimeDynamicalSystem, x, t) -> Any

Computes the divergence of the drift field $f(x)$ at state x. For time- dependent systems, the time can be specified as a keyword argument t (by default t=0).

covariance_matrix(ds::CoupledSDEs)

Returns the covariance matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.

See also diffusion_matrix.

diffusion_matrix(ds::CoupledSDEs)

Returns the diffusion matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.

Note: The diffusion matrix $Σ$ is the square root of the noise covariance matrix $Q$ (see covariance_matrix), defined via the Cholesky decomposition $Q = Σ Σ^\top$.

noise_process(sys::CoupledSDEs) -> Any

Fetches the stochastic process $\mathcal{N}$ specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.

noise_strength(sys::CoupledSDEs) -> Any

Effective noise strength $\sigma_{\mathrm{eff}} = \sqrt{\mathrm{tr}(\mathbf{\Sigma})/D}$ of sys, where $\mathbf{\Sigma} =$ covariance_matrix(sys) and $D =$ dimension(sys). Together with normalize_covariance! it satisfies $\sigma_{\mathrm{eff}}^2 \cdot \mathbf{Q}_{\mathrm{can}} = \mathbf{\Sigma}$.

For an SDE built as $\mathrm{d}\mathbf{x} = \mathbf{b}\,\mathrm{d}t + \sigma\sqrt{\mathbf{Q}}\,\mathrm{d}\mathbf{W}$:

• Isotropic covariance = c·I: $\sigma_{\mathrm{eff}} = \sigma\sqrt{c}$. The default covariance = I ($c=1$) recovers the construction-time noise_strength keyword exactly.
• Anisotropic Q: $\sigma_{\mathrm{eff}} = \sigma\sqrt{\mathrm{tr}(\mathbf{Q})/D}$, the per-direction average noise amplitude. Equals $\sigma$ whenever $\mathrm{tr}(\mathbf{Q}) = D$.
• Custom g (not of the form $\sigma\sqrt{\mathbf{Q}}$): still well-defined as the trace-based effective σ of the resulting $\mathbf{\Sigma}$.

See Large deviation theory for the role of $\sigma_{\mathrm{eff}}$ in the action functionals.

normalize_covariance!(covariance) -> Any

Returns $\mathbf{Q}\cdot D/\mathrm{tr}(\mathbf{Q})$, the trace-normalized covariance ($\mathrm{tr} = D$, average eigenvalue 1).

For an SDE built as $\mathrm{d}\mathbf{x} = \mathbf{b}\,\mathrm{d}t + \sigma\sqrt{\mathbf{Q}}\,\mathrm{d}\mathbf{W}$ the SDE only fixes the product $\sigma^2\mathbf{Q} =$ covariance_matrix(sys); this function picks the canonical pair $(\sigma_{\mathrm{eff}}, \mathbf{Q}_{\mathrm{can}})$ defined by $\mathrm{tr}(\mathbf{Q}_{\mathrm{can}}) = D$ (see noise_strength for the matching $\sigma_{\mathrm{eff}}$). Used internally by fw_action, om_action, and geometric_action. For diagonal positive $\mathbf{Q}$ this coincides numerically with dividing by $L_1(\mathbf{Q})/D$.

See Large deviation theory for the rotation-invariance and conversion factors.

deterministic_orbit(
    sys::CoupledSDEs,
    T;
    ...
) -> Tuple{Any, Any}
deterministic_orbit(
    sys::CoupledSDEs,
    T,
    init;
    diffeq,
    kwargs...
) -> Tuple{Any, Any}

Simulates the deterministic (noise-free) dynamics of CoupledSDEs sys in time for a duration T, starting at initial condition init.

This function is equivalent to calling trajectory on the deterministic part of the CoupledSDEs (with noise_strength=0). It works with the ODE solvers used for CoupledODEs.

Keyword arguments

• diffeq=(alg=Tsit5(), abstol = 1e-6, reltol = 1e-6): ODE solver settings (see CoupledODEs)
• kwargs...: keyword arguments passed to trajectory

For more info, see ODEProblem. For stochastic integration, see trajectory.

transition(
    sys::CoupledSDEs,
    x_i,
    x_f;
    radii,
    tmax,
    radius_directions,
    cut_start,
    seed,
    kwargs...
) -> Tuple{Any, Any, Bool}

Generates a sample transition from point x_i to point x_f.

This function simulates sys in time, starting from initial condition x_i, until entering a ball of given radius around x_f. By default the solver is initialized with the seed stored on sys; pass seed to override it.

Keyword arguments

• radii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively
• tmax=1e3: maximum time until the simulation stops even x_f has not been reached
• radius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.
• cut_start=true: if false, returns the whole trajectory up to the transition
• seed=nothing: if given, override the solver seed for this call (else use sys's seed)
• kwargs...: keyword arguments passed to CommonSolve.solve. These take precedence over the options stored in sys.diffeq, which are forwarded by default.

Output

[path, times, success]

• path (Matrix): transition path (size [dim × N], where N is the number of time points)
• times (Vector): time values (since start of simulation) of the path points (size N)
• success (bool): if true, a transition occurred (i.e. the ball around x_f has been reached), else false

See also transitions, trajectory.

transitions(
    sys::CoupledSDEs,
    x_i,
    x_f;
    ...
) -> CriticalTransitions.TransitionEnsemble{_A, _B, EnsembleSolution{T, N, S}} where {_A, _B, T, N, S}
transitions(
    sys::CoupledSDEs,
    x_i,
    x_f,
    N::Int64;
    radii,
    tmax,
    Nmax,
    cut_start,
    radius_directions,
    show_progress,
    EnsembleAlg,
    seed,
    kwargs...
) -> CriticalTransitions.TransitionEnsemble{_A, _B, EnsembleSolution{T, N, S}} where {_A, _B, T, N, S}

Generates an ensemble of N transition samples of sys from point x_i to point x_f. The transitions is by default simulated using threading. To sample the transitions in serial, GPU or Distributed enverionment, pass the desired SciMLBase.EnsembleAlgorithm to the EnsembleAlg algorithm.

Keyword arguments

• radii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively
• Nmax: number of attempts before the algorithm stops even if less than N transitions occurred.
• tmax=1e3: maximum time when the simulation stops even x_f has not been reached
• radius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.
• cut_start=true: if false, returns the whole trajectory up to the transition
• show_progress=true: shows a progress bar with respect to Nmax
• seed=nothing: master seed for per-simulation reseeding. If nothing, a fresh random master seed is drawn each call (so independent calls give independent ensembles even when the same sys is reused). Pass an integer for reproducibility.
• kwargs...: keyword arguments passed to CommonSolve.solve. These take precedence over the options stored in sys.diffeq, which are forwarded by default.

See also transition.

Returns a TransitionEnsemble object.

struct TransitionEnsemble{SSS, T, ES}

Ensemble of transition paths between two points in a state space.

Fields

• paths::Vector: paths sampled from the transition process

• times::Array{Vector{T}, 1} where T: coresponsing times of the paths

• stats::CriticalTransitions.TransitionStatistics: statistics of the ensemble

• sciml_ensemble::Any: original ensemble solution of the SciML Ensemble problem

Constructors

TransitionEnsemble(sim, success_rate)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/TransitionEnsemble.jl:55.

struct TransitionStatistics{T}

Statistics of the ensemble of transition paths between two points in a state space.

Fields

• success_rate::Any: success rate of the transition process

• residence_time::Any: mean residence time of the transition process

• transition_time::Any: mean transition time of the transition process

• rareness::Any: rareness of the transition process

Constructors

TransitionStatistics(sim, success_rate)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/TransitionEnsemble.jl:23.

fw_action(sys::CoupledSDEs, path, time) -> Any

Freidlin-Wentzell action of path (a D × N matrix with D = dimension(sys)) with time points time (length N) for the drift b = dynamic_rule(sys):

\[S_T[\phi] \;=\; \tfrac{1}{2}\int_0^T \big\| \dot\phi - \mathbf{b}(\phi)\big\|_{\mathbf{Q}(\phi)}^2 \,\mathrm{d}t\]

where $\|a\|_\mathbf{Q}^2 = \langle a, \mathbf{Q}^{-1} a\rangle$ (see anorm), $T$ is the total time of the path, and $\mathbf{Q}(x)$ is the trace-normalized diffusion tensor a(x) / s with s = tr(a(u₀))/D pinned at current_state(ds); for state-dependent diffusions s depends on the chosen u₀.

geometric_action(sys::CoupledSDEs, path) -> Any
geometric_action(sys::CoupledSDEs, path, arclength) -> Any

Geometric Freidlin-Wentzell action of path (a D × N matrix) with the given arclength, for the drift b = dynamic_rule(sys). Equals $\inf_T S_T[\varphi]$ of the time-parameterized fw_action on the instanton:

\[\bar S[\varphi] \;=\; \int_0^L \Big( \|\varphi'\|_\mathbf{Q}\,\|\mathbf{b}(\varphi)\|_\mathbf{Q} \;-\; \langle \varphi', \mathbf{b}(\varphi)\rangle_\mathbf{Q} \Big)\,\mathrm{d}s\]

where $s$ is the arclength coordinate, $L$ the arclength, and $\|a\|_\mathbf{Q}^2 = \langle a, b\rangle_\mathbf{Q} = a^\top \mathbf{Q}^{-1} b$ (see anorm). $\mathbf{Q}$ is the trace-normalized covariance_matrix(sys). As with fw_action, the returned value is independent of the noise_strength keyword and rotation-invariant.

geometric_action(b::Function, path, arclength=1.0; A=nothing)

Geometric Freidlin-Wentzell action for a drift field b and a discrete path.

This is a drift-only convenience overload that uses the same integrand as geometric_action(sys::CoupledSDEs, ...), but requires an explicit inverse covariance (A = Q^{-1}) if you want anything other than the identity metric.

If A === nothing, it defaults to the identity metric.

The path must be a (D × N) matrix with points in columns.

om_action(sys::CoupledSDEs, path, time, noise_strength)

Onsager-Machlup action of path (a D × N matrix) with time points time (length N) for the drift b = dynamic_rule(sys), at the given noise_strength σ:

\[S^\sigma_T[\phi] \;=\; \tfrac{1}{2}\int_0^T \Big(\big\|\dot\phi - \mathbf{b}(\phi)\big\|_\mathbf{Q}^2 \;+\; \sigma^2 \,\nabla\cdot \mathbf{b}(\phi)\Big)\,\mathrm{d}t\]

where $\|a\|_\mathbf{Q}^2 = \langle a, \mathbf{Q}^{-1} a\rangle$ (see anorm), $T$ is the total time, and $\mathbf{Q}$ is the trace-normalized covariance_matrix(sys). The first term is exactly fw_action and is independent of the noise_strength keyword under the trace-normalize convention. The second term is the Onsager-Machlup finite-noise correction parameterized by the explicit noise_strength argument: pass the σ at which you want $-\log P[\phi]$ evaluated (typically the same noise_strength you used at construction). As $\sigma \to 0$, om_action → fw_action.

action(
    sys::CoupledSDEs,
    path::AbstractMatrix,
    time,
    functional::AbstractString;
    noise_strength
) -> Any

Computes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.

• functional = "FW": Returns the Freidlin-Wentzell action (fw_action)
• functional = "OM": Returns the Onsager-Machlup action (om_action)

minimize_action(sys::ContinuousTimeDynamicalSystem, x_i, x_f, T::Real, optimizer=Optimisers.Adam(); kwargs...)

Minimizes an action functional to obtain a minimum action path (instanton) between an initial state x_i and final state x_f in phase space.

This algorithm uses the Optimization.jl package to minimize the specified action functional (either fw_action or om_action) for the system sys over paths connecting x_i to x_f in time T.

The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is $\Delta t = T/N$. To set an initial path different from a straight line, see the multiple dispatch method

minimize_action(sys::ContinuousTimeDynamicalSystem, init::AbstractMatrix, T::Real, optimizer=Optimisers.Adam(); kwargs...).

Returns a MinimumActionPath object containing the optimized path and the action value.

The optional positional argument optimizer selects the Optimization.jl solver. It defaults to Optimisers.Adam().

Keyword arguments

• functional = "FW": type of action functional to minimize. Defaults to fw_action, alternative: "OM" for om_action
• npoints = 100: number of path points to use for the discretization of the path
• noise_strength = nothing: noise strength for the action functional. Specify only if functional = "OM"
• ad_type = OptimizationBase.AutoFiniteDiff(): type of automatic differentiation to use (see Optimization.jl)
• maxiters = 100: maximum number of iterations before the algorithm stops
• abstol::Real=NaN: absolute tolerance of action gradient to determine convergence
• reltol::Real=NaN: relative tolerance of action gradient to determine convergence
• verbose = false: whether to print Optimization information during the run
• show_progress = false: whether to print a progress bar

minimize_action(sys::ContinuousTimeDynamicalSystem, init::AbstractMatrix, T::Real, optimizer=Optimisers.Adam(); kwargs...)

Minimizes the specified action functional to obtain a minimum action path (instanton) between fixed end points given a system sys and total path time T.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is $\Delta t = T/N$.

The optional positional argument optimizer selects the Optimization.jl solver. It defaults to Optimisers.Adam().

minimize_geometric_action(sys::FreidlinWentzellHamiltonian, x_initial, optimizer = GeometricGradient(; stepsize = 1e3); kwargs...)

Runs the simplified geometric Minimum Action Method (sgMAM, (Grafke et al., 2017)) on the FreidlinWentzellHamiltonian sys, returning the instanton (minimizer of the Freidlin-Wentzell action) that connects the two endpoints of x_initial.

sgMAM works directly in the doubled $(x, p)$ phase space and minimizes the symplectic action $\int p \cdot \mathrm{d}x$ under the zero-energy constraint $H(x, p) = 0$. Here $p = a(x)^{-1}(\dot x - b(x))$ is the conjugate momentum, i.e. the (dimensionless) noise force that the Wiener process must supply to drive the trajectory over the barrier.

Arguments

• sys::FreidlinWentzellHamiltonian: the Hamiltonian, typically obtained from a CoupledSDEs via FreidlinWentzellHamiltonian(ds).
• x_initial::AbstractMatrix{T} (or StateSpaceSet): initial guess for the instanton, shape D × Nt with state points in columns. Endpoints x_initial[:, 1] and x_initial[:, end] are held fixed; the interior is reparameterized to uniform arclength on the first iteration.
• optimizer: step-size control, either
  • GeometricGradient (default; backtracking projected gradient), or
  • AdaptiveGeometricGradient (multi-phase probe variant; more robust on underdamped systems at higher per-iteration cost).

Keyword arguments

• maxiters::Int = 1000: maximum outer iterations.
• abstol::Real = NaN, reltol::Real = NaN: convergence on the absolute / relative action change between accepted iterations. NaN disables that criterion.
• verbose::Bool = false: log convergence and backtracking events.
• show_progress::Bool = false: show a ProgressMeter bar.

Returns a MinimumActionPath

minimize_geometric_action(sys::CoupledSDEs, x_i, x_f, optimizer=GeometricGradient(); kwargs...)
minimize_geometric_action(sys::CoupledSDEs, init::AbstractMatrix, optimizer=GeometricGradient(); kwargs...)

Computes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of Optimization.jl or the original formulation by Heymann and Vanden-Eijnden (2008) (projected gradient descent with optional backtracking).

The minimizer is computed for system sys over all paths from x_i to x_f. To set an initial path different from a straight line, see the multiple-dispatch method minimize_geometric_action(sys::CoupledSDEs, init::AbstractMatrix, optimizer; kwargs...).

Returns a MinimumActionPath.

Keyword arguments

• npoints::Int = 100: number of discretization points (endpoint form only)
• maxiters::Int = 100: maximum optimizer iterations
• abstol::Real = NaN, reltol::Real = NaN: action-change convergence criteria
• ad_type = OptimizationBase.AutoFiniteDiff() (Optimization.jl path only)
• verbose::Bool = false, show_progress::Bool = true

string_method(
    sys::Union{Function, FreidlinWentzellHamiltonian},
    x_initial::AbstractMatrix;
    stepsize,
    integrator,
    maxiters,
    show_progress
) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}

Compute the string method for a given system using the method of E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

References

• (E et al., 2007)
• StringMethod.jl

Keyword arguments

• stepsize::Real=1e-1: Step size for the evolution step. Default: 0.1
• integrator=Euler(): SciML integrator algorithm (e.g. Euler(), Tsit5())
• maxiters::Int=1000: Maximum number of iterations for path convergence
• show_progress::Bool=false: Whether to display a progress meter during computation

Returns

A MinimumActionPath containing the converged path and action value.

For sys::CoupledSDEs the .action field is computed as the geometric Freidlin-Wentzell action via geometric_action. For drift-only systems (e.g. sys::Function) the same geometric action is computed assuming an identity noise covariance.

string_method(
    sys::ContinuousTimeDynamicalSystem,
    init;
    kwargs...
) -> MinimumActionPath{_A, _B, _B1, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _B1}

Compute the string method for a given system using the method of E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

References

• (E et al., 2007)

Keyword arguments

• stepsize::Real=1e-1: Step size for the evolution step. Default: 0.1
• integrator=Euler(): SciML integrator algorithm (e.g. Euler(), Tsit5())
• maxiters::Int=1000: Maximum number of iterations for path convergence
• show_progress::Bool=false: Whether to display a progress meter during computation

Returns

A MinimumActionPath containing the converged path and action value.

For sys::CoupledSDEs, the action is computed using the system's noise covariance.

struct MinimumActionPath{D, T<:Real, V, Phis, Ahis, Lambda, GPV, PV}

The minimum action path between two points in a D-dimensional phase space.

Fields

• path::StateSpaceSet{D, T, V} where {D, T<:Real, V}: The path matrix.

• action::Real: The action value associated to the path.

• path_history::Any: The history of action of the paths in the optimisation algorithm (optional).

• action_history::Any: The history of action of the paths in the optimisation algorithm (optional).

• λ::Any: The Lagrange multiplier parameter for the minimum action path (optional).

• generalized_momentum::Any: The generalized momentum of the phase space variables (optional).

• path_velocity::Any: The path velocity (optional).

Constructors

MinimumActionPath(
    path,
    action;
    path_history,
    action_history,
    λ,
    generalized_momentum,
    path_velocity
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.

struct GeometricGradient{T<:Real} <: CriticalTransitions.GMAMOptimizer

Optimizer configuration for the (s)gMAM projected-gradient update with built-in backtracking step-size control.

By default, backtracking is enabled (max_backtracks=10): each iteration tries the current step size and, if the action increases, shrinks it by shrink and retries up to max_backtracks times. On accepted steps the step size grows by grow for the next iteration. This makes the algorithm insensitive to the initial step size choice: small starting values are grown automatically, and large values are safely reduced by backtracking. In practice, prefer a large initial stepsize — a rejected step costs only a few extra action evaluations, whereas a step size that is too small slows every iteration throughout the entire run. Set max_backtracks=0 to disable backtracking and use a fixed step size.

Fields

• stepsize::Real: Initial step size for the projected gradient update.

• shrink::Real: Step-size shrink factor on rejected steps (backtracking).

• grow::Real: Step-size growth factor after accepted steps (backtracking).

• max_backtracks::Int64: Maximum number of backtracking attempts per iteration.

• stepsize_min::Real: Lower clamp for step size.

• stepsize_max::Real: Upper clamp for step size.

Keyword constructors

GeometricGradient(
    stepsize,
    shrink,
    grow,
    max_backtracks,
    stepsize_min,
    stepsize_max
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/methods.jl:89.

GeometricGradient(
;
    stepsize,
    shrink,
    grow,
    max_backtracks,
    stepsize_min,
    stepsize_max
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/methods.jl:103.

struct AdaptiveGeometricGradient{T<:Real} <: CriticalTransitions.GMAMOptimizer

Optimizer configuration for the adaptive multi-phase step-size selection variant of the (s)gMAM projected-gradient update.

Unlike GeometricGradient, which adapts the step size from one iteration to the next, this optimizer compares two probes of probe_length iterations each, both starting from the same path: one at the current step size ϵ, one at ϵ * shrink. Whichever probe yields the lower action after the full window is kept, and the step size is updated accordingly (shrunk if the small probe wins, grown otherwise).

This is significantly more robust than single-iteration backtracking on underdamped systems where a large step size makes immediate progress but over-smooths the path, and a small step size accumulates to a better final result. By comparing over many iterations rather than a single update, the cumulative effect of over-smoothing becomes detectable.

The cost is roughly 2× per probe window relative to a fixed-step run.

Fields

• stepsize::Real: Initial step size for the projected gradient update.

• probe_length::Int64: Number of inner iterations per probe window.

• shrink::Real: Step-size shrink factor (also defines the small-probe step ϵ * shrink).

• grow::Real: Step-size growth factor applied when the larger probe wins.

• stepsize_min::Real: Lower clamp for step size.

• stepsize_max::Real: Upper clamp for step size.

Keyword constructors

AdaptiveGeometricGradient(
    stepsize,
    probe_length,
    shrink,
    grow,
    stepsize_min,
    stepsize_max
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/methods.jl:29.

AdaptiveGeometricGradient(
;
    stepsize,
    probe_length,
    shrink,
    grow,
    stepsize_min,
    stepsize_max
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/methods.jl:43.

struct MultipleShooting{S, J, T} <: CriticalTransitions.GMAMOptimizer

Multiple-shooting BVP optimizer for the Freidlin-Wentzell instanton on a FreidlinWentzellHamiltonian. Integrates the arclength-reparameterized Hamilton equations

\[\frac{\mathrm{d}\varphi}{\mathrm{d}s} = \alpha\,H_p,\qquad \frac{\mathrm{d}p}{\mathrm{d}s} = -\alpha\,H_x,\qquad \alpha = L / \|H_p\|\]

on s ∈ [0, 1] with path length L a Newton unknown. Boundary states are parameterized by the unstable / stable eigenvectors of the Hamiltonian Jacobian at each fixed-point endpoint. The heteroclinic must not cross a drift fixed point in its interior; through-saddle problems must be user-split into attractor → saddle legs and the actions summed.

Fields

• nshoots::Int64: Number of shooting segments (≥ 2).

• ode_solver::Any: ODE integrator passed to the per-segment solve.

• nlsolve::Any: Override for the Newton solver; nothing builds a default sparse-AD NewtonRaphson.

• abstol::Any: Newton absolute tolerance.

• reltol::Any: Newton relative tolerance.

• eps_lin::Any: Anchor magnitude: distance of y_0 from (xa, 0) in 2D-space.

• invariant_tol::Any: Warning threshold on max_s |H(φ(s), p(s))| along the converged path.

• maxiters::Int64: Maximum Newton iterations.

quasipotential(sys, grid, attractor; band_radius, near_source_layers,
               verbose, show_progress) -> QuasiPotential{D}

Compute the Freidlin-Wentzell quasipotential field U_A(x) from attractor using the Ordered Line Integral Method (Dahiya and Cameron 2018). The state dimension D is taken from sys::CoupledSDEs{IIP, D, I, P} and must match grid::CartesianGrid{D}. A warning is emitted for D > 4.

Keyword arguments

• band_radius::Int = default_K(grid): accepted-band radius in grid cells.
• near_source_layers::Int = 3: size of the analytic CARE seed box; 0 disables analytic seeding.
• regularization::Real = default_regularization(grid): for a system with a single noiseless coordinate (rank-1 diffusion, e.g. momentum-only noise in a second-order Langevin or van der Pol oscillator) the trace-normalized diffusion is singular and has no metric. This amount is added to that coordinate's diagonal (a structural vanishing-viscosity perturbation that leaves the noisy block exact) to make it invertible; the degenerate Freidlin-Wentzell quasipotential is recovered as it tends to 0. The default scales like 1 / minimum(grid.nbox) (about 0.04 at 80 cells), so it shrinks under refinement; smaller values are more accurate but stiffer. The escape sheet and saddle barrier carry negligible bias, the non-escape region carries the regularization bias. Ignored for full-rank systems; two or more noiseless coordinates (sub-Riemannian) or a non-coordinate-aligned null space are rejected.
• verbose::Bool = false
• show_progress::Bool = true

See also: QuasiPotential, BackRef.

struct QuasiPotential{D, T<:AbstractFloat}

Discrete quasipotential field U_A(x) computed by quasipotential. U is +Inf on cells the sweep did not reach. back_pointer encodes the winning sub-cell predecessor per cell (zeroed for the source and unreached cells).

Fields

• U::Array{T, D} where {D, T<:AbstractFloat}

• back_pointer::Array{CriticalTransitions.BackRef{D}, D} where D

• source::CartesianIndex

• grid::CartesianGrid

struct BackRef{D}

Sub-cell back-reference used to reconstruct instantons from a QuasiPotential grid. For one-point (vertex) updates v1 == v0 and s == NaN32. For two-point (edge) updates the predecessor lies at (1 - s) * v0 + s * v1 with s ∈ [0, 1].

Fields

• v0::CartesianIndex

• v1::CartesianIndex

• s::Float32

struct DiffusionGenerator{D, BC<:Tuple, T<:AbstractFloat}

The discretised generator of an autonomous Itô diffusion on a CartesianGrid, built by the Scharfetter–Gummel exponential- fitting finite-volume scheme.

The matrix Q generates the transition semigroup exp(t Q) of the discretised process. With probability densities and observables represented as column vectors (the convention used throughout this module), it propagates observables forward (du/dt = Q u, the discrete backward Kolmogorov) and densities forward via the transpose (dρ/dt = Qᵀ ρ, the discrete Fokker–Planck, exposed as fokker_planck_operator). Equivalently, Q[i, j] for i ≠ j is the transition rate from cell i to cell j.

The same matrix has several traditional names:

Use m_matrix for the M-matrix view (PDE convention), rate_matrix for the CTMC view, or fokker_planck_operator for the forward Kolmogorov / FP operator (Qᵀ).

Fields

• grid::CartesianGrid: The Cartesian grid the generator lives on.

• Q::SparseArrays.SparseMatrixCSC{T, Int64} where T<:AbstractFloat: Sparse rate matrix (off-diagonals = transition rates ≥ 0).

• bc::Tuple: Per-axis boundary condition (one BoundaryCondition per axis).

Construction

DiffusionGenerator(sys::CoupledSDEs, grid::CartesianGrid; bc=Reflecting())

Discretises the SDE's backward Kolmogorov operator on grid. Drift b(x) is read via drift; diffusion is read once from covariance_matrix(sys) and must be diagonal (axis-aligned noise).

bc controls how the outer grid faces are treated and is a singleton BoundaryCondition instance applied to every axis, or a D-tuple for per-axis control:

• Reflecting() (default) — outer faces are omitted; chain bounces off the edges.
• Periodic() — outer faces wrap to the opposite end of the same axis.
• Absorbing() — boundary cells leak probability through their outer faces. Q is a sub-generator (row sums on the boundary become negative).

rate_matrix(
    generator::DiffusionGenerator
) -> SparseArrays.SparseMatrixCSC{T, Int64} where T<:AbstractFloat

Rate matrix of gen: returns gen.Q, also called the discrete backward Kolmogorov / generator.

The off-diagonals are transition rates between adjacent cells, the diagonal is minus the escape rate. For the PDE M-matrix sign convention see m_matrix; for the forward Kolmogorov / Fokker–Planck operator (acts on densities), see fokker_planck_operator.

m_matrix(
    generator::DiffusionGenerator
) -> SparseArrays.SparseMatrixCSC{Tv, Int64} where Tv

M-matrix of the generator. Equivalent to -rate_matrix(gen).

fokker_planck_operator(
    generator::DiffusionGenerator
) -> SparseArrays.SparseMatrixCSC{Tv, Int64} where Tv

Discrete Fokker–Planck operator (forward Kolmogorov operator) of generator.

It acts on probability densities: the discrete Fokker–Planck equation is

dρ/dt = fokker_planck_operator(gen) * ρ ,

and the invariant density solves fokker_planck_operator(gen) * ρ = 0.

struct CartesianGrid{D, T<:AbstractFloat}

A D-dimensional uniform axis-aligned Cartesian grid with coordinates stored in the floating-point type T (default Float64).

Fields

• nbox::Tuple{Vararg{Int64, D}} where D: Number of cells along each axis.

• h::Tuple{Vararg{T, D}} where {D, T<:AbstractFloat}: Cell width along each axis.

• centers::Tuple{Vararg{LinRange{T, Int64}, D}} where {D, T<:AbstractFloat}: Cell-center coordinates along each axis (length N_k).

• edges::Tuple{Vararg{LinRange{T, Int64}, D}} where {D, T<:AbstractFloat}: Cell-edge coordinates along each axis (length N_k + 1).

Construction

CartesianGrid((u_lo, u_hi, N_u), (v_lo, v_hi, N_v), ...)            # Float64
CartesianGrid{Double64}((u_lo, u_hi, N_u), (v_lo, v_hi, N_v), ...)  # extended

Each axis is (lo, hi, N) and creates N equal-width cells. The storage type T propagates into DiffusionGenerator and the resulting rate matrix.

Abstract supertype of boundary conditions accepted by DiffusionGenerator. Concrete subtypes: Reflecting, Periodic, Absorbing.

No-flux Neumann boundary: outer faces are omitted from the stencil. The chain bounces off the grid edges, mass is preserved.

Periodic boundary: outer faces wrap to the opposite end of the same axis. Use for systems on a torus / ring (angle variables).

Absorbing boundary: cells on the outer boundary leak probability through their outer faces (escape rate from the cell-center drift). Row sums on boundary cells become negative; the resulting Q is a sub-generator.

stationary_distribution(
    generator::DiffusionGenerator{D, BC, T};
    ...
) -> Any
stationary_distribution(
    generator::DiffusionGenerator{D, BC, T},
    alg;
    verbose,
    multimodal_tol,
    kwargs...
) -> Any

Invariant probability density of the discretised process generated by generator: solves the discrete adjoint stationary equation Qᵀ ρ = 0 augmented with the normalisation ∫ ρ dV = 1, where Q = generator.Q is the rate matrix.

Returns a length-ncells(generator) vector.

Algorithm (alg)

• A LinearSolve.SciMLLinearSolveAlgorithm — pinned linear solve with the given algorithm. Default is UMFPACKFactorization() (sparse direct LU). Pass e.g. KrylovJL_GMRES() for an iterative solver. Extra kwargs flow to LinearSolve.solve.
• DenseEigen() — LinearAlgebra.eigen on Matrix(Qᵀ), picking the eigenvector at λ ≈ 0. Bounded by O(N²) storage.
• KrylovKitSolver() — shift-invert Arnoldi via KrylovKit.eigsolve near σ = 0. Use for large sparse problems; pass inner_alg = … to control the LinearSolve algorithm used for the inner shift-invert solve. Extra kwargs flow to KrylovKit.eigsolve.

Diagnostic

The result is always validated post-hoc:

1. Sign check: entries below -16ε · max|ρ| are flagged.
2. Residual check: ‖Qᵀ ρ‖ / ‖ρ‖ should be near machine ε.

Pass verbose = true to additionally run the multi-modality probe:

• Linear-solve path: re-pins at several spread-out cells and checks the answer is invariant. Adds three full LU solves.
• Eigensolver path (DenseEigen, KrylovKitSolver): inspects the spectral gap above the kernel by computing λ₂, λ₃ of Qᵀ and flagging |λ₂| / |λ₃| < multimodal_tol. Adds one extra eigenmodes call (a single LU factorisation reused across Arnoldi iterations).

In both cases a positive result reveals ≥ 2 near-zero eigenvalues within machine precision (multi-basin metastability): the returned vector is one quasi-stationary mode rather than the global invariant, and you should switch to quasi_stationary_distribution per basin. The probe is off by default.

A single combined warning is emitted with concrete suggestions if any check fires.

quasi_stationary_distribution(
    generator::DiffusionGenerator,
    basin;
    ...
) -> Tuple{Vector{T} where T<:AbstractFloat, Any}
quasi_stationary_distribution(
    generator::DiffusionGenerator,
    basin,
    alg::Union{DenseEigen, KrylovKitSolver};
    kwargs...
) -> Tuple{Vector{T} where T<:AbstractFloat, Any}

Quasi-stationary distribution (QSD) of a metastable basin.

Returns (ρ_QSD, λ_exit) where ρ_QSD is the length-ncells(gen) QSD vector (with the generator's float element type, zero on cells outside basin, normalised to ∫ ρ_QSD dV = 1 on basin) and λ_exit > 0 is the corresponding exit rate. The mean lifetime in the basin, conditional on staying in it long enough to reach quasi-equilibrium, is 1 / λ_exit.

basin accepts a predicate x -> Bool, a BitVector of length ncells(gen), or a vector of linear cell indices.

Theory

The QSD is the principal left eigenvector of the basin sub-generator Q_S (the rows and columns of Q indexed by the basin S):

\[Q_S^\top \rho_S = -\lambda_{\text{exit}}\, \rho_S ,\qquad \rho_S \ge 0 ,\qquad \int \rho_S\, dV = 1 ,\]

with -λ_exit the eigenvalue of Q_S^T closest to zero (it is real and strictly negative because mass leaks across the basin boundary). Equivalently, λ_exit is the smallest-magnitude eigenvalue of -Q_S. For a single metastable set inside an otherwise-reflecting domain, the QSD is the long-time limit of the distribution conditioned on not yet having left S, and -D log ρ_QSD is the local quasipotential of S at noise strength D << 1.

This routine is the standard tool for visualising the quasipotential well of a metastable basin at small noise, where the global stationary_distribution becomes ill-conditioned.

Algorithm (alg)

• KrylovKitSolver() (default) — shift-invert Arnoldi via KrylovKit.eigsolve. Kwargs flow to KrylovKit.eigsolve; pass inner_alg = … to control the LinearSolve algorithm used for the inner shift-invert solve.
• DenseEigen() — LinearAlgebra.eigen on Matrix(Q_Sᵀ). Use on small basins where the dense factorisation fits in memory.

Passing a LinearSolve.SciMLLinearSolveAlgorithm is rejected: the QSD eigenvalue −λ_exit is unknown a priori, so a single linear solve is insufficient.

mean_first_passage_time(
    generator::DiffusionGenerator{D, BC, T},
    target;
    alg,
    kwargs...
) -> Any

Mean first-passage time τ[i] from cell i to the target set: the expected time for a trajectory of the process generated by generator to first hit target starting at i. Solves Q τ = -1 on the free cells with Dirichlet condition τ = 0 on target.

target accepts a predicate, a BitVector, or a vector of linear cell indices.

Default alg is UMFPACKFactorization() (sparse direct LU). Pass any SciMLLinearSolveAlgorithm for an iterative solver.

See also first_passage_variance.

first_passage_variance(
    generator::DiffusionGenerator{D, BC, T},
    target;
    alg,
    kwargs...
) -> Any

Variance Var[τ | X_0 = i] = E[τ²] − E[τ]² of the first-passage time to target, as a function of starting cell i.

The first moment τ_1 = E[τ] solves Q τ_1 = -1 with τ_1|_target = 0 (mean_first_passage_time); the second moment τ_2 = E[τ²] solves Q τ_2 = -2 τ_1 with τ_2|_target = 0. The variance follows from Var[τ] = τ_2 − τ_1².

Returns a length-ncells(generator) vector of variances, all ≥ 0.

Default alg is UMFPACKFactorization() (sparse direct LU). Pass any SciMLLinearSolveAlgorithm for an iterative solver.

See also mean_first_passage_time.

eigenmodes(gen::DiffusionGenerator; ...) -> Tuple{Any, Any}
eigenmodes(
    gen::DiffusionGenerator,
    k::Integer;
    ...
) -> Tuple{Any, Any}
eigenmodes(
    gen::DiffusionGenerator,
    k::Integer,
    alg::Union{DenseEigen, KrylovKitSolver};
    kwargs...
) -> Tuple{Any, Any}

Slowest-decaying eigenmodes of the generator gen. Returns (λ, V) where λ is a vector of Complex{T} eigenvalues (with T = floattype(gen)) with largest real part (least negative — slowest decay), sorted descending, and V collects the corresponding right eigenvectors of Q columnwise (Q * V[:, i] = λ[i] * V[:, i]).

For mass-preserving boundary conditions (Reflecting / Periodic) the eigenvalue closest to zero is exactly 0, with eigenvector the constant function (since Q * 1 = 0); the corresponding left eigenvector is the invariant density. For Absorbing boundaries Q is a sub-generator and all eigenvalues lie strictly in the open left half-plane (no zero mode). The remaining eigenvalues lie in the open left half-plane, and -1 / Re(λ_k) gives the metastable timescale of the k-th mode. The slow non-trivial right eigenvectors are the canonical reaction coordinates / metastable-state indicators used in Markov state model decomposition (e.g. PCCA+).

Algorithm (alg)

• KrylovKitSolver() (default) — shift-invert Arnoldi via KrylovKit.eigsolve near σ = 0. The right choice for large sparse generators. Kwargs flow to KrylovKit.eigsolve; pass inner_alg (a SciMLLinearSolveAlgorithm, defaults to UMFPACKFactorization()) to override the LinearSolve algorithm for the inner shift-invert solve.
• DenseEigen() — LinearAlgebra.eigen on Matrix(gen.Q), returning all eigenpairs sliced to the first k. Fine up to ~5000 cells.

Passing a LinearSolve.SciMLLinearSolveAlgorithm is rejected: the target eigenvalues are unknown a priori, so a single linear solve is insufficient.

propagate_density(
    gen::DiffusionGenerator{D, BC, S},
    T::Real,
    ρ_0::AbstractVector;
    Δt,
    Ttr,
    tol,
    m,
    adaptive
) -> Tuple{Any, Any}

Time-evolved probability density. Solves the discrete Fokker–Planck equation dρ/dt = Qᵀ ρ (where Qᵀ =fokker_planck_operator(gen)) starting from ρ_0, and records ρ on a uniform time grid.

Returns (ρs, t) where t = Ttr:Δt:(Ttr + T) and ρs is a Matrix (with the generator's float element type) of size (ncells(gen), length(t)) whose i-th column is ρ(t[i]).

ρs, t = propagate_density(gen, T, ρ_0)                      # snapshot at 0 and T
ρs, t = propagate_density(gen, T, ρ_0; Δt = 0.1)            # snapshots every 0.1
ρs, t = propagate_density(gen, T, ρ_0; Δt = 0.1, Ttr = 5)   # skip transient

Arguments

• gen — diffusion generator (provides Qᵀ).
• T — total time over which the trajectory is recorded.
• ρ_0 — initial density (length-ncells(gen) vector).

Keyword arguments

• Δt = T — sampling interval. Default returns just the initial and final snapshots.
• Ttr = 0 — transient time skipped before recording starts.
• tol = 1e-7 — Krylov local-error tolerance.
• m = 30 — maximum Krylov subspace dimension before restart.
• adaptive = true — use adaptive Krylov substepping. Disable only for debugging or to force fixed-step Krylov of dimension m.

For mass-preserving boundary conditions (Reflecting / Periodic) the total mass is preserved across the trajectory; for Absorbing boundaries it decays as probability leaks through. Internal substepping is automatic — tol/m/adaptive control the Krylov-subspace matrix-exponential algorithm in ExponentialUtilities, which handles stiffness without ever forming exp(t·Qᵀ) explicitly.

struct DenseEigen

Dense eigensolver backend: LinearAlgebra.eigen on Matrix(Q).

Returns every eigenpair; fine up to a few thousand cells, becomes infeasible above that because of the dense O(N²) storage. Always applicable.

struct KrylovKitSolver

KrylovKit.eigsolve with shift-invert through LinearSolve.jl.

The workhorse for large sparse problems and for any analysis where the target eigenvalue is not known a priori. Tuning kwargs (tol, krylovdim, maxiter, …) flow to KrylovKit.eigsolve; pass inner_alg to override the LinearSolve algorithm used for the inner shift-invert solve.

ncells(grid::CartesianGrid) -> Int64

Total number of cells in grid.

cell_volume(grid::CartesianGrid) -> Any

Volume of one cell (product of axis spacings).

cell_center(
    grid::CartesianGrid{D, T},
    I::CartesianIndex{D}
) -> Any

Center coordinates of cell I of grid.

ball(
    center,
    radius::Real
) -> CriticalTransitions.var"#187#188"{_A, <:Real} where _A

Return a predicate for the open ball with given center and radius.

cuboid(lo, hi) -> CriticalTransitions.var"#189#191"

Predicate matching points inside an axis-aligned box: returns the closure x -> all(lo .≤ x .≤ hi). lo and hi may be tuples, SVectors, or any indexable types of the right length.

sublevel(
    f,
    c::Real
) -> CriticalTransitions.var"#193#194"{_A, <:Real} where _A

Predicate for the sublevel set {x : f(x) < c}. Useful for masking by a potential, energy, or any scalar-valued function on state space.

using CriticalTransitions: sublevel
U(x) = 0.25 * (x[1]^2 - 1)^2 + 0.5 * x[2]^2
basin_low_energy = sublevel(U, 0.1)

reshape_to_grid(
    v::AbstractVector,
    gen::DiffusionGenerator
) -> Any

Reshape a per-cell vector (length ncells(grid)) back to a grid-shaped Array of size grid.nbox. Useful for plotting and for spatial operations on grid quantities.

using CriticalTransitions: reshape_to_grid
ρ = stationary_distribution(gen)
ρ_grid = reshape_to_grid(ρ, gen)
heatmap(grid.centers[1], grid.centers[2], ρ_grid)

Both DiffusionGenerator and CartesianGrid are accepted as the second argument.

struct ReactiveTransition{D, BC, T<:AbstractFloat}

Transition-path-theory result for transitions from set A to set B.

ReactiveTransition stores the generator, set masks, invariant density, forward and backward committors, discrete adjoint generator, and reactive rate. It represents the stationary ensemble of trajectories that most recently left A and will hit B before returning to A.

Fields

• generator::DiffusionGenerator{D, BC} where {D, BC}: The diffusion generator the analysis was built on.

• A_mask::BitVector: Boolean mask selecting cells in set A (length ncells(grid)).

• B_mask::BitVector: Boolean mask selecting cells in set B (length ncells(grid)).

• ρ::Vector{T} where T<:AbstractFloat: Invariant probability density (dot(ρ, weights) = 1).

• qplus::Vector{T} where T<:AbstractFloat: Forward committor q⁺ (probability of reaching B before A).

• qminus::Vector{T} where T<:AbstractFloat: Backward committor q⁻.

• Qadj::SparseArrays.SparseMatrixCSC{T, Int64} where T<:AbstractFloat: Discrete adjoint generator w.r.t. ρ; used by reactive_current_{reversible,irreversible}.

• rate::AbstractFloat: A → B reactive transition rate (cached).

• physical_reverse::Bool: True if qminus came from a user-supplied physical reverse generator; false if from the discrete adjoint.

Construction

ReactiveTransition(gen::DiffusionGenerator, A, B; alg=UMFPACKFactorization(), reverse=nothing)
ReactiveTransition(sys::CoupledSDEs, grid, A, B; bc=Reflecting(), reverse=nothing)

The first form takes a pre-built DiffusionGenerator — useful for sweeping multiple (A, B) pairs without rediscretising. The second is a convenience that builds the generator from sys internally.

A and B accept:

• a predicate x::SVector{D} -> Bool evaluated at each cell center,
• a BitVector of length ncells(grid), or
• a vector of linear cell indices.

If reverse is provided (a CoupledSDEs for the high-level form, a DiffusionGenerator for the generator form), the backward committor is computed on that physical reverse-drift generator instead of the discrete adjoint.

res = ReactiveTransition(gen, x -> x[1] < -1, x -> x[1] > 1)
k = reactive_rate(res)
J_nodes, J_faces = reactive_current(res)

Default alg is UMFPACKFactorization() (sparse direct LU) for all internal solves. Pass any SciMLLinearSolveAlgorithm (e.g. KrylovJL_GMRES()) for an iterative solver — useful for large grids where direct LU runs out of memory. Further kwargs flow to LinearSolve.solve.

See also forward_committor, backward_committor, reactive_density, and reactive_current.

forward_committor(
    generator::DiffusionGenerator{D, BC, T},
    A,
    B;
    alg,
    kwargs...
) -> Any

Forward committor q⁺[i]: probability that a trajectory of the process generated by generator started at cell i reaches B before A. Solves the homogeneous linear system Q q⁺ = 0 with Dirichlet conditions q⁺ = 0 on A and q⁺ = 1 on B.

A and B accept a predicate x -> Bool, a BitVector of length ncells(generator), or a vector of linear cell indices.

Default alg is UMFPACKFactorization() (sparse direct LU). Pass any SciMLLinearSolveAlgorithm (e.g. KrylovJL_GMRES()) for an iterative solver; further kwargs flow to LinearSolve.solve.

See also backward_committor and ReactiveTransition.

backward_committor(
    generator::DiffusionGenerator{D, BC, T},
    A,
    B;
    alg,
    reverse,
    kwargs...
) -> Any

Backward committor q⁻[i]: probability that a trajectory in cell i last visited A (rather than B). For reversible processes q⁻ = 1 - q⁺.

By default q⁻ is solved on the discrete adjoint of generator.Q (time- reversal w.r.t. the invariant measure — this requires solving the stationary distribution as well). Pass reverse::DiffusionGenerator to instead solve on a separately-supplied physical reverse-drift generator; useful as a cross-check or for non-reversible systems where the exact reverse SDE is known.

Default alg is UMFPACKFactorization() (sparse direct LU). Pass any SciMLLinearSolveAlgorithm for an iterative solver; further kwargs flow to LinearSolve.solve.

See also forward_committor, stationary_distribution, and ReactiveTransition.

reactive_rate(r::ReactiveTransition) -> AbstractFloat

A → B reactive transition rate $k_{AB}$: events per unit time, the rate at which probability leaves A along reactive trajectories.

In CTMC form,

k_{AB} = ∑_{i ∈ A, j ∉ A} ρ[i] · v · q⁻[i] · Q[i, j] · q⁺[j] ,

with v = cell_volume(grid).

reactive_density(r::ReactiveTransition) -> Any

Per-cell reactive density ρ[i] * q⁺[i] * q⁻[i]. This is a density per unit volume. For the integrated reactive probability use probability_reactive.

reactive_current(r::ReactiveTransition) -> Tuple{Any, Any}

Reactive probability current as (J_nodes, J_faces):

• J_faces :: NTuple{D, Array{Float64, D}} — net signed flux across each axis-aligned face, divided by the corresponding axis spacing. Element k has shape nbox along all axes except axis k, where it has nbox[k] - 1 entries for reflecting / absorbing boundary conditions (one per interior face), or nbox[k] entries for periodic boundary conditions (the extra slot is the wraparound face from cell nbox[k] back to cell 1).
• J_nodes :: NTuple{D, Array{Float64, D}} — face fluxes averaged onto cell centers along each axis; element k has shape nbox.

The face flux from n to m is μ[n] · q⁻[n] · Q[n, m] · q⁺[m], where μ = ρ · v and Q is the rate matrix.

For non-equilibrium systems the current splits into a reversible (gradient-flow) part and an irreversible (cyclic) part — see reactive_current_reversible and reactive_current_irreversible to split gradient-flow and cyclic components.

reactive_current_reversible(
    r::ReactiveTransition
) -> Tuple{Any, Any}

Reversible (gradient-flow) part of the reactive current. Built from the symmetric part of the generator, Q_sym = (Q + Q̃)/2, where Q̃ is the discrete adjoint with respect to the invariant density. For a reversible system this equals the full reactive_current.

reactive_current_irreversible(
    r::ReactiveTransition
) -> Tuple{Any, Any}

Irreversible (cyclic / time-asymmetric) part of the reactive current. Built from the antisymmetric part of the generator, Q_antisym = (Q − Q̃)/2. Vanishes for systems satisfying detailed balance; for non-equilibrium systems it exposes the circulating component of the reactive flow.

Returns (J_nodes, J_faces) with the same shapes as reactive_current.

probability_reactive(r::ReactiveTransition{D, BC, T}) -> Any

Total reactive probability integrated over the grid. This is the stationary probability that the process is currently on a reactive segment from A to B.

See also reactive_density.

probability_last_A(r::ReactiveTransition{D, BC, T}) -> Any

Total probability that the most recent metastable visit was set A.

The complementary probability is the chance that the most recent visit was B.

RateSystem <: ContinuousTimeDynamicalSystem
RateSystem(ds::ContinuousTimeDynamicalSystem, forcing_profile::AbstractDict; kw...)

Construct a non-autonomous dynamical system by applying time-dependent parametric forcing protocols to an underlying autonomous continuous-time dynamical system ds. The unforced system parameters are copied and are always considered the starting parameter values irrespectively of what happens with ds afterwards.

forcing_profile is a Dict mapping parameter indices (anything valid for set_parameter!) to ForcingProfile instances. Each entry defines the functional form of how the corresponding parameter evolves over time. The profile is then combined with how long, and how large, each parameter forcing is, using the keywords below.

Keyword arguments

• forcing_start_time = initial_time(ds): start time(s) for the parameter ramping. You may supply an AbstractDict mapping keys to start times, or a scalar value which will be applied to all forcing profiles.
• forcing_duration = one(initial_time(ds)): duration of the parameter ramping (in system time units). Can be an AbstractDict mapping keys to durations or a scalar applied to all forcing profiles.
• forcing_scale = 1.0: amplitude multiplicative factor of the forcing profile. Can be an AbstractDict mapping keys to scales or a scalar applied to all forcing profiles.
• reverse = false: whether to reverse the forcing after the first forcing interval. Can be an AbstractDict mapping keys to booleans or a scalar boolean applied to all forcing profiles. If true, forcing is reversed over a second interval with the same duration.
• t0 = initial_time(ds): initial time of the RateSystem.

Description

The profile interval defines the domain of the forcing function; when applied to the system the profile is rescaled in system time units using the configured start and duration values - this allow changing the rate of the parameter ramping. Before a forcingprofile's start time the parameter equals its initial autonomous value; during the forcing interval it follows the rescaled profile (multiplied by the corresponding `forcingscalefactor). Ifreverse = true, the same profile is traversed backwards over a second interval of equal duration, after which the parameter returns to its initial autonomous value. Ifreverse = false`, after the interval the parameter is frozen at its final value.

To modify a rate system after it has been created, use set_forcing_duration!, set_forcing_scale!, set_forcing_start!, set_forcing_reverse!, and to obtain time dependent parameters use parameters, parameter.

Single parameter

You can use the convenience syntax

RateSystem(ds::ContinuousTimeDynamicalSystem, fp::ForcingProfile, pidx; kw...)

if you are only varying a single parameter with index pidx.

ForcingProfile(profile::Function, interval)

Time-dependent forcing profile $p(t)$ describing the evolution of a parameter over a domain interval = (start, end). Used to define a parametric forcing when constructing a non-autonomous RateSystem.

Arguments

• profile: function $p(t)$ from $ℝ → ℝ$ describing the parameter change
• interval: domain interval over which the profile is considered

Note: The units of $t$ are arbitrary; the forcing profile is rescaled to system time units.

Convenience functions

• ForcingProfile(:linear): Create a linear ramp from 0 to 1.
• ForcingProfile(:tanh): Create a hyperbolic tangent ramping from 0 to 1 with interval (-3, 3).

unforced_system(rs::RateSystem, t) -> Any

Returns an autonomous dynamical system corresponding to the frozen system of the non-autonomous RateSystem rs at time t.

parameters(rs::RateSystem, t)

Return the parameter container of a RateSystem at time t (in system time units). This function returns a non-aliased copy but also modifies the parameter state of rs to be on time t.

parameter(rs::RateSystem, t, pkey)

Returns the parameter with key pkey of a RateSystem at time t (in system time units).

set_forcing_start!(
    rs::RateSystem,
    start_time;
    pkey
) -> RateSystem

Sets the start time (forcing_start_time) of the forcing protocol applied to the RateSystem rs. If the optional keyword pkey is provided, only the start time for that parameter is changed; otherwise all parameters are updated.

set_forcing_duration!(
    rs::RateSystem,
    duration;
    pkey
) -> RateSystem

Sets the duration (forcing_duration) of the forcing protocol applied to the RateSystem rs. If the optional keyword pkey is provided, only the duration for that parameter is changed; otherwise all parameters are updated.

set_forcing_scale!(rs::RateSystem, scale; pkey=nothing)

Sets the amplitude (forcing_scale) of the forcing of the RateSystem to scale. If the keyword pkey is provided, only the corresponding parameter's forcing scale is changed; otherwise the forcing scale of all parameters is updated to scale.

set_forcing_reverse!(
    rs::RateSystem,
    value;
    pkey
) -> RateSystem

Set whether the rate forcing is also reversed for the RateSystem rs. If the optional keyword pkey is provided, only the reversing for that parameter is changed; otherwise all forcings are updated to the given value.

rate_track_return_tip(rs::RateSystem, Δts, Δps, mapper, ics; kw...)

Utilize the global_continuation functionality of Attractors.jl to calculate a rate track-return-tip diagram for rs for a variety of forcing duration and forcing scales Δts, Δps with the rate system starting always at u0. Return:

1. a matrix of size length(Δps) × length(Δts) encoding the type of rate tipping behavior.
2. attractors_cont, the attractors of the global continuation of the unforced system.

Keyword arguments

• distance = Centroid(): Distance function used when (1) matching attractors, and (2) mapping the end state of a rate simulation to its closest attractor through the AttractorsViaProximity.
• proximity_kw: Keywords propagated to AttractorsViaProximity, for mapping the rate system state to the unforced attractors at the middle and end (reverse) of forcing. Do not provide a distance as this is used from the above keyword.
• decide_rate_outcome: A three input function f(track_id, return_id, start_id) that returns an integer representing a possible rate tipping case of a simulation with a specific Δt, Δp. It inputs the attractor ID the system converges after the forwards rate forcing, the one it converges after the forcing is reverse, and the ID that u0 converges at the starting parameter of the system.

Description

This function formalizes and generalizes the concept of tracking, returning, or tipping, in rate forced systems introduced in (Ritchie et al., 2023). To achieve this, it uses global continuation using mapper, ics. The function will run a whole global_continuation run over the parameter range prange = p0 .+ Δps to establish the unfrozen system attractors and assign unique IDs to them throughout prange. If you instead want to run the global continuation yourself, which allows you to change Δts, Δps without re-running the global continuation, then simply do:

pcurve = unforced_pcurve(rs, Δps)
matcher = MatchBySSSetDistance(distance)
ascm = AttractorsSeedContinueMatch(mapper, matcher)
_, attractors_cont = global_continuation(ascm, pcurve, ics)
rate_track_return_tip(rs, Δts, Δps, attractors_cont; distance, kw...)

After the global continunation, the function will perform a rate simulation with duration Δt and scale Δp for all combinations. For each, it will then assign an integer corresponding to the type of rate-dependent behaviour as in (Ritchie et al., 2023). By default, these integers are ∈ (1, 2, 3) and mean:

1. Tracking always.
2. Return but not tracking (safe overshoot).
3. Tipping always (failure to track).

You can however provide a custom function that may have more involved decision logic via the keyword decide_rate_outcome. To find the unforced system attractor IDs at the middle and and of the rate simulation an AttractorsViaProximity is used.

Notes

This function applies the same duration and scaling to all parameters that are forced in rs, and enforces all forcings to be reversed as well. The profiles of each parameter are individual though. If any profile does not have the reverse = true option, an error is thrown.