Re-implement common interface bindings #114

PerezHz · 2020-10-14T04:32:04Z

This PR proposes to re-implement the common interface bindings to the DifferentialEquations.jl ecosystem essentially by adding TaylorMethod as an ODE solver method of OrdinaryDiffEq.jl in the sense TaylorAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm. So the proposed approach is to extend initialize! and perform_step! methods in OrdinaryDiffEq.jl, while handling step size control via corresponding method extensions of stepsize_controller! and step_accept_controller!. By doing this, calling solve yields the same results as a direct call to taylorinteg with comparable performance, while allowing the use of a larger set of keywords with TaylorMethod, including saveat, tstops and even callback (currently, discrete callbacks seem to work okay while continuous callbacks do not).

On the TaylorIntegration.jl side, I extended jetcoeffs! and stepsize, as well as TaylorSeries.evaluate! to allow use of high-dimension arrays. Following this strategy, current tests in master pass, and also added some more tests. In particular, this PR might help to find a way out of #103 and #109.

There are some issues remaining: besides the continuous callbacks (which I'm currently trying to sort out) I'm not completely sure how to handle the FSAL stuff, but since it is used in 3rd-order Hermite interpolation in OrdinaryDiffEq, I added FSAL-related updates in initialize! and perform_step!, although I'm not sure if I did it right. Also, I left most of the keywords in the warning list warnkeywords, although a better approach might now be the inverse: to ditch the keyword warning list altogether (or at least most of it?), and handle any issues that will arise in the future.

Surely there are lots of things that can be improved, so feedback is appreciated!

cc: @ChrisRackauckas @lbenet

coveralls · 2020-10-14T04:40:55Z

Coverage decreased (-0.1%) to 94.253% when pulling c3eebe9 on jp/update-diffeq-interface into dbbd504 on master.

PerezHz · 2020-10-25T07:36:25Z

Updates: I found a way, thanks to @SebastianM-C's work in #108, to finally handle DynamicalODEProblems with TaylorIntegration.jl. A thing to note is that use of @taylorize for speedup is currently not allowed for DynamicalODEProblem. Anyway, with the latest changes this PR fixes #109. We can always take care of performance for DynamicalODEProblem, via use of @taylorize, in a future PR.

lbenet · 2020-10-25T16:30:28Z

Thanks for the great effort!

Please ping me, and perhaps @SebastianM-C too, whenever you think this is ready for review.

PerezHz · 2020-10-26T00:47:24Z

Please ping me, and perhaps @SebastianM-C too, whenever you think this is ready for review.

The latest update is that I think I found a way to make continuous and vector callbacks work (discrete callbacks were already working), which I hope will help users of TaylorIntegration.jl through the common interface.

I'm still not 100% percent sure about the FSAL thing, but things appear to be working well, so... 🤷‍♂️. The keyword warning list perhaps needs to be updated in some way as well, but here's where some feedback might help to know what's the best way to proceed going forward. I think we should add support, respectively, for in-place format in scalar problems and out-of-place format in array problems, but we might tackle that in another PR. Finally, in the tests some funny warnings started to appear:

WARNING: Method definition vec(Number) in module FiniteDiff at
C:\Users\appveyor\.julia\packages\FiniteDiff\EBys0\src\jacobians.jl:128 
overwritten in module DiffEqDiffTools at
C:\Users\appveyor\.julia\packages\DiffEqDiffTools\uD0fb\src\jacobians.jl:114.

but that seems to be related to FiniteDiff and DiffEqDiffTools, since we're not overwriting any vec methods here. EDIT: the error appears only on julia 1.0 (see e.g., here and here)

Otherwise, this PR is ready for review!

cc: @lbenet @SebastianM-C

SebastianM-C

Thanks for adding support for DynamicalODEProblems. I was thinking if a more general approach where the conversion to ODEProblem is no longer required would be possible.

SebastianM-C · 2020-10-28T00:37:01Z

src/common.jl

+    sizeu = size(prob.u0)
+    # This if block handles DynamicalODEProblems
+    # credit to @SebastianM-C for coming up with this solution
+    if prob.f isa DynamicalODEFunction


From what I understand, this special casing is needed because the jetcoeffs! function calls the user function with xaux which is an Array. I was wondering if it would be possible to just use the DiffEqFunction interface (i.e. call prob.f with f(du,u,p,t) and u with the same type as prob.u0), but I'm not sure I understand how the jetcoeffs! code works and if the above suggestion would make sense.

Thank you for your review and your suggestions! You're right, it'd be nice if we didn't have to convert everything to an ODEProblem... after reading your comment here I gave it a shot, and I think just got to overcome this conversion, at least for the in-place case (I should push that commit in a couple of minutes). The out-of-place version won't work as of now, since TaylorIntegration.jl currently does not support oop for prob.u0 <: AbstractArray. But once we add that, perhaps in another PR, we'll be able to get rid of most of the ifs and elses, while avoiding the conversion from DynamicalODEProblem to ODEProblem.

PerezHz · 2020-10-28T19:24:31Z

Just pushed a commit which fixed an error in the casing for out-of-place DynamicalODEProblem (the returned value from f wasn't being passed correctly) and added some tests.

PerezHz · 2020-10-28T19:28:13Z

@SebastianM-C I was thinking... is there a way to convert an out-of-place DynamicalODEProblem to an in-place DynamicalODEProblem? Since we now can handle in-place DynamicalODEProblem if somehow we were able to do that conversion, we could avoid the conversion to ODEProblem for out-of-place DynamicalODEProblem...

PerezHz · 2020-10-28T19:41:52Z

I was trying stuff like

if prob.f isa DynamicalODEFunction && !isinplace
    f1! = (dv, v, u, p, t) -> (dv .= prob.f.f1(v, u, p, t); 0)
    f2! = (du, v, u, p, t) -> (du .= prob.f.f2(v, u, p, t); 0)
    prob.f = DynamicalODEFunction{true}(f1!, f2!)
    # ...

which is similar to what we do for out-of-place ODEProblem, but then I get an error:

ERROR: MethodError: Cannot `convert` an object of type
  DynamicalODEFunction{true,ODEFunction{true,typeof(iip_ṗ),LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}},ODEFunction{true,typeof(iip_q̇),LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}},LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}} to an object of type
  DynamicalODEFunction{false,ODEFunction{false,typeof(oop_ṗ),LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}},ODEFunction{false,typeof(oop_q̇),LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}},LinearAlgebra.UniformScaling{Bool{}},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{},Nothing{}}
Closest candidates are:
  convert(::Type{T}, ::T) where T at essentials.jl:171
  DynamicalODEFunction{false,ODEFunction{false,typeof(oop_ṗ),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{false,typeof(oop_q̇),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) where {iip, F1, F2, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, TCV} at /Users/Jorge/.julia/packages/DiffEqBase/3iigH/src/diffeqfunction.jl:64
Stacktrace:
 [1] setproperty!(::ODEProblem{RecursiveArrayTools.ArrayPartition{Float64,Tuple{SArray{Tuple{2},Float64,1,2},SArray{Tuple{2},Float64,1,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,DynamicalODEFunction{false,ODEFunction{false,typeof(oop_ṗ),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},ODEFunction{false,typeof(oop_q̇),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Symbol, ::Function) at ./Base.jl:34
 [2] top-level scope at REPL[35]:1

In the tests, if I try

oop_prob.f = iip_prob.f

the same error shows up.

SebastianM-C

This is a much better solution as it preserves the problem type and avoids converting the user input to an Array. This should improve both performance and generality/composability.

test/common.jl

src/common.jl

SebastianM-C

The translation of out-of-place problems to in-place ones is fundamentally incompatible with immutable types. From what I understand from the DiffEq code, the OOP code path is made specifically for immutable types.

SebastianM-C · 2020-10-31T22:30:53Z

src/common.jl

+    f = prob.f
+    if !isinplace && typeof(prob.u0) <: AbstractArray
+        if prob.f isa DynamicalODEFunction
+            f1! = (dv, v, u, p, t) -> (dv .= prob.f.f1(v, u, p, t); 0)


This is what causes test to fail in the OOP case with immutable types.

Good catch, thanks! I'm proposing, as a temporary solution for this specific case, to convert the initial condition prob.u0 from an immutable to a mutable array type (see below)

PerezHz · 2020-11-02T03:05:36Z

Ok, so I think I found a compromise: I found a way to make things to work with SVectors when dealing with oop DynamicalODEProblem by converting specifically for this case the initial condition prob.u0 from an immutable array type to a mutable one (it's highly non-elegant, but does the job). I fully agree that this does not solve the incompatibility of the oop->iip translation with immutable arrays, but solving that deeper issue requires some changes under the hood of TaylorIntegration.jl itself, which I believe we should think through and tackle in another PR.

lbenet

Thanks for the tremendous work that there is in this PR! I have left mostly trivial remarks to clarify a couple of points.

LGTM and I am in favor of merging this, though maybe we may wait for @SebastianM-C's opinion.

.gitignore

lbenet · 2020-11-17T00:10:30Z

src/common.jl

+import OrdinaryDiffEq: OrdinaryDiffEqAdaptiveAlgorithm,
+OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache,
+alg_order, alg_cache, initialize!, perform_step!, @unpack,
+@cache, stepsize_controller!, step_accept_controller!
+
+# TODO: check which keywords work fine
+const warnkeywords = (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep,
+:save_end, :initialize_save, :adaptive, :dt, :reltol, :dtmax,
+:dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2,
+:qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor,
+:isoutofdomain, :unstable_check,
+:calck, :progress, :timeseries_steps, :dense)


Can you include some indentation; makes nicer reading the code.

Can you open an issue to have present the TO DO indicated here?

Thanks, added indentation and also opened the issue corresponding to this TODO

lbenet · 2020-11-17T00:12:30Z

src/common.jl

+struct _TaylorMethod <: TaylorAlgorithm
+    order::Int
+    parse_eqs::Bool
+end
+
+_TaylorMethod(order; parse_eqs=true) = _TaylorMethod(order, parse_eqs) # set `parse_eqs` to `true` by default


Naive question: why do we need two (essentially identical) structs, instead of only one? Perhaps by adjusting some constructors, everything works with only one...

Good catch! Yes, it was redundant to have those two structures, so I removed one of them

lbenet · 2020-11-17T00:14:46Z

src/common.jl

+# interpolation may be helpful (since Taylor interpolation is memory-intensive),
+# so let's set isfsal to true for now, which is the default, so no method
+# definition needed here, so we in principle should support FSAL
+# isfsal(::TaylorMethod) = true # see discussion above


Yet, this line is commented... Is it needed? We can certainly have it there commented, but then clarifying a bit would be very enlightening.

It was more of a personal reminder when trying to figure the internals of DiffEqBase and OrdinaryDiffEq; I've now removed these commented lines

lbenet · 2020-11-17T00:16:56Z

src/common.jl

+# @show macroexpand(@__MODULE__, :(@cache struct TaylorMethodCache{uType,rateType,tTType,uTType} <: OrdinaryDiffEqMutableCache
+#   u::uType
+#   uprev::uType
+#   tmp::uType
+#   k::rateType
+#   fsalfirst::rateType
+#   tT::tTType
+#   uT::uTType
+#   duT::uTType
+#   uauxT::uTType
+#   parse_eqs::Bool
+# end) )


Do we need this? Seems it is a left-over of something you were testing...

😅 sorry, thanks! Yes, this was a leftover, it's now gone

lbenet · 2020-11-17T00:25:34Z

src/common.jl

+        if typeof(cache) <: OrdinaryDiffEqMutableCache
+        rtmp = similar(u, eltype(eltype(k)))
+        f(rtmp,uprev,p,t)
+        copyat_or_push!(k,1,rtmp)
+        f(rtmp,u,p,t+dt)
+        copyat_or_push!(k,2,rtmp)
+        else
+        copyat_or_push!(k,1,f(uprev,p,t))
+        copyat_or_push!(k,2,f(u,p,t+dt))


Can you include some indentation here?

Thanks, done!

lbenet · 2020-11-17T00:26:36Z

src/common.jl

+# This function was modified from TaylorSeries.jl; MIT-licensed
+# evaluate! overload to handle AbstractArray
+import TaylorSeries: evaluate!
+function evaluate!(x::AbstractArray{Taylor1{T},N}, δt::S,
+        x0::AbstractArray{T,N}) where {T<:Number, S<:Number, N}

-    tstops_internal, saveat_internal, d_discontinuities_internal
+    # @assert length(x) == length(x0)
+    @inbounds for i in eachindex(x, x0)
+        x0[i] = evaluate( x[i], δt )
+    end
+    nothing
 end


Why not adding it directly to TaylorSeries?

This is the only point which I haven't addressed yet; I will open the corresponding PR in TaylorSeries, and it will perhaps require a new TaylorSeries release

Update: just opened the corresponding PR. See: JuliaDiff/TaylorSeries.jl#252

lbenet · 2020-11-17T00:28:15Z

test/common.jl

@@ -1,16 +1,24 @@
-using TaylorIntegration, Test, DiffEqBase
+using Test, DiffEqBase, TaylorIntegration


Does the order matter here? If so, i think it is worth adding a couple of sentences to the docs.

The order makes no difference, I just forgot to leave that line untouched. I've now gone back to the original order, so that this line doesn't show up when doing diffs in the future

lbenet · 2020-11-17T00:34:46Z

src/common.jl

-      end
+    # DiffEqBase.solve(prob, _alg, args...; kwargs...)
+    integrator = DiffEqBase.__init(_prob, _alg, args...; kwargs...)
+    integrator.dt = stepsize(integrator.cache.uT, integrator.opts.abstol) # override handle_dt! setting of initial dt


Previously, there was a default value set for abstol; I cannot see it now. So the question is if the user must now to specify it, or is it set somewhere else?

Good catch! If the user doesn't specify abstol it is set automatically by DiffEqBase.__init to abstol=1e-6

PerezHz · 2020-11-17T04:59:33Z

Thanks for the tremendous work that there is in this PR! I have left mostly trivial remarks to clarify a couple of points.

Thanks @lbenet for your review! I addressed all your comments except one, which will have to wait a bit while we add the new evaluate!method to TaylorSeries. In the meantime, please let me know if there's something that you think can be further improved.

LGTM and I am in favor of merging this, though maybe we may wait for @SebastianM-C's opinion.

I agree we should wait to hear Sebastian's opinion.

…ty warnkeywords, add verbose kwarg

…od for high-dim arrays

[ci skip]

[skip ci]

[ci skip]

via overloading DiffEqBase.addsteps! for ::TaylorMethodCache

[skip ci]

…itial condition to mutable array

[skip ci]

lbenet · 2020-12-04T15:01:05Z

Tests are passing, except in nightly; I am just including the (first) problem spotted by travis here to address it later.

lbenet · 2020-12-04T15:01:15Z

Shall I go ahead and merge?

PerezHz · 2020-12-05T05:54:09Z

Shall I go ahead and merge?

Sorry for the late reply, sure, I agree that we should merge!

lbenet · 2020-12-05T19:23:07Z

Merged! Thanks a lot!

PerezHz mentioned this pull request Oct 22, 2020

Compatibility with DynamicalODEProblem #109

Closed

SebastianM-C reviewed Oct 28, 2020

View reviewed changes

SebastianM-C reviewed Oct 31, 2020

View reviewed changes

test/common.jl Outdated Show resolved Hide resolved

src/common.jl Show resolved Hide resolved

SebastianM-C reviewed Oct 31, 2020

View reviewed changes

lbenet approved these changes Nov 17, 2020

View reviewed changes

PerezHz mentioned this pull request Nov 17, 2020

Check keyword warning list in common interface #116

Open

PerezHz mentioned this pull request Nov 17, 2020

Add evaluate! method for higher-dim AbstractArray JuliaDiff/TaylorSeries.jl#252

Merged

SebastianM-C approved these changes Nov 17, 2020

View reviewed changes

PerezHz added 13 commits December 3, 2020 12:24

Re-add OrdinaryDiffEq and remove DataStructures from Project

a431c02

Update Project.toml

09bbd9e

WIP: re-design common interface with DiffEq ecosystem

d7caaee

Add support for parse_eqs kwarg

73c6480

Add stepsize control for TaylorMethod, fix TaylorMethodCache, add emp…

cc16e9f

…ty warnkeywords, add verbose kwarg

Support high-dimensional arrays in jetcoeffs!, stepsize

b0a95fc

Fix stepping methods, re-add keyword warning list, add evaluate! meth…

387188f

…od for high-dim arrays

Update common interface tests

4b43d06

Add time-dependent scalar ODE in common interface; fix other tests

fe2f4cd

[ci skip]

Fix error message

0479369

[skip ci]

Remove unused imports

b21a83f

[skip ci]

Last minute fixes

3e83d1f

Update .gitignore

967bac9

PerezHz added 22 commits December 3, 2020 12:24

Add vector continuous callback tests in common interface

cf09f50

Update docs

219c3ab

[ci skip]

Allow AbstractArray{...,N} in parsed jetcoeffs!

828f9c8

Update docs

e64f60c

Simplify oop if else block

63f287a

Add StaticArrays to test deps

21a21d4

Add support for DynamicalODE in common interface

b8dd3ee

Add DynamicalODE tests from #108 (thanks to @SebastianM-C)

82686af

Remove timed integrations

176a778

Fix continuous and vector callbacks

894857f

via overloading DiffEqBase.addsteps! for ::TaylorMethodCache

Add comments

3254334

[skip ci]

In-place DynamicalODEProblem: don't convert to ODEProblem

d319f32

Fix oop DynamicalODEProblem and add some tests

badafd6

Fix oop ODEProblem and DynamicalODEProblem

23a3c77

Add/update tests

80440eb

DynamicalODEProblem: when using SVector in oop problems, convert in…

88d7049

…itial condition to mutable array

Add comments

da7714b

[skip ci]

Update .gitignore

0b5a55b

Update test/common.jl

dc053fe

Address review

e544fa4

Remove show

86a5bbe

Update TaylorSeries version

c3eebe9

PerezHz force-pushed the jp/update-diffeq-interface branch from 9e2b0e2 to c3eebe9 Compare December 3, 2020 18:24

lbenet merged commit 6cb5399 into master Dec 5, 2020

lbenet mentioned this pull request Jan 22, 2021

Compatibility with DynamicalODEProblem #108

Closed

lbenet deleted the jp/update-diffeq-interface branch April 26, 2022 22:31

		@@ -1,16 +1,24 @@
		using TaylorIntegration, Test, DiffEqBase
		using Test, DiffEqBase, TaylorIntegration

Re-implement common interface bindings #114

Re-implement common interface bindings #114

Conversation

PerezHz commented Oct 14, 2020 • edited Loading

coveralls commented Oct 14, 2020 • edited Loading

PerezHz commented Oct 25, 2020

lbenet commented Oct 25, 2020

PerezHz commented Oct 26, 2020 • edited Loading

SebastianM-C left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

PerezHz commented Oct 28, 2020

PerezHz commented Oct 28, 2020 • edited Loading

PerezHz commented Oct 28, 2020

SebastianM-C left a comment

Choose a reason for hiding this comment

SebastianM-C left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

PerezHz commented Nov 2, 2020 • edited Loading

lbenet left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

PerezHz commented Nov 17, 2020

lbenet commented Dec 4, 2020

lbenet commented Dec 4, 2020

PerezHz commented Dec 5, 2020

lbenet commented Dec 5, 2020

PerezHz commented Oct 14, 2020 •

edited

Loading

coveralls commented Oct 14, 2020 •

edited

Loading

PerezHz commented Oct 26, 2020 •

edited

Loading

PerezHz commented Oct 28, 2020 •

edited

Loading

PerezHz commented Nov 2, 2020 •

edited

Loading