convert from Taylor model to SparsePolynomialZonotope

src/Initialization/init_TaylorModels.jl

@@ -0,0 +1,8 @@
+using .TaylorModels: domain
+
+# check that a vector of Taylor models has the [-1, 1] domain
+function _has_normalized_domain(vTM)
+    return all(p -> all(==(IA.interval(-1, 1)), domain(p)), vTM)
+end
+
+eval(load_taylormodels_convert_polynomial_zonotope())

src/LazySets.jl

@@ -6,7 +6,6 @@ module LazySets
 using LinearAlgebra, RecipesBase, Reexport, Requires, SparseArrays
 import GLPK, JuMP, Random, ReachabilityBase
 import IntervalArithmetic as IA
-
 using IntervalArithmetic: mince
 import IntervalArithmetic: radius, ⊂
 using LinearAlgebra: checksquare
@@ -15,13 +14,17 @@ using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle, randperm
 import RecipesBase: apply_recipe
 import SparseArrays: permute
 
+@static if VERSION < v"1.9"
+    stack(vecs) = hcat(vecs...)
+end
+
 export Arrays
 export ×, normalize, ⊂,
        subtypes
 
-# ==============
+# ================
 # ReachabilityBase
-# ==============
+# ================
 
 using ReachabilityBase.Assertions: @assert
 import ReachabilityBase.Assertions: activate_assertions, deactivate_assertions

src/convert.jl

@@ -1149,3 +1149,62 @@ end
 function convert(::Type{VPolytope}, T::Tetrahedron)
     return VPolytope(T.vertices)
 end
+
+function load_taylormodels_convert_polynomial_zonotope()
+    return quote
+        using .TaylorModels: TaylorModelN, polynomial, remainder, constant_term
+        using .TaylorModels.TaylorSeries: coeff_table
+
+        # implements Proposition 3.1.12 in thesis
+        function convert(::Type{SparsePolynomialZonotope},
+                         vTM::Vector{<:TaylorModelN{r,N}}) where {r,N}
+            @assert _has_normalized_domain(vTM) "normalized domain (-1, 1) required"
+
+            # upper bound on the number of terms/columns (ignores duplicates)
+            # - 1 per iteration because we treat the constant terms separately
+            num_coeffs = 0
+            for TMi in vTM
+                sum(TMi -> length(polynomial(TMi).coeffs) - 1, vTM)
+            end
+
+            n = length(vTM)
+            c = Vector{N}(undef, n)
+            Gs = Vector{Vector{N}}()
+            GI_diag = Vector{N}(undef, n)
+            Es = Vector{Vector{Int}}()
+
+            total_columns = 0
+            @inbounds for (i, TMi) in enumerate(vTM)
+                pol = polynomial(TMi)
+                rem = remainder(TMi)
+                c[i] = IA.mid(rem) + constant_term(pol)
+                GI_diag[i] = radius(rem)
+                for (order, term_order) in enumerate(pol.coeffs)
+                    # we skip the first (= constant) term
+                    if order == 1 || iszero(term_order)
+                        continue
+                    end
+                    for (k, coeff_k) in enumerate(term_order.coeffs)
+                        if iszero(coeff_k)
+                            continue
+                        end
+                        Ej = coeff_table[order][k]
+                        j = findfirst(e -> e == Ej, Es)
+                        if isnothing(j)
+                            total_columns += 1
+                            j = total_columns
+                            push!(Es, Ej)
+                            push!(Gs, zeros(N, n))
+                        end
+                        Gs[j][i] += coeff_k
+                    end
+                end
+            end
+            G = stack(Gs)
+            GI = remove_zero_columns(diagm(GI_diag))
+            E = stack(Es)
+            idx = uniqueID(n)
+            return SparsePolynomialZonotope(c, G, GI, E, idx)
+        end
+    end
+end  # quote / load_taylormodels_convert_polynomial_zonotope

src/init.jl

@@ -13,6 +13,7 @@ function __init__()
     # @require StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" include("Initialization/init_StaticArraysCore.jl")
     @require RangeEnclosures = "1b4d18b6-9e5d-11e9-236c-f792b01831f8" include("Initialization/init_RangeEnclosures.jl")
     @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("Initialization/init_Symbolics.jl")
+    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("Initialization/init_TaylorModels.jl")
     @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("Initialization/init_WriteVTK.jl")
 
     return nothing

test/Sets/SparsePolynomialZonotope.jl

@@ -160,4 +160,27 @@ let
 
     @test expmat(Pred) == [1 0 1 2; 0 1 1 0]
     @test indexvector(Pred) == [1, 2]
+
+    # conversion from Taylor model
+    x₁, x₂, x₃ = set_variables(Float64, ["x₁", "x₂", "x₃"]; order=3)
+    dom1 = IA.interval(N(-1), N(1))
+    dom = dom1 × dom1 × dom1
+    x0 = IA.IntervalBox(IA.mid.(dom)...)
+    rem = IA.interval(N(0), N(0))
+    p₁ = 33 + 2x₁ + 3x₂ + 4x₃ + 5x₁^2 + 6x₂ * x₃ + 7x₃^2 + 8x₁ * x₂ * x₃
+    p₂ = x₃ - x₁
+    vTM = [TaylorModels.TaylorModelN(pi, rem, x0, dom) for pi in [p₁, p₂]]
+    PZ = convert(SparsePolynomialZonotope, vTM)
+    # the following tests check for equality (but equivalence should be tested)
+    @test PZ.c == N[33, 0]
+    @test PZ.G == N[2 3 4 5 6 7 8;
+                    -1 0 1 0 0 0 0]
+    @test PZ.GI == Matrix{N}(undef, 2, 0)
+    @test PZ.E == [1 0 0 2 0 0 1;
+                   0 1 0 0 1 0 1;
+                   0 0 1 0 1 2 1]
+    # interestingly, the zonotope approximations are equivalent
+    Zt = overapproximate(vTM, Zonotope)
+    Zp = overapproximate(PZ, Zonotope)
+    @test isequivalent(Zt, Zp)
 end