`save_idxs` with a scalar in a lazy interpolation algorithm

[SebastianM-C]

Using a scalar save_idxs with a lazy interpolation algorithm breaks the plot recipe for the solution object.

using OrdinaryDiffEq
using DiffEqCallbacks
using Plots

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

p = (10, 28, 8/3)
u0 = [0,10.,0]

tspan = (0., T+Ttr)
prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob, Vern9(), save_idxs=1)
plot(sol)

gives:

MethodError: no method matching eachindex(::Float64)
Closest candidates are:
  eachindex(::Cmd) at C:\Users\sebastian\.julia\v0.6\Compat\src\Compat.jl:409
  eachindex(::IndexLinear, ::AbstractArray) at abstractarray.jl:820
  eachindex(::IndexLinear, ::AbstractArray, ::AbstractArray...) at abstractarray.jl:822
  ...

Stacktrace:
 [1] ode_addsteps! at C:\Users\sebastian\.julia\v0.6\OrdinaryDiffEq\src\dense\verner_addsteps.jl:749 [inlined]
 [2] ode_addsteps! at C:\Users\sebastian\.julia\v0.6\OrdinaryDiffEq\src\dense\verner_addsteps.jl:665 [inlined] (repeats 3 times)
 [3] ode_interpolation(::Array{Float64,1}, ::Function, ::Void, ::Type{T} where T, ::Tuple{Int64,Int64,Float64}) at C:\Users\sebastian\.julia\v0.6\OrdinaryDiffEq\src\dense\generic_dense.jl:147
 [4] macro expansion at C:\Users\sebastian\.julia\v0.6\DiffEqBase\src\solutions\solution_interface.jl:99 [inlined]
 [5] apply_recipe(::Dict{Symbol,Any}, ::DiffEqBase.ODESolution{Float64,2,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,Tuple{Int64,Int64,Float64},#lorenz,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Vern9,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}}) at C:\Users\sebastian\.julia\v0.6\RecipesBase\src\RecipesBase.jl:291
 [6] _process_userrecipes(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Float64,2,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,Tuple{Int64,Int64,Float64},#lorenz,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Vern9,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}}}) at C:\Users\sebastian\.julia\v0.6\Plots\src\pipeline.jl:81
 [7] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Float64,2,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,Tuple{Int64,Int64,Float64},#lorenz,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Vern9,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}}}) at C:\Users\sebastian\.julia\v0.6\Plots\src\plot.jl:179
 [8] plot(::DiffEqBase.ODESolution{Float64,2,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,Tuple{Int64,Int64,Float64},#lorenz,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Vern9,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}}) at C:\Users\sebastian\.julia\v0.6\Plots\src\plot.jl:52

[ChrisRackauckas]

This and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/421 are the same issue. It's due to the loop in ode_addsteps! doing an eachindex loop using the saved k values, but this assumes all k are saved. It should just view to the right subset to build the subset of values to save. This should get fixed during the v0.7 changes to broadcasting.

[ChrisRackauckas]

@YingboMa did this get fixed by #716 ?

[YingboMa]

This issue wasn't fixed by 716. Here is a MWE and a workaround.

julia> using OrdinaryDiffEq
[ Info: Precompiling OrdinaryDiffEq [1dea7af3-3e70-54e6-95c3-0bf5283fa5ed]

julia> prob = ODEProblem((du,u,p,t)->du.=u, [0.1], (0., 10.0));

julia> sol = solve(prob, Vern9(), save_idxs=1);

julia> sol([0.1, 0.2])
ERROR: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type Float64
Closest candidates are:
  convert(::Type{T<:Number}, ::T<:Number) where T<:Number at number.jl:6
  convert(::Type{T<:Number}, ::Number) where T<:Number at number.jl:7
  convert(::Type{T<:Number}, ::Base.TwicePrecision) where T<:Number at twiceprecision.jl:250
  ...
Stacktrace:
 [1] push!(::Array{Float64,1}, ::Array{Float64,1}) at ./array.jl:853
 [2] copyat_or_push!(::Array{Float64,1}, ::Int64, ::Array{Float64,1}, ::Type{Val{true}}) at /home/scheme/.julia/packages/RecursiveArrayTools/peaIl/src/utils.jl:68
 [3] copyat_or_push!(::Array{Float64,1}, ::Int64, ::Array{Float64,1}) at /home/scheme/.julia/packages/RecursiveArrayTools/peaIl/src/utils.jl:58
 [4] addsteps!(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64, ::ODEFunction{true,getfield(Main, Symbol("##3#4")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Nothing, ::OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}, ::Bool, ::Bool, ::Bool) at /home/scheme/.julia/dev/OrdinaryDiffEq/src/dense/verner_addsteps.jl:232
 [5] addsteps!(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64, ::Function, ::Nothing, ::OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}) at /home/scheme/.julia/dev/OrdinaryDiffEq/src/dense/verner_addsteps.jl:180
 [6] ode_interpolation(::Array{Float64,1}, ::Function, ::Nothing, ::Type, ::Nothing, ::Symbol) at /home/scheme/.julia/dev/OrdinaryDiffEq/src/dense/generic_dense.jl:155
 [7] #call#8 at /home/scheme/.julia/dev/OrdinaryDiffEq/src/interp_func.jl:73 [inlined]
 [8] ODESolution at /home/scheme/.julia/dev/DiffEqBase/src/solutions/ode_solutions.jl:15 [inlined] (repeats 2 times)
 [9] top-level scope at none:0

julia> sol = solve(prob, Vern9(), save_idxs=1:1); # workaround

julia> sol([0.1, 0.2])
t: 2-element Array{Float64,1}:
 0.1
 0.2
u: 2-element Array{Array{Float64,1},1}:
 [0.110517]
 [0.12214]

[brainsMAKER]

@YingboMa Working on other issues I revisited this one and I find the workaround isn't fail safe.

Using the original problem :

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

p = (10, 28, 8/3)
u0 = [0,10.,0]
tspan = (0., 10.)
prob = ODEProblem(lorenz, u0, tspan, `p)`

I get inconsistent results :

julia> sol = solve(prob, Vern6());

julia> sol([0.1,0.2])
t: 2-element Vector{Float64}:
 0.1
 0.2
u: 2-element Vector{Vector{Float64}}:
 [8.924800440492561, 19.684304487138352, 6.028399107817964]
 [20.055499831507554, 24.827537624341137, 39.9791212199783]

julia> sol = solve(prob, Vern6(), save_idxs=1:1);

julia> sol([0.1,0.2])
t: 2-element Vector{Float64}:
 0.1
 0.2
u: 2-element Vector{Vector{Float64}}:
 [18.266342386514168]
 [26.839956722239855]

I wonder if this may be because of a bad @inbounds.