Update to DiffEq5 (#78)

* update type parameters of ODEIntegrators

* trajectory works for OrdinaryDiffEq

* all tests pass on DiffEq5

* Switch default solver to SimpleATsit5

* use .alg instead

* use [:alg] god damn it

Author: Datseris

REQUIRE

@@ -1,7 +1,7 @@
 julia 0.7.0-beta2
 StaticArrays 0.8
-DiffEqBase 4.10
-OrdinaryDiffEq 4.8
+DiffEqBase 5.0
 ForwardDiff 0.8
 DelayEmbeddings
 Reexport 0.1
+SimpleDiffEq 0.1

src/continuous.jl

@@ -1,14 +1,14 @@
-using OrdinaryDiffEq, StaticArrays
-import OrdinaryDiffEq: ODEIntegrator, ODEProblem
-using DiffEqBase: __init, ODEFunction
+using DiffEqBase, StaticArrays
+using DiffEqBase: __init, ODEFunction, AbstractODEIntegrator
 
 export CDS_KWARGS
 #####################################################################################
 #                                    Defaults                                       #
 #####################################################################################
-const DEFAULT_SOLVER = Vern9()
-const DEFAULT_DIFFEQ_KWARGS = (abstol = 1e-9,
-reltol = 1e-9, maxiters = typemax(Int))
+using SimpleDiffEq: SimpleATsit5
+const DEFAULT_SOLVER = SimpleATsit5()
+const DEFAULT_DIFFEQ_KWARGS = (abstol = 1e-6,
+reltol = 1e-6, maxiters = typemax(Int))
 
 const CDS_KWARGS = (alg = DEFAULT_SOLVER, DEFAULT_DIFFEQ_KWARGS...)
 
@@ -22,26 +22,23 @@ function ContinuousDynamicalSystem(prob::ODEProblem, args...)
            t0 = prob.tspan[1])
 end
 
-function ODEProblem(ds::CDS{IIP}, tspan, args...) where {IIP}
-    # when stable, do ODEFunction(ds.f; jac = ds.jacobian)
-    return ODEProblem{IIP}(ds.f, ds.u0, tspan, args...)
+function DiffEqBase.ODEProblem(ds::CDS{IIP}, tspan, args...) where {IIP}
+    return ODEProblem{IIP}(ODEFunction(ds.f; jac = ds.jacobian), ds.u0, tspan, args...)
 end
 
 #####################################################################################
 #                                 Integrators                                       #
 #####################################################################################
-stateeltype(::ODEIntegrator{Alg, S}) where {Alg, S} = eltype(S)
-stateeltype(::ODEIntegrator{Alg, S}) where {
-    Alg, S<:Vector{<:AbstractArray{T}}} where {T} = T
+stateeltype(::AbstractODEIntegrator{A, IIP, S}) where {A, IIP, S} = eltype(S)
+stateeltype(::AbstractODEIntegrator{A, IIP, S}) where {
+    A, IIP, S<:Vector{<:AbstractArray{T}}} where {T} = T
 
 function integrator(ds::CDS{iip}, u0 = ds.u0;
     tfinal = Inf, diffeq...) where {iip}
 
     u = safe_state_type(Val{iip}(), u0)
     prob = ODEProblem{iip}(ds.f, u, (ds.t0, typeof(ds.t0)(tfinal)), ds.p)
 
-    (haskey(diffeq, :saveat) && tfinal == Inf) && error("Infinite solving!")
-
     solver = _get_solver(diffeq)
     integ = __init(prob, solver; DEFAULT_DIFFEQ_KWARGS...,
                    save_everystep = false, diffeq...)
@@ -108,12 +105,12 @@ function create_parallel(ds::CDS{true}, states)
     return paralleleom, st
 end
 
-const STIFFSOLVERS = (ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2,
-GenericImplicitEuler,
-GenericTrapezoid, SDIRK2, Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, Kvaerno4,
-KenCarp4, Kvaerno5, KenCarp5, Rosenbrock23,
-Rosenbrock32, ROS3P, Rodas3, RosShamp4, Veldd4, Velds4, GRK4T,
-GRK4A, Ros4LStab, Rodas4, Rodas42, Rodas4P)
+# const STIFFSOLVERS = (ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2,
+# GenericImplicitEuler,
+# GenericTrapezoid, SDIRK2, Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, Kvaerno4,
+# KenCarp4, Kvaerno5, KenCarp5, Rosenbrock23,
+# Rosenbrock32, ROS3P, Rodas3, RosShamp4, Veldd4, Velds4, GRK4T,
+# GRK4A, Ros4LStab, Rodas4, Rodas42, Rodas4P)
 
 function parallel_integrator(ds::CDS, states; diffeq...)
     peom, st = create_parallel(ds, states)
@@ -134,61 +131,74 @@ function trajectory(ds::ContinuousDynamicalSystem, T, u = ds.u0;
 
     t0 = ds.t0
     tvec = (t0+Ttr):dt:(T+t0+Ttr)
-    integ = integrator(ds, u; tfinal = t0 + Ttr + T, diffeq..., saveat = tvec)
-    solve!(integ)
-    return Dataset(integ.sol.u)
+    sol = Vector{SVector{dimension(ds), stateeltype(ds)}}(undef, length(tvec))
+    integ = integrator(ds, u; dt = dt, tfinal = tvec[end]+2dt, diffeq...)
+    step!(integ, Ttr)
+    for (i, t) in enumerate(tvec)
+        while t > integ.t
+            step!(integ)
+        end
+        if integ.tprev ≤ t ≤ integ.t
+            sol[i] = integ(t)
+        else
+            error("should be integ.tprev ≤ t ≤ integ.t")
+        end
+    end
+    return Dataset(sol)
 end
 
 #####################################################################################
 #                                    Get States                                     #
 #####################################################################################
-get_state(integ::ODEIntegrator{Alg, S}) where {Alg, S<:AbstractVector} = integ.u
-get_state(integ::ODEIntegrator{Alg, S}) where {Alg, S<:AbstractMatrix} =
-integ.u[:, 1]
-get_state(integ::ODEIntegrator{Alg, S}) where {Alg, S<:Vector{<:AbstractVector}} =
+get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:AbstractVector} =
+    integ.u
+get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:AbstractMatrix} =
+    integ.u[:, 1]
+get_state(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:Vector{<:AbstractVector}} =
     integ.u[1]
-get_state(integ::ODEIntegrator{Alg, S}, k::Int) where {
-    Alg, S<:Vector{<:AbstractVector}} = integ.u[k]
-get_state(integ::ODEIntegrator{Alg, S}, k::Int) where {Alg, S<:AbstractMatrix} =
+get_state(integ::AbstractODEIntegrator{Alg, IIP, S}, k::Int) where {Alg, IIP, S<:Vector{<:AbstractVector}} =
+    integ.u[k]
+get_state(integ::AbstractODEIntegrator{Alg, IIP, S}, k::Int) where {Alg, IIP, S<:AbstractMatrix} =
     integ.u[:, k]
 
 function set_state!(
-    integ::ODEIntegrator{Alg, S}, u::AbstractVector, k::Int = 1
-    ) where {Alg, S<:Vector{<:AbstractVector}}
+    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector, k::Int = 1
+    ) where {Alg, IIP, S<:Vector{<:AbstractVector}}
     integ.u[k] = u
     u_modified!(integ, true)
 end
 function set_state!(
-    integ::ODEIntegrator{Alg, S}, u::AbstractVector) where {Alg, S<:Matrix}
+    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector
+    ) where {Alg, IIP, S<:Matrix}
     integ.u[:, 1] .= u
     u_modified!(integ, true)
 end
 function set_state!(
-    integ::ODEIntegrator{Alg, S}, u::AbstractVector
-    ) where {Alg, S<:SMatrix{D, K}} where {D, K}
+    integ::AbstractODEIntegrator{Alg, IIP, S}, u::AbstractVector
+    ) where {Alg, IIP, S<:SMatrix{D, K}} where {D, K}
     integ.u = hcat(SVector{D}(u), integ.u[:, SVector{K-1}(2:K...)])
     u_modified!(integ, true)
 end
 
-get_deviations(integ::ODEIntegrator{Alg, S}) where {Alg, S<:Matrix} =
+get_deviations(integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:Matrix} =
     @view integ.u[:, 2:end]
 
 
 @generated function get_deviations(
-    integ::ODEIntegrator{Alg, S}) where {Alg, S<:SMatrix{D,K}} where {D,K}
+    integ::AbstractODEIntegrator{Alg, IIP, S}) where {Alg, IIP, S<:SMatrix{D,K}} where {D,K}
     gens = [:($k) for k=2:K]
     quote
         sind = SVector{$(K-1)}($(gens...))
         integ.u[:, sind]
     end
 end
 
-set_deviations!(integ::ODEIntegrator{Alg, S}, Q) where {Alg, S<:Matrix} =
+set_deviations!(integ::AbstractODEIntegrator{Alg, IIP, S}, Q) where {Alg, IIP, S<:Matrix} =
     (integ.u[:, 2:end] .= Q; u_modified!(integ, true))
-set_deviations!(integ::ODEIntegrator{Alg, S}, Q) where {Alg, S<:SMatrix} =
+set_deviations!(integ::AbstractODEIntegrator{Alg, IIP, S}, Q) where {Alg, IIP, S<:SMatrix} =
     (integ.u = hcat(integ.u[:,1], Q); u_modified!(integ, true))
 
-function DiffEqBase.reinit!(integ::ODEIntegrator, u0::AbstractVector,
+function DiffEqBase.reinit!(integ::AbstractODEIntegrator, u0::AbstractVector,
     Q0::AbstractMatrix; kwargs...)
 
     set_state!(integ, u0)

src/discrete.jl

@@ -176,7 +176,7 @@ end
 #####################################################################################
 #                            Creating Integrators                                   #
 #####################################################################################
-function integrator(ds::DDS{IIP, S, D, F, P}, u0 = ds.u0) where {IIP, S, D, F, P}
+function integrator(ds::DDS{IIP, S, D, F, P}, u0 = ds.u0; kwargs...) where {IIP, S, D, F, P}
     return MinimalDiscreteIntegrator{IIP, S, D, F, P}(
     ds.f, S(u0), ds.t0, S(copy(u0)), ds.p, ds.t0)
 end
@@ -187,7 +187,7 @@ function tangent_integrator(ds::DDS, k::Int = dimension(ds); kwargs...)
 end
 
 function tangent_integrator(ds::DDS{IIP, S, D, F, P, JAC, JM}, Q0::AbstractMatrix;
-    u0 = ds.u0) where {IIP, S, D, F, P, JAC, JM}
+    u0 = ds.u0, kwargs...) where {IIP, S, D, F, P, JAC, JM}
 
     Q = safe_matrix_type(Val{IIP}(), Q0)
     s = safe_state_type(Val{IIP}(), u0)
@@ -204,7 +204,7 @@ end
 # Auto-diffed in-place version
 function tangent_integrator(ds::DDS{true, S, D, F, P, JAC, JM, true},
     Q0::AbstractMatrix;
-    u0 = ds.u0) where {S, D, F, P, JAC, JM}
+    u0 = ds.u0, kwargs...) where {S, D, F, P, JAC, JM}
 
     t0 = ds.t0
     R = D + length(Q0)
@@ -226,7 +226,7 @@ function tangent_integrator(ds::DDS{true, S, D, F, P, JAC, JM, true},
     return MDI{true, SS, R, TF, P}(tangentf, s, t0, deepcopy(s), ds.p, t0)
 end
 
-function parallel_integrator(ds::DDS{IIP, S, D}, states) where {IIP, S, D}
+function parallel_integrator(ds::DDS{IIP, S, D}, states; kwargs...) where {IIP, S, D}
     peom, st = create_parallel(ds, states)
     F = typeof(peom); X = typeof(st); P = typeof(ds.p)
     return MDI{true, X, D, F, P}(peom, st, ds.t0, deepcopy(st), ds.p, ds.t0)

src/dynamicalsystem.jl

@@ -1,5 +1,4 @@
 using LinearAlgebra, DiffEqBase, ForwardDiff, StaticArrays
-import DiffEqBase: isinplace, reinit!
 using SparseArrays
 
 export dimension, get_state, DynamicalSystem
@@ -140,7 +139,7 @@ const DDS = DiscreteDynamicalSystem
 systemtype(::DDS) = "discrete"
 
 
-isinplace(::DS{IIP}) where {IIP} = IIP
+DiffEqBase.isinplace(::DS{IIP}) where {IIP} = IIP
 statetype(::DS{IIP, S}) where {IIP, S} = S
 stateeltype(::DS{IIP, S}) where {IIP, S} = eltype(S)
 isautodiff(::DS{IIP, S, D, F, P, JAC, JM, IAD}) where
@@ -289,7 +288,7 @@ end
 #                                    Jacobians                                        #
 #######################################################################################
 function create_jacobian(
-    f::F, ::Val{IIP}, s::S, p::P, t::T, ::Val{D}) where {F, IIP, S, P, T, D}
+    @nospecialize(f::F), ::Val{IIP}, s::S, p::P, t::T, ::Val{D}) where {F, IIP, S, P, T, D}
     if IIP
         dum = deepcopy(s)
         cfg = ForwardDiff.JacobianConfig(
@@ -339,8 +338,9 @@ ds.jacobian(u, ds.p, t)
 #######################################################################################
 #                                 Tanget Dynamics                                     #
 #######################################################################################
+
 # IIP Tangent Space dynamics
-function create_tangent(f::F, jacobian::JAC, J::JM,
+function create_tangent(@nospecialize(f::F), @nospecialize(jacobian::JAC), J::JM,
     ::Val{true}, ::Val{k}) where {F, JAC, JM, k}
     J = deepcopy(J)
     tangentf = (du, u, p, t) -> begin
@@ -352,7 +352,6 @@ function create_tangent(f::F, jacobian::JAC, J::JM,
     end
     return tangentf
 end
-
 # for the case of autodiffed systems, a specialized version is created
 # so that f! is not called twice in ForwardDiff
 function create_tangent_iad(f::F, J::JM, u, p, t, ::Val{k}) where {F, JM, k}
@@ -374,7 +373,6 @@ function create_tangent_iad(f::F, J::JM, u, p, t, ::Val{k}) where {F, JM, k}
 end
 
 
-
 # OOP Tangent Space dynamics
 function create_tangent(f::F, jacobian::JAC, J::JM,
     ::Val{false}, ::Val{k}) where {F, JAC, JM, k}

test/REQUIRE

@@ -1 +1,2 @@
 DiffEqCallbacks 2.1
+OrdinaryDiffEq 5.0

test/dynsys_tangent.jl

@@ -1,15 +1,18 @@
 using DynamicalSystemsBase
-using StaticArrays, LinearAlgebra
+using OrdinaryDiffEq, SimpleDiffEq, StaticArrays, LinearAlgebra
 using Test
 
+algs = (Vern9(), Tsit5(), SimpleATsit5())
+
 using DynamicalSystemsBase: CDS, DDS, DS
 using DynamicalSystemsBase.Systems: hoop, hoop_jac, hiip, hiip_jac
 using DynamicalSystemsBase.Systems: loop, loop_jac, liip, liip_jac
 using DiffEqBase
 
 println("\nTesting tangent dynamics...")
 
-let
+for alg in algs
+
 u0 = [0, 10.0, 0]
 p = [10, 28, 8/3]
 u0h = ones(2)
@@ -23,7 +26,7 @@ PARAMS = [p, ph]
 # minimalistic lyapunov
 function lyapunov_iip(ds::DS, k)
     D = dimension(ds)
-    tode = tangent_integrator(ds, orthonormal(D,k))
+    tode = tangent_integrator(ds, orthonormal(D,k); alg = alg)
     λ = zeros(k)
     for t in 1:1000
         while tode.t < t
@@ -41,7 +44,7 @@ function lyapunov_iip(ds::DS, k)
 end
 function lyapunov_oop(ds::DS, k)
     D = dimension(ds)
-    tode = tangent_integrator(ds, orthonormal(D,k))
+    tode = tangent_integrator(ds, orthonormal(D,k); alg = alg)
     λ = zeros(k)
     ws_idx = SVector{k, Int}(collect(2:k+1))
     for t in 1:1000
@@ -60,7 +63,7 @@ function lyapunov_oop(ds::DS, k)
 end
 
 for i in 1:8
-    @testset "combination $i" begin
+    @testset "$alg combination $i" begin
         sysindx = i < 5 ? 1 : 2
         if i < 5
             if isodd(i)
@@ -84,7 +87,7 @@ for i in 1:8
         end
 
         if i < 5
-            @test 0.8 < λ[1] < 0.9
+            @test 0.8 < λ[1] < 0.92
         else
             @test 0.4 < λ[1] < 0.45
         end

test/dynsys_types.jl

@@ -1,11 +1,14 @@
 using DynamicalSystemsBase
-using Test, StaticArrays, LinearAlgebra
+using Test, StaticArrays, LinearAlgebra, OrdinaryDiffEq, SimpleDiffEq
 using DynamicalSystemsBase: CDS, DDS
 using DynamicalSystemsBase.Systems: hoop, hoop_jac, hiip, hiip_jac
 using DynamicalSystemsBase.Systems: loop, loop_jac, liip, liip_jac
 
 println("\nTesting dynamical systems...")
-let
+
+algs = (Vern9(), Tsit5(), SimpleATsit5())
+
+for alg in algs
 u0 = [0, 10.0, 0]
 p = [10, 28, 8/3]
 u0h = ones(2)
@@ -16,7 +19,7 @@ INITCOD = [u0, u0h]
 PARAMS = [p, ph]
 
 for i in 1:8
-    @testset "combination $i" begin
+    @testset "$alg combination $i" begin
         sysindx = i < 5 ? 1 : 2
         if i < 5
             if isodd(i)
@@ -39,12 +42,12 @@ for i in 1:8
         J = jacobian(ds)
         @test typeof(J) <: (iip ? Matrix : SMatrix)
 
-        tinteg = tangent_integrator(ds, orthonormal(dimension(ds), dimension(ds)))
+        tinteg = tangent_integrator(ds, orthonormal(dimension(ds), dimension(ds)); alg = alg)
         tuprev = deepcopy(get_state(tinteg))
         step!(tinteg)
         @test tuprev != get_state(tinteg)
 
-        integ = integrator(ds)
+        integ = integrator(ds; alg = alg)
         uprev = deepcopy(get_state(integ))
         step!(integ)
         @test uprev != get_state(integ)
@@ -56,14 +59,14 @@ for i in 1:8
                 step!(integ)
             end
 
-            @test get_state(tinteg) ≈ integ(tt)
+            @test get_state(tinteg) ≈ integ(tt) atol = 1e-2
         else
             @test get_state(tinteg) == get_state(integ)
         end
 
         # Currently the in-place version does not work from DiffEq's side:
         if i > 2
-            pinteg = parallel_integrator(ds, [INITCOD[sysindx], INITCOD[sysindx]])
+            pinteg = parallel_integrator(ds, [INITCOD[sysindx], INITCOD[sysindx]]; alg = alg)
             puprev = deepcopy(get_state(pinteg))
             step!(pinteg)
             @test get_state(pinteg, 1) == get_state(pinteg, 2) == get_state(pinteg)
@@ -83,7 +86,7 @@ for i in 1:8
                 @test get_state(pinteg) == get_state(integ)
             end
         else
-            pinteg = parallel_integrator(ds, [INITCOD[sysindx], INITCOD[sysindx]])
+            pinteg = parallel_integrator(ds, [INITCOD[sysindx], INITCOD[sysindx]]; alg = alg)
             puprev = deepcopy(get_state(pinteg))
             step!(pinteg)
             @test get_state(pinteg, 1) == get_state(pinteg, 2) == get_state(pinteg)
@@ -92,7 +95,7 @@ for i in 1:8
             while integ.t < tt
                 step!(integ)
             end
-            @test get_state(pinteg, 1) ≈ integ(tt)
+            @test get_state(pinteg, 1) ≈ integ(tt) atol = 1e-2
         end
     end
 end