Upgrade TaylorSeries, TaylorIntegration and TaylorModels

Project.toml

@@ -44,7 +44,7 @@ Reexport = "0.2, 1"
 Requires = "0.5, 1"
 SparseArrays = "<0.0.1, 1.6"
 StaticArrays = "0.12, 1"
-TaylorIntegration = "0.9 - 0.10, =0.10.0"  # new versions require updates
-TaylorModels = "0.6"
-TaylorSeries = "=0.12"  # new versions require updates
+TaylorIntegration = "0.15"
+TaylorModels = "0.7"
+TaylorSeries = "0.17"
 julia = "1.6"

src/Algorithms/TMJets/TMJets.jl

@@ -0,0 +1,96 @@
+"""
+    TMJets1{N, DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost
+
+Validated integration using Taylor models.
+
+### Fields
+
+- `orderQ`       -- (optional, default: `2`) order of the Taylor models for jet transport variables
+- `orderT`       -- (optional, default: `8`) order of the Taylor model in time
+- `abstol`       -- (optional, default: `1e-10`) absolute tolerance
+- `maxsteps`     -- (optional, default: `2000`) maximum number of steps in the
+                     validated integration ``x' = f(x)``
+- `adaptive`     -- (optional, default: `true`) if `true`, try decreasing the absolute
+                     tolerance each time step validation fails, until `min_abs_tol` is reached
+- `minabstol`    -- (optional, default: `1e-29`) minimum absolute tolerance for the adaptive algorithm
+- `disjointness` -- (optional, default: `ZonotopeEnclosure()`) defines the method to
+                     perform the disjointness check between the taylor model flowpipe and the invariant
+
+### Notes
+
+The argument `disjointness` allows to control how are disjointness checks
+computed, in the case where the invariant is not universal. In particular,
+`ZonotopeEnclosure()` pre-processes the taylor model with a zonotopic
+overapproximation, then performs the disjointness check with that zonotope and the
+invariant. For other options, see the documentation of `AbstractDisjointnessMethod`.
+
+This algorithm is an adaptation of the implementation in `TaylorModels.jl`
+(see copyright license in the file `reach.jl` of the current folder). The package
+`TaylorIntegration.jl` is used for jet-transport of ODEs using the Taylor method,
+and `TaylorSeries.jl` is used to work with truncated Taylor series.
+"""
+@with_kw struct TMJets1{N,DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost
+    orderQ::Int = DEFAULT_ORDER_Q_TMJETS
+    orderT::Int = DEFAULT_ORDER_T_TMJETS
+    abstol::N = DEFAULT_ABS_TOL_TMJETS
+    maxsteps::Int = DEFAULT_MAX_STEPS_TMJETS
+    adaptive::Bool = true
+    minabstol::N = Float64(TaylorModels._DEF_MINABSTOL)
+    absorb::Bool = false
+    disjointness::DM = ZonotopeEnclosure()
+end
+
+numtype(::TMJets1{N}) where {N} = N
+rsetrep(::TMJets1{N}) where {N} = TaylorModelReachSet{N}
+
+
+"""
+    TMJets1{N, DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost
+
+Validated integration using Taylor models.
+
+### Fields
+
+- `orderQ`       -- (optional, default: `2`) order of the Taylor models for jet transport variables
+- `orderT`       -- (optional, default: `8`) order of the Taylor model in time
+- `abstol`       -- (optional, default: `1e-10`) absolute tolerance
+- `maxsteps`     -- (optional, default: `2000`) maximum number of steps in the
+                     validated integration ``x' = f(x)``
+- `adaptive`     -- (optional, default: `true`) if `true`, try decreasing the absolute
+                     tolerance each time step validation fails, until `min_abs_tol` is reached
+- `minabstol`    -- (optional, default: `1e-29`) minimum absolute tolerance for the adaptive algorithm
+- `disjointness` -- (optional, default: `ZonotopeEnclosure()`) defines the method to
+                     perform the disjointness check between the taylor model flowpipe and the invariant
+
+### Notes
+
+The argument `disjointness` allows to control how are disjointness checks
+computed, in the case where the invariant is not universal. In particular,
+`ZonotopeEnclosure()` pre-processes the taylor model with a zonotopic
+overapproximation, then performs the disjointness check with that zonotope and the
+invariant. For other options, see the documentation of `AbstractDisjointnessMethod`.
+
+This algorithm is an adaptation of the implementation in `TaylorModels.jl`
+(see copyright license in the file `reach.jl` of the current folder). The package
+`TaylorIntegration.jl` is used for jet-transport of ODEs using the Taylor method,
+and `TaylorSeries.jl` is used to work with truncated Taylor series.
+"""
+@with_kw struct TMJets2{N,DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost
+    orderQ::Int = DEFAULT_ORDER_Q_TMJETS
+    orderT::Int = DEFAULT_ORDER_T_TMJETS
+    abstol::N = DEFAULT_ABS_TOL_TMJETS
+    maxsteps::Int = DEFAULT_MAX_STEPS_TMJETS
+    absorb::Bool = false
+    adaptive::Bool = true
+    minabstol::N = Float64(TaylorModels._DEF_MINABSTOL)
+    validatesteps::Int = 30
+    ε::N = 1e-10
+    δ::N = 1e-6
+    absorb_steps::Int = 3
+    disjointness::DM = ZonotopeEnclosure()
+end
+
+numtype(::TMJets2{N}) where {N} = N
+rsetrep(::TMJets2{N}) where {N} = TaylorModelReachSet{N}
+
+include("post.jl")

src/Algorithms/TMJets/common.jl

@@ -15,9 +15,9 @@ const DEFAULT_ORDER_Q_TMJETS = 2
 """
     TMJets
 
-The algorithm TMJets defaults to `TMJets21b`.
+The algorithm TMJets defaults to `TMJets2`.
 """
-const TMJets = TMJets21b
+const TMJets = TMJets2
 
 # =======================================================================
 # Initialization functions to prepare the input for validated integration

src/Algorithms/TMJets/post.jl

@@ -0,0 +1,149 @@
+function post(alg::TMJets1{N}, ivp::IVP{<:AbstractContinuousSystem}, timespan;
+        Δt0::TimeInterval=zeroI,
+        kwargs...) where {N}
+    @unpack orderQ, orderT, abstol, maxsteps, adaptive, minabstol, absorb, disjointness = alg
+
+    # initial time and final time
+    t0 = tstart(timespan)
+    T = tend(timespan)
+
+    # vector field
+    if islinear(ivp) || isaffine(ivp) # TODO: refactor with inplace_field!
+        f! = inplace_field!(ivp)
+    else
+        f! = VectorField(ivp)
+    end
+
+    # algorithm optional arguments
+    parse_eqs = get(kwargs, :parse_eqs, false)
+    params = get(kwargs, :params, nothing)
+    check_property = get(kwargs, :check_property, (t, x) -> true)
+
+    n = statedim(ivp)
+    ivp_norm = _normalize(ivp)
+    X = stateset(ivp_norm)
+
+    # fix the working variables and maximum order in the global
+    # parameters struct (_params_TaylorN_)
+    set_variables("x"; numvars=n, order=2 * orderQ)
+
+    # initial set
+    X0 = initial_state(ivp_norm)
+    X0tm = _initialize(X0, orderQ, orderT)
+
+    # optionally absorb initial remainder
+    shrink_wrapping = get(kwargs, :shrink_wrapping, true)
+    if shrink_wrapping && isa(X0tm, TaylorModelReachSet)
+        if !all(iszero, remainder(X0tm))
+            X0tm = _shrink_wrapping(X0tm)
+        end
+    end
+
+    # call external solver
+    TMSol = TaylorModels.validated_integ(f!, X0tm, t0, T, orderQ, orderT,
+                                abstol, params;
+                                maxsteps=maxsteps,
+                                parse_eqs=parse_eqs,
+                                adaptive=adaptive,
+                                minabstol=minabstol,
+                                absorb=absorb,
+                                check_property=check_property)
+    tv = TMSol.time
+    xv = TMSol.fp
+    xTM1v = TMSol.xTM
+
+    # preallocate flowpipe
+    F = Vector{TaylorModelReachSet{N}}()
+    sizehint!(F, maxsteps)
+
+    # loop over reach-sets (the first reach-set at the initial time point is ignored)
+    @inbounds for i in 2:length(tv)
+
+        # create Taylor model reach-set
+        δt = TimeInterval(tv[i] + domain(xTM1v[1, i]))
+        Ri = TaylorModelReachSet(xTM1v[:, i], δt + Δt0)
+
+        # check intersection with invariant
+        _is_intersection_empty(Ri, X, disjointness) && break
+
+        push!(F, Ri)
+    end
+
+    ext = Dict{Symbol,Any}(:tv => tv, :xv => xv, :xTM1v => xTM1v)
+    return Flowpipe(F, ext)
+end
+
+function post(alg::TMJets2{N}, ivp::IVP{<:AbstractContinuousSystem}, timespan;
+        Δt0::TimeInterval=zeroI,
+        kwargs...) where {N}
+    @unpack orderQ, orderT, abstol, maxsteps, absorb, adaptive, minabstol,
+    validatesteps, ε, δ, absorb_steps, disjointness = alg
+
+    # initial time and final time
+    t0 = tstart(timespan)
+    T = tend(timespan)
+
+    # vector field
+    if islinear(ivp) || isaffine(ivp) # TODO: refactor with inplace_field!
+        f! = inplace_field!(ivp)
+    else
+        f! = VectorField(ivp)
+    end
+
+    # algorithm optional arguments
+    parse_eqs = get(kwargs, :parse_eqs, false)
+    params = get(kwargs, :params, nothing)
+
+    n = statedim(ivp)
+    ivp_norm = _normalize(ivp)
+    X = stateset(ivp_norm)
+
+    # fix the working variables and maximum order in the global
+    # parameters struct (_params_TaylorN_)
+    set_variables("x"; numvars=n, order=2 * orderQ)
+
+    # initial set
+    X0 = initial_state(ivp_norm)
+    X0tm = _initialize(X0, orderQ, orderT)
+
+    # optionally absorb initial remainder
+    shrink_wrapping = get(kwargs, :shrink_wrapping, true)
+    if shrink_wrapping && isa(X0tm, TaylorModelReachSet)
+        if !all(iszero, remainder(X0tm))
+            X0tm = _shrink_wrapping(X0tm)
+        end
+    end
+
+    # call external solver
+    TMSol = TaylorModels.validated_integ2(f!, X0tm, t0, T, orderQ, orderT,
+                                    abstol, params;
+                                    parse_eqs=parse_eqs,
+                                    maxsteps=maxsteps,
+                                    absorb=absorb,
+                                    adaptive=adaptive,
+                                    minabstol=minabstol,
+                                    validatesteps=validatesteps,
+                                    ε=ε, δ=ε, absorb_steps=absorb_steps)
+    tv = TMSol.time
+    xv = TMSol.fp
+    xTM1v = TMSol.xTM
+
+    # build flowpipe
+    F = Vector{TaylorModelReachSet{N}}()
+    sizehint!(F, maxsteps)
+
+    # loop over reach-sets (the first reach-set at the initial time point is ignored)
+    @inbounds for i in 2:length(tv)
+
+        # create Taylor model reach-set
+        δt = TimeInterval(tv[i] + domain(xTM1v[1, i]))
+        Ri = TaylorModelReachSet(xTM1v[:, i], δt + Δt0)
+
+        # check intersection with invariant
+        _is_intersection_empty(Ri, X, disjointness) && break
+
+        push!(F, Ri)
+    end
+    ext = Dict{Symbol,Any}(:tv => tv, :xv => xv, :xTM1v => xTM1v)
+    return Flowpipe(F, ext)
+end

src/Initialization/exports.jl

@@ -18,9 +18,8 @@ export
       ORBIT,
       QINT,
       TMJets,
-      TMJets20,
-      TMJets21a,
-      TMJets21b,
+      TMJets1,
+      TMJets2,
       VREP,
 
 # Flowpipes

src/ReachabilityAnalysis.jl

@@ -49,10 +49,8 @@ include("Algorithms/FLOWSTAR/FLOWSTAR.jl")
 
 # Nonlinear
 include("Algorithms/CARLIN/CARLIN.jl")
-include("Algorithms/TMJets/TMJets21b/TMJets21b.jl")
+include("Algorithms/TMJets/TMJets.jl")
 include("Algorithms/TMJets/common.jl")
-include("Algorithms/TMJets/TMJets21a/TMJets21a.jl")
-include("Algorithms/TMJets/TMJets20/TMJets20.jl")
 include("Algorithms/QINT/QINT.jl")
 
 # ===========================================================

test/Project.toml

@@ -8,6 +8,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
+ReachabilityAnalysis = "1e97bd63-91d1-579d-8e8d-501d2b57c93f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
@@ -21,10 +22,10 @@ CDDLib = "0.5 - 0.9"
 DifferentialEquations = "6.17, 7"
 Expokit = "0.2"
 Flowstar = "0.2.4"
-IntervalArithmetic = "0.16 - 0.20, =0.20.9"  # new versions require updates and are incompatible with dependencies
+IntervalArithmetic = "0.16 - 0.20, =0.20.9"
 LazySets = "2.3"
 Polyhedra = "0.5 - 0.7"
 StaticArrays = "0.12, 1"
 Suppressor = "0.2"
 Symbolics = "4 - 5"
-TaylorModels = "0.6"
+TaylorModels = "0.7"

test/algorithms/TMJets.jl

@@ -1,4 +1,4 @@
-@testset "TMJets algorithm (TMJets21b)" begin
+@testset "TMJets algorithm (TMJets2)" begin
     prob, tspan = vanderpol()
 
     # default algorithm for nonlinear systems
@@ -33,7 +33,7 @@
     @test flowpipe(sols) isa MixedFlowpipe
 end
 
-@testset "TMJets algorithm (TMJets21b): linear IVPs" begin
+@testset "TMJets algorithm (TMJets2): linear IVPs" begin
     prob, dt = exponential_1d()
     sol = solve(prob; tspan=dt, alg=TMJets())
     @test sol.alg isa TMJets
@@ -65,10 +65,10 @@ end
     #    @test sol.alg isa TMJets
 end
 
-@testset "TMJets algorithm (TMJets20): linear IVPs" begin
+@testset "TMJets algorithm (TMJets1): linear IVPs" begin
     prob, dt = exponential_1d()
-    sol = solve(prob; tspan=dt, alg=TMJets20())
-    @test sol.alg isa TMJets20
+    sol = solve(prob; tspan=dt, alg=TMJets1())
+    @test sol.alg isa TMJets1
 
     # getter functions for a taylor model reach-set
     R = sol[1]
@@ -80,22 +80,22 @@ end
 
     # test intersection with invariant
     prob, dt = exponential_1d(; invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
-    sol_inv = solve(prob; tspan=dt, alg=TMJets20())
+    sol_inv = solve(prob; tspan=dt, alg=TMJets1())
     @test [0.3] ∈ overapproximate(sol_inv[end], Zonotope)
     m = length(sol_inv)
     # check that the following reach-set escapes the invariant
     @test [0.3] ∉ overapproximate(sol[m + 1], Zonotope)
 end
 
-@testset "1D Burgers equation (TMJets21b)" begin
+@testset "1D Burgers equation (TMJets2)" begin
     L0 = 1.0 # domain length
     U0 = 1.0 # Re = 20.
     x = range(-0.5 * L0, 0.5 * L0; length=4)
     # Initial velocity
     X0 = Singleton(-U0 * sin.(2 * π / L0 * x))
     # IVP definition
     prob = @ivp(x' = burgers!(x), dim = 4, x(0) ∈ X0)
-    sol = solve(prob; tspan=(0.0, 1.0), alg=TMJets())
+    sol = solve(prob; tspan=(0.0, 1.0), alg=TMJets2())
     @test dim(sol) == 4
 end
 

test/models/generated/Brusselator.jl

@@ -18,7 +18,7 @@ prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2)
 
 T = 18.0;
 
-alg = TMJets20(; orderT=6, orderQ=2)
+alg = TMJets1(; orderT=6, orderQ=2)
 sol = solve(prob; T=T, alg=alg);
 
 U0(r) = BallInf([1.0, 1.0], r);

test/models/generated/DuffingOscillator.jl

@@ -22,4 +22,4 @@ prob = @ivp(x' = duffing!(x), x(0) ∈ X0, dim:2)
 
 T = 2 * pi / ω;
 
-sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets21a());
+sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets1());

test/models/generated/SEIR.jl

@@ -28,5 +28,5 @@ p = [α, β, γ]
 X0 = IntervalBox(vcat(x₀, p));
 prob = @ivp(x' = seir!(x), dim:7, x(0) ∈ X0);
 
-sol = solve(prob; T=200.0, alg=TMJets21a(; orderT=7, orderQ=1))
+sol = solve(prob; T=200.0, alg=TMJets1(; orderT=7, orderQ=1))
 solz = overapproximate(sol, Zonotope);