Support in-place interpolation of symbolic idxs

[hersle]

An attempt at fixing SciML/OrdinaryDiffEq.jl#2562 and making sure the documented in-place interpolation API actually works consistently with MTK, too.

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

[ChrisRackauckas]

This isn't the best solution since it allocates. But one step at a time, at least making this work should come first.

[ChrisRackauckas]

I think the tests are showing this is leading some issues / ambiguities?

[hersle]

Yeah, I got started on this and left it hanging a little. Will see if I can get back on it in not too long.

[hersle]

Can you help me understand how the downstream testing works?

For example, IntegrationTest / OrdinaryDiffEq.jl/Interface/1 (pull_request) passes when I dev OrdinaryDiffEq SciMLBase locally (both at or ahead of the latest master branches):

julia> Pkg.status()
  ...
  [1dea7af3] OrdinaryDiffEq v6.95.0 `C:\Users\herma\.julia\dev\OrdinaryDiffEq`      
  [0bca4576] SciMLBase v2.85.0 `C:\Users\herma\.julia\dev\SciMLBase`
  ...

julia> include(raw"C:\Users\herma\.julia\dev\OrdinaryDiffEq\test\interface\inplace_interpolation.jl");
Test Summary:                      | Pass  Total  Time
ODESolution interpolation (normal) |    4      4  0.0s
  1D                               |    2      2  0.0s
  2D                               |    2      2  0.0s
Test Summary:                         | Pass  Total  Time
ODESolution interpolation (not dense) |    4      4  0.0s
  1D                                  |    2      2  0.0s
  2D                                  |    2      2  0.0s
Test Summary:                          | Pass  Total  Time
ODESolution interpolation (fixed-step) |    4      4  0.0s
  1D                                   |    2      2  0.0s
  2D                                   |    2      2  0.0s
Test Summary:                                  | Pass  Total  Time
ODESolution interpolation (CompositeAlgorithm) |    4      4  0.0s
  1D                                           |    2      2  0.0s
  2D                                           |    2      2  0.0s

But the Github run fails. I have not checked the other tests, though.

[ChrisRackauckas]

It grabs this package https://github.com/SciML/SciMLBase.jl/actions/runs/14602553998/job/41007048826?pr=988#step:6:177 devs the downstream https://github.com/SciML/SciMLBase.jl/actions/runs/14602553998/job/41007048826?pr=988#step:6:626 and then runs. Rebase, I did a lot of releases with minor test fixes today so maybe if it's just rerun on the latest master things will be better? If that doesn't work, make sure you dev'd OrdinaryDiffEqCore? And if you check that, I can dive in and see if I can figure out what's going on.

[hersle]

Okay, thank you. Does that mean it's equivalent to an environment where SciMLBase and OrdinaryDiffEq are both deved and on the tip of master? I just rebased onto the latest SciMLBase to be sure, and that test still passes locally on my computer. I don't think the precise state of OrdinaryDiffEq is relevant anyway, since I'm just running this test file.

Can you try to run tests once more, unless it sounds like I'm misunderstanding something? 😅

[hersle]

I was missing an in-place dispatch for t::Number, so more tests should pass now...

I'm still confused about this MethodError vs BoundsError, though. It might not be a problem with versioning, but simply that my computer and Github give different results. It passes locally with MethodError, e.g.:

julia> sol_ODE(out_VF, tt; idxs = 1:1)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  convert(::Type{T}, ::Union{Static.StaticBool{N}, Static.StaticFloat64{N}, Static.StaticInt{N}} where N) where T<:Number
   @ Static C:\Users\herma\.julia\packages\Static\SeEGr\src\Static.jl:427
  convert(::Type{T}, ::MultivariatePolynomials.AbstractPolynomialLike) where T<:Number
   @ MultivariatePolynomials C:\Users\herma\.julia\packages\MultivariatePolynomials\iTfZE\src\conversion.jl:96
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number
   @ Base multidimensional.jl:136
  ...

Stacktrace:
  [1] setindex!(A::Memory{Float64}, x::Vector{Float64}, i1::Int64)
    @ Base .\genericmemory.jl:243
  [2] unsafe_copyto!(dest::Memory{Float64}, doffs::Int64, src::Memory{Vector{Float64}}, soffs::Int64, n::Int64)
    @ Base .\genericmemory.jl:153
  [3] unsafe_copyto!
    @ .\genericmemory.jl:133 [inlined]
  [4] _copyto_impl!
    @ .\array.jl:308 [inlined]
  [5] copyto!
    @ .\array.jl:294 [inlined]
  [6] copyto!
    @ .\array.jl:319 [inlined]
  [7] copyto!
    @ .\broadcast.jl:966 [inlined]
  [8] copyto!
    @ .\broadcast.jl:925 [inlined]
  [9] materialize!(::Base.Broadcast.DefaultArrayStyle{…}, dest::Vector{…}, bc::Base.Broadcast.Broadcasted{…})
    @ Base.Broadcast .\broadcast.jl:883
 [10] materialize!
    @ .\broadcast.jl:880 [inlined]
 [11] (::ODESolution{…})(v::Vector{…}, t::StepRangeLen{…}, ::Type{…}, idxs::UnitRange{…}, continuity::Symbol)
    @ SciMLBase C:\Users\herma\.julia\dev\SciMLBase\src\solutions\ode_solutions.jl:409
 [12] #_#538
    @ C:\Users\herma\.julia\dev\SciMLBase\src\solutions\ode_solutions.jl:229 [inlined]
 [13] AbstractODESolution
    @ C:\Users\herma\.julia\dev\SciMLBase\src\solutions\ode_solutions.jl:224 [inlined]
 [14] top-level scope
    @ REPL[106]:1
Some type information was truncated. Use `show(err)` to see complete types.

But the Github run somehow gives BoundsError, which also sounds reasonable, since the test is testing an incompatible output array on purpose. Is @test_throws Union{MethodError, BoundsError} OK, or should we get to the bottom of this?

[ChrisRackauckas]

I'm still confused about this MethodError vs BoundsError, though. It might not be a problem with versioning, but simply that my computer and Github give different results. It passes locally with MethodError, e.g.:

That's because tests are run with --bounds-check=true, so the difference of force enabling bounds check is making that a bounds check error on the indexing of a scalar.

[ChrisRackauckas]

Check a=1.0; @inbounds a[1] in your REPL with --bounds-check=false/true and you'll see the difference.

[ChrisRackauckas]

Is @test_throws Union{MethodError, BoundsError} OK

yes

[ChrisRackauckas]

The issue is that the symbolic dispatch is kept open, since there's more than one issymbolic type. We should union that in the future, but it needs JuliaSymbolics/SymbolicUtils.jl#737

[ChrisRackauckas]

Suggested change

	v .= map(eachindex(t)) do ti
	sol.interp(u, t[ti], nothing, deriv, p, continuity)
	if eltype(v) <: Number
	v[ti] = sol.interp(t[ti], nothing, deriv, p, continuity)
	else
	for ti in eachindex(t)
	sol.interp(v[ti], t[ti], nothing, deriv, p, continuity)
	end
	end

[hersle]

That's because tests are run with --bounds-check=true, so the difference of force enabling bounds check is making that a bounds check error on the indexing of a scalar.

Oh, well that explains it. Thanks.

Let's see how this goes...