Update to OrdinaryDiffEq split version

Project.toml

@@ -17,7 +17,10 @@ Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
 Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 PSDMatrices = "fe68d972-6fd8-4755-bdf0-97d4c54cefdc"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -60,7 +63,10 @@ Kronecker = "0.5.4"
 LinearAlgebra = "1"
 MatrixEquations = "2"
 Octavian = "0.3.17"
-OrdinaryDiffEq = "6.52"
+OrdinaryDiffEqCore = "1.3"
+OrdinaryDiffEqDifferentiation = "1.1"
+OrdinaryDiffEqRosenbrock = "1.1"
+OrdinaryDiffEqVerner = "1.1"
 PSDMatrices = "0.5.0"
 PrecompileTools = "1"
 Printf = "1"
@@ -77,7 +83,7 @@ TaylorIntegration = "0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 TaylorSeries = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18"
 Test = "1"
 ToeplitzMatrices = "0.7, 0.8"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/ProbNumDiffEq.jl

@@ -17,7 +17,8 @@ using Reexport
 @reexport using DiffEqBase
 import SciMLBase
 import SciMLBase: interpret_vars, getsyms, remake
-using OrdinaryDiffEq
+using OrdinaryDiffEqCore, OrdinaryDiffEqDifferentiation, OrdinaryDiffEqVerner,
+    OrdinaryDiffEqRosenbrock
 using ToeplitzMatrices
 using FastBroadcast
 using StaticArrayInterface
@@ -67,7 +68,7 @@ include("blocksofdiagonals.jl")
 include("covariance_structure.jl")
 export IsometricKroneckerCovariance, DenseCovariance, BlockDiagonalCovariance
 
-abstract type AbstractODEFilterCache <: OrdinaryDiffEq.OrdinaryDiffEqCache end
+abstract type AbstractODEFilterCache <: OrdinaryDiffEqCore.OrdinaryDiffEqCache end
 
 include("gaussians.jl")
 export Gaussian

src/alg_utils.jl

@@ -1,47 +1,47 @@
 ############################################################################################
-# For the equivalent parts in OrdinaryDiffEq.jl, see:
-# https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/alg_utils.jl
+# For the equivalent parts in OrdinaryDiffEqCore.jl, see:
+# https://github.com/SciML/OrdinaryDiffEqCore.jl/blob/master/src/alg_utils.jl
 ############################################################################################
 
-OrdinaryDiffEq._alg_autodiff(::AbstractEK) = Val{true}()
-OrdinaryDiffEq.standardtag(::AbstractEK) = false
-OrdinaryDiffEq.concrete_jac(::AbstractEK) = nothing
+OrdinaryDiffEqDifferentiation._alg_autodiff(::AbstractEK) = Val{true}()
+OrdinaryDiffEqDifferentiation.standardtag(::AbstractEK) = false
+OrdinaryDiffEqDifferentiation.concrete_jac(::AbstractEK) = nothing
 
 @inline DiffEqBase.get_tmp_cache(integ, alg::AbstractEK, cache::AbstractODEFilterCache) =
     (cache.tmp, cache.atmp)
-OrdinaryDiffEq.isfsal(::AbstractEK) = false
+OrdinaryDiffEqCore.isfsal(::AbstractEK) = false
 
 for ALG in [:EK1, :DiagonalEK1]
-    @eval OrdinaryDiffEq._alg_autodiff(::$ALG{CS,AD}) where {CS,AD} = Val{AD}()
-    @eval OrdinaryDiffEq.alg_difftype(::$ALG{CS,AD,DiffType}) where {CS,AD,DiffType} =
+    @eval OrdinaryDiffEqDifferentiation._alg_autodiff(::$ALG{CS,AD}) where {CS,AD} = Val{AD}()
+    @eval OrdinaryDiffEqDifferentiation.alg_difftype(::$ALG{CS,AD,DiffType}) where {CS,AD,DiffType} =
         DiffType
-    @eval OrdinaryDiffEq.standardtag(::$ALG{CS,AD,DiffType,ST}) where {CS,AD,DiffType,ST} =
+    @eval OrdinaryDiffEqDifferentiation.standardtag(::$ALG{CS,AD,DiffType,ST}) where {CS,AD,DiffType,ST} =
         ST
-    @eval OrdinaryDiffEq.concrete_jac(
+    @eval OrdinaryDiffEqDifferentiation.concrete_jac(
         ::$ALG{CS,AD,DiffType,ST,CJ},
     ) where {CS,AD,DiffType,ST,CJ} = CJ
-    @eval OrdinaryDiffEq.get_chunksize(::$ALG{CS}) where {CS} = Val(CS)
-    @eval OrdinaryDiffEq.isimplicit(::$ALG) = true
+    @eval OrdinaryDiffEqDifferentiation.get_chunksize(::$ALG{CS}) where {CS} = Val(CS)
+    @eval OrdinaryDiffEqCore.isimplicit(::$ALG) = true
 end
 
 ############################################
 # Step size control
-OrdinaryDiffEq.isadaptive(::AbstractEK) = true
-OrdinaryDiffEq.alg_order(alg::AbstractEK) = num_derivatives(alg.prior)
-# OrdinaryDiffEq.alg_adaptive_order(alg::AbstractEK) =
+OrdinaryDiffEqCore.isadaptive(::AbstractEK) = true
+OrdinaryDiffEqCore.alg_order(alg::AbstractEK) = num_derivatives(alg.prior)
+# OrdinaryDiffEqCore.alg_adaptive_order(alg::AbstractEK) =
 
 # PI control is the default!
-OrdinaryDiffEq.isstandard(::AbstractEK) = false # proportional
-OrdinaryDiffEq.ispredictive(::AbstractEK) = false # not sure, maybe Gustafsson acceleration?
+OrdinaryDiffEqCore.isstandard(::AbstractEK) = false # proportional
+OrdinaryDiffEqCore.ispredictive(::AbstractEK) = false # not sure, maybe Gustafsson acceleration?
 
-# OrdinaryDiffEq.qmin_default(alg::AbstractEK) =
-# OrdinaryDiffEq.qmax_default(alg::AbstractEK) =
-# OrdinaryDiffEq.beta2_default(alg::AbstractEK) = 2 // (5(OrdinaryDiffEq.alg_order(alg) + 1))
-# OrdinaryDiffEq.beta1_default(alg::AbstractEK, beta2) = 7 // (10(OrdinaryDiffEq.alg_order(alg) + 1))
-# OrdinaryDiffEq.gamma_default(alg::AbstractEK) =
+# OrdinaryDiffEqCore.qmin_default(alg::AbstractEK) =
+# OrdinaryDiffEqCore.qmax_default(alg::AbstractEK) =
+# OrdinaryDiffEqCore.beta2_default(alg::AbstractEK) = 2 // (5(OrdinaryDiffEqCore.alg_order(alg) + 1))
+# OrdinaryDiffEqCore.beta1_default(alg::AbstractEK, beta2) = 7 // (10(OrdinaryDiffEqCore.alg_order(alg) + 1))
+# OrdinaryDiffEqCore.gamma_default(alg::AbstractEK) =
 
-# OrdinaryDiffEq.uses_uprev(alg::, adaptive::Bool) = adaptive
-OrdinaryDiffEq.is_mass_matrix_alg(::AbstractEK) = true
+# OrdinaryDiffEqCore.uses_uprev(alg::, adaptive::Bool) = adaptive
+OrdinaryDiffEqCore.is_mass_matrix_alg(::AbstractEK) = true
 
 SciMLBase.isautodifferentiable(::AbstractEK) = true
 SciMLBase.allows_arbitrary_number_types(::AbstractEK) = true

src/algorithms.jl

@@ -1,7 +1,7 @@
 ########################################################################################
 # Algorithm
 ########################################################################################
-abstract type AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm end
+abstract type AbstractEK <: OrdinaryDiffEqCore.OrdinaryDiffEqAdaptiveAlgorithm end
 
 function ekargcheck(
     alg;
@@ -352,9 +352,9 @@ function DiffEqBase.prepare_alg(
     p,
     prob,
 ) where {T}
-    # See OrdinaryDiffEq.jl: ./src/alg_utils.jl (where this is copied from).
+    # See OrdinaryDiffEqCore.jl: ./src/alg_utils.jl (where this is copied from).
     # In the future we might want to make EK1 an OrdinaryDiffEqAdaptiveImplicitAlgorithm and
-    # use the prepare_alg from OrdinaryDiffEq; but right now, we do not use `linsolve` which
+    # use the prepare_alg from OrdinaryDiffEqCore; but right now, we do not use `linsolve` which
     # is a requirement.
 
     if (isbitstype(T) && sizeof(T) > 24) || (

src/caches.jl

@@ -72,7 +72,9 @@ mutable struct EKCache{
     jac_config::JC
 end
 
-function OrdinaryDiffEq.alg_cache(
+OrdinaryDiffEqCore.get_fsalfirstlast(cache::AbstractODEFilterCache, rate_prototype) = (nothing, nothing)
+
+function OrdinaryDiffEqCore.alg_cache(
     alg::AbstractEK,
     u,
     rate_prototype,
@@ -221,8 +223,8 @@ function OrdinaryDiffEq.alg_cache(
     du1 = similar(rate_prototype)
     dw1 = zero(u)
     atmp = similar(u, uEltypeNoUnits)
-    if OrdinaryDiffEq.isimplicit(alg)
-        jac_config = OrdinaryDiffEq.build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
+    if OrdinaryDiffEqCore.isimplicit(alg)
+        jac_config = OrdinaryDiffEqDifferentiation.build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)
     else
         jac_config = nothing
     end
@@ -251,5 +253,5 @@ function OrdinaryDiffEq.alg_cache(
     )
 end
 
-get_uf(f, t, p, ::Val{true}) = OrdinaryDiffEq.UJacobianWrapper(f, t, p)
-get_uf(f, t, p, ::Val{false}) = OrdinaryDiffEq.UDerivativeWrapper(f, t, p)
+get_uf(f, t, p, ::Val{true}) = OrdinaryDiffEqDifferentiation.UJacobianWrapper(f, t, p)
+get_uf(f, t, p, ::Val{false}) = OrdinaryDiffEqDifferentiation.UDerivativeWrapper(f, t, p)

src/derivative_utils.jl

@@ -7,12 +7,12 @@ function calc_H!(H, integ, cache)
     elseif integ.alg isa EK1
         calc_H_EK0!(H, integ, cache)
         # @assert integ.u == @view x_pred.μ[1:(q+1):end]
-        OrdinaryDiffEq.calc_J!(ddu, integ, cache, true)
+        OrdinaryDiffEqDifferentiation.calc_J!(ddu, integ, cache, true)
         _ddu = size(ddu, 2) != d ? view(ddu, 1:d, :) : ddu
         _matmul!(H, _ddu, cache.SolProj, -1.0, 1.0)
     elseif integ.alg isa DiagonalEK1
         calc_H_EK0!(H, integ, cache)
-        OrdinaryDiffEq.calc_J!(ddu, integ, cache, true)
+        OrdinaryDiffEqDifferentiation.calc_J!(ddu, integ, cache, true)
         _ddu = size(ddu, 2) != d ? view(ddu, 1:d, :) : ddu
         ddu_diag = if size(ddu, 2) == d
             # the normal case: just extract the diagonal

src/initialization/classicsolverinit.jl

@@ -47,7 +47,7 @@ function initial_update!(integ, cache, ::ClassicSolverInit)
     # Compute the other parts with classic solvers
     t0 = integ.sol.prob.tspan[1]
     dt =
-        10 * OrdinaryDiffEq.ode_determine_initdt(
+        10 * OrdinaryDiffEqCore.ode_determine_initdt(
             u,
             t,
             1,

src/initialization/common.jl

@@ -65,7 +65,7 @@ ForwardDiffInit() = begin
 end
 
 """
-    ClassicSolverInit(; alg=OrdinaryDiffEq.Tsit5(), init_on_ddu=false)
+    ClassicSolverInit(; alg=OrdinaryDiffEqCore.Tsit5(), init_on_ddu=false)
 
 Initialization via regression on a few steps of a classic ODE solver.
 
@@ -80,7 +80,7 @@ optionally the second derivative can also be set via automatic differentiation b
 `init_on_ddu=true`.
 
 # Arguments
-- `alg`: The solver to be used. Can be any solver from OrdinaryDiffEq.jl.
+- `alg`: The solver to be used. Can be any solver from OrdinaryDiffEqCore.jl.
 - `init_on_ddu`: If `true`, the second derivative is also initialized exactly via
   automatic differentiation with ForwardDiff.jl.
 

src/integrator_utils.jl

@@ -1,14 +1,14 @@
 """
-    OrdinaryDiffEq.postamble!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})
+    OrdinaryDiffEqCore.postamble!(integ::OrdinaryDiffEqCore.ODEIntegrator{<:AbstractEK})
 
-ProbNumDiffEq.jl-specific implementation of OrdinaryDiffEq.jl's `postamble!`.
+ProbNumDiffEq.jl-specific implementation of OrdinaryDiffEqCore.jl's `postamble!`.
 
-In addition to calling `OrdinaryDiffEq._postamble!(integ)`, calibrate the diffusion and
+In addition to calling `OrdinaryDiffEqCore._postamble!(integ)`, calibrate the diffusion and
 smooth the solution.
 """
-function OrdinaryDiffEq.postamble!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})
-    # OrdinaryDiffEq.jl-related calls:
-    OrdinaryDiffEq._postamble!(integ)
+function OrdinaryDiffEqCore.postamble!(integ::OrdinaryDiffEqCore.ODEIntegrator{<:AbstractEK})
+    # OrdinaryDiffEqCore.jl-related calls:
+    OrdinaryDiffEqCore._postamble!(integ)
     copyat_or_push!(integ.sol.k, integ.saveiter_dense, integ.k)
     pn_solution_endpoint_match_cur_integrator!(integ)
 
@@ -118,56 +118,56 @@ function smooth_solution!(integ)
     return nothing
 end
 
-"Inspired by `OrdinaryDiffEq.solution_match_cur_integrator!`"
+"Inspired by `OrdinaryDiffEqCore.solution_match_cur_integrator!`"
 function pn_solution_endpoint_match_cur_integrator!(integ)
     if integ.opts.save_end
         if integ.alg.smooth
-            OrdinaryDiffEq.copyat_or_push!(
+            OrdinaryDiffEqCore.copyat_or_push!(
                 integ.sol.x_filt,
                 integ.saveiter_dense,
                 integ.cache.x,
             )
         end
 
-        OrdinaryDiffEq.copyat_or_push!(
+        OrdinaryDiffEqCore.copyat_or_push!(
             integ.sol.pu,
             integ.saveiter,
             _gaussian_mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x),
         )
     end
 end
 
-"Extends `OrdinaryDiffEq._savevalues!` to save ProbNumDiffEq.jl-specific things."
+"Extends `OrdinaryDiffEqCore._savevalues!` to save ProbNumDiffEq.jl-specific things."
 function DiffEqBase.savevalues!(
-    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
+    integ::OrdinaryDiffEqCore.ODEIntegrator{<:AbstractEK},
     force_save=false,
     reduce_size=true,
 )
 
-    # Do whatever OrdinaryDiffEq would do
-    out = OrdinaryDiffEq._savevalues!(integ, force_save, reduce_size)
+    # Do whatever OrdinaryDiffEqCore would do
+    out = OrdinaryDiffEqCore._savevalues!(integ, force_save, reduce_size)
 
     # Save our custom stuff that we need for the posterior
     if integ.opts.save_everystep
         i = integ.saveiter
-        OrdinaryDiffEq.copyat_or_push!(integ.sol.diffusions, i, integ.cache.local_diffusion)
-        OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, i, integ.cache.x)
+        OrdinaryDiffEqCore.copyat_or_push!(integ.sol.diffusions, i, integ.cache.local_diffusion)
+        OrdinaryDiffEqCore.copyat_or_push!(integ.sol.x_filt, i, integ.cache.x)
         _gaussian_mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x)
-        OrdinaryDiffEq.copyat_or_push!(integ.sol.pu, i, integ.cache.pu_tmp)
+        OrdinaryDiffEqCore.copyat_or_push!(integ.sol.pu, i, integ.cache.pu_tmp)
 
         if integ.alg.smooth
-            OrdinaryDiffEq.copyat_or_push!(
+            OrdinaryDiffEqCore.copyat_or_push!(
                 integ.sol.backward_kernels, i, integ.cache.backward_kernel)
         end
     end
 
     return out
 end
 
-function OrdinaryDiffEq.update_uprev!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})
-    @assert !OrdinaryDiffEq.alg_extrapolates(integ.alg)
+function OrdinaryDiffEqCore.update_uprev!(integ::OrdinaryDiffEqCore.ODEIntegrator{<:AbstractEK})
+    @assert !OrdinaryDiffEqCore.alg_extrapolates(integ.alg)
     @assert isinplace(integ.sol.prob)
-    @assert !(integ.alg isa OrdinaryDiffEq.DAEAlgorithm)
+    @assert !(integ.alg isa OrdinaryDiffEqCore.DAEAlgorithm)
 
     recursivecopy!(integ.uprev, integ.u)
     recursivecopy!(integ.cache.xprev, integ.cache.x)

src/perform_step.jl

@@ -1,5 +1,5 @@
-# Called in the OrdinaryDiffEQ.__init; All `OrdinaryDiffEqAlgorithm`s have one
-function OrdinaryDiffEq.initialize!(integ::OrdinaryDiffEq.ODEIntegrator, cache::EKCache)
+# Called in the OrdinaryDiffEqCore.__init; All `OrdinaryDiffEqAlgorithm`s have one
+function OrdinaryDiffEqCore.initialize!(integ::OrdinaryDiffEqCore.ODEIntegrator, cache::EKCache)
     check_secondorderode(integ)
     check_densesmooth(integ)
     check_saveiter(integ)
@@ -13,15 +13,15 @@ function OrdinaryDiffEq.initialize!(integ::OrdinaryDiffEq.ODEIntegrator, cache::
     copy!(integ.cache.xprev, integ.cache.x)
 
     # These are necessary since the solution object is not 100% initialized by default
-    OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, integ.saveiter, cache.x)
+    OrdinaryDiffEqCore.copyat_or_push!(integ.sol.x_filt, integ.saveiter, cache.x)
     initial_pu = _gaussian_mul!(cache.pu_tmp, cache.SolProj, cache.x)
-    OrdinaryDiffEq.copyat_or_push!(integ.sol.pu, integ.saveiter, initial_pu)
+    OrdinaryDiffEqCore.copyat_or_push!(integ.sol.pu, integ.saveiter, initial_pu)
 
     return nothing
 end
 
 function make_new_transitions(integ, cache, repeat_step)::Bool
-    # Similar to OrdinaryDiffEq.do_newJ
+    # Similar to OrdinaryDiffEqCore.do_newJ
     if integ.iter <= 1
         return true
     elseif cache.prior isa IOUP && cache.prior.update_rate_parameter
@@ -59,10 +59,10 @@ Basically consists of the following steps
 - Kalman update step
 - (optional) Update the global diffusion MLE
 
-As in OrdinaryDiffEq.jl, this step is not necessarily successful!
-For that functionality, use `OrdinaryDiffEq.step!(integ)`.
+As in OrdinaryDiffEqCore.jl, this step is not necessarily successful!
+For that functionality, use `OrdinaryDiffEqCore.step!(integ)`.
 """
-function OrdinaryDiffEq.perform_step!(integ, cache::EKCache, repeat_step=false)
+function OrdinaryDiffEqCore.perform_step!(integ, cache::EKCache, repeat_step=false)
     @unpack t, dt = integ
     @unpack d = integ.cache
     @unpack xprev, x_pred, u_pred, x_filt, err_tmp = integ.cache
@@ -73,7 +73,7 @@ function OrdinaryDiffEq.perform_step!(integ, cache::EKCache, repeat_step=false)
     if make_new_transitions(integ, cache, repeat_step)
         # Rosenbrock-style update of the IOUP rate parameter
         if cache.prior isa IOUP && cache.prior.update_rate_parameter
-            OrdinaryDiffEq.calc_J!(cache.prior.rate_parameter, integ, cache, false)
+            OrdinaryDiffEqDifferentiation.calc_J!(cache.prior.rate_parameter, integ, cache, false)
         end
 
         make_transition_matrices!(cache, cache.prior, dt)
@@ -142,7 +142,7 @@ In addition, compute the measurement mean (`z`) and the measurement function Jac
 Results are saved into `integ.cache.du`, `integ.cache.ddu`, `integ.cache.measurement.μ`
 and `integ.cache.H`.
 Jacobians are computed either with the supplied `f.jac`, or via automatic differentiation,
-as in OrdinaryDiffEq.jl.
+as in OrdinaryDiffEqCore.jl.
 """
 function evaluate_ode!(integ, x_pred, t)
     @unpack f, p, dt = integ

src/precompile.jl

@@ -1,5 +1,5 @@
 # Inspired by:
-# https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.21.0/src/OrdinaryDiffEq.jl#L195-L221
+# https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.21.0/src/OrdinaryDiffEqCore.jl#L195-L221
 
 import PrecompileTools
 

src/solution.jl

@@ -74,7 +74,7 @@ function DiffEqBase.solution_new_retcode(sol::ProbODESolution{T,N}, retcode) whe
     )
 end
 
-# Used to build the initial empty solution in OrdinaryDiffEq.__init
+# Used to build the initial empty solution in OrdinaryDiffEqCore.__init
 function DiffEqBase.build_solution(
     prob::DiffEqBase.AbstractODEProblem,
     alg::AbstractEK,
@@ -87,7 +87,7 @@ function DiffEqBase.build_solution(
     kwargs...,
 )
     # By making an actual cache, interpolation can be written very closely to the solver
-    cache = OrdinaryDiffEq.alg_cache(
+    cache = OrdinaryDiffEqCore.alg_cache(
         alg,
         prob.u0,
         recursivecopy(prob.u0),