Handle OrdinaryDiffEq v7 / DiffEqBase v7 / OrdinaryDiffEqCore v4 breaking changes (#410)

* Handle OrdinaryDiffEq v7 / DiffEqBase v7 / OrdinaryDiffEqCore v4 breaking changes

Bump compats to the new major versions (matching CompatHelper PRs #407, #408, #409)
and gate the OrdinaryDiffEqCore internals we extend so the package continues to build
against both OrdinaryDiffEqCore v3 and v4.

v4 removed a handful of internal hooks that ProbNumDiffEq relied on:

- `OrdinaryDiffEqCore._default_dae_init!(integrator, prob, x, alg)` was deleted
  along with the per-solver DAE-init dispatch; initialization now routes through
  `_initialize_dae!(integrator, prob, alg::DefaultInit, x)` and the default
  algorithm switched from `BrownFullBasicInit` to `CheckInit`. Restore the prior
  behaviour for `AbstractEK` by overriding `_initialize_dae!` for our integrator
  type on v4.
- `OrdinaryDiffEqCore.isstandard` / `ispredictive` traits were removed; PI control
  is now the default via `default_controller`, so the trait overrides are only
  needed on v3.
- `integrator.EEst` was removed from the `ODEIntegrator` struct and replaced with
  `OrdinaryDiffEqCore.get_EEst` / `set_EEst!` accessors backed by
  `integrator.controller_cache`. Wrap both call sites in `perform_step!` with
  helpers that pick the right API per version.
- `OrdinaryDiffEqCore._process_AD_choice` is gone on v4 (the `chunk_size` /
  `diff_type` / `standardtag` solver kwargs were removed in favour of
  `AutoForwardDiff(chunksize=…)` / `AutoFiniteDiff(fdtype=…)`). Fall back to an
  inlined equivalent so `EK1` / `DiagonalEK1` keep accepting the legacy kwargs
  on both versions.

See the OrdinaryDiffEq v7 NEWS.md for the full list of breaking changes.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

* JuliaFormatter

* Fix v4/v7 stack CI failures: scope reduction, RAT v4, Bool autodiff, DiffEqDevTools v4

Four distinct failures surfaced once CI started resolving the v4/v7 stack
(which became possible after the gating deps registered in General).

1. UndefVarError: RadauIIA5. OrdinaryDiffEq v7's package scope reduction
   means `using OrdinaryDiffEq` no longer re-exports the full solver suite.
   Add `using OrdinaryDiffEqFIRK: RadauIIA5` in test/stiff_problem.jl and
   test/mass_matrix.jl, plus the dep UUID in test/Project.toml.

2. length(sol) returned the wrong value at test/solution.jl:29,96.
   RecursiveArrayTools v4 changed length(sol) from "number of timesteps"
   to prod(size(sol)) since AbstractVectorOfArray now subtypes AbstractArray.
   Switch to length(sol.u) per the v7 NEWS migration.

3. MethodError: prepare_ADType(::Bool, ...) when EK1 or DiagonalEK1 was
   constructed with the legacy autodiff=true or autodiff=false kwarg.
   The v4 fallback for _process_AD_choice in src/algorithms.jl only handled
   AutoForwardDiff and AutoFiniteDiff; the v3 OrdinaryDiffEqCore version
   it replaced also converted Bool to an ADType. Mirror that conversion so
   the Bool no longer propagates as the AD type parameter on the algorithm.

4. ArgumentError: Passing a Bool for verbose, raised by DiffEqDevTools v3's
   WorkPrecisionSet on test/diffeqdevtools.jl:79. DiffEqDevTools v4 removed
   the offending kwarg. Widen Project.toml compat to "2.46, 3, 4" so the
   resolver picks v4 on the v4/v7 stack.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

* Switch default DAE init to CheckInit; drop verbose=false from broken WP test

Two changes:

1. src/algorithms.jl: replace the AbstractEK DAE-init override with one that
   routes through CheckInit on both stacks. Matches the OrdinaryDiffEq v7 /
   DifferentialEquations v8 default (errors on inconsistent ICs instead of
   silently correcting them). Users who want the old auto-fix behavior must
   now pass initializealg=BrownFullBasicInit() explicitly. On v4 the per-alg
   hook is gone and CheckInit is already the default, so the override is only
   needed on v3. Verified locally against test/mass_matrix.jl on the v4/v7
   stack: both UniformScaling and Robertson mass-matrix tests pass with
   consistent ICs.

2. test/diffeqdevtools.jl: drop verbose=false from the two WorkPrecisionSet
   calls in the "is broken" testset. DiffEqDevTools v3 forwards verbose as a
   Bool which OrdinaryDiffEq v7 rejects (DEVerbosity / SciMLLogging required).
   DiffEqDevTools v4 removed the kwarg, but v4.0.0 is yanked in General, so
   we drop the kwarg test-side. On v3/v6 this only removes a deprecation
   warning the test was already eating.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

* Address code review feedback: fix AutoFiniteDiff regression, trim comments

Two functional changes in src/algorithms.jl:

- Drop the stray em-dash from the `_process_AD_choice` overview comment to
  comply with the project's no-em-dashes style.

- Guard the `AutoFiniteDiff` branch of the v4 `_process_AD_choice` fallback
  to only rebuild the ADType when the user passed a non-default legacy
  `diff_type` kwarg. The previous version unconditionally rebuilt with
  `fdtype=diff_type`, which silently dropped a user's pre-configured
  `AutoFiniteDiff(fdtype=Val(:central))` (or similar) when they didn't also
  set the legacy `diff_type`. Mirrors v3's behavior.

Comment trim across src/algorithms.jl and src/alg_utils.jl: the
`_alg_autodiff` block, the `_process_AD_choice` overview, the Bool-branch
inline comment, and the DAE-init block are all condensed without losing
the migration rationale. Saves ~10 lines in the compat-shim region.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

---------

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Nathanael Bosch <nathanael.bosch@epfl.ch>
Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: 4 people

Project.toml

@@ -53,9 +53,9 @@ RecipesBaseExt = "RecipesBase"
 ADTypes = "1.16"
 ArrayAllocators = "0.3"
 BlockArrays = "1"
-DiffEqBase = "6.194"
+DiffEqBase = "6.194, 7"
 DiffEqCallbacks = "4"
-DiffEqDevTools = "2.46"
+DiffEqDevTools = "2.46, 3, 4"
 DocStringExtensions = "0.9.4"
 FastBroadcast = "0.3.5, 1"
 FillArrays = "1.9"
@@ -66,10 +66,10 @@ Kronecker = "0.5.4"
 LinearAlgebra = "1"
 MatrixEquations = "2"
 Octavian = "0.3.29"
-OrdinaryDiffEqCore = "3"
-OrdinaryDiffEqDifferentiation = "2"
-OrdinaryDiffEqRosenbrock = "1.23"
-OrdinaryDiffEqVerner = "1.10"
+OrdinaryDiffEqCore = "3, 4"
+OrdinaryDiffEqDifferentiation = "2, 3"
+OrdinaryDiffEqRosenbrock = "1.23, 2"
+OrdinaryDiffEqVerner = "1.10, 2"
 PSDMatrices = "0.5"
 PrecompileTools = "1.2.1"
 Printf = "1"

src/alg_utils.jl

@@ -3,7 +3,14 @@
 # https://github.com/SciML/OrdinaryDiffEqCore.jl/blob/master/src/alg_utils.jl
 ############################################################################################
 
-OrdinaryDiffEqDifferentiation._alg_autodiff(::AbstractEK) = Val{true}()
+# v3 reaches autodiff settings via the private `_alg_autodiff`; v4 promoted
+# `alg_autodiff` to the public entry point and deleted `_alg_autodiff`. Define both:
+# the `_alg_autodiff` override is gated so v4 doesn't see it, and `alg_autodiff` is
+# harmlessly redefined on v3 (an existing function gets a new method).
+@static if isdefined(OrdinaryDiffEqDifferentiation, :_alg_autodiff)
+    OrdinaryDiffEqDifferentiation._alg_autodiff(::AbstractEK) = Val{true}()
+end
+OrdinaryDiffEqDifferentiation.alg_autodiff(::AbstractEK) = ADTypes.AutoForwardDiff()
 OrdinaryDiffEqDifferentiation.standardtag(::AbstractEK) = false
 OrdinaryDiffEqDifferentiation.concrete_jac(::AbstractEK) = nothing
 
@@ -12,8 +19,11 @@ OrdinaryDiffEqDifferentiation.concrete_jac(::AbstractEK) = nothing
 OrdinaryDiffEqCore.isfsal(::AbstractEK) = false
 
 for ALG in [:EK1, :DiagonalEK1]
-    @eval OrdinaryDiffEqDifferentiation._alg_autodiff(alg::$ALG{CS,AD}) where {CS,AD} =
-        alg.autodiff
+    @static if isdefined(OrdinaryDiffEqDifferentiation, :_alg_autodiff)
+        @eval OrdinaryDiffEqDifferentiation._alg_autodiff(alg::$ALG{CS,AD}) where {CS,AD} =
+            alg.autodiff
+    end
+    @eval OrdinaryDiffEqDifferentiation.alg_autodiff(alg::$ALG) = alg.autodiff
     @eval OrdinaryDiffEqDifferentiation.alg_difftype(
         ::$ALG{CS,AD,DiffType},
     ) where {CS,AD,DiffType} =
@@ -35,9 +45,13 @@ OrdinaryDiffEqCore.isadaptive(::AbstractEK) = true
 OrdinaryDiffEqCore.alg_order(alg::AbstractEK) = num_derivatives(alg.prior)
 # OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractEK) =
 
-# PI control is the default!
-OrdinaryDiffEqCore.isstandard(::AbstractEK) = false # proportional
-OrdinaryDiffEqCore.ispredictive(::AbstractEK) = false # not sure, maybe Gustafsson acceleration?
+# PI control is the default. On OrdinaryDiffEqCore v3 we have to explicitly set the
+# `isstandard`/`ispredictive` traits to false; on v4 those traits were removed and the
+# default is PI already (via `default_controller`), so we only opt in on the old version.
+@static if isdefined(OrdinaryDiffEqCore, :isstandard)
+    OrdinaryDiffEqCore.isstandard(::AbstractEK) = false # proportional
+    OrdinaryDiffEqCore.ispredictive(::AbstractEK) = false # not sure, maybe Gustafsson acceleration?
+end
 
 # OrdinaryDiffEqCore.qmin_default(alg::AbstractEK) =
 # OrdinaryDiffEqCore.qmax_default(alg::AbstractEK) =

src/algorithms.jl

@@ -3,6 +3,49 @@
 ########################################################################################
 abstract type AbstractEK <: OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end
 
+# `_process_AD_choice` normalizes the legacy `(autodiff, chunk_size, diff_type)`
+# kwargs into an `ADTypes.AbstractADType`. v3's OrdinaryDiffEqCore exports this
+# helper; v4 (OrdinaryDiffEq v7) removed it along with the legacy kwargs. We keep
+# accepting the kwargs on EK1/DiagonalEK1 by delegating to OrdinaryDiffEqCore on v3
+# and inlining an equivalent shim on v4.
+@static if isdefined(OrdinaryDiffEqCore, :_process_AD_choice)
+    _process_AD_choice(autodiff, chunk_size, diff_type) =
+        OrdinaryDiffEqCore._process_AD_choice(autodiff, chunk_size, diff_type)
+else
+    function _process_AD_choice(autodiff, chunk_size, diff_type)
+        ad = autodiff
+        if ad isa Bool
+            # Mirror v3: `true` → AutoForwardDiff(chunksize=chunk_size), `false` →
+            # AutoFiniteDiff(fdtype=diff_type, dir=1). Without this a Bool would
+            # propagate as EK1's AD type parameter and break `AbstractADType` dispatch.
+            if ad
+                cs_int = _unwrap_val(chunk_size)
+                cs = cs_int == 0 ? nothing : cs_int
+                ad = ADTypes.AutoForwardDiff(; chunksize=cs)
+            else
+                fdtype = diff_type isa Val ? diff_type : Val(diff_type)
+                ad = ADTypes.AutoFiniteDiff(; fdtype=fdtype, dir=1)
+            end
+        elseif ad isa ADTypes.AutoForwardDiff
+            cs_int = _unwrap_val(chunk_size)
+            # If the user set an explicit chunk_size kwarg, fold it into the ADType.
+            # Otherwise defer to whatever chunksize is already on `ad`.
+            if cs_int != 0
+                ad = ADTypes.AutoForwardDiff(; chunksize=cs_int, tag=ad.tag)
+            end
+        elseif ad isa ADTypes.AutoFiniteDiff
+            # Only override the ADType's fdtype if the user passed a non-default
+            # legacy `diff_type` kwarg; otherwise preserve their pre-configured
+            # AutoFiniteDiff (matching v3's behavior).
+            fdtype = diff_type isa Val ? diff_type : Val(diff_type)
+            if fdtype !== Val(:forward)
+                ad = ADTypes.AutoFiniteDiff(; fdtype=fdtype)
+            end
+        end
+        return (ad, chunk_size, diff_type)
+    end
+end
+
 # Tell SciMLBase that EK1/DiagonalEK1 use ForwardDiff on the model function,
 # so that FunctionWrappersWrapper registers Dual-compatible wrappers.
 function SciMLBase.forwarddiffs_model(alg::AbstractEK)
@@ -13,19 +56,21 @@ function SciMLBase.forwarddiffs_model(alg::AbstractEK)
     return false
 end
 
-# Newer OrdinaryDiffEqCore dispatches _default_dae_init! by algorithm type, but only
-# extends it for its own implicit algorithm types in OrdinaryDiffEqNonlinearSolve.
-# Replicate the same behavior here so that DAE initialization (BrownFullBasicInit)
-# works for singular mass matrices, as it did before the refactor.
-function OrdinaryDiffEqCore._default_dae_init!(integrator, prob, x, alg::AbstractEK)
-    initializealg = DiffEqBase.BrownFullBasicInit(integrator.opts.abstol)
-    if applicable(OrdinaryDiffEqCore._initialize_dae!, integrator, prob, initializealg, x)
-        OrdinaryDiffEqCore._initialize_dae!(integrator, prob, initializealg, x)
-    else
-        error(
-            "`OrdinaryDiffEqNonlinearSolve` is not loaded, which is required for " *
-            "DAE initialization with singular mass matrices. " *
-            "To fix this, do `using OrdinaryDiffEqNonlinearSolve` or `using OrdinaryDiffEq`.",
+# DAE-initialization hook.
+#
+# v3's `_default_dae_init!` is only extended for OrdinaryDiffEq's own implicit
+# algorithms (in OrdinaryDiffEqNonlinearSolve), so without an override the EK
+# solvers hit a MethodError on DAE / singular-mass-matrix problems. Route through
+# `CheckInit` to match the OrdinaryDiffEq v7 default; users who want the previous
+# silent-fix behavior pass `initializealg = BrownFullBasicInit()` explicitly.
+# v4 dropped the hook and `CheckInit` is already its `DefaultInit` target.
+@static if isdefined(OrdinaryDiffEqCore, :_default_dae_init!)
+    function OrdinaryDiffEqCore._default_dae_init!(integrator, prob, x, alg::AbstractEK)
+        return OrdinaryDiffEqCore._initialize_dae!(
+            integrator,
+            prob,
+            DiffEqBase.CheckInit(),
+            x,
         )
     end
 end
@@ -226,7 +271,7 @@ struct EK1{CS,AD,DiffType,ST,CJ,PT,DT,IT,RT,CF} <: AbstractEK
     ) where {PT,DT,IT,RT,CF} = begin
         ekargcheck(EK1; diffusionmodel, pn_observation_noise, covariance_factorization)
         AD_choice, chunk_size, diff_type =
-            OrdinaryDiffEqCore._process_AD_choice(autodiff, chunk_size, diff_type)
+            _process_AD_choice(autodiff, chunk_size, diff_type)
         new{
             _unwrap_val(chunk_size),
             typeof(AD_choice),
@@ -323,7 +368,7 @@ struct DiagonalEK1{CS,AD,DiffType,ST,CJ,PT,DT,IT,RT,CF} <: AbstractEK
     ) where {PT,DT,IT,RT,CF} = begin
         ekargcheck(DiagonalEK1; diffusionmodel, pn_observation_noise, covariance_factorization)
         AD_choice, chunk_size, diff_type =
-            OrdinaryDiffEqCore._process_AD_choice(autodiff, chunk_size, diff_type)
+            _process_AD_choice(autodiff, chunk_size, diff_type)
         new{
             _unwrap_val(chunk_size),
             typeof(AD_choice),

src/perform_step.jl

@@ -1,3 +1,14 @@
+# On OrdinaryDiffEqCore v4+ the local error estimate `EEst` lives on `integrator.controller_cache`
+# and must be accessed through `get_EEst`/`set_EEst!`. On v3 it was a mutable field on the
+# integrator struct itself. These helpers paper over that difference.
+@static if isdefined(OrdinaryDiffEqCore, :set_EEst!)
+    _set_EEst!(integ, val) = OrdinaryDiffEqCore.set_EEst!(integ, val)
+    _get_EEst(integ) = OrdinaryDiffEqCore.get_EEst(integ)
+else
+    _set_EEst!(integ, val) = (integ.EEst = val)
+    _get_EEst(integ) = integ.EEst
+end
+
 # Called in the OrdinaryDiffEqCore.__init; All `OrdinaryDiffEqAlgorithm`s have one
 function OrdinaryDiffEqCore.initialize!(
     integ::OrdinaryDiffEqCore.ODEIntegrator,
@@ -100,8 +111,9 @@ function OrdinaryDiffEqCore.perform_step!(integ, cache::EKCache, repeat_step=fal
         cache.local_diffusion = estimate_local_diffusion(cache.diffusionmodel, integ)
     end
     if integ.opts.adaptive
-        integ.EEst = compute_scaled_error_estimate!(integ, cache)
-        if integ.EEst >= one(integ.EEst)
+        _set_EEst!(integ, compute_scaled_error_estimate!(integ, cache))
+        _EEst = _get_EEst(integ)
+        if _EEst >= one(_EEst)
             return
         end
     end

test/Project.toml

@@ -14,6 +14,7 @@ LinearRegression = "92481ed7-9fb7-40fd-80f2-46fd0f076581"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

test/diffeqdevtools.jl

@@ -91,7 +91,6 @@ end
         numruns=2,
         maxiters=Int(1e7),
         timeseries_errors=false,
-        verbose=false,
     )
     @test_nowarn WorkPrecisionSet(
         prob, abstols, reltols, setups;
@@ -101,6 +100,5 @@ end
         numruns=2,
         maxiters=Int(1e7),
         timeseries_errors=false,
-        verbose=false,
     )
 end

test/mass_matrix.jl

@@ -1,5 +1,6 @@
 using ProbNumDiffEq
 using OrdinaryDiffEq
+using OrdinaryDiffEqFIRK: RadauIIA5
 using LinearAlgebra
 using Test
 

test/solution.jl

@@ -26,7 +26,7 @@ using ODEProblemLibrary: prob_ode_lotkavolterra
         @testset "Alg $Alg" for Alg in (EK0, EK1)
             sol = solve(prob, Alg())
 
-            @test length(sol) > 2
+            @test length(sol.u) > 2
             @test length(sol.t) == length(sol.u)
             @test length(prob.u0) == length(sol.u[end])
 
@@ -93,7 +93,7 @@ using ODEProblemLibrary: prob_ode_lotkavolterra
                     @test samples isa Array
 
                     m, n, o = size(samples)
-                    @test m == length(sol)
+                    @test m == length(sol.u)
                     @test n == length(sol.u[1])
                     @test o == n_samples
 

test/stiff_problem.jl

@@ -1,6 +1,7 @@
 using ProbNumDiffEq
 using Test
 using OrdinaryDiffEq
+using OrdinaryDiffEqFIRK: RadauIIA5
 using ODEProblemLibrary: prob_ode_vanderpol_stiff
 
 prob = prob_ode_vanderpol_stiff