Merge pull request #3132 from ChrisRackauckas-Claude/unify-sde-integrator-methods

Unify SDE integrator methods via field-gated conditionals

Author: ChrisRackauckas

lib/OrdinaryDiffEqCore/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqCore"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "3.19.0"
+version = "3.20.0"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

lib/OrdinaryDiffEqCore/src/dense/generic_dense.jl

@@ -1321,12 +1321,17 @@ function interpolation_differential_vars(differential_vars, y₀, idxs)
     end
 end
 
-# If no dispatch found, assume Hermite
+# If no dispatch found, assume Hermite (or linear when k is empty, e.g. SDE)
 function _ode_interpolant(
         Θ, dt, y₀, y₁, k, cache, idxs, T::Type{Val{TI}}, differential_vars
     ) where {TI}
     TI > 3 && throw(DerivativeOrderNotPossibleError())
 
+    # Linear fallback when no dense output vectors (e.g. SDE)
+    if isempty(k)
+        return linear_interpolant(Θ, dt, y₀, y₁, idxs, T)
+    end
+
     differential_vars = interpolation_differential_vars(differential_vars, y₀, idxs)
     return hermite_interpolant(
         Θ, dt, y₀, y₁, k, Val{cache isa OrdinaryDiffEqMutableCache},
@@ -1339,6 +1344,10 @@ function _ode_interpolant!(
     ) where {TI}
     TI > 3 && throw(DerivativeOrderNotPossibleError())
 
+    if isempty(k)
+        return linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, T)
+    end
+
     differential_vars = interpolation_differential_vars(differential_vars, y₀, idxs)
     return hermite_interpolant!(out, Θ, dt, y₀, y₁, k, idxs, T, differential_vars)
 end

lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl

@@ -17,11 +17,21 @@ function _change_t_via_interpolation!(
         else
             integrator(integrator.u, t)
         end
+        # SDE path: reject noise to rewind W/P, update sqdt
+        W = _get_W(integrator)
+        if !isnothing(W)
+            reject_noise!(W, t - integrator.tprev, integrator.u, integrator.p)
+            reject_noise!(_get_P(integrator), t - integrator.tprev, integrator.u, integrator.p)
+        end
         integrator.t = t
         integrator.dt = integrator.t - integrator.tprev
-        SciMLBase.reeval_internals_due_to_modification!(
-            integrator; callback_initializealg = reinitialize_alg
-        )
+        if isnothing(W)
+            SciMLBase.reeval_internals_due_to_modification!(
+                integrator; callback_initializealg = reinitialize_alg
+            )
+        else
+            integrator.sqdt = sqrt(abs(integrator.dt))
+        end
         if T
             solution_endpoint_match_cur_integrator!(integrator)
         end
@@ -374,6 +384,15 @@ function SciMLBase.set_rng!(integrator::ODEIntegrator, rng)
         )
     end
     integrator.rng = rng
+    # Sync framework-constructed noise processes (SDE only, no-op when W/P are nothing)
+    W = _get_W(integrator)
+    if !isnothing(W) && integrator.noise === nothing
+        W.rng = rng
+    end
+    P = _get_P(integrator)
+    if !isnothing(P)
+        P.rng = rng
+    end
     return nothing
 end
 
@@ -435,8 +454,10 @@ function SciMLBase.reinit!(
 
     tType = typeof(integrator.t)
     tspan = (tType(t0), tType(tf))
-    reinit_tstops!(tType, integrator.opts.tstops, tstops, d_discontinuities, tspan;
-        p = parameter_values(integrator))
+    reinit_tstops!(
+        tType, integrator.opts.tstops, tstops, d_discontinuities, tspan;
+        p = parameter_values(integrator)
+    )
     reinit_saveat!(tType, integrator.opts.saveat, saveat, tspan)
     reinit_d_discontinuities!(tType, integrator.opts.d_discontinuities, d_discontinuities, tspan)
     if erase_sol
@@ -447,7 +468,9 @@ function SciMLBase.reinit!(
         end
         resize!(integrator.sol.u, resize_start)
         resize!(integrator.sol.t, resize_start)
-        resize!(integrator.sol.k, resize_start)
+        if _has_ks(integrator)
+            resize!(integrator.sol.k, resize_start)
+        end
 
         if integrator.opts.save_start || (!isempty(saveat) && saveat[1] == tType(t0))
             copyat_or_push!(integrator.sol.t, 1, t0)
@@ -461,7 +484,7 @@ function SciMLBase.reinit!(
         if integrator.sol.u_analytic !== nothing
             resize!(integrator.sol.u_analytic, 0)
         end
-        if integrator.alg isa OrdinaryDiffEqCompositeAlgorithm
+        if is_composite_algorithm(integrator.alg)
             resize!(integrator.sol.alg_choice, resize_start)
         end
         integrator.saveiter = resize_start
@@ -505,17 +528,25 @@ function SciMLBase.reinit!(
     if reinit_retcode
         integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.Default)
     end
+
+    reinit_noise!(_get_W(integrator), integrator.dt)
     return nothing
 end
 
-function SciMLBase.auto_dt_reset!(integrator::ODEIntegrator)
-    integrator.dt = ode_determine_initdt(
+# Extensible initdt hook: ODE defaults to ode_determine_initdt.
+# SDE extends for SDE algorithm types to call sde_determine_initdt.
+function _determine_initdt(integrator)
+    return ode_determine_initdt(
         integrator.u, integrator.t,
         integrator.tdir, integrator.opts.dtmax,
         integrator.opts.abstol, integrator.opts.reltol,
         integrator.opts.internalnorm, integrator.sol.prob,
         integrator
     )
+end
+
+function SciMLBase.auto_dt_reset!(integrator::ODEIntegrator)
+    integrator.dt = _determine_initdt(integrator)
     integrator.dtpropose = integrator.dt
     return increment_nf!(integrator.stats, 2)
 end

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

@@ -5,6 +5,7 @@ reject_noise!(::Nothing, args...) = nothing
 save_noise!(::Nothing) = nothing
 noise_curt(::Nothing) = nothing
 is_noise_saveable(::Nothing) = false
+reinit_noise!(::Nothing, dt) = nothing
 
 # Noise field accessors — safe for any integrator type.
 # ODEIntegrator has W/P/sqdt; other integrators (DDEIntegrator) don't.