Hello,
I was running into an issue with indexing symbolic arrays from a solution and was wondering if anyone had some insight. Below is a simplified example:
using ModelingToolkit, StochasticDiffEq
using Plots
@variables t1 x(t1)[1:10]
@parameters μ σ
Dx = Differential(t1)
x_eq = Symbolics.scalarize(Dx.(x) .~ 1)
x_noise_eq = Symbolics.scalarize(x .* σ)
@named sde_sys = SDESystem(x_eq,
x_noise_eq, t1,
[x...],
[μ, σ]; tspan=(0, 10.0))
sde_prob = SDEProblem(sde_sys, Symbolics.scalarize(x .=> 1), (0.0, 10.0),
[μ => 1, σ => 0.2])
sde_sol = solve(sde_prob, SOSRI())
sde_sol(0; idxs=x[1:2]) # leads to an error
plot(sde_sol; idxs=x[1:2]) # does not lead to an error
Does anyone know why the plot line works but the direct indexing does not?
Thanks!
Edit:
This does the trick for retrieving the integer indices in the above and appears to also work with multi-D arrays on a separate example I was working on:
indices = []
for i in 1:2
push!(indices,SciMLBase.sym_to_index(x[i], SciMLBase.getsyms(sde_sol)))
end
sde_sol(0.5; idxs=indices)
It’s not an issue with x[1:2], you get the same issue without indexing:
julia> sde_sol(0; idxs=x)
ERROR: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
[1] checkbounds
@ .\abstractarray.jl:668 [inlined]
[2] view
@ .\subarray.jl:177 [inlined]
[3] maybeview
@ .\views.jl:146 [inlined]
[4] dotview
@ .\broadcast.jl:1201 [inlined]
[5] sde_interpolant(Θ::Float64, dt::Float64, u0::Vector{Float64}, u1::Vector{Float64}, idxs::Symbolics.Arr{Num, 1}, deriv::Type{Val{0}})
@ StochasticDiffEq C:\Users\User\.julia\packages\StochasticDiffEq\rZj4R\src\dense.jl:32
[6] sde_interpolation(tval::Int64, id::StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, idxs::Symbolics.Arr{Num, 1}, deriv::Type, p::Vector{Float64}, continuity::Symbol)
@ StochasticDiffEq C:\Users\User\.julia\packages\StochasticDiffEq\rZj4R\src\dense.jl:169
[7] LinearInterpolationData
@ C:\Users\User\.julia\packages\StochasticDiffEq\rZj4R\src\interp_func.jl:7 [inlined]
[8] #_#420
@ C:\Users\User\.julia\packages\SciMLBase\o0qGB\src\solutions\rode_solutions.jl:66 [inlined]
[9] top-level scope
@ Untitled-1:16
julia> plot(sde_sol; idxs=x)
When you put in this form of idxs into sde_sol(0; idxs = x), the underlying interpolant in sde_sol is trying to use it, but I guess it hasn’t been implemented to interpret it correctly (you can see how it interprets it by looking at the functions in the stacktrace).
It works for the plot call because it gets interpreted correctly in the call which you can see here SciMLBase.jl/solution_interface.jl at 9e5f84dd748f462266d52927f1b5dd7202e40064 · SciML/SciMLBase.jl · GitHub, e.g. with SciMLBase.interpret_vars(x[1:2], sde_sol, SciMLBase.getsyms(sde_sol)).
The issue causing the error is an expression like checkbounds(A, x) for an array A when doing A[x] in StochasticDiffEq.sde_interpolant(0.0,0.0,zeros(10),zeros(10),x,Val{0}) when evaluating the interpolant (the arrays here are just for example). e.g.
julia> @inbounds checkbounds([], x)
ERROR: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
[1] checkbounds(A::Vector{Any}, I::Symbolics.Arr{Num, 1})
@ Base .\abstractarray.jl:668
[2] top-level scope
@ REPL[82]:1
I guess there could be a better way to handle this array access by converting somewhere along the way.
Sorry. I should have been clearer. My issue is not with indexing as x[1:2] but generally trying to index a solution’s results by using the variable names and indices as opposed to figuring out the integer slices for each variable. The issue is one of both readability (since I think it is clearer to index by variable names) and for working with multi-D arrays where it can be a bit harder to tell how Julia flattens it out.
Thanks for pointing to SciMLBase docs! I think the sym_to_index does the trick for what I am looking for here
I think this is a case we haven’t handled well yet. In the near future it should be automatic but not all of the interpolation cases seem to be handled with SDEs.