From db7d707e032329ccbc6b3d2665b33e3a9642cf85 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 24 Apr 2026 01:36:15 -0400 Subject: [PATCH 1/3] Migrate to SciMLBase v3 / RecursiveArrayTools v4 / DiffEqBase v7 RecursiveArrayTools v4 makes `AbstractVectorOfArray <: AbstractArray`, which changes the semantics of several operations on `ODESolution`, `DAESolution`, and `EnsembleSolution`: - `length(sol)` is now `prod(size(sol))` (total scalar count) rather than the number of timesteps. Replaced with `length(sol.u)` in the timestep-counting loops in `L2Loss` and in the ensemble retcode check in `LogLikeLoss`. - Iterating `sol` on an `AbstractEnsembleSolution` now yields scalar elements, not trajectories. Replaced `for s in sol` with `for s in sol.u` in both `LogLikeLoss(::AbstractEnsembleSolution)` and the `AbstractSciMLSolution` fallback branch. - For the steady-state `AbstractNoTimeSolution` path, moved to iterating `sol.u` (the final-state vector) explicitly rather than via `length(sol)` / `sol[i]`, which are now scalar AbstractArray operations on v4. 2-argument indexing `sol[j, i]` is unchanged: on both RAT v3 and v4 it resolves to `sol.u[i][j]` (component j at timestep i), so the interior loops do not need to change. Widen `Project.toml` compat: - `DiffEqBase = "6, 7"` (was "6") - `SciMLBase = "1.69, 2, 3"` (was "1.69, 2, 3.1") - Bump package version to 2.5.0. Test-suite fix: `result.minimizer` on a SciMLBase `OptimizationSolution` was removed in v3; replaced with `result.u` in `regularization_test.jl`. Other `result.minimizer` references in the test suite come from `Optim.optimize(...)` (not `solve(optprob, ...)`), which still returns an `Optim.OptimizationResults` with `.minimizer` and is unaffected. Supersedes and extends Dependabot PR #292, which only widened the `DiffEqBase` compat without the accompanying code migration. Co-Authored-By: Claude Opus 4.7 (1M context) Co-Authored-By: Chris Rackauckas --- Project.toml | 6 ++-- src/cost_functions.jl | 39 +++++++++++++++-------- test/tests_on_odes/regularization_test.jl | 4 +-- 3 files changed, 31 insertions(+), 18 deletions(-) diff --git a/Project.toml b/Project.toml index c3d1e59..2cd9c41 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "DiffEqParamEstim" uuid = "1130ab10-4a5a-5621-a13d-e4788d82bd4c" -version = "2.4.0" +version = "2.5.0" authors = ["Chris Rackauckas "] [deps] @@ -21,13 +21,13 @@ StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" Calculus = "0.5" CommonSolve = "0.2.6" Dierckx = "0.4, 0.5" -DiffEqBase = "6" +DiffEqBase = "6, 7" Distributions = "0.25" ForwardDiff = "0.10" PenaltyFunctions = "0.1, 0.2, 0.3" PreallocationTools = "0.2, 0.3, 0.4, 1.0" RecursiveArrayTools = "1.0, 2.0, 3, 4" -SciMLBase = "1.69, 2, 3.1" +SciMLBase = "1.69, 2, 3" SciMLSensitivity = "7" Statistics = "1.10" StatsAPI = "1.8.0" diff --git a/src/cost_functions.jl b/src/cost_functions.jl index 8edfd88..9de49ee 100644 --- a/src/cost_functions.jl +++ b/src/cost_functions.jl @@ -48,16 +48,20 @@ function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution) sumsq = 0.0 + # `sol.u` is the final-state vector for a NoTimeSolution (steady-state / + # nonlinear). Iterating `sol.u` is the forward-compatible form under both + # RAT v3 and v4. + solu = sol.u if weight === nothing - @inbounds for i in 1:length(sol) - sumsq += (data[i] - sol[i])^2 + @inbounds for i in 1:length(solu) + sumsq += (data[i] - solu[i])^2 end else - @inbounds for i in 1:length(sol) + @inbounds for i in 1:length(solu) if weight isa Real - sumsq = sumsq + ((data[i] - sol[i])^2) * weight + sumsq = sumsq + ((data[i] - solu[i])^2) * weight else - sumsq = sumsq + ((data[i] - sol[i])^2) * weight[i] + sumsq = sumsq + ((data[i] - solu[i])^2) * weight[i] end end end @@ -80,8 +84,13 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) sumsq = 0.0 + # Under RAT v4, `AbstractVectorOfArray <: AbstractArray`, so `length(sol)` + # is the total scalar element count, not the number of timesteps. Iterate + # `sol.u` explicitly for the forward-compatible form. 2-argument indexing + # `sol[j, i]` keeps the same (component, time) meaning under both v3 and v4. + nsteps = length(sol.u) if weight === nothing - @inbounds for i in 1:length(sol) + @inbounds for i in 1:nsteps for j in 1:length(sol.u[i]) sumsq += (data[j, i] - sol[j, i])^2 end @@ -108,7 +117,7 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) end end else - @inbounds for i in 1:length(sol) + @inbounds for i in 1:nsteps if weight isa Real for j in 1:length(sol.u[i]) sumsq = sumsq + ((data[j, i] - sol[j, i])^2) * weight @@ -195,7 +204,8 @@ end function (f::LogLikeLoss)(sol::SciMLBase.AbstractSciMLSolution) distributions = f.data_distributions if sol isa DiffEqBase.AbstractEnsembleSolution - failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol) + # `sol.u` is the vector of trajectories under both RAT v3 and v4. + failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) end @@ -244,7 +254,10 @@ end function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution) distributions = f.data_distributions - failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol) + # Under RAT v4, `AbstractEnsembleSolution <: AbstractArray` and iterating + # `sol` yields scalar elements, not trajectories. Iterate `sol.u` for + # trajectory-wise access in a way that is forward-compatible with v3/v4. + failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) failure && return Inf ll = 0.0 if eltype(distributions) <: UnivariateDistribution @@ -252,7 +265,7 @@ function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution) # i is the number of time points # j is the size of the system # corresponds to distributions[i,j] - vals = [s[i, j] for s in sol] + vals = [s[i, j] for s in sol.u] ll -= loglikelihood(distributions[i, j], vals) end else @@ -260,7 +273,7 @@ function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution) # i is the number of time points # j is the size of the system # corresponds to distributions[i,j] - vals = [s[i, j] for i in 1:length(sol.u[1].u[1]), s in sol] + vals = [s[i, j] for i in 1:length(sol.u[1].u[1]), s in sol.u] ll -= loglikelihood(distributions[j], vals) end end @@ -271,12 +284,12 @@ function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution) if eltype(distributions) <: UnivariateDistribution for j in 2:length(f.t), i in 1:length(sol.u[1].u[1]) - vals = [s[i, j] - s[i, j - 1] for s in sol] + vals = [s[i, j] - s[i, j - 1] for s in sol.u] fdll -= logpdf(distributions[j - 1, i], vals)[1] end else for j in 2:length(f.t) - vals = [s[i, j] - s[i, j - 1] for i in 1:length(sol.u[1].u[1]), s in sol] + vals = [s[i, j] - s[i, j - 1] for i in 1:length(sol.u[1].u[1]), s in sol.u] fdll -= logpdf(distributions[j - 1], vals)[1] end end diff --git a/test/tests_on_odes/regularization_test.jl b/test/tests_on_odes/regularization_test.jl index 54f2e4b..039672f 100644 --- a/test/tests_on_odes/regularization_test.jl +++ b/test/tests_on_odes/regularization_test.jl @@ -36,8 +36,8 @@ result = solve(optprob, Optim.BFGS()) optprob = Optimization.OptimizationProblem(cost_function_2, [1.2, 2.7]) result = solve(optprob, Optim.BFGS()) -@test result.minimizer ≈ [1.5; 3.0] atol = 3.0e-1 +@test result.u ≈ [1.5; 3.0] atol = 3.0e-1 optprob = Optimization.OptimizationProblem(cost_function_3, [1.3, 0.8, 2.8, 1.2]) result = solve(optprob, Optim.BFGS()) -@test result.minimizer ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 +@test result.u ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 From ac0c3cf3eea70b53865df18b0bdf67645cd5c7ad Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 24 Apr 2026 01:36:59 -0400 Subject: [PATCH 2/3] Fix remaining ensemble iteration sites in L2Loss Two `L2Loss` branches (the `AbstractNoTimeSolution` and `AbstractSciMLSolution` methods) still iterated `for s in sol` in the `AbstractEnsembleSolution` retcode-check. Under RAT v4 that iterates scalar elements, not trajectories. Switched both to `for s in sol.u`, matching the fix already applied in `LogLikeLoss`. Co-Authored-By: Claude Opus 4.7 (1M context) Co-Authored-By: Chris Rackauckas --- src/cost_functions.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/cost_functions.jl b/src/cost_functions.jl index 9de49ee..a334388 100644 --- a/src/cost_functions.jl +++ b/src/cost_functions.jl @@ -40,7 +40,8 @@ function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution) dudt = f.dudt if sol isa DiffEqBase.AbstractEnsembleSolution - failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol) + # `sol.u` is the vector of trajectories under both RAT v3 and v4. + failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) end @@ -76,7 +77,8 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) dudt = f.dudt if sol isa DiffEqBase.AbstractEnsembleSolution - failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol) + # `sol.u` is the vector of trajectories under both RAT v3 and v4. + failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) end From ec5c26ecec3a70125bdc2b0030c4390751450543 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 24 Apr 2026 05:59:02 +0000 Subject: [PATCH 3/3] Apply suggestions from code review Co-authored-by: Christopher Rackauckas --- src/cost_functions.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/cost_functions.jl b/src/cost_functions.jl index a334388..e329893 100644 --- a/src/cost_functions.jl +++ b/src/cost_functions.jl @@ -40,7 +40,6 @@ function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution) dudt = f.dudt if sol isa DiffEqBase.AbstractEnsembleSolution - # `sol.u` is the vector of trajectories under both RAT v3 and v4. failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) @@ -49,9 +48,6 @@ function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution) sumsq = 0.0 - # `sol.u` is the final-state vector for a NoTimeSolution (steady-state / - # nonlinear). Iterating `sol.u` is the forward-compatible form under both - # RAT v3 and v4. solu = sol.u if weight === nothing @inbounds for i in 1:length(solu) @@ -77,7 +73,6 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) dudt = f.dudt if sol isa DiffEqBase.AbstractEnsembleSolution - # `sol.u` is the vector of trajectories under both RAT v3 and v4. failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) @@ -86,10 +81,6 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) sumsq = 0.0 - # Under RAT v4, `AbstractVectorOfArray <: AbstractArray`, so `length(sol)` - # is the total scalar element count, not the number of timesteps. Iterate - # `sol.u` explicitly for the forward-compatible form. 2-argument indexing - # `sol[j, i]` keeps the same (component, time) meaning under both v3 and v4. nsteps = length(sol.u) if weight === nothing @inbounds for i in 1:nsteps @@ -206,7 +197,6 @@ end function (f::LogLikeLoss)(sol::SciMLBase.AbstractSciMLSolution) distributions = f.data_distributions if sol isa DiffEqBase.AbstractEnsembleSolution - # `sol.u` is the vector of trajectories under both RAT v3 and v4. failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) else failure = !SciMLBase.successful_retcode(sol.retcode) @@ -256,9 +246,6 @@ end function (f::LogLikeLoss)(sol::DiffEqBase.AbstractEnsembleSolution) distributions = f.data_distributions - # Under RAT v4, `AbstractEnsembleSolution <: AbstractArray` and iterating - # `sol` yields scalar elements, not trajectories. Iterate `sol.u` for - # trajectory-wise access in a way that is forward-compatible with v3/v4. failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol.u) failure && return Inf ll = 0.0