Taylor series

[andreapasquale94]

This is to introduce TaylorSeries.jl interoperability.

[PerezHz]

Many thanks for adding this!

julia> using Ephemerides, TaylorSeries

julia> eph_spk = EphemerisProvider("de440.bsp")
 1-kernel EphemerisProvider:
 "de440.bsp"

julia> time = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> pos = ephem_vector3(eph_spk, 3, 399, time)
 3-element StaticArraysCore.SVector{3, Taylor1{Float64}} with indices SOneTo(3):
  3543.2122880199786 - 0.007819282437845182 t - 1.1070930264252339e-8 t² + 8.950613857694201e-15 t³ + 3.761528551396154e-21 t⁴ - 1.9973114109467468e-27 t⁵ + 𝒪(t⁶)
     3240.7653939449465 + 0.008093354620577769 t - 9.94896857688314e-9 t² - 7.43180726250546e-15 t³ + 5.169442616501638e-21 t⁴ + 3.7299514157257047e-28 t⁵ + 𝒪(t⁶)
  924.6896922398621 + 0.0036612834090198375 t - 2.8165693209656376e-9 t² - 3.499002547324463e-15 t³ + 1.6203643125542688e-21 t⁴ + 2.988340338110004e-28 t⁵ + 𝒪(t⁶)

I think this is working in principle and can be really useful to do computations over the JPL DE440/441 ephemerides in Julia. The only caveat would be that, to get the best performance, we might want to add some ephem_vector3 dispatches (plus sibling methods) so that allocations are minimal.

[andreapasquale94]

I just created the minimal working interface and have not optimized it for performance. This was an easy update as it doesn't imply any work on Ephemerides :) We can for sure discuss improvements.

Can TaylorSeries work in-place? In this case, we can use the following interface, which is tailored for in-place evaluations: https://github.com/JuliaSpaceMissionDesign/Ephemerides.jl/blob/9d0766d9323f62433cc8a24fdad342ea2652dfb5/src/interfaces.jl#L45 and implement some overloads for the siblings' methods.

[PerezHz]

This was an easy update as it doesn't imply any work on Ephemerides :)

If it was an easy plug-and-play then this is a good sign of the Ephemerides.jl interface design!

Can TaylorSeries work in-place? In this case, we can use the following interface, which is tailored for in-place evaluations:

It can in the sense that variables of type <:TaylorSeries.AbstractSeries (i.e., Taylor1, etc.) can be updated in-place. By reading the code for jEph.ephem_compute! it looks like internally it calls ephem_vectorX (where X=3, etc...), right? Are there in-place versions of ephem_vectorX functions?

[MicheleCeresoli]

It can in the sense that variables of type <:TaylorSeries.AbstractSeries (i.e., Taylor1, etc.) can be updated in-place. By reading the code for jEph.ephem_compute! it looks like internally it calls ephem_vectorX (where X=3, etc...), right? Are there in-place versions of ephem_vectorX functions?

Hi there! Not at the moment. In the current implementation all the ephem_vectorX functions (after the DAF and SPK segment that contain the desired data are identified) are dispatched to the corresponding spk_vectorX functions which have dedicated implementations for each SPK type and return a static array of dimension 3/6/9 or 12.

In this regard, what would be needed to improve the compatibility and performance of TaylorSeries with Ephemerides?

[MicheleCeresoli]

@PerezHz I've noticed that this extension will not be compatible with all versions of TaylorSeries. For some compatibility reasons with my registry, the highest TaylorSeries version I was able to install was the v0.13. However, the code you reported failed on that one because it did not have an implementation for Base.isless between a Float64 and Taylor1:

julia> t = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> 1.0 < t
ERROR: MethodError: no method matching isless(::Float64, ::Taylor1{Float64})

Do you by any chance recall when were those features added? This way we can properly limit the minimum supported version of TaylorSeries in the extension

[PerezHz]

If I'm not mistaken I think those features were added in v0.15 though perhaps it'd be better if @lbenet can confirm this.

[PerezHz]

In this regard, what would be needed to improve the compatibility and performance of TaylorSeries with Ephemerides?

We would have to take a deeper look at the code, but I think a way forward could be to add in-place ephem_vectorX! functions which are called in ephem_compute! and add TaylorSeries dispatches on those. But as this PR is now I think this gives us a very good starting point to play around; we can improve performance for more intensive computations down the road.

[lbenet]

If I'm not mistaken I think those features were added in v0.15 though perhaps it'd be better if @lbenet can confirm this.

Indeed, the minimum required version of TaylorSeries for < and > to work is 0.15.0. Once said this, I would strongly recommend using the last version, which is 0.18.0. It includes many improvements which allow using some inplace functions that do improve allcations.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 0% with 16 lines in your changes missing coverage. Please review.

Project coverage is 97.21%. Comparing base (cd43ddc) to head (ea3224e).

Files	Patch %	Lines
ext/TaylorSeriesExt.jl	0.00%	16 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## dev #48 +/- ## ========================================== - Coverage 97.93% 97.21% -0.72% ========================================== Files 26 27 +1 Lines 2466 2481 +15 ========================================== - Hits 2415 2412 -3 - Misses 51 69 +18 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[andreapasquale94]

We would have to take a deeper look at the code, but I think a way forward could be to add in-place ephem_vectorX! functions which are called in ephem_compute! and add TaylorSeries dispatches on those. But as this PR is now I think this gives us a very good starting point to play around; we can improve performance for more intensive computations down the road.

Ok, I think that this would anyway involve some work on Ephemerides.jl.

In short, speaking about the most relevant format (type 2/3), we now basically read the data from the binary file and create an interpolation cache. In this cache, we store both the binary data and some intermediate results that are associated with the independent variable (so, on Taylor1 in this case). So I would foresee at least the following steps on the way forward:

1. To have a look at how to extend that cache for generic independent variable types -> this would imply an Ephemerides.jl change. On this, me and @MicheleCeresoli we can for sure support.
2. Computation optimizations (mainly on the interpolators). @PerezHz, @lbenet What kind of improvements do you think would be possible here? You can use the following as a reference: https://github.com/JuliaSpaceMissionDesign/Ephemerides.jl/blob/9d0766d9323f62433cc8a24fdad342ea2652dfb5/src/interp/chebyshev.jl#L17

[PerezHz]

About point 1, the place to look would be InterpCache then? Many thanks in advance for offering support!

About point 2, in the chebyshev definition, line 20: https://github.com/JuliaSpaceMissionDesign/Ephemerides.jl/blob/9d0766d9323f62433cc8a24fdad342ea2652dfb5/src/interp/chebyshev.jl#L20

Is the eltype of Tn <:Real, say e.g. Float64 in the most general case? If this is the case, then I would think of defining a dispatch of chebyshev for t::Taylor1, i.e. function chebyshev(cache::InterpCache, cₖ, t::Taylor1, idx::Int, N::Int, ibuff=1), such that internally it uses TaylorSeries mutating functions.

[MicheleCeresoli]

@andreapasquale94 I think this pull request is ready to be merged as soon as the last test errors are fixed. I also had to add an additional function check_linktime which is dispatched based on the type of the time variable and checks whether a given SPK contains or not data for the requested epoch.

For the errors, I would just like a small confirmation from @PerezHz, assuming I want to compute the value of a function and its derivative at time t = 54, is this code correct?

julia> t = Taylor1(5)
1.0 t + 𝒪(t⁶)

julia> x = fcn(t + 54)

julia> evaluate(x)

julia> evaluate(differentiate(x))

If so, the current errors can easily be solved by reducing the test tolerances.

[MicheleCeresoli]

About point 1, the place to look would be InterpCache then? Many thanks in advance for offering support!

The InterpCache structure exists for two main reasons. First of all it provides an easy way to store and handle the buffers needed to interpolate without allocations the Chebyshev\Lagrange\Hermite polynomials. In this regard, this structure only holds the buffers to temporarily store the intermediate results. All the other SPK-related information (e.g., the polynomial coefficients), are stored in an SPK-segment specific cache, which in turns store the InterpCache.

The second reason is to introduce compatibility with ForwardDiff. Indeed, you can see the DiffCache type is used, which allows to store and retrieve the Float64 or ForwardDiff.Dual buffer type depending on the input argument.

Is the eltype of Tn <:Real, say e.g. Float64 in the most general case? If this is the case, then I would think of defining a dispatch of chebyshev for t::Taylor1, i.e. function chebyshev(cache::InterpCache, cₖ, t::Taylor1, idx::Int, N::Int, ibuff=1), such that internally it uses TaylorSeries mutating functions.

To conclude the previous sentence, the eltype of Tn can vary depending on the type of the time argument. If time is a standard number (e.g., Int orFloat64), then the elements of Tn will be of type Float64. However, if time is a ForwardDiff.Dual, then the eltype of Tn becomes something in the form of: ForwardDiff.Dual{ForwardDiff.Tag{var"#32#52", Int64}, Float64, 1}

I think it might be worth to spend a little bit more time on thinking whether a better design of this cache could more easily support different differentiation backends because, please correct me if I'm wrong, but I don't see many use cases in which one would use both TaylorSeries and ForwardDiff.

[PerezHz]

To answer your question @MicheleCeresoli if you have a function fcn and you want to compute the value of the function and its first derivative, then yes, doing an evaluate of a Taylor1 without a second argument gives you the result. Just for completeness it's probably worth mentioning that with the code you wrote, you can get the function value with x[0] which is a short-hand notation for the 0-th order Taylor coefficient of x at the point of expansion (in this case t0 = 54). Similarly, you can get the first derivative simply with x[1].

[PerezHz]

I think it might be worth to spend a little bit more time on thinking whether a better design of this cache could more easily support different differentiation backends because, please correct me if I'm wrong, but I don't see many use cases in which one would use both TaylorSeries and ForwardDiff.

I mean you probably could find a way to do it, but yeah I don't think it's a very common use case right now to combine TaylorSeries with ForwardDiff (and if it turns out to be that case can be tackled later on).

[PerezHz]

Hey, just found out that the current dispatches don't work with TaylorSeries mixtures (but can be included easily):

dq = get_variables()
time = Taylor1([dq[1], one(dq[1])])
ephem_vector3(eph_spk, 3, 399, time) # throws error

I'll post a suggestion in a bit