WIP: Automatic Domain Splitting

[LuEdRaMo]

At present, TaylorIntegration (TI) supports jet transport (JT) "naively". This means that we assume that polynomials converge both with respect to time and JT variables. However, we don't enforce it. The time step is a control mechanism in the time domain, but there is no analog in the JT space. As a result, there are examples where the coefficients grow too large, causing the polynomial to no longer represent the solution accurately.

To solve this issue, this PR introduces a new taylorinteg method that implements a technique called Automatic Domain Splitting (ADS), proposed by Wittig, et. al. (2015). The idea is simple: we begin with a domain where we can trust JT works (ADSDomain). As we propagate, whenever an error metric exceeds a tolerance (stol), the current polynomial splits into two polynomials, each representing half of the original domain. The new taylorinteg method returns a customized binary tree data structure (ADSBinaryNode).

For instance, Wittig, et. al. (2015) studied a 2D Kepler problem. In the animation below, the left panel depicts a Monte Carlo propagation of a small box of initial conditions. On the other hand, the middle panel shows the "naive" JT propagation that we would perform in the current TI version. Finally, the right panel shows the equivalent ADS integration. It is noticeable that the naive JT fails to represent the solution accurately after the second revolution.

[lbenet]

Thanks a lot @LuEdRaMo for this initiative! Please let us know when shall we review your PR. Incidentally, the video was not uploaded properly.

[LuEdRaMo]

Thanks a lot @LuEdRaMo for this initiative! Please let us know when shall we review your PR. Incidentally, the video was not uploaded properly.

I've changed the video format, it should play properly now.

[lbenet]

Thanks a lot @LuEdRaMo, very nice work! I've left some comments, mostly questions, and perhaps a couple of suggestions, among other, adding a section on the documentation, and perhaps for improving marginally the performance.

I see this PR is progressing nicely, and while it is not necesary for merging this PR, I still pose the question: are you planing adding in this PR the possibility of using common interface to run it from DifferentialEquations? (If not, perhaps it is a good idea to open an issue about this for the future.)

[lbenet]

I think this test should be incorporated in the documentation too, with some very basic explanation.

[LuEdRaMo]

Thanks a lot @LuEdRaMo, very nice work! I've left some comments, mostly questions, and perhaps a couple of suggestions, among other, adding a section on the documentation, and perhaps for improving marginally the performance.

Thanks @lbenet for your suggestions.

I see this PR is progressing nicely, and while it is not necesary for merging this PR, I still pose the question: are you planing adding in this PR the possibility of using common interface to run it from DifferentialEquations? (If not, perhaps it is a good idea to open an issue about this for the future.)

I'm not familiar with the DifferentialEquations interface, so I think we should open an issue and leave that for a future PR.

[lbenet]

#192 has been merged (and the process of registering the new version is on the way), so I think you should rebase and try to adapt this PR to those changes.

[LuEdRaMo]

@lbenet @PerezHz I've adapted this PR to current master. There are still some things pending, like adding more tests and documenting the ADS integrator, but before that I'd like to hear if you have any comments.

[lbenet]

Thanks a lot @LuEdRaMo; I'll try to work on this this afternoon...

[lbenet]

Thanks a lot @LuEdRaMo, this is truly a very nice new addition!

I have left some comments, mostly related to improving allocations, which in my opinion are partially related to using some recursive functions (aside the exponential growing number of the computed regions). There are also few questions, which are purely related to my ignorance on using AbstractTrees.

Also, as you pointed out, test coverage should be improved.

[lbenet]

Suggested change

	if isnothing(node.left) && isnothing(node.right)
	()
	elseif isnothing(node.left) && !isnothing(node.right)
	(node.right,)
	elseif !isnothing(node.left) && isnothing(node.right)
	(node.left,)
	else
	(node.left, node.right)
	if isnothing(node.left) && isnothing(node.right)
	return (nothing, nothing)
	elseif isnothing(node.left) && !isnothing(node.right)
	return (nothing, node.right,)
	elseif !isnothing(node.left) && isnothing(node.right)
	return (node.left, nothing)
	else
	return (node.left, node.right)

[lbenet]

The suggestion above is related to type stability; writing it as proposed, I think you get the return type Tuple{Union{Nothing, ADSTaylorSolution}, Union{Nothing, ADSTaylorSolution}}. Otherwise, you get either Tuple{}, Tuple{ADSTaylorSolution} or Tuple{ADSTaylorSolution,ADSTaylorSolution}.

Not sure if it really matters, but I rather mention it...

[LuEdRaMo]

This function is used to iterate a tree, in this case ADSTaylorSolution, so we cannot return a ::Nothing since nothing cannot be converted to ADSTaylorSolution.

[lbenet]

But then we do have some type unstability. My point is to somehow avoid it.

[PerezHz]

Just thought that it's perhaps worth mentioning that the current proposal is following both the AbstractTrees.jl interface and AFAIU the implementation here follows closely the binary tree example provided in AbstractTrees.jl itself (see e.g. here). This is only to say that as it is now children is in principle as type unstable as the example provided by them in their documentation. Maybe the compiler is somehow able to resolve this safely? @LuEdRaMo what is the output of @code_warntype children(...) called on, say, a leaf node and a non-leaf node?

[LuEdRaMo]

The output of @code_warntype is the same for all nodes:

Arguments
  #self#::Core.Const(AbstractTrees.children)
  node::ADSTaylorSolution{Float64, 2, 4}
Body::Tuple{Vararg{Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}}}
1 ── %1  = Base.getproperty(node, :left)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %2  = TaylorIntegration.isnothing(%1)::Bool
└───       goto #4 if not %2
2 ── %4  = Base.getproperty(node, :right)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %5  = TaylorIntegration.isnothing(%4)::Bool
└───       goto #4 if not %5
3 ── %7  = ()::Core.Const(())
└───       return %7
4 ┄─ %9  = Base.getproperty(node, :left)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %10 = TaylorIntegration.isnothing(%9)::Bool
└───       goto #7 if not %10
5 ── %12 = Base.getproperty(node, :right)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %13 = TaylorIntegration.isnothing(%12)::Bool
│    %14 = !%13::Bool
└───       goto #7 if not %14
6 ── %16 = Base.getproperty(node, :right)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %17 = Core.tuple(%16)::Tuple{Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}}
└───       return %17
7 ┄─ %19 = Base.getproperty(node, :left)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %20 = TaylorIntegration.isnothing(%19)::Bool
│    %21 = !%20::Bool
└───       goto #10 if not %21
8 ── %23 = Base.getproperty(node, :right)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %24 = TaylorIntegration.isnothing(%23)::Bool
└───       goto #10 if not %24
9 ── %26 = Base.getproperty(node, :left)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %27 = Core.tuple(%26)::Tuple{Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}}
└───       return %27
10 ┄ %29 = Base.getproperty(node, :left)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %30 = Base.getproperty(node, :right)::Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}
│    %31 = Core.tuple(%29, %30)::Tuple{Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}, Union{Nothing, ADSTaylorSolution{Float64, 2, 4}}}
└───       return %31

[LuEdRaMo]

AFAIK the AbstractTrees API does not require children to return a Tuple, so I tried returning a Vector{ADSTaylorSolution{T, N, M}}, it works and is type stable.

[lbenet]

If the returned tuple yields no @code_warning problems, I think we can go with it; otherwise, returning a Vector seems the way to go!

Thanks a lot for taking a look into this!

[PerezHz]

I think I would go for the tuple, otherwise, the performance we gain from type stability we might lose by allocating a vector each time children is called. Thanks for checking this!

[lbenet]

Did you check/notice any improvements (after using for loops instead of recursive functions) in the allocations or time?

[PerezHz]

Besides what @lbenet already pointed out I only left a couple of comments. Also I noticed ADSDomain is not there anymore and there was a comment I had left there before about a constructor but now I guess this is what they call a no-code solution 😄 Then just out of curiosity: why did you decide to get rid of ADSDomain? Otherwise I think this is looking pretty solid to go ahead with adding more tests and improving docs!

[PerezHz]

Is it clear why in the end a separate RetAlloc is needed for each split? If it's just some aliasing occurring in the arguments of jetcoeffs! maybe there's a way to work around it. Or maybe it's something else; either way I think it can help to understand where this is coming from.

[lbenet]

My naive guess is that you want some sort of "thread safety" in the tree context, which I'm unable to formulate properly...

[LuEdRaMo]

My guess is that it is not a "thread safety" issue since splits are processed sequentially. When I started doing ADS integrations in NEOs I noticed unexpected results that were solved by introducing one RetAlloc per split. So, there is something in NEOs dynamical functions that tests here do not notice and affect the RetAlloc object.

[LuEdRaMo]

Then just out of curiosity: why did you decide to get rid of ADSDomain?

That change was motivated by one of your comments:

I'm thinking it can make sense to add the reference initial condition q0 as a field to ADSDomain? Then, an ADS integration can be a call to taylorinteg on an initial condition <: ADSDomain, in a similar way as JT is a dispatch on Taylor1{TaylorN{T}}.

I decided to concentrate all necessary parts of ADS in a single struct, then taylorinteg can use the ADS integrator via dispatch.

[PerezHz]

Ok I see, thanks! I'm all in favor of condensing things as much as we can and avoid repeating ourselves, and I do see the value of keeping all ADS related things in one struct; I'm just wondering at the design level (or abstraction level, if you will) whether in this case it makes sense to actually put together everything in one place. What I mean is, for me the notion of an ADSDomain as an initial condition for an ADS integration, is different from an ADSTaylorSolution as the output of an ADS integration. So I would either try to keep ADSDomain with the minimal things we need, or keep the dispatch on ADSTaylorSolution well documented.