From 6cb539955805c5e41c3d9891aee88f3860c81f84 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Sat, 5 Dec 2020 13:22:56 -0600 Subject: [PATCH] Re-implement common interface bindings (#114) * Re-add OrdinaryDiffEq and remove DataStructures from Project * Update Project.toml * WIP: re-design common interface with DiffEq ecosystem * Add support for `parse_eqs` kwarg * Add stepsize control for TaylorMethod, fix TaylorMethodCache, add empty warnkeywords, add verbose kwarg * Support high-dimensional arrays in jetcoeffs!, stepsize * Fix stepping methods, re-add keyword warning list, add evaluate! method for high-dim arrays * Update common interface tests * Add time-dependent scalar ODE in common interface; fix other tests [ci skip] * Fix error message [skip ci] * Remove unused imports [skip ci] * Last minute fixes * Update .gitignore * Remove unused lines in tests * Update tests * Update tests comments * Minor fixes * Add working version of continuous callback test * Add update_jetcoeffs_cache! for callback handling * Add vector continuous callback tests in common interface * Update docs [ci skip] * Allow AbstractArray{...,N} in parsed jetcoeffs! * Update docs * Simplify oop if else block * Add StaticArrays to test deps * Add support for DynamicalODE in common interface * Add DynamicalODE tests from #108 (thanks to @SebastianM-C) * Remove timed integrations * Fix continuous and vector callbacks via overloading DiffEqBase.addsteps! for ::TaylorMethodCache * Add comments [skip ci] * In-place DynamicalODEProblem: don't convert to ODEProblem * Fix oop DynamicalODEProblem and add some tests * Fix oop ODEProblem and DynamicalODEProblem * Add/update tests * DynamicalODEProblem: when using `SVector` in oop problems, convert initial condition to mutable array * Add comments [skip ci] * Update .gitignore * Update test/common.jl * Address review * Remove show * Update TaylorSeries version --- .gitignore | 1 + Project.toml | 13 +- docs/src/common.md | 8 +- docs/src/taylorize.md | 4 +- src/common.jl | 352 +++++++++++++++++++++++------------------- src/explicitode.jl | 8 +- src/parse_eqs.jl | 4 +- test/common.jl | 204 +++++++++++++++++++++++- test/taylorize.jl | 6 +- 9 files changed, 419 insertions(+), 181 deletions(-) diff --git a/.gitignore b/.gitignore index 3df4c353f..23aa59dca 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ examples/.ipynb_checkpoints examples/JuliaCon2017/.ipynb_checkpoints/ docs/build/ docs/site/ +.vscode/ diff --git a/Project.toml b/Project.toml index d4865fd5c..c26ad85ee 100644 --- a/Project.toml +++ b/Project.toml @@ -4,24 +4,28 @@ repo = "https://github.com/PerezHz/TaylorIntegration.jl.git" version = "0.8.5" [deps] -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Espresso = "6912e4f1-e036-58b0-9138-08d1e6358ea9" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" [compat] -DataStructures = "0.17.2" DiffEqBase = "6" Elliptic = "0.5, 1.0" Espresso = "0.6.1" +OrdinaryDiffEq = "5" +RecursiveArrayTools = "2.2" Reexport = "0.2" Requires = "0.5.2, 1.0" -TaylorSeries = "0.10.2" +StaticArrays = "0.12.5" +TaylorSeries = "0.10.9" julia = "1" [extras] @@ -29,8 +33,9 @@ Elliptic = "b305315f-e792-5b7a-8f41-49f472929428" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["TaylorSeries", "LinearAlgebra", "Test", "Elliptic", "InteractiveUtils", "Pkg"] +test = ["TaylorSeries", "LinearAlgebra", "Test", "Elliptic", "InteractiveUtils", "Pkg", "StaticArrays"] diff --git a/docs/src/common.md b/docs/src/common.md index 241e4e5be..d24da4407 100644 --- a/docs/src/common.md +++ b/docs/src/common.md @@ -13,8 +13,10 @@ aspects, we do not rely on them simply because using properly `@taylorize` improves the performance. !!! note - Currently, the only keyword argument supported by `DiffEqBase.solve` that is - implemented in `TaylorIntegration.jl` is `:saveat`. The keyword argument + Currently, the only keyword arguments supported by `DiffEqBase.solve` that + are implemented in `TaylorIntegration.jl` are `:saveat` and `:tstops`. There + is also experimental support for `:callback`, both discrete and continuous; + some examples may be found in `test/common.jl`. The keyword argument `:parse_eqs` is available in order to control the use of methods defined via [`@taylorize`](@ref). @@ -270,4 +272,4 @@ containing the differential equations. [1] Murray, Carl D., Stanley F. Dermott. Solar System dynamics. Cambridge University Press, 1999. -[2] [DiffEqBenchmarks.jl/DynamicalODE](https://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.ipynb) +[2] [SciMLBenchmarks.jl/DynamicalODE](https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html) diff --git a/docs/src/taylorize.md b/docs/src/taylorize.md index a0c04fe57..0c3cf6883 100644 --- a/docs/src/taylorize.md +++ b/docs/src/taylorize.md @@ -188,8 +188,8 @@ list some limitations and advices. - The macro allows to use array declarations through `Array`, but other ways (e.g. `similar`) are not yet implemented. -- Avoid using variables prefixed by an underscore, in particular `_T` and `_S`; - using them may lead to name collisions with some internal variables. +- Avoid using variables prefixed by an underscore, in particular `_T`, `_S` and + `_N`; using them may lead to name collisions with some internal variables. - Broadcasting is not recognized by `@taylorize`. diff --git a/src/common.jl b/src/common.jl index bce2789d5..2cb3d86e5 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,194 +1,232 @@ -using DiffEqBase, DataStructures +using DiffEqBase, OrdinaryDiffEq +using StaticArrays: SVector, SizedArray +using RecursiveArrayTools: ArrayPartition -import DiffEqBase: ODEProblem, solve +import DiffEqBase: ODEProblem, solve, ODE_DEFAULT_NORM, @.., addsteps! -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, - :maxiters, :isoutofdomain, :unstable_check, - :calck, :progress, :timeseries_steps, :tstops, :dense) +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) global warnlist = Set(warnkeywords) -abstract type TaylorAlgorithm <: DiffEqBase.DEAlgorithm end +abstract type TaylorAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm end 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 + +alg_order(alg::TaylorMethod) = alg.order + TaylorMethod() = error("Maximum order must be specified for the Taylor method") export TaylorMethod -function DiffEqBase.solve( - prob::DiffEqBase.AbstractODEProblem{uType,tupType,isinplace}, - alg::AlgType, - timeseries=[],ts=[],ks=[]; - verbose=true, - saveat=eltype(tupType)[], - abstol = 1e-6, - save_start = isempty(saveat) || saveat isa Number || prob.tspan[1] in saveat, - save_end = isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat, - timeseries_errors=true, maxiters = 1_000_000, - callback=nothing, kwargs...) where - {uType, tupType, isinplace, AlgType <: TaylorAlgorithm} - - tType = eltype(tupType) +# overload DiffEqBase.ODE_DEFAULT_NORM for Taylor1 arrays +ODE_DEFAULT_NORM(x::AbstractArray{Taylor1{T}, N},y) where {T<:Number, N} = norm(x, Inf) + +### cache stuff +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::Ref{Bool} +end - if verbose - warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist) - warned && warn_compat() - end +full_cache(c::TaylorMethodCache) = begin + tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs) +end - if haskey(prob.kwargs, :callback) || haskey(kwargs, :callback) || !isa(callback, Nothing) - error("TaylorIntegration is not compatible with callbacks.") - end +struct TaylorMethodConstantCache{uTType} <: OrdinaryDiffEqConstantCache + uT::uTType + parse_eqs::Ref{Bool} +end - sizeu = size(prob.u0) - f = prob.f.f +function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, + uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, + calck,::Val{true}) + tT = Taylor1(typeof(t), alg.order) + tT[0] = t + uT = Taylor1.(u, tT.order) + duT = zero.(Taylor1.(u, tT.order)) + uauxT = similar(uT) + TaylorMethodCache( + u, + uprev, + similar(u), + zero(rate_prototype), + zero(rate_prototype), + tT, + uT, + duT, + uauxT, + Ref(alg.parse_eqs) + ) +end - tspan = prob.tspan +alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits, + uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) = TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(alg.parse_eqs)) + +function initialize!(integrator, c::TaylorMethodConstantCache) + @unpack u, t, f, p = integrator + tT = Taylor1(typeof(t), integrator.alg.order) + tT[0] = t + c.uT .= Taylor1(u, tT.order) + c.parse_eqs.x = _determine_parsing!(c.parse_eqs.x, f, tT, c.uT, p) + __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p) + # FSAL stuff + integrator.kshortsize = 2 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + integrator.destats.nf += 1 + # Avoid undefined entries if k is an array of arrays + integrator.fsallast = zero(integrator.fsalfirst) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast +end - _, saveat_vec_, _ = tstop_saveat_disc_handling(1,saveat,1,tspan) - saveat_vec = sort(saveat_vec_.valtree) +function perform_step!(integrator,cache::TaylorMethodConstantCache) + @unpack u, t, dt, f, p = integrator + tT = Taylor1(typeof(t), integrator.alg.order) + tT[0] = t+dt + u = evaluate(cache.uT, dt) + cache.uT[0] = u + __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p) + k = f(u, p, t+dt) # For the interpolation, needs k at the updated point + integrator.destats.nf += 1 + integrator.fsallast = k + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + integrator.u = u +end - if !isempty(saveat) && save_end && saveat_vec[end] != tspan[2] - push!(saveat_vec,tspan[2]) - end - if !isempty(saveat) && saveat_vec[1] != tspan[1] - prepend!(saveat_vec,tspan[1]) - end +function initialize!(integrator, cache::TaylorMethodCache) + @unpack u, t, f, p = integrator + @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs = cache + parse_eqs.x = _determine_parsing!(parse_eqs.x, f, tT, uT, duT, p) + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + # FSAL for interpolation + integrator.fsalfirst = fsalfirst + integrator.fsallast = k + integrator.kshortsize = 1 + resize!(integrator.k, integrator.kshortsize) + integrator.k[1] = integrator.fsalfirst + # integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) + integrator.fsalfirst = constant_term.(duT) + integrator.destats.nf += 1 +end - if !isinplace && typeof(prob.u0) <: Vector{Float64} - f! = (du, u, p, t) -> (du .= f(u, p, t); 0) - if isempty(saveat_vec) - t, vectimeseries = taylorinteg(f!, prob.u0, prob.tspan[1], prob.tspan[2], - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - else - t = saveat_vec - vectimeseries = taylorinteg(f!, prob.u0, t, - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - end +function perform_step!(integrator, cache::TaylorMethodCache) + @unpack t, dt, u, f, p = integrator + @unpack k, tT, uT, duT, uauxT, parse_eqs = cache + evaluate!(uT, dt, u) + tT[0] = t+dt + for i in eachindex(u) + @inbounds uT[i][0] = u[i] + duT[i].coeffs .= zero(duT[i][0]) + end + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + k = constant_term.(duT) # For the interpolation, needs k at the updated point + integrator.destats.nf += 1 +end - elseif !isinplace && typeof(prob.u0) <: AbstractArray - f! = (du, u, p, t) -> (du .= vec(f(reshape(u, sizeu), p, t)); 0) - if isempty(saveat_vec) - t, vectimeseries = taylorinteg(f!, prob.u0, prob.tspan[1], prob.tspan[2], - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - else - t = saveat_vec - vectimeseries = taylorinteg(f!, prob.u0, t, - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - end +stepsize_controller!(integrator,alg::TaylorMethod) = stepsize(integrator.cache.uT, integrator.opts.abstol) +step_accept_controller!(integrator, alg::TaylorMethod, q) = q - elseif isinplace && typeof(prob.u0) <: AbstractArray{eltype(uType),2} - f! = (du, u, p, t) -> ( - dd = reshape(du, sizeu); uu = reshape(u, sizeu); - f(dd, uu, p, t); u = vec(uu); du = vec(dd); 0) - u0 = vec(prob.u0) - if isempty(saveat_vec) - t, vectimeseries = taylorinteg(f!, u0, prob.tspan[1], prob.tspan[2], - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - else - t = saveat_vec - vectimeseries = taylorinteg(f!, prob.u0, t, - alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false) - end - else - if haskey(kwargs, :parse_eqs) - parse_eqs = kwargs[:parse_eqs] - else - parse_eqs = true - end - if isempty(saveat_vec) - t, vectimeseries = taylorinteg(f, prob.u0, prob.tspan[1], - prob.tspan[2], alg.order, abstol, prob.p, maxsteps=maxiters, - parse_eqs=parse_eqs) - else - t = saveat_vec - vectimeseries = taylorinteg(f, prob.u0, t, - alg.order, abstol, prob.p, maxsteps=maxiters, - parse_eqs=parse_eqs) - end - end +function DiffEqBase.solve( + prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace}, + alg::TaylorMethod, args...; + verbose=true, + kwargs...) where + {uType, tupType, isinplace} - if save_start - start_idx = 1 - _t = t - else - start_idx = 2 - _t = t[2:end] + if verbose + warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist) + warned && warn_compat() end - if typeof(prob.u0) <: AbstractArray - _timeseries = Vector{uType}(undef, 0) - for i=start_idx:size(vectimeseries, 1) - push!(_timeseries, reshape(view(vectimeseries, i, :, )', sizeu)) - end + 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) + f2! = (du, v, u, p, t) -> (du .= prob.f.f2(v, u, p, t); 0) + _alg = TaylorMethod(alg.order, parse_eqs = false) + ### workaround use of `SVector` with oop `DynamicalODEProblem` + ### TODO: add proper support for oop problems with arrays + if eltype(prob.u0.x) <: SVector + _u0 = ArrayPartition(SizedArray{Tuple{length(prob.u0.x[1])}}.(prob.u0.x)) + else + _u0 = prob.u0 + end + _prob = DynamicalODEProblem(f1!, f2!, _u0.x[1], _u0.x[2], prob.tspan, prob.p; prob.kwargs...) + else + f! = (du, u, p, t) -> (du .= f(u, p, t); 0) + _alg = TaylorMethod(alg.order, parse_eqs = false) + _prob = ODEProblem(f!, prob.u0, prob.tspan, prob.p; prob.kwargs...) + end + elseif haskey(kwargs, :parse_eqs) + _alg = TaylorMethod(alg.order, parse_eqs = kwargs[:parse_eqs]) + _prob = prob else - _timeseries = vec(vectimeseries)[start_idx:end] + _alg = TaylorMethod(alg.order) + _prob = prob end - DiffEqBase.build_solution(prob, alg, _t, _timeseries, - timeseries_errors = timeseries_errors, retcode = :Success) + # 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 + DiffEqBase.solve!(integrator) + integrator.sol end -# This function is part of OrdinaryDiffEq.jl; MIT-licensed -# https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/8b5af68ee4126d3772fe690a7128468f1334b4d6/src/solve.jl#L411 -function tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan) - t0, tf = tspan - tType = eltype(tspan) - tdir = sign(tf - t0) - - tdir_t0 = tdir * t0 - tdir_tf = tdir * tf - - # time stops - tstops_internal = BinaryMinHeap{tType}() - if isempty(d_discontinuities) && isempty(tstops) # TODO: Specialize more - push!(tstops_internal, tdir_tf) - else - for t in tstops - tdir_t = tdir * t - tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t) - end - - for t in d_discontinuities - tdir_t = tdir * t - tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t) - end - - push!(tstops_internal, tdir_tf) +# used in continuous callbacks and related methods to update Taylor expansions cache +function update_jetcoeffs_cache!(u,f,p,cache::TaylorMethodCache) + @unpack tT, uT, duT, uauxT, parse_eqs = cache + for i in eachindex(u) + @inbounds uT[i][0] = u[i] + duT[i].coeffs .= zero(duT[i][0]) end + __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p) + return nothing +end - # saving time points - saveat_internal = BinaryMinHeap{tType}() - if typeof(saveat) <: Number - if (t0:saveat:tf)[end] == tf - for t in (t0 + saveat):saveat:tf - push!(saveat_internal, tdir * t) - end - else - for t in (t0 + saveat):saveat:(tf - saveat) - push!(saveat_internal, tdir * t) +# This function was modified from OrdinaryDiffEq.jl; MIT-licensed +# DiffEqBase.addsteps! overload for ::TaylorMethodCache to handle continuous +# and vector callbacks with TaylorIntegration.jl via the common interface +function DiffEqBase.addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache, + always_calc_begin = false, allow_calc_end = true,force_calc_end = false) + if length(k)<2 || always_calc_begin + 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)) end - end - elseif !isempty(saveat) - for t in saveat - tdir_t = tdir * t - tdir_t0 < tdir_t < tdir_tf && push!(saveat_internal, tdir_t) - end end - - # discontinuities - d_discontinuities_internal = BinaryMinHeap{tType}() - sizehint!(d_discontinuities_internal.valtree, length(d_discontinuities)) - for t in d_discontinuities - push!(d_discontinuities_internal, tdir * t) - end - - tstops_internal, saveat_internal, d_discontinuities_internal + update_jetcoeffs_cache!(u,f,p,cache) + nothing end diff --git a/src/explicitode.jl b/src/explicitode.jl index 73fb51272..96a889327 100644 --- a/src/explicitode.jl +++ b/src/explicitode.jl @@ -61,8 +61,8 @@ computes recursively the high-order derivates back into `x`. """ function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T}, - x::AbstractVector{Taylor1{U}}, dx::AbstractVector{Taylor1{U}}, - xaux::AbstractVector{Taylor1{U}}, params) where {T<:Real, U<:Number} + x::AbstractArray{Taylor1{U}, N}, dx::AbstractArray{Taylor1{U}, N}, + xaux::AbstractArray{Taylor1{U}, N}, params) where {T<:Real, U<:Number, N} order = x[1].order for ord in 0:order-1 ordnext = ord+1 @@ -132,8 +132,8 @@ function stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number} return h::R end -function stepsize(q::AbstractArray{Taylor1{U},1}, epsilon::T) where - {T<:Real, U<:Number} +function stepsize(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where + {T<:Real, U<:Number, N} R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T) h = convert(R, Inf) for i in eachindex(q) diff --git a/src/parse_eqs.jl b/src/parse_eqs.jl index 8b6b390d1..086cd016d 100644 --- a/src/parse_eqs.jl +++ b/src/parse_eqs.jl @@ -17,8 +17,8 @@ end) ); const _HEAD_PARSEDFN_VECTOR = sanitize(:( -function TaylorIntegration.jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, __x::AbstractVector{Taylor1{_S}}, - __dx::AbstractVector{Taylor1{_S}}, __params) where {_T<:Real, _S<:Number} +function TaylorIntegration.jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, __x::AbstractArray{Taylor1{_S}, _N}, + __dx::AbstractArray{Taylor1{_S}, _N}, __params) where {_T<:Real, _S<:Number, _N} order = __tT.order nothing diff --git a/test/common.jl b/test/common.jl index bc58f49cc..ea6ac638f 100644 --- a/test/common.jl +++ b/test/common.jl @@ -1,16 +1,24 @@ using TaylorIntegration, Test, DiffEqBase using LinearAlgebra: norm +using StaticArrays @testset "Testing `common.jl`" begin f(u,p,t) = u + g(u,p,t) = cos(t) @testset "Test integration of ODE with numbers in common interface" begin u0 = 0.5 tspan = (0.0, 1.0) prob = ODEProblem(f, u0, tspan) sol = solve(prob, TaylorMethod(50), abstol=1e-20) - - @test sol[end] - 0.5*exp(1) < 1e-12 + @test TaylorIntegration.alg_order(TaylorMethod(50)) == 50 + @test TaylorIntegration.alg_order(sol.alg) == TaylorIntegration.alg_order(TaylorMethod(50)) + @test abs(sol[end] - u0*exp(1)) < 1e-12 + u0 = 0.0 + tspan = (0.0, 11pi) + prob = ODEProblem(g, u0, tspan) + sol = solve(prob, TaylorMethod(50), abstol=1e-20) + @test abs(sol[end] - sin(sol.t[end])) < 1e-12 end f!(du, u, p, t) = (du .= u) @@ -19,8 +27,13 @@ using LinearAlgebra: norm tspan = (0.0, 1.0) prob = ODEProblem(f!, u0, tspan) sol = solve(prob, TaylorMethod(50), abstol=1e-20) - @test norm(sol[end] - u0.*exp(1)) < 1e-12 + f_oop(u, p, t) = u + prob_oop = ODEProblem(f_oop, u0, tspan) + sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20) + @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12 + @test sol.t == sol_oop.t + @test sol.u == sol_oop.u end tspan = (0.0,5.0) @@ -63,6 +76,138 @@ using LinearAlgebra: norm end end + function harmosc!(dx, x, p, t) + dx[1] = x[2] + dx[2] = - x[1] + return nothing + end + tspan = (0.0, 10pi) + abstol=1e-20 # 1e-16 + order = 25 # Taylor expansion order wrt time + u0 = [1.0; 0.0] + prob = ODEProblem(harmosc!, u0, tspan) + @testset "Test consistency with taylorinteg" begin + sol = solve(prob, TaylorMethod(order), abstol=abstol) + tv1, xv1 = taylorinteg(harmosc!, u0, tspan[1], tspan[2], order, abstol) + @test sol.t == tv1 + @test xv1[end,:] == sol[end] + end + + @testset "Test discrete callback in common interface" begin + # discrete callback: example taken from DifferentialEquations.jl docs: + # https://diffeq.sciml.ai/dev/features/callback_functions/#Using-Callbacks + t_cb = 1.0pi + condition(u,t,integrator) = t == t_cb + affect!(integrator) = integrator.u[1] += 0.1 + cb = DiscreteCallback(condition,affect!) + sol = solve(prob, TaylorMethod(order), abstol=abstol, tstops=[t_cb], callback=cb) + @test sol.t[4] == t_cb + @test sol.t[4] == sol.t[5] + @test sol[4][1] + 0.1 == sol[5][1] + end + + @testset "Test continuous callback in common interface" begin + # continuous callback: example taken from DifferentialEquations.jl docs: + # https://diffeq.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball + @taylorize function f(du,u,p,t) + local g_acc = p + du[1] = u[2] + du[2] = -g_acc + zero(u[2]) + end + function condition(u,t,integrator) + u[1] + end + function affect!(integrator) + integrator.u[2] = -integrator.u[2] + return nothing + end + cb = ContinuousCallback(condition,affect!) + u0 = [50.0,0.0] + tspan = (0.0,15.0) + p = 9.8 + prob = ODEProblem(f,u0,tspan,p) + sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + tb = sqrt(2*50/9.8) # bounce time + @test abs(tb - sol.t[9]) < 1e-14 + @test sol.t[9] == sol.t[10] + @test sol[9][1] == sol[10][1] + @test sol[9][2] == -sol[10][2] # check that callback was applied correctly (1st bounce) + @test abs(3tb - sol.t[19]) < 1e-14 + @test sol.t[19] == sol.t[20] + @test sol[19][1] == sol[20][1] + @test sol[19][2] == -sol[20][2] # check that callback was applied correctly (2nd bounce) + end + + @testset "Test vector continuous callback in common interface" begin + # vector continuous callback: example taken from DifferentialEquations.jl docs: + # https://diffeq.sciml.ai/dev/features/callback_functions/#VectorContinuousCallback-Example + @taylorize function f(du,u,p,t) + local g_acc = p + du[1] = u[2] + du[2] = -g_acc + zero(u[2]) + du[3] = u[4] + du[4] = zero(u[4]) + end + + function condition(out,u,t,integrator) # Event when event_f(u,t) == 0 + out[1] = u[1] + out[2] = (u[3] - 10.0)u[3] + end + + function affect!(integrator, idx) + if idx == 1 + integrator.u[2] = -0.9integrator.u[2] + elseif idx == 2 + integrator.u[4] = -0.9integrator.u[4] + end + end + + cb = VectorContinuousCallback(condition, affect!, 2) + + u0 = [50.0, 0.0, 0.0, 2.0] + tspan = (0.0, 15.0) + p = 9.8 + prob = ODEProblem(f, u0, tspan, p) + sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb) + tb = sqrt(2*50/9.8) # bounce time + @test abs(tb - sol.t[8]) < 1e-14 + @test sol.t[8] == sol.t[9] + @test sol[9][1] == sol[8][1] + @test sol[9][2] == -0.9sol[8][2] + @test sol[9][3] == sol[8][3] + @test sol[9][4] == sol[8][4] + @test sol.t[13] == sol.t[14] + @test sol[14][1] == sol[13][1] + @test sol[14][2] == sol[13][2] + @test sol[14][3] == sol[13][3] + @test sol[14][4] == -0.9sol[13][4] + end + + @testset "Test parsed jetcoeffs! method in common interface" begin + @taylorize function integ_vec(dx, x, p, t) + local λ = p[1] + dx[1] = cos(t) + dx[2] = -λ*sin(t) + return dx + end + @test (@isdefined integ_vec) + x0 = [0.0, 1.0] + tspan = (0.0, pi) + prob = ODEProblem(integ_vec, x0, tspan, [1.0]) + sol1 = solve(prob, TaylorMethod(order), abstol=abstol, parse_eqs=false) + sol2 = solve(prob, TaylorMethod(order), abstol=abstol) # parse_eqs=true + @test length(sol1.t) == length(sol2.t) + @test sol1.t == sol2.t + @test sol1.u == sol2.u + tv, xv = taylorinteg(integ_vec, x0, tspan[1], tspan[2], order, abstol, [1.0]) + @test sol1.t == tv + @test sol1[1,:] == xv[:,1] + @test sol1[2,:] == xv[:,2] + @test transpose(sol1[:,:]) == xv[:,:] + @test norm(sol1[end][1] - sin(sol1.t[end]), Inf) < 1.0e-15 + @test norm(sol1[end][2] - cos(sol1.t[end]), Inf) < 1.0e-15 + end + @testset "Test throwing errors in common interface" begin u0 = rand(4, 2) tspan = (0.0, 1.0) @@ -71,9 +216,56 @@ using LinearAlgebra: norm # `order` is not specified @test_throws ErrorException solve(prob, TaylorMethod(), abstol=1e-20) + end + + ### DynamicalODEProblem tests (see #108, #109) + @testset "Test integration of DynamicalODEPoblem" begin + function iip_q̇(dq,p,q,params,t) + dq[1] = p[1] + dq[2] = p[2] + end + + function iip_ṗ(dp,p,q,params,t) + dp[1] = -q[1] * (1 + 2q[2]) + dp[2] = -q[2] - (q[1]^2 - q[2]^2) + end + + iip_q0 = [0.1, 0.] + iip_p0 = [0., 0.5] + + + function oop_q̇(p, q, params, t) + p + end + + function oop_ṗ(p, q, params, t) + dp1 = -q[1] * (1 + 2q[2]) + dp2 = -q[2] - (q[1]^2 - q[2]^2) + @SVector [dp1, dp2] + end + + oop_q0 = @SVector [0.1, 0.] + oop_p0 = @SVector [0., 0.5] + + T(p) = 1//2 * (p[1]^2 + p[2]^2) + V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2]- 2//3 * q[2]^3) + H(p,q, params) = T(p) + V(q) + + E = H(iip_p0, iip_q0, nothing) - # Using a `callback` - prob2 = ODEProblem(f!, u0, tspan, callback=nothing) - @test_throws ErrorException solve(prob2, TaylorMethod(10), abstol=1e-20) + energy_err(sol) = maximum(i->H([sol[1,i], sol[2,i]], [sol[3,i], sol[4,i]], nothing)-E, 1:length(sol.u)) + + iip_prob = DynamicalODEProblem(iip_ṗ, iip_q̇, iip_p0, iip_q0, (0., 100.)) + oop_prob = DynamicalODEProblem(oop_ṗ, oop_q̇, oop_p0, oop_q0, (0., 100.)) + + sol1 = solve(iip_prob, TaylorMethod(50), abstol=1e-20) + @test energy_err(sol1) < 1e-10 + + sol2 = solve(oop_prob, TaylorMethod(50), abstol=1e-20) + @test energy_err(sol2) < 1e-10 + + @test sol1.t == sol2.t + @test sol1[:,:] == sol2[:,:] end + end diff --git a/test/taylorize.jl b/test/taylorize.jl index 19aa877eb..ccac77681 100644 --- a/test/taylorize.jl +++ b/test/taylorize.jl @@ -1070,9 +1070,9 @@ import Pkg # Ignore declarations in the function @test newex.args[1] == :( TaylorIntegration.jetcoeffs!(::Val{f!}, t::Taylor1{_T}, - q::AbstractVector{Taylor1{_S}}, - dq::AbstractVector{Taylor1{_S}}, p) where - {_T <: Real, _S <: Number}) + q::AbstractArray{Taylor1{_S}, _N}, + dq::AbstractArray{Taylor1{_S}, _N}, p) where + {_T <: Real, _S <: Number, _N}) # Include not recognized functions as they appear @test newex.args[2].args[2] == :(aa = my_simple_function(q, p, t))