Compatibility with DynamicalODEProblem

[SebastianM-C]

This PR tries to make the minimal amount of changes required with compatibility with DynamicalODEProblems. Unfortunately, the explicit Array conversions like

TaylorIntegration.jl/src/explicitode.jl

Line 634 in f3575ee

q0_ = convert(Array{typeof(elq0)}, q0)

make this complicated.

With this PR, I tried to convert the DynamicalODEFunctions to something that can work with the plain Array representation. The next issue was that the creation of a Val{f}

TaylorIntegration.jl/src/explicitode.jl

Line 256 in f3575ee

!isempty(methodswith(Val{f}, TaylorIntegration.jetcoeffs!))

would fail with the above created closure. I tried replicating the problem in the REPL, but the above worked and returned an empty method list, so I disable the check to continue.

The last problem was due to the fact that the result timeseries is created using uType, but this is no longer valid because the data representation was converted to Array as mentioned above. As a hack I used the eltype of uType to create a Vector{Vector{eltype(uType)}} structure.

With all these hacks I can successfully integrate a DynamicalODEProblem (in the in-place form).

cc: @ChrisRackauckas

[SebastianM-C]

I checked the integration results and it seems that something strange happens (and I get obviously wrong results), so clearly my current implementation affects the correctness of the algorithm. I'm not sure exactly where is the mistake, I will investigate further.

[lbenet]

First of all, thanks a lot for your pull request and for addressing this.

Aside from the problem you have due to array conversions, which we should discuss, there are a couple of things that I think are worth mentioning/clarifying.

I am not too familiar with DynamicalODEProblems, AFAIK they are related to Hamiltonian systems or other special second-order ODEs. The current implementation assumes that the ODEs are provided as first-order ODEs. It is not difficult to work out the actual algorithmic changes to implement a pure Taylor method for the second-order ODEs, but this may be the slow path due to many extensions that would be required with respect to the API.

Another option that may be interesting, if I understand properly the documentation, is to transform the partitioned ODE mentioned in the docs, to a system of first-order ODEs. I am thinking in something like defining a function f(dy, y, params, t) which can be passed to TaylorIntegration.jl, where dy = [du, dv], y= [u, v], u and v being the dependent variables, and du and dv can be computed from f1 and f2.

Does this make any sense?

[SebastianM-C]

My first try was to use f = prob.f which would automatically do that conversion and things would work, provided that the array type doesn't change. See
https://github.com/SciML/DiffEqBase.jl/blob/50fe37da2c1f60b213236de81cbb3adb65ad000c/src/diffeqfunction.jl#L256-L262

But this needs a compatible data representation for y in your example. That's why I hacked a function that would work on Arrays. Having AbstractArray support would solve this issue.

[ChrisRackauckas]

Another option that may be interesting, if I understand properly the documentation, is to transform the partitioned ODE mentioned in the docs, to a system of first-order ODEs. I am thinking in something like defining a function f(dy, y, params, t) which can be passed to TaylorIntegration.jl, where dy = [du, dv], y= [u, v], u and v being the dependent variables, and du and dv can be computed from f1 and f2.

This PR is transforming it from a partitioned ODE on an ArrayPartition to a first-order ODE on an ArrayPartition. By default, it transforms into the first order ODEs through its dispatch on ArrayPartition.

[lbenet]

This PR tries to make the minimal amount of changes required with compatibility with DynamicalODEProblems. Unfortunately, the explicit Array conversions like

TaylorIntegration.jl/src/explicitode.jl

Line 634 in f3575ee

q0_ = convert(Array{typeof(elq0)}, q0)

make this complicated.

I understand the problems that these lines may bring; trying to recall, we introduced that motivated by the fact that we dealt with Arrays, and that way, we could specify the actual eltype needed, to have everything type stable.

Once said that, please feel free to open-up whatever you may need to allow using AbstractArrays. We would obviously require that whatever change allows as much as possible back-ward compatibility, which means no to change other code which runs with the present version. We do use the current implementation in research.

[lbenet]

The next issue was that the creation of a Val{f}

TaylorIntegration.jl/src/explicitode.jl

Line 256 in f3575ee

!isempty(methodswith(Val{f}, TaylorIntegration.jetcoeffs!))

would fail with the above created closure. I tried replicating the problem in the REPL, but the above worked and returned an empty method list, so I disable the check to continue.

Val{f} is actually required to exploit @taylorize, which is some very low level way of optimizing allocations and operations when using TaylorIntegration. While it would be great to make this more user friendly, I would suggest to avoid this which is non-trivial (at least to me), for the time being at least, and get around this by fixing parse_eqs to false in common.jl for the DynamicalODEProblems; this is actually done for some cases when using the common interoperability.

[lbenet]

The last problem was due to the fact that the result timeseries is created using uType, but this is no longer valid because the data representation was converted to Array as mentioned above. As a hack I used the eltype of uType to create a Vector{Vector{eltype(uType)}} structure.

With all these hacks I can successfully integrate a DynamicalODEProblem (in the in-place form).

This is truly great and I am look very much forward to have this in the production version!!!

[lbenet]

Could you please add some tests, or a minimal working example, so we can help in making everything work together?

[SebastianM-C]

I added some tests based on the Henon-Heiles problem.
The in-place form gives strange results, but the out of place one seems fine. I'm not sure why.

I still have to fix the Val{f} thing. Could you expand a bit more on how that should work?

I'm not sure why, but there are a lot of errors in other tests.

[lbenet]

Thanks! With the examples in hand, I'll take a look on them and provide some feedback.

[lbenet]

Closing, since I think this has been covered by #114