Developer API

AbstractJ

Supertype for the objects that evaluate the iterated stochastic integrals J needed by higher order SDE solvers.

All stochastic iterated integrals are written in the Stratonovich sense, as indicated by the J. An AbstractJ is produced from a problem and an algorithm by get_Jalg and is then evaluated with get_iterated_I (out-of-place) or get_iterated_I! (in-place, storing the result in the object's J field).

Subtypes: AbstractJDiagonal, AbstractJCommute, and the Lévy area algorithms of StochasticDiffEqLevyArea together with IteratedIntegralAlgorithm_iip.

AbstractJCommute <: AbstractJ

Iterated-integral evaluators for commutative noise, where the Lévy area terms cancel and the integrals reduce to the outer product 1/2 ΔW ΔWᵀ.

Subtypes: JCommute_oop, JCommute_iip.

AbstractJDiagonal <: AbstractJ

Iterated-integral evaluators for diagonal (or scalar) noise, where the integrals reduce to 1/2 ΔWᵢ² componentwise and no Lévy area approximation is required.

Subtypes: JDiagonal_oop, JDiagonal_iip.

AutoAlgSwitch(nonstiffalg, stiffalg; kwargs...) -> StochasticCompositeAlgorithm

Build a stiffness-switching composite of nonstiffalg and stiffalg.

This is the constructor behind the Auto* SDE solvers (for example AutoSOSRI2): it wraps the two algorithms in a StochasticCompositeAlgorithm whose choice function is an AutoSwitch. All keyword arguments are forwarded to AutoSwitch.

alg = AutoAlgSwitch(SOSRI(), ImplicitEM(); maxstiffstep = 5)

AutoSwitch(
    nonstiffalg, stiffalg; maxstiffstep = 10, maxnonstiffstep = 3,
    nonstifftol = 9 // 10, stifftol = 9 // 10, dtfac = 2,
    stiffalgfirst = false, switch_max = 5
)

Choice function that switches a StochasticCompositeAlgorithm between a nonstiff and a stiff method based on an online stiffness estimate.

Each step the estimate abs(eigen_est * dt / alg_stability_size(alg)) is compared against a tolerance. Consecutive verdicts are counted, and the algorithm is only switched once the count exceeds the corresponding threshold, which prevents thrashing on a borderline problem.

Keyword Arguments

• maxstiffstep: Number of consecutive stiff verdicts required before switching to stiffalg (default: 10).
• maxnonstiffstep: Number of consecutive nonstiff verdicts required before switching back to nonstiffalg (default: 3).
• nonstifftol, stifftol: Stiffness-ratio tolerances used while the nonstiff and the stiff algorithm is active, respectively (default: 9//10 for both).
• dtfac: Factor applied to dt at a switch — dt is multiplied by it when moving to the stiff method and divided by it when moving back (default: 2).
• stiffalgfirst: Start with stiffalg instead of nonstiffalg (default: false).
• switch_max: Number of successive nonstiff verdicts after which the error check is re-enabled (default: 5).

Use AutoAlgSwitch to build the composite algorithm directly.

DiffCache(u)
DiffCache(u, nlsolve)
DiffCache(u, ::Type{Val{chunk_size}})
DiffCache(T, size, ::Type{Val{chunk_size}})

Dual-buffered cache holding both a plain array and a ForwardDiff.Dual array of the same shape.

An in-place function that must work both on ordinary numbers and on duals cannot use a single preallocated buffer, because the element types differ. DiffCache stores one buffer of each and SciMLBase.get_du picks the right one from the element type at the call site.

The chunk size defaults to ForwardDiff.pickchunksize(length(u)), or is taken from the nonlinear solver via get_chunksize when one is supplied.

DiffEqNLSolveTag

ForwardDiff tag type used for the Jacobians of the nonlinear solves inside the SDE integrators.

Tagging the duals keeps derivatives taken by the integrator distinguishable from derivatives a user takes through solve, which is what allows nested differentiation to work (see the ForwardDiff.jl documentation on tags).

IICommutative()

Iterated-integral strategy that assumes the noise is commutative.

Under commutative noise the Lévy area terms cancel and the iterated integrals reduce to 1/2 ΔW ΔWᵀ, which is cheap and exact for that problem class. Selecting this for a genuinely non-commutative problem silently loses the method's strong order, so use IILevyArea unless commutativity is known to hold.

IIFNLSolveFunc(f)

Wrapper holding the residual function f of an IIF nonlinear solve.

NLSOLVEJL_SETUP returns one of these from its Val{:init} call so that the in-place residual f(resid, u) can later be turned into a NonlinearProblem without recapturing the surrounding cache.

IILevyArea()

Iterated-integral strategy that approximates the Lévy area numerically.

This is the general-purpose choice: the number of terms is selected automatically from the step size and the requested accuracy, and an appropriate simulation algorithm is picked by StochasticDiffEqLevyArea.optimal_algorithm. It is more expensive than IICommutative but does not assume anything about the noise structure.

Ihat2(...)

Approximation of the second-order multiple stochastic integral used by the weak order 2 schemes.

Only the generic function lives here; the methods are defined in StochasticDiffEqWeak alongside the cache types they dispatch on, because their argument lists differ between the DRI/RI (Rößler) and RDI families.

IteratedIntegralAlgorithm_iip(ΔW, levyalg)

In-place wrapper around a StochasticDiffEqLevyArea iterated-integral algorithm.

StochasticDiffEqLevyArea's algorithms are allocating; this wrapper pairs one of them (levyalg) with a preallocated length(ΔW) × length(ΔW) buffer J, so that get_iterated_I! can serve in-place solvers with general non-commutative noise.

IteratedIntegralApprox

Supertype for the strategies used to approximate the iterated stochastic integrals (Lévy areas) that higher order SDE solvers need for non-diagonal noise.

An algorithm carries its choice in the ii_approx field, and get_Jalg turns that choice plus the problem's noise structure into the concrete AbstractJ object that computes the integrals.

Subtypes: IICommutative, IILevyArea.

JCommute_iip(ΔW)

In-place iterated-integral evaluator for commutative noise, preallocated to a length(ΔW) × length(ΔW) matrix.

get_iterated_I! writes 1/2 * vec(ΔW) * vec(ΔW)' into the J field.

JCommute_oop()

Out-of-place iterated-integral evaluator for commutative noise.

get_iterated_I returns a freshly allocated 1/2 .* vec(ΔW) .* vec(ΔW)'.

JDiagonal_iip(ΔW)

In-place iterated-integral evaluator for diagonal noise, preallocated to the shape of the Brownian increment ΔW.

get_iterated_I! writes 1/2 * ΔW^2 into the J field.

JDiagonal_oop()

Out-of-place iterated-integral evaluator for diagonal noise.

get_iterated_I returns a freshly allocated 1/2 .* ΔW .* ΔW.

NLSOLVEJL_SETUP(; autodiff = AutoForwardDiff())

Nonlinear-solver setup used by the IIF (implicit integrating factor) methods.

The name is historical — the solve is no longer performed by NLsolve.jl but by SimpleTrustRegion from SimpleNonlinearSolve.jl. A setup object is callable in two ways: setup(Val{:init}, f, u0_prototype) wraps f in an IIFNLSolveFunc, and setup(wrapped_f, u0) solves f(resid, u) = 0 starting from u0 and returns the root.

Keyword Arguments

• autodiff: ADTypes.jl backend used for the Jacobian of the inner solve (default: AutoForwardDiff()).

SDEAlgTypes

Union of the algorithm supertypes that the SDE integrator loop handles, i.e. StochasticDiffEqAlgorithm and StochasticDiffEqRODEAlgorithm.

Used to constrain SDEIntegrator and to write methods that apply to both SDE and RODE algorithms at once.

SDEIntegrator

The integrator type for SDE and RODE problems: an ODEIntegrator specialized to an algorithm in SDEAlgTypes.

SDE integration reuses OrdinaryDiffEqCore's integrator machinery and adds the noise process in the W field (and P for jump problems), so this alias is what SDE-specific methods — resizing, noise cache handling, the interpolation and step interface — dispatch on.

SDEOptions

Alias for OrdinaryDiffEqCore.DEOptions, the container holding the solver options (tolerances, callbacks, save settings) of a running integrator.

Kept as its own name for downstream packages such as StochasticDelayDiffEq that were written against the pre-migration StochasticDiffEq.SDEOptions.

StochasticCompositeCache(caches, choice_function, current)

Cache of a StochasticCompositeAlgorithm.

Holds one member cache per member algorithm in caches, the choice_function that selects among them, and current, the index chosen for the step being taken. unwrap_alg and get_current_alg_order read current so that perform_step! and the adaptivity see the member that is actually running.

TauLeapingDrift{C, R, RateCache, IIP}(c, rate, rate_cache)

Callable that presents the drift of a tau-leaping jump problem as an ordinary ODE right-hand side.

The drift is ν ⋅ a(u) written in terms of the RegularJump pieces: rate is the propensity function a and c is the stoichiometry function that applies the jump counts. The wrapper evaluates c(u, p, t, rate(u, p, t), nothing), which is exactly the function the implicit tau-leaping methods hand to their nonlinear solver.

The IIP type parameter selects the calling convention: out-of-place instances are called as drift(u, p, t), in-place ones as drift(du, u, p, t) and use rate_cache as scratch space for the propensities.

_resolve_rng(rng, seed, prob) -> (rng, seed, rng_provided)

Resolve the RNG and seed for an SDE/RODE integration from the user-provided rng and seed kwargs plus the problem's stored seed.

_sde_init(prob, alg; kwargs...) -> SDEIntegrator

Build the SDEIntegrator for prob and alg.

This is the body of SciMLBase.__init for SDE, RODE, and jump problems: it resolves the options and the RNG, allocates the solution object, the algorithm cache (alg_cache) and the noise process, and returns the integrator positioned at the initial condition. solve is init followed by solve!, so every documented solve keyword is accepted here.

It is exposed separately from __init so that downstream packages which build on the SDE integrator (for example StochasticDelayDiffEq) can construct one without going through SciMLBase.__init dispatch.

_z_prototype(alg, rand_prototype, iip::Bool) -> rand_prototype2

Compute the Z process prototype for algorithms that need an extra Brownian process. Default: use rand_prototype itself as Z prototype (same shape as W). Solver subpackages override this for algorithms with special Z requirements (e.g., PL1WM, RKMilGeneral, W2Ito1).

_z_prototype(alg, rand_prototype, iip::Bool, dt) -> rand_prototype2

Step-size aware form of _z_prototype, used when the size of the Z process depends on the step size — as it does for Lévy area truncations, where the number of retained series terms is chosen from the accuracy the step size demands.

dt is the initial step size as passed to solve, which is zero when the caller left it to be determined automatically. Since the Z prototype must be built before the initial step size is known, an algorithm that sizes from dt needs a fallback for that case.

Defaults to the three-argument form, so an override only has to be added by algorithms that actually need the step size.

addat_noise!(integrator, cache, idxs) -> nothing

Insert new components at positions idxs in the noise process of integrator.

The counterpart of deleteat_noise!, called from addat! on the integrator. Space is made in every noise buffer and in the cached increments of the interpolation stacks, and the new slots are filled with freshly sampled increments through fill_new_noise_caches!.

alg_cache(
    alg, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype,
    jump_rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits},
    ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{iip}}, verbose
)

Construct the per-algorithm cache used by perform_step!.

This is the main extension point for SDE solver subpackages: each algorithm defines a cache type (usually with the @cache macro, which also generates the full_cache, rand_cache, and ratenoise_cache accessors used for resize!) and one alg_cache method that allocates it. The Val{iip} argument selects the in-place or out-of-place variant, and the *NoUnits type arguments give the element types the error estimates should be computed in.

The fallback method defined here errors, pointing at the solver subpackage that needs to be loaded for alg.

alg_can_repeat_jac(alg) -> Bool

Whether alg may reuse a Jacobian across a rejected-and-repeated step.

When true the integrator keeps the factorized W matrix after a step rejection instead of recomputing it, which is the common case. Algorithms whose stage structure changes between attempts must override this to false.

alg_compatible(prob, alg) -> Bool

Whether alg can solve prob.

solve consults this trait before setting up an integrator and errors with a descriptive message when it returns false. Solver subpackages add methods for the problem classes their algorithm supports, for example a Milstein-type method that requires diagonal noise:

alg_compatible(prob::SciMLBase.AbstractSDEProblem, alg::RKMil) = is_diagonal_noise(prob)

For a StochasticCompositeAlgorithm the result is the maximum over its members, and for a JumpProblem compatibility additionally requires that the algorithm supports_regular_jumps whenever the problem carries a RegularJump.

alg_control_rate(alg) -> Bool

Whether the adaptive step-size controller of alg acts on the jump/leap rate rather than on the usual solution error estimate.

Tau-leaping methods choose dt from how much the propensities are allowed to change over a step, so they override this to true and the integrator routes step-size control through the leap-specific controller instead of the standard EEst path.

alg_mass_matrix_compatible(alg) -> Bool

Whether alg can solve a problem with a non-identity mass matrix.

solve errors when a mass matrix is supplied to an algorithm for which this is false. The implicit theta-method solvers accept a mass matrix only in the configurations for which the discretization stays consistent (symplectic, or theta == 1) and throw a descriptive error otherwise.

alg_needs_extra_process(alg) -> Bool

Whether alg needs a second noise process ΔZ alongside the Brownian increment ΔW.

Higher order schemes (for example the Rößler SRA/SRI families) require an auxiliary independent process to approximate the extra stochastic integrals appearing in their order conditions. When this trait is true the integrator allocates and resizes the Z process in addition to W, so solver subpackages must set it for any algorithm that reads integrator.W.dZ.

For a StochasticCompositeAlgorithm the result is the maximum over its members, so the extra process is present if any member needs it.

alg_stability_size(alg) -> Real

Radius of the (deterministic) stability region of alg along the negative real axis.

This is the scale used by the automatic stiffness detection of AutoSwitch: the stiffness ratio is abs(eigen_est * dt / alg_stability_size(alg)), so an algorithm participating in a stiffness-switching composite must report a nonzero value. The default of 0 means "not characterized", and is overridden per algorithm in the solver subpackages.

calc_threepoint_random(sq3dt, quantile, dW_scaled)

Three-point discrete random variable replacing the Brownian increment in a weak (moment-matching) scheme.

Returns -sq3dt, 0, or +sq3dt depending on where the standardized increment dW_scaled falls relative to quantile. Unlike the two-point variable of calc_twopoint_random this matches the fourth moment of the Gaussian as well, which weak order 2 schemes require.

sq3dt is sqrt(3dt), and quantile is the (negative) standard-normal quantile that gives the zero outcome its correct probability.

calc_threepoint_random!(_dW, sq3dt, quantile, dW_scaled) -> nothing

In-place form of calc_threepoint_random, writing the discrete increments into _dW.

calc_twopoint_random(sqdt, dW)

Two-point discrete random variable replacing the Brownian increment in a weak (moment-matching) scheme.

Returns ±sqdt, taking the sign of dW, so the result has the same mean and variance as the Brownian increment while taking only two values. Weak-order methods may substitute such a variable because they only need the moments of the noise to match, and sampling two points is cheaper than sampling a Gaussian.

sqdt is sqrt(dt). See calc_threepoint_random for the variant that also matches the fourth moment.

calc_twopoint_random!(_dW, sqdt, dW) -> nothing

In-place form of calc_twopoint_random, writing the discrete increments into _dW.

concrete_prob(prob) -> prob

The underlying differential equation problem carried by prob.

For an SDEProblem, RODEProblem, or DiscreteProblem this is prob itself; for a JumpProblem it is the wrapped prob.prob. Setup code that needs to inspect the problem's f, u0, or tspan goes through this so it works uniformly with and without a jump wrapper.

deleteat_noise!(integrator, cache, idxs) -> nothing

Delete the components idxs from the noise process of integrator.

The counterpart of addat_noise!, called from deleteat! on the integrator: the same components are removed from every noise buffer and from the cached increments in the interpolation stacks, keeping the process dimension in step with the shrunken state.

delta_default(alg) -> Real

Default value of the delta option for alg.

delta is the SDE-specific mixing parameter of the adaptive error estimate: it weighs the drift (deterministic) error contribution against the diffusion (stochastic) one when the two are combined into a single EEst. A value of 1 uses the drift estimate at full strength.

Solver subpackages override this for algorithms whose published adaptivity scheme prescribes a different weighting.

determine_chunksize(u, alg) -> Int
determine_chunksize(u, CS) -> Int

ForwardDiff chunk size to use when differentiating with respect to a state like u.

A nonzero explicitly configured chunk size CS (or the one reported by get_chunksize for alg) is used as-is; 0 means "not configured" and falls back to ForwardDiff.pickchunksize(length(u)).

fill_new_noise_caches!(integrator, c, scaling_factor, idxs) -> nothing

Sample fresh noise increments into positions idxs of a cached noise entry c.

c is one entry of the noise process's interpolation stacks (W.S₁/W.S₂), a tuple whose first element is the step's scaling factor and whose remaining elements hold the cached ΔW and, when alg_needs_extra_process holds, ΔZ values. The new values are drawn from the process's own distribution so that the extended path has the correct law.

Used by resize_noise! and addat_noise! whenever the state grows.

get_Jalg(ΔW, dt, prob, alg) -> AbstractJ

Select the iterated-integral evaluator to use for prob with alg.

The choice combines the algorithm's ii_approx field — an IteratedIntegralApprox — with the problem's noise structure and whether it is in-place:

• IILevyArea uses a diagonal evaluator when the noise is diagonal or scalar, and otherwise the StochasticDiffEqLevyArea algorithm chosen by optimal_algorithm(length(ΔW), dt) (wrapped in IteratedIntegralAlgorithm_iip for in-place problems).
• IICommutative uses a diagonal evaluator for diagonal/scalar noise and a commutative evaluator otherwise.
• Any other ii_approx is taken to be a concrete Lévy area algorithm and is used directly.

A solver subpackage can override the choice for a specific algorithm by adding a method, e.g. get_Jalg(ΔW, dt, prob, alg::MySolver) = MronRoe().

get_chunksize(x) -> Int

ForwardDiff chunk size configured on x, or 0 when x does not configure one.

0 is the "unset" sentinel that makes determine_chunksize fall back to ForwardDiff.pickchunksize.

get_current_alg_order(alg, cache) -> Real

Strong order of the algorithm that is currently being stepped with.

For a plain algorithm this is just alg_order(alg). For a StochasticCompositeAlgorithm it is the order of the member algorithm selected by cache.current, so that adaptivity uses the order of the method that actually took the step rather than the order of the composite as a whole.

get_iterated_I(dt, dW, dZ, alg, p = nothing, c = 1, γ = 1 // 1)

Evaluate and return the Stratonovich iterated stochastic integrals for the step, out of place.

Arguments

• dt: The step size.
• dW, dZ: The Brownian increment and, where the method needs one, the auxiliary increment for the step.
• alg: The AbstractJ evaluator, as produced by get_Jalg.
• p: Number of terms in the Lévy area series. nothing (the default) selects it automatically from the requested accuracy ε = c * dt^(γ + 1/2).
• c, γ: Constant and exponent of that accuracy target, normally the constant and strong order of the calling solver.

See get_iterated_I! for the in-place form.

get_iterated_I!(dt, dW, dZ, alg, p = nothing, c = 1, γ = 1 // 1)

In-place form of get_iterated_I: compute the Stratonovich iterated stochastic integrals for the step and store them in alg.J, returning nothing.

is_split_step(alg) -> Bool

Whether alg is a split-step method, i.e. one that advances the drift to an intermediate state before applying the diffusion rather than adding both contributions to the same base point.

The split-step solvers (ISSEM, ISSEulerHeun) override this to true; the integrator uses it to select the matching residual and error-estimate paths.

resize_noise!(integrator, cache, bot_idx, i) -> nothing

Resize the noise process of integrator to length i after the state was resized.

All of the noise process's buffers (dW, dWtilde, dWtmp, curW, the dZ family when alg_needs_extra_process holds, and the cached increments in the interpolation stacks S₁/S₂) are grown or shrunk together. New entries from bot_idx up to i are filled with freshly sampled increments via fill_new_noise_caches! so the process stays consistent with the already generated path.

This is called from the resize! integrator interface; see deleteat_noise! and addat_noise! for the index-wise variants.

supports_regular_jumps(alg) -> Bool

Whether alg can integrate a JumpProblem that carries a RegularJump.

Regular jumps are leaped over with a Poisson count per step rather than simulated one event at a time, which only some SDE integrators implement. The default is false; EM and ImplicitEM override it to true in their subpackages. alg_compatible uses this trait to reject a RegularJump problem paired with an algorithm that cannot handle it.

unwrap_alg(integrator, is_nlsolve) -> alg

The concrete algorithm that integrator should use for the current step.

For a plain algorithm this is integrator.alg itself. For a StochasticCompositeAlgorithm it is the currently selected member: the stiff member when is_nlsolve is true under AutoSwitch-style stiffness switching, and the member indicated by integrator.cache.current otherwise.

perform_step! implementations call this instead of reading integrator.alg directly so that they work unchanged inside a composite algorithm.

checkSRAOrder(RosslerSRI)

Determines whether the order conditions are met via the tableaus of the SRA method.

checkSRIOrder(RosslerSRI)

Determines whether the order conditions are met via the tableaus of the SRI method.

constructExplicitSKenCarp()

Constructs the tableau type for the explicit part of SKenCarp as a RosslerSRA tableau.

constructSKenCarp()

Constructs the tableau type for the implicit SKenCarp method as a RosslerSRA tableau.

constructSOSRA()

Constructs the taleau type for the SOSRA method.

constructSOSRA2()

Constructs the taleau type for the SOSRA method.

constructSRA1()

Constructs the taleau type for the SRA1 method.

constructSRA2()

Constructs the taleau type for the SRA2 method.

constructSRA3()

Constructs the taleau type for the SRA3 method.

constructSRIOpt1()

Opti6-12-11-10-01-47

constructSRIOpt2()

Opti6-12-11-10-01-47

constructSRIW1()

Constructs the tableau type for the SRIW1 method.

constructSRIW2()

Constructs the tableau type for the SRIW1 method.

du_cache(cache::SRACache)

Return derivative and noise-work buffers stored by an SRACache. This is a developer extension hook for the SRA cache implementation.

Returns

A tuple of internal work buffers. Its order and contents are a solver-internal contract and must not be used by application code.

u_cache(cache::SRACache)

Return additional mutable state buffers for an SRACache. Solver implementations use this developer hook when resizing or inspecting a high-order SRA cache.

Returns

An empty tuple for SRACache; user code must not rely on cache-hook tuple layouts.

user_cache(cache::SRACache)

Return the user-state buffers stored by an SRACache. This developer hook is used by solver infrastructure; application code should not depend on its tuple layout.

Returns

A tuple containing the current and previous user states plus internal temporaries. Its order and contents are only a solver-extension contract.