Context:
In my personal project I am working with running ensemble problems of SDEs. To improve performance, I used SVectors wherever possible. Depending on my needs, I used componentwise_vectors_timepoint to interpolate my solution to various timepoints. Notably, the t values I provided to the function were also SVectors, as this way the returned EnsembleSolutions u values are also SVectors, simplifying the rest of my code. Recently I decided to expand the scope to ODEs as well, and was surprised to find that the same process failed when the problem is instead an ODEProblem.
MWE:
using StaticArrays
using DifferentialEquations
using DifferentialEquations.EnsembleAnalysis
SDEprob = SDEProblem((u, p, t) -> 1.01u, (u, p, t) -> 1.01u, 0.5, (0.0, 1.0))
SDEensemble_prob = EnsembleProblem(SDEprob)
prob = ODEProblem((u, p, t) -> 1.01u, 0.5, (0.0, 1.0))
ensemble_prob = EnsembleProblem(prob)
SDEysim = solve(SDEensemble_prob, SOSRI(), trajectories = 1)
ysim = solve(ensemble_prob, Tsit5(), trajectories = 1)
t = SVector{11}(0:0.1:1) # t values we want to interpolate to, converted to SVector
componentwise_vectors_timepoint(SDEysim, t)[1] # Returns the interpolated u as an SVector
componentwise_vectors_timepoint(ysim, t)[1] # MethodError due to using setindex! on t
The SDE timepoint interpolation returns the following:
t: 11-element SVector{11, Float64} with indices SOneTo(11):
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
u: 11-element SVector{11, Float64} with indices SOneTo(11):
0.5
0.5546109594698383
0.8321574072368957
0.6965533677964221
0.6559293124311955
0.4941314584280177
0.4658099736486656
0.32344010904393883
0.28964405720750147
0.31472125794773376
0.2209161971944798
The error that results from the ODE interpolation is the following:
ERROR: LoadError: setindex!(::SVector{11, Float64}, value, ::Int) is not defined.
Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] setindex!(a::SVector{11, Float64}, value::Float64, i::Int64)
@ StaticArrays ~/.julia/packages/StaticArrays/9Yt0H/src/indexing.jl:3
[3] macro expansion
@ ~/.julia/packages/StaticArrays/9Yt0H/src/indexing.jl:199 [inlined]
[4] _setindex!
@ ~/.julia/packages/StaticArrays/9Yt0H/src/indexing.jl:192 [inlined]
[5] setindex!
@ ~/.julia/packages/StaticArrays/9Yt0H/src/indexing.jl:167 [inlined]
[6] invpermute!(v::SVector{11, Float64}, p::SVector{11, Int64})
@ Base ./combinatorics.jl:237
[7] ode_interpolation(tvals::SVector{…}, id::OrdinaryDiffEqCore.InterpolationData{…}, idxs::Nothing, deriv::Type{…}, p::SciMLBase.NullParameters, continuity::Symbol)
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/NnA60/src/dense/generic_dense.jl:574
[8] (::OrdinaryDiffEqCore.InterpolationData{…})(tvals::SVector{…}, idxs::Nothing, deriv::Type, p::SciMLBase.NullParameters, continuity::Symbol)
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/NnA60/src/interp_func.jl:46
[9] (::ODESolution{…})(t::SVector{…}, ::Type{…}, idxs::Nothing, continuity::Symbol)
@ SciMLBase ~/.julia/packages/SciMLBase/fqdZV/src/solutions/ode_solutions.jl:242
[10] (::ODESolution{…})(t::SVector{…}, ::Type{…}; idxs::Nothing, continuity::Symbol)
@ SciMLBase ~/.julia/packages/SciMLBase/fqdZV/src/solutions/ode_solutions.jl:224
[11] (::ODESolution{Float64, 1, Vector{…}, Nothing, Nothing, Vector{…}, Vector{…}, Nothing, ODEProblem{…}, Tsit5{…}, OrdinaryDiffEqCore.InterpolationData{…}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing})(t::SVector{11, Float64}, ::Type{Val{…}})
@ SciMLBase ~/.julia/packages/SciMLBase/fqdZV/src/solutions/ode_solutions.jl:219
[12] (::SciMLBase.EnsembleAnalysis.var"#3#4"{…})(sol::ODESolution{…})
@ SciMLBase.EnsembleAnalysis ./none:0
[13] iterate
@ ./generator.jl:48 [inlined]
[14] componentwise_vectors_timepoint(sim::EnsembleSolution{Float64, 2, Vector{ODESolution{…}}}, t::SVector{11, Float64})
@ SciMLBase.EnsembleAnalysis ~/.julia/packages/SciMLBase/fqdZV/src/ensemble/ensemble_analysis.jl:17
[15] top-level scope
@ ./toyensemble.jl:18
[16] include(fname::String)
@ Main ./sysimg.jl:38
[17] top-level scope
@ REPL[3]:1
in expression starting at ./toyensemble.jl:18
Some type information was truncated. Use `show(err)` to see complete types.
I realize that here t does not need to be an SVector, but the odd behavior makes me wonder if this is a bug.