setindex!(::SVector{3, Float64}, value, ::Int) is not defined when trying to use a non-default solver

[jarroyoe]

Describe the bug I'm trying to run AttractorsViaRecurrences using the Vern9 solver, but when I run it, it gives me the following error:

ERROR: setindex!(::SVector{3, Float64}, value, ::Int) is not defined.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] setindex!(a::SVector{3, Float64}, value::Float64, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/0T5rI/src/indexing.jl:3
  [3] copyto_unaliased!(deststyle::IndexLinear, dest::SVector{3, Float64}, srcstyle::IndexLinear, src::SVector{3, Float64})
    @ Base ./abstractarray.jl:970
  [4] copyto!(dest::SVector{3, Float64}, src::SVector{3, Float64})
    @ Base ./abstractarray.jl:950
  [5] recursivecopy!(b::SVector{3, Float64}, a::SVector{3, Float64})
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/DPGwv/src/utils.jl:30
  [6] copyat_or_push!
    @ ~/.julia/packages/RecursiveArrayTools/DPGwv/src/utils.jl:117 [inlined]
  [7] copyat_or_push!
    @ ~/.julia/packages/RecursiveArrayTools/DPGwv/src/utils.jl:110 [inlined]
  [8] reinit!(integrator::OrdinaryDiffEq.ODEIntegrator{Vern9, false, SVector{3, Float64}, Nothing, Float64, Nothing, Float64, Float64, Float64, Float64, Vector{SVector{3, Float64}}, ODESolution{Float64, 2, Vector{SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, Nothing, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern9, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{SVector{3, Float64}}, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, SVector{3, Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, u0::SVector{3, Float64}; t0::Float64, tf::Float64, erase_sol::Bool, tstops::Tuple{}, saveat::Tuple{}, d_discontinuities::Tuple{}, reset_dt::Bool, reinit_callbacks::Bool, initialize_save::Bool, reinit_cache::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/QXAKd/src/integrators/integrator_interface.jl:287
  [9] reinit!
    @ ~/.julia/packages/OrdinaryDiffEq/QXAKd/src/integrators/integrator_interface.jl:249 [inlined]
 [10] automatic_Δt_basins(integ::OrdinaryDiffEq.ODEIntegrator{Vern9, false, SVector{3, Float64}, Nothing, Float64, Nothing, Float64, Float64, Float64, Float64, Vector{SVector{3, Float64}}, ODESolution{Float64, 2, Vector{SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, Nothing, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern9, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{SVector{3, Float64}}, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, SVector{3, Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, grid::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}; N::Int64)
    @ ChaosTools ~/.julia/packages/ChaosTools/PHPDF/src/basins/attractor_mapping_recurrences.jl:228
 [11] automatic_Δt_basins(integ::OrdinaryDiffEq.ODEIntegrator{Vern9, false, SVector{3, Float64}, Nothing, Float64, Nothing, Float64, Float64, Float64, Float64, Vector{SVector{3, Float64}}, ODESolution{Float64, 2, Vector{SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, Nothing, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern9, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{SVector{3, Float64}}, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(eq), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern9ConstantCache{OrdinaryDiffEq.Vern9Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, SVector{3, Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, grid::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}})
    @ ChaosTools ~/.julia/packages/ChaosTools/PHPDF/src/basins/attractor_mapping_recurrences.jl:216
 [12] basininfo_and_integ(ds::ContinuousDynamicalSystem{false, SVector{3, Float64}, 3, typeof(eq), Nothing, DynamicalSystemsBase.var"#9#15"{typeof(eq)}, SMatrix{3, 3, Float64, 9}, true}, grid::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, Δt::Nothing, diffeq::NamedTuple{(:alg, :reltol, :abstol), Tuple{Vern9, Float64, Float64}})
    @ ChaosTools ~/.julia/packages/ChaosTools/PHPDF/src/basins/attractor_mapping_recurrences.jl:165
 [13] AttractorsViaRecurrences(ds::ContinuousDynamicalSystem{false, SVector{3, Float64}, 3, typeof(eq), Nothing, DynamicalSystemsBase.var"#9#15"{typeof(eq)}, SMatrix{3, 3, Float64, 9}, true}, grid::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}; Δt::Nothing, diffeq::NamedTuple{(:alg, :reltol, :abstol), Tuple{Vern9, Float64, Float64}}, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ChaosTools ~/.julia/packages/ChaosTools/PHPDF/src/basins/attractor_mapping_recurrences.jl:79

The weird thing is that other people who have access to this code can run it on their computers without issues. If I don't specify the algorithm, this error doesn't occur, but I need to use a better solver than Tsit5. I've tried with other solvers and still get the same error.

Minimal Working Example

using DynamicalSystems, OrdinaryDiffEq
function eq(N,p,t)
    dN1 = ((N[2]*N[3]) / (1.0 + N[2]^2))
    dN2 = ((N[3]*N[2]) / (1.0 + N[3]^2)) - ((N[2]*N[3]) / (1.0 + N[2]^2))
    dN3 = N[3]*(1.0 - N[3]) - ((N[3]*N[2]) / (1.0 + N[3]^2))

    return SVector{3}(dN1,dN2,dN3)
end

system = ContinuousDynamicalSystem(eq,[0.0,0.0,0.0],nothing)

xg = yg = zg = range(0,100; length = 100)
## Find basins of attraction
mapper = AttractorsViaRecurrences(system,(xg,yg,zg);diffeq = (alg = Vern9(),reltol=1e-16, abstol=1e-16))
basins, attractors = basins_of_attraction(mapper)

[Datseris]

Thanks, this is a bug in OrdinaryDiffEq.jl. Can you please use this MWE and open an issue there?

using OrdinaryDiffEq, StaticArrays

function eq(N,p,t)
    dN1 = ((N[2]*N[3]) / (1.0 + N[2]^2))
    dN2 = ((N[3]*N[2]) / (1.0 + N[3]^2)) - ((N[2]*N[3]) / (1.0 + N[2]^2))
    dN3 = N[3]*(1.0 - N[3]) - ((N[3]*N[2]) / (1.0 + N[3]^2))
    return SVector{3}(dN1,dN2,dN3)
end

system = ODEProblem(eq,@SVector([0.0,0.0,0.0]), (0.0, 100.0), nothing)
integ = init(system, Vern9(); reltol = 1e-16, abstol = 1e-16)
step!(integ) # works
reinit!(integ, [0.0, 0, 0]) # doesn't
reinit!(integ, @SVector([0.0, 0, 0])) # also doesn't

[Datseris]

The reason this works in other colleague's computres is that they haven't updated. This piece of code 100% worked in previous package versions so some kind of code change in DiffEq broke reinit!.

[Datseris]

Did it here https://github.com/SciML/OrdinaryDiffEq.jl/issues/1688