Introducing TaylorSolution

[PerezHz]

This PR introduces the TaylorSolution struct, a return type for taylorinteg output. In summary, whereas we now currently do:

t, x = taylorinteg(...)

the t and the x are now wrapped in a TaylorSolution, such that

sol = taylorinteg(...)

where sol isa TaylorSolution and we can retrieve the output time vector and solution array from sol.t, sol.x. Dense output is also now supported via the dense = true keyword, such that a solution sol can be called at time t as sol(t) to return the value of the solution at t. Here, "interpolation" occurs by simply evaluating the Taylor series polynomial at the appropriate time, thanks to time-step control ensuring the Taylor polynomial can be expanded at any time within the time-step size. In the proposed behavior, when dense = false, simply no dense output is produced. Besides some reworking of the TaylorIntegration internals, the DiffEq interface has also been updated, such that dense output is now also supported (we get dense output almost for free in the common interface just by adding the appropriate dispatches). This PR depends on JuliaDiff/TaylorSeries.jl#355 to work properly. (UPDATE: this PR no longer requires JuliaDiff/TaylorSeries.jl#355 to work properly).

A quick example how things work:

using TaylorIntegration

order = 25
abstol = 1e-20

@taylorize function kepler1!(dq, q, p, t)
    local μ = -1.0
    r_p2 = ((q[1]^2)+(q[2]^2))
    r_p3d2 = r_p2^1.5
    dq[1] = q[3]
    dq[2] = q[4]
    newtonianCoeff = μ / r_p3d2
    dq[3] = q[1] * newtonianCoeff
    dq[4] = q[2] * newtonianCoeff
    return nothing
end

q0 = [0.2, 0.0, 0.0, 3.0]

sol = taylorinteg(kepler1!, q0, 0.0, 10pi, order, abstol, dense=true, maxsteps=1000)
sol.t # vector of time-steps
sol.x # solution array
sol(1.547) # solution evaluated at t = 1.547

Points to solve before merging:

• Adapting TaylorSolution for root-finding and Lyapunov spectrum taylorinteg specializations (still in the works). Update: a working version is now in place.
• Agreeing on how we can add a TaylorSolution struct for ADS integrations (WIP: Automatic Domain Splitting #178).
• Update docs.
• Fix TaylorSeries v0.17.7 breaks some tests in TaylorIntegration JuliaDiff/TaylorSeries.jl#359

[lbenet]

... Dense output is also now supported via the dense = true keyword, such that a solution sol can be called at time t as sol(t) to return the value of the solution at t. Here, "interpolation" occurs by simply evaluating the Taylor series polynomial at the appropriate time, thanks to time-step control ensuring the Taylor polynomial can be expanded at any time within the time-step size. In the proposed behavior, when dense = false, simply no dense output is produced...

I haven't begun reviewing this PR. Yet, from the comment motivating this PR, I think the evaluation sol(t) should only be performed for dense output (so an @assert may be needed), or we must save all intermediate Taylor polynomials also for dense=false. Take for instance the current implementation with (time) ranges: we do not return all intermediate steps, which is fine, because the user may not care about them; to be consistent with this, sol(t) should return an error. (Note that with dense=true, even for this example with time ranges the output should be truly the dense one.)

[PerezHz]

I haven't begun reviewing this PR. Yet, from the comment motivating this PR, I think the evaluation sol(t) should only be performed for dense output (so an @Assert may be needed), or we must save all intermediate Taylor polynomials also for dense=false. Take for instance the current implementation with (time) ranges: we do not return all intermediate steps, which is fine, because the user may not care about them; to be consistent with this, sol(t) should return an error. (Note that with dense=true, even for this example with time ranges the output should be truly the dense one.)

Thanks for the comment! I agree this PR still has to be polished a bit. Currently there is already in place a mechanism to save the polynomial coefficients only when dense=false; a similar mechanism can be put in place to return an error when sol(t) is called dense=false. Otherwise, similar to what is done in the DiffEq ecosystem, a linear interpolation can be returned when dense=false, so that sol(t) still has a meaning when dense=false, but without saving the whole polynomial.

[lbenet]

In an offline discussion with @LuEdRaMo, we discussed how easy can it be to adapt this proposal with the binary-tree structure involved in ADS; see #178. The crutial aspect are the fields x (and perhaps p too) of TaylorSolution, which are abstract vectors.

The point I want to rise is, if TaylorSolution can be taylored so it also accepts theADSBinaryNode, or should it (at some point, either in the implementation of #178 or as an outer constructor) convert an ADSBinaryNode to a vector of vectors? Or should ADSBinaryNode be a subtype of vector?

Thoughs are welcome...

[lbenet]

LGTM! Thanks a lot! Have you had chance to look up the performance impact of this PR?

[lbenet]

I just left 4 suggestions (remove @time from some tests). Aside from that, I think we should simply bump the minor version, merge and register

[PerezHz]

I just left 4 suggestions (remove @time from some tests).

Thanks, I added those to run some checks; I've removed them now.

Aside from that, I think we should simply bump the minor version, merge and register

Just as a way to help users, I would go actually for bumping the patch version with deprecations (I went ahead and added some in src/deprecated.jl as well as a depwarn in the TaylorSolution inner constructor). The caveat is that, since there is a warning and we test the warning logs, tests will fail. I would afterwards propose release the new minor version (v0.16) with the deprecations removed (I can take care of both releases).

[lbenet]

Just as a way to help users, I would go actually for bumping the patch version with deprecations (I went ahead and added some in src/deprecated.jl as well as a depwarn in the TaylorSolution inner constructor). The caveat is that, since there is a warning and we test the warning logs, tests will fail. I would afterwards propose release the new minor version (v0.16) with the deprecations removed (I can take care of both releases).

As you prefer; I though about the minor version because it is a (quite) breaking release, which is not backward compatible (in the sense of the returned types/content of, say, taylorinteg...

[lbenet]

I'm trying now to have some very naive benchmarks...

[lbenet]

I run (two) extremely naive benchmarks, based on the time reported by the tests for jet transport and taylorize (second run of includeing the tested file), comparing this PR (with d360531) and v0.15.3 (main). It seems there is not a large performance hit due to this PR, and perhaps even better performance in taylorize tests, but these statements are based on very naive benchmarks...

[PerezHz]

I'm also running some naive timings on 3000 steps of the Kepler problem with 2nd-order jet transport:

varorder = 2 #the order of the variational expansion
p = set_variables("ξ", numvars=4, order=varorder) #TaylorN setup
q0 = [0.2, 0.0, 0.0, 3.0]
q0TN = q0 + p # JT initial condition
t0 = 0.0
tf = (2π)*20

@taylorize function kepler1!(dq, q, p, t)
    local μ = -1.0
    r_p2 = ( (q[1]^2) + (q[2]^2) )
    r_p3d2 = r_p2^1.5
    dq[1] = q[3]
    dq[2] = q[4]
    newtonianCoeff = μ / r_p3d2
    dq[3] = q[1] * newtonianCoeff
    dq[4] = q[2] * newtonianCoeff
    return nothing
end

On v0.15.2:

@time taylorinteg(kepler1!, q0TN, t0, tf, 25, 1e-20, Val(true), maxsteps=3000);
# 4.353918 seconds (18.52 M allocations: 1.403 GiB, 4.02% gc time)

On main (v0.15.3):

@time taylorinteg(kepler1!, q0TN, t0, tf, 25, 1e-20, Val(true), maxsteps=3000);
# 4.295033 seconds (13.29 M allocations: 1013.443 MiB, 4.12% gc time)

On jp/solution:

@time taylorinteg(kepler1!, q0TN, t0, tf, 25, 1e-20, maxsteps=3000, dense=true);
# 2.970461 seconds (12.05 M allocations: 989.286 MiB, 8.68% gc time)

Does this look more or less similar to what you're seeing?

[lbenet]

Yes, they seem consistent wrt to the elapsed time; I'm not checking allocations, but I can guess that they reproduce your results. So, aside from the depration message, I am in favor of merging and registering the new version! Thanks truly a lot!

[PerezHz]

As you prefer; I though about the minor version because it is a (quite) breaking release, which is not backward compatible (in the sense of the returned types/content of, say, taylorinteg...

Ah, you're right, I've not managed to catch all the deprecations so indeed if we release a patch version like this dependant code will break. The clean way would be to catch all deprecated behavior, but this is something I overlooked at the beginning of this PR and to do it now would be a bit of a hassle. So indeed perhaps it'd be better to go straight for the minor version, without deprecations (and I'll add a breaking changes section to the release notes).

[lbenet]

So, let's go for the minor version patch then!!