BoundsError when solving ODE wirh TaylorMethod

[SebastianGuadalupe]

using DifferentialEquations
using TaylorIntegration
@taylorize function benchmark7!(du, u, p, t)
    du[1] = u[3]^3 - u[2] + u[4]
    du[2] = u[3]
    du[3] = 2
    du[4] = u[4]
    return nothing
end
tspan = (0.0, 10.0)
q0 = [0.4, 0.5, 0.3, 0]
prob = ODEProblem(benchmark7!, q0, tspan)
sol = solve(prob, TaylorMethod(10), abstol=1e-20);

UndefVarError: false not defined

Stacktrace:
 [1] macro expansion at ./logging.jl:322 [inlined]
 [2] _determine_parsing!(::Bool, ::Function, ::Taylor1{Float64}, ::Array{Taylor1{Float64},1}, ::Array{Taylor1{Float64},1}, ::DiffEqBase.NullParameters) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:276
 [3] taylorinteg(::Function, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Float64, ::DiffEqBase.NullParameters; maxsteps::Int64, parse_eqs::Bool) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/explicitode.jl:421
 [4] solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(benchmark7!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::TaylorMethod, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}; verbose::Bool, saveat::Array{Float64,1}, abstol::Float64, save_start::Bool, save_end::Bool, timeseries_errors::Bool, maxiters::Int64, callback::Nothing, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/sguadalupe/.julia/packages/TaylorIntegration/KHxcX/src/common.jl:107
 [5] top-level scope at In[15]:1

There is a workaround using 2.0+zero(u[1]) instead of only 2 for du[3]

[lbenet]

Thanks for reporting!

There are two issues here: The error, which is due to this line (and maybe others). This is solved by #105.

The second one, which is the one you are able to solve and in fact triggers the former error, is related to the declaration of u and du, which are Array{Taylor1{T}, 1}. The assignment du[3] = 2 (if not properly parsed by the macro) tries to change the type of du. This one is related to the way @taylorize parses the code; I'll take a look on this later.

[lbenet]

The first issue is solved by #105; I'm tagging a new version. The other one needs a bit more time...

[lbenet]

TL;DR: Not using the macro still throws an error. I am reluctant to make the macro work in the case the original function does not work directly (which is equivalent to use the option parse_eqs=false).

First, let me explain the (second) problem. I begin stating that@taylorize creates an internal method (TaylorIntegration.jetcoeffs!) which is supposed to optimize the integration in several aspects, but also "evaluates" the actual function so it exists in the actual scope. If by some reason the optimization does not work, it uses the "standard" integration procedure (with a warning message), which uses directly the function that has the equations of motion; in your MWE the function is benchmark7!. This is then equivalent to include the option parse_eqs=false in your solve statement, or to simply avoid the use of @taylorize.

If you do that, you get the following error:

julia> sol = solve(prob, TaylorMethod(10), abstol=1e-20, parse_eqs=false);
ERROR: BoundsError: attempt to access 1-element Array{Float64,1} at index [2]
Stacktrace:
 [1] getindex at ./array.jl:788 [inlined]
 [2] getindex at /Users/benet/.julia/packages/TaylorSeries/qn7NV/src/auxiliary.jl:80 [inlined]
 [3] jetcoeffs!(::typeof(benchmark7!), ::Taylor1{Float64}, ::Array{Taylor1{Float64},1}, ::Array{Taylor1{Float64},1}, ::Array{Taylor1{Float64},1}, ::DiffEqBase.NullParameters) at /Users/benet/.julia/dev/TaylorIntegration/src/explicitode.jl:82
...

That error is related to your actual function benchmark7!, in particular to du[3] = 2. That line, in a way, imposes that du[3] is of type Int (or Float64 after some conversion). In the best case, a conversion is performed but it may not work. The error is thrown because the function jetcoeffs! expects that du[3] to be a Taylor1{T}, whose components can then be called as du[3][i], but du[3] is a Float64.

Your solution du[3] = 2 + zero(u[1]) allows du[3] to be properly converted into a Taylor1{Float64}, and then everything works fine, with and without @taylorize. (This is what we actually do in some tests!)

It is possible to make that the macro @taylorize parses (under the hood) that equation as du[3] = 2 + zero(u[1]). Then, using the macro would work properly. However, I am reluctant to do that, since that would break the integration with parse_eqs=false, which should always work.

[mforets]

Thanks for the write-up! Even for parse_eqs=false, i would expect that Julia compiles a specific method for benchmark7!(du, u, p, t) and that we know the types of the arguments. Wouldn't that be enough for promotion in each cordinate equation du[i]? In other words, i fail to see why are du[i] so uncoupled if in fact, they are components of the same vector.

[mforets]

In simpler terms, my mental model is (let alone changing Float -> Taylor1{T})):

julia> x = [1.0, 2.0, 3.0]

julia> x[1] = 10 # implicit promotion

julia> x
3-element Array{Float64,1}:
 10.0
  2.0
  3.0

[lbenet]

I reopened the issue until thins is fully clarified

I agree with your mental view @mforets , and to some extend, that actually happens. Internally (assume parse_eqs=false) benchmark7!(du, u, p, t) is called as follows:

julia> t = Taylor1(5) # 5th order for illustration
 1.0 t + 𝒪(t⁶)

julia> u = q0 .+ zero(t)
4-element Array{Taylor1{Float64},1}:
  0.4 + 𝒪(t⁶)
  0.5 + 𝒪(t⁶)
  0.3 + 𝒪(t⁶)
  0.0 + 𝒪(t⁶)

julia> du = similar(u);

julia> benchmark7!(du, u, nothing, t)  # returns `nothing` but updates `du`.

julia> du
4-element Array{Taylor1{Float64},1}:
  - 0.473 + 𝒪(t⁶)
      0.3 + 𝒪(t⁶)
      2.0 + 𝒪(t¹)
      0.0 + 𝒪(t⁶)

Note that du[3] is (correctly!) a Taylor1{Float64}, but of order 0 which is incorrect, and therefore has no 1st component, which is the reason for the error thrown:

julia> du[3][1]
ERROR: BoundsError: attempt to access 1-element Array{Float64,1} at index [2]
Stacktrace:
 [1] getindex at ./array.jl:788 [inlined]
 [2] getindex(::Taylor1{Float64}, ::Int64) at /Users/benet/.julia/packages/TaylorSeries/qn7NV/src/auxiliary.jl:80
 [3] top-level scope at REPL[20]:1

While the conversion works (Float64 -> Taylor1{Float64}) it fails to produce the correct order: in general you do not know the order that the user would like to have. This is solved if you use zero(x[1]) or zero(t), since x[1](or t) has the correct order.

[mforets]

Ok, it makes sense, thanks for the explanation.

in general you do not know the order that the user would like to have

True; i see that taylor series order is not known at compile time. Adding the order as a type parameter is possible, though perhaps not very convenient because you have to recompile a new type for each operations that produces a taylor series of different order.

Otoh, one way out could be to define a new Taylor1Vector type, which contains the vector and stores a shared order, that behaves like

julia> t = Taylor1(5) # 5th order for illustration
 1.0 t + 𝒪(t⁶)

julia> u = q0 .+ zero(t) |> Taylor1Vector
4-element Taylor1Vector{Taylor1{Float64}, Array{Taylor1{Float64},1}:
  0.4 + 𝒪(t⁶)
  0.5 + 𝒪(t⁶)
  0.3 + 𝒪(t⁶)
  0.0 + 𝒪(t⁶)

# the order of u is 5

julia> du = similar(u);

julia> benchmark7!(du, u, nothing, t)  # returns `nothing` but updates `du`.

julia> du # if you update du[i] it will promote to the order of the Taylor1Vector
4-element Array{Taylor1{Float64},1}:
  - 0.473 + 𝒪(t⁶)
      0.3 + 𝒪(t⁶)
      2.0 + 𝒪(t⁶)
      0.0 + 𝒪(t⁶)

[mforets]

i mean writing du[3] = 2 + zero(u[1]) is just fine; the "problem" is that newcomers hit an unexpected error by using du[3] = 2.

[lbenet]

@mforets I like the idea you propose of defining a Taylor1Vector, actually a lot!

cc @PerezHz

[PerezHz]

@mforets I like the idea you propose of defining a Taylor1Vector, actually a lot!

Looks like a slick solution! It'd be nice to try it! Should that be implemented here, or rather in TaylorSeries.jl?

[lbenet]

Maybe TaylorSeries is a better place...