Diffusion operator

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.

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.

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.

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.

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.

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.

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.

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.

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

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

cuboid(lo, hi) -> CriticalTransitions.var"#174#176"

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"#178#179"{_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.