UndefRefError: access to undefined reference

[blas-ko]

I'm trying to model a differential equation (in polar coordenates) that has a periodic orbit in r=a

using TaylorSeries, TaylorIntegration
periodic_orbit(t,x) = [x[1](x[1]-a),b]

with a and b constants.

Integrating this using

taylorinteg(periodic_orbit,x0,tv,15,1.0e-15)

with

x0 = [1.0 ξ₁ + 𝒪(‖x‖¹⁶) , 0.005 + 1.0 ξ₂ + 𝒪(‖x‖¹⁶)]
tv = t0:Δt:t

outputs the following error

UndefRefError: access to undefined reference

 in jetcoeffs!(::#per_orb_bad, ::Float64, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:54
 in taylorstep!(::Function, ::Float64, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:147
 in #taylorinteg#10(::Int64, ::Function, ::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:412
 in taylorinteg(::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/<something>/.julia/v0.5/TaylorIntegration/src/jettransport.jl:392

if instead

periodic_orbit2(t,x) = [x[1](x[1]-a),b*x[2]]

is used, the integrator runs as it should.

I believe this has something to do with not making any operations for the second element of the array of the periodic_orbit function, but honestly I have no idea why his wouldn't work.

[lbenet]

I'm not sure if this answers the question, but I think your problem is related to the way you type the equations.

If I doing the following:

julia> a = 1; b=-1;

julia> periodic_orbit(t,x) = [x[1](x[1]-a),b]
periodic_orbit (generic function with 1 method)

julia> periodic_orbit2(t,x) = [x[1](x[1]-a),b*x[2]]
periodic_orbit2 (generic function with 1 method)

I get errors:

julia> periodic_orbit(0.0, [0.0, 0.0])
ERROR: MethodError: objects of type Float64 are not callable
 in periodic_orbit(::Float64, ::Array{Float64,1}) at ./REPL[13]:1

julia> periodic_orbit(0.0, [0.0, 0.0])
ERROR: MethodError: objects of type Float64 are not callable
 in periodic_orbit(::Float64, ::Array{Float64,1}) at ./REPL[13]:1

I do not get the errors if I define the equations of motion for example as: periodic_orbit(t,x) = [ x[1]*(x[1]-a), b ], where you have to note the * before the parenthesis. Maybe this is the problem.

[PerezHz]

Hi, @blas-ko! Glad you're interested in TaylorIntegration! In Julia 0.5, if I try to run

using TaylorSeries, TaylorIntegration
 a = 1; b=-1;
periodic_orbit2(t,x) = [x[1]*(x[1]-a),b]
x0 = [0.01, 0.005]
p = set_variables("ξ", numvars=2, order=15)
x0TN = x0+p
t0 = 0.0; Δt = 0.1; t = 1.0; tr = t0:Δt:t

@time xv = taylorinteg(periodic_orbit2,x0TN,tr,15,1.0e-15);

Then I get the following error:

BoundsError: attempt to access 1-element Array{TaylorSeries.TaylorN{Float64},1} at index [2]

 in getindex(::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64) at ./array.jl:386
 in jetcoeffs!(::#periodic_orbit2, ::Float64, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}) at /Users/Jorge/TaylorIntegration.jl/src/jettransport.jl:54
 in taylorstep!(::Function, ::Float64, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/jettransport.jl:147
 in #taylorinteg#14(::Int64, ::Function, ::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/jettransport.jl:412
 in taylorinteg(::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::FloatRange{Float64}, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/jettransport.jl:392

But if instead I try to run

@time xv = taylorinteg(periodic_orbit3,x0TN,tr,15,1.0e-15);

where

periodic_orbit3(t,x) = [x[1]*(x[1]-a),b*one(x[1])]

Then everything seems to work fine (note the *one(x[1]) in the second element). I realized this is related to the following: if we define x0T1N = [Taylor1([x0TN[1]], 15),Taylor1([x0TN[2]], 15)] (something similar to what taylorinteg does internally), then periodic_orbit2(t0, x0T1N) returns

2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
   - 0.0099 - 0.98 ξ₁ + 1.0 ξ₁² + 𝒪(‖x‖¹⁶) + 𝒪(t¹⁶)
                            - 1.0 + 𝒪(‖x‖¹) + 𝒪(t¹)

whereas running periodic_orbit3(t0, x0T1N) gives

2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
   - 0.0099 - 0.98 ξ₁ + 1.0 ξ₁² + 𝒪(‖x‖¹⁶) + 𝒪(t¹⁶)
                          - 1.0 + 𝒪(‖x‖¹⁶) + 𝒪(t¹⁶)

the difference being in the order of the second Taylor1 polynomial. So jetcoeffs! tries to access the higher-order coefficients of - 1.0 + 𝒪(‖x‖¹) + 𝒪(t¹) but it can't, since it is a zero-th order Taylor1 polynomial. So the fastest workaround for the time being would be to multiply b by one(x[1]), as I'm proposing in periodic_orbit3. This happens only for jet transport methods (i.e., periodic_orbit2 works fine with taylorinteg when the initial condition is a Vector{Float64}).

[lbenet]

Does the problem persist?

[PerezHz]

Hi, @blas-ko! Any update on this issue? 😄

[blas-ko]

Hey! Sorry for the long long delay on answering.
As you say, @PerezHz, it has to do on how the order on the polynomials are managed. As the different array elements have different orders, the integrator doesn't know what to do when jetcoeffs! takes them.

A solution may be a kind of promotion in the order of the polynomials in the array, where it takes the higher order so none of the internal functions have problems with this. I don't know if this may the integrator lose some generality, but if it doesn't I can do a PR with this little change.

[lbenet]

@blas-ko Sorry for the long silence, and thanks again for opening this issue; yesterday @PerezHz reminded me about it.

What I'd like to point out here is that, indeed, things are messed up due to some ambiguous behavior related to default promotion of a simple number to a TaylorN or to a Taylor1.

The following is the simplest case I could construct: xdot(t,x) = constant; below I use 3.0.

julia> using TaylorSeries, TaylorIntegration

julia> tv = 0.0:0.25:1.0 ; 

julia> xdot(t,x) = 3.0;

julia> xv = taylorinteg(xdot, 0.0, tv, 20, 1.0e-20)
5-element Array{Float64,1}:
 0.0     
 0.785268
 1.57054 
 2.3558  
 3.14107 

Notice that there is no error thrown in this case, but the result is simply wrong, and it is not related with the order or the absolute tolerance used. The correct result should be (x(t) = 3*t) the vector [0.0, 0.75, 1.5, 2.25, 3.0].

The problem is, as @PerezHz points out, related to the promotion of a constant to a Taylor1 for this case, which happens internally, probably inside of jetcoeffs!. A way of getting around this is by defining the equations of motion as xdot(t,x) = 3.0 + zero(x), or by multiplying by one(x), or using t instead of x.

I think the problem comes from TaylorSeries.jl, where zero(Taylor1{Float64}) yields an order 1 Taylor1 polynomial, and something similar occurs for TaylorN.

@blas-ko, please go ahead and send us your PR. Thanks again!

[PerezHz]

Hi, @blas-ko! Do you have yet an idea of an implementation to solve this issue?

[blas-ko]

I've discussed with @lbenet offline and it seems that #31 solves this issue when @taylorize_ode is used.
In case that this macro is not used, the problem arises when updating dx in the function. Look at the following:

julia> using TaylorSeries, TaylorIntegration;

julia> a,b = 1., 0.5;
julia> function periodic_orbit1!(t,x,dx)
           dx[1] = x[1]*(x[1]-a) 
           dx[2] = b
           nothing
       end;
julia> ξ1, ξ2 = set_variables("ξ", numvars=2, order=15);
julia> q0 =  [ξ1, ξ2 + 0.05];
julia> U, dof = eltype(q0), length(q0);
julia> dx = Array{Taylor1{U}}(dof) #note this is a Taylor1{TaylorN{U}}
2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
 #undef
 #undef
julia> periodic_orbit1!(0.0,q0,dx)
julia> dx
2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
   - 1.0 ξ₁ + 1.0 ξ₁² + 𝒪(‖x‖¹⁶) + 𝒪(t¹)
                   0.5 + 𝒪(‖x‖¹) + 𝒪(t¹)

Note that, by construction of dx, it is indeed a Taylor1{TaylorN{Float64}}} but, as it never sees the order of q0, it never knows what order does dx[2] should have; it just assigns dx[2] = b #with typeof(b) == Float64.

The only things I can think about are:

• Use TaylorSeries.fixorder to set a common order if any of the variables doesn't match. I don't know how to efficiently this could be generalized to more than 2 variables, anyway.
• Require an extra order argument in eqs_diffs! for internal use, but I think this isn't elegant at all.

My suggestion is to use @taylorize_ode always, as it's much faster anyways.

[PerezHz]

Well, I definitely agree with you @blas-ko, @taylorize_ode is the way to go! Still, while trying to find a way around this, I noticed that part of the problem could be solved. In TaylorSeries.jl src/conversion.jl, lines 77 and 78, if we define

convert{T<:Number}(::Type{TaylorN{T}}, b::T) =
    TaylorN( [HomogeneousPolynomial([b], 0)], get_order())

instead of what is currently defined:

convert{T<:Number}(::Type{TaylorN{T}}, b::T) =
    TaylorN( [HomogeneousPolynomial([b], 0)], 0) #note the 0 instead of get_order()

then the TaylorN-related part of the problem is solved:

julia> using TaylorSeries, TaylorIntegration;
INFO: Recompiling stale cache file /Users/Jorge/.julia/lib/v0.6/TaylorSeries.ji for module TaylorSeries.

julia> a,b = 1., 0.5;

julia> function periodic_orbit1!(t,x,dx)
           dx[1] = x[1]*(x[1]-a) 
           dx[2] = b
           nothing
       end;

julia> ξ1, ξ2 = set_variables("ξ", numvars=2, order=15);

julia> q0 =  [ξ1, ξ2 + 0.05]
2-element Array{TaylorSeries.TaylorN{Float64},1}:
         1.0 ξ₁ + 𝒪(‖x‖¹⁶)
  0.05 + 1.0 ξ₂ + 𝒪(‖x‖¹⁶)

julia> U, dof = eltype(q0), length(q0);

julia> dx = Array{Taylor1{U}}(dof)
2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
 #undef
 #undef

julia> periodic_orbit1!(0.0,q0,dx)

julia> dx
2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
   - 1.0 ξ₁ + 1.0 ξ₁² + 𝒪(‖x‖¹⁶) + 𝒪(t¹)
                  0.5 + 𝒪(‖x‖¹⁶) + 𝒪(t¹) #note the 16-th order TaylorN here

At least now the order of the TaylorN has the correct value. @lbenet: Incidentally, I executed TaylorSeries.jl tests on my machine with this change, and everything runs smoothly! 👌 I think this gives more consistency to conversions from constants to TaylorNs, what do you think?

Now, returning to the original problem, unfortunately, this is not enough to solve it:

julia> @time tv, xv = taylorinteg(periodic_orbit1!, q0, 0.0, 1.0, 28, 1e-20);
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] jetcoeffs!(::#periodic_orbit1!, ::Float64, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{Float64,1}) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:85
 [2] taylorstep!(::Function, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}, ::Float64, ::Float64, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Int64, ::Float64, ::Array{Float64,1}) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:189
 [3] #taylorinteg#4(::Int64, ::Function, ::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Float64, ::Float64, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:335
 [4] taylorinteg(::Function, ::Array{TaylorSeries.TaylorN{Float64},1}, ::Float64, ::Float64, ::Int64, ::Float64) at /Users/Jorge/TaylorIntegration.jl/src/explicitode.jl:313

I also realized that evaluating periodic_orbit1!(0.0,q0,dx) eventually calls convert(Taylor1{TaylorN{Float64}}, b), and I think that's where the core of the problem is... Running @which convert(Taylor1{TaylorN{Float64}}, b) tells us that the convert method we are dealing with is in line 26 of src/conversion.jl in TaylorSeries.jl:

convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::S) = Taylor1([convert(T,b)], 0)

In order to prove that this indeed the method that is called, I tried changing the 0 to a 5 and got crazy things like

julia> periodic_orbit1!(0.0,q0,dx)

julia> dx
2-element Array{TaylorSeries.Taylor1{TaylorSeries.TaylorN{Float64}},1}:
   - 1.0 ξ₁ + 1.0 ξ₁² + 𝒪(‖x‖¹⁶) + 𝒪(t¹)
                  0.5 + 𝒪(‖x‖¹⁶) + 𝒪(t⁶) #Note the t^6!!!

This tells us that the problem indeed has to do with this convert method

I also tried overloading setindex! (that's another story...), e.g. doing fill!(dx, ...) (so that dx is filled with Taylor1s of suitable order), and then defining:

setindex!{T<:Number}(a::Array{Taylor1{TaylorN{T}},1}, x::T, n::Int) = a[n] = Taylor1( TaylorN( [HomogeneousPolynomial([x], 0)], get_order()), a[n].order)

before calling periodic_orbit1!(0.0,q0,dx) but haven't been able to make it work yet...

So, yeah, it'd be tricky to turn this around since, as has been said before, as it is written now, the code has no way to tell what the "good" Taylor1 order is for the constant b before assigning it to dx[2]... A possibility would be remove this issue from the list in #32, document this "feature" 😉 , and focus ourselves on #31... what do you guys think?

[lbenet]

A short note to simply state that @taylorize solves the problem; see the tests at test/taylorize.jl for a concrete example and a way also to get around this.

@PerezHz and @blas-ko, do you agree to close this?

[PerezHz]

Thanks for the ping, @lbenet! I agree to close this

[blas-ko]

Great!! I agree also :)
Thanks !