diff --git a/Project.toml b/Project.toml index f54b2f60d..3ec471bb2 100644 --- a/Project.toml +++ b/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" diff --git a/src/ProbNumDiffEq.jl b/src/ProbNumDiffEq.jl index 1bd5cce8f..ed144c39c 100644 --- a/src/ProbNumDiffEq.jl +++ b/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 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4fbbddd1b..3094060e4 100644 --- a/src/alg_utils.jl +++ b/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 diff --git a/src/algorithms.jl b/src/algorithms.jl index ba535e334..479b26148 100644 --- a/src/algorithms.jl +++ b/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) || ( diff --git a/src/caches.jl b/src/caches.jl index 8cd96353f..5bb019bf4 100644 --- a/src/caches.jl +++ b/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) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index fe6136522..c13d2ad64 100644 --- a/src/derivative_utils.jl +++ b/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 diff --git a/src/initialization/classicsolverinit.jl b/src/initialization/classicsolverinit.jl index 123ddadee..dfc1dbb8a 100644 --- a/src/initialization/classicsolverinit.jl +++ b/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, diff --git a/src/initialization/common.jl b/src/initialization/common.jl index 0d6909053..08a382c4c 100644 --- a/src/initialization/common.jl +++ b/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. diff --git a/src/integrator_utils.jl b/src/integrator_utils.jl index b27e0b60c..0e0e4cfc8 100644 --- a/src/integrator_utils.jl +++ b/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,18 +118,18 @@ 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), @@ -137,26 +137,26 @@ function pn_solution_endpoint_match_cur_integrator!(integ) 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 @@ -164,10 +164,10 @@ function DiffEqBase.savevalues!( 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) diff --git a/src/perform_step.jl b/src/perform_step.jl index f27cc6de3..3a791b585 100644 --- a/src/perform_step.jl +++ b/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 diff --git a/src/precompile.jl b/src/precompile.jl index ab48649dc..1a1fea9e4 100644 --- a/src/precompile.jl +++ b/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 diff --git a/src/solution.jl b/src/solution.jl index c85792666..d4562a52f 100644 --- a/src/solution.jl +++ b/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),