Merge pull request #3447 from JuliaReach/schillic/convert2

`convert` from SparsePolynomialZonotope to Taylor model

docs/src/lib/sets/SparsePolynomialZonotope.md

@@ -16,6 +16,7 @@ uniqueID(::Int)
 ngens_dep(::SparsePolynomialZonotope)
 ngens_indep(::SparsePolynomialZonotope)
 nparams(::SparsePolynomialZonotope)
+polynomial_order(::SparsePolynomialZonotope)
 order(::SparsePolynomialZonotope)
 linear_map(::AbstractMatrix, ::SparsePolynomialZonotope)
 translate(::SparsePolynomialZonotope, ::AbstractVector)

src/Initialization/init_IntervalArithmetic.jl

@@ -11,6 +11,12 @@ function vertices_list(H::IA.IntervalBox)
     return vertices_list(convert(Hyperrectangle, H))
 end
 
+const zero_itv(N) = IA.interval(zero(N))
+const sym_itv(N) = IA.interval(-one(N), one(N))
+
+zero_box(n::Int, N=Float64) = IA.IntervalBox(zero_itv(N), n)
+sym_box(n::Int, N=Float64) = IA.IntervalBox(sym_itv(N), n)
+
 """
     fast_interval_pow(a::IntervalArithmetic.Interval, n::Int)
 

src/Initialization/init_TaylorModels.jl

@@ -0,0 +1,11 @@
+using .TaylorModels: domain
+using .TaylorModels.TaylorSeries: numtype
+
+_eltype_TM(TMi) = numtype(polynomial(TMi))
+
+# check that a vector of Taylor models has the [-1, 1] domain
+function _has_normalized_domain(vTM::Vector)
+    return all(TMi -> all(==(sym_itv(_eltype_TM(TMi))), domain(TMi)), vTM)
+end
+
+eval(load_taylormodels_convert_polynomial_zonotope())

src/LazySets.jl

@@ -3,7 +3,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
@@ -12,13 +11,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
-# ==============
+# ================
 
 import ReachabilityBase.Assertions
 using ReachabilityBase.Assertions: @assert

src/Sets/SparsePolynomialZonotope.jl

@@ -1,5 +1,5 @@
 export SparsePolynomialZonotope, expmat, nparams, ngens_dep, ngens_indep,
-       genmat_dep, genmat_indep, indexvector, translate,
+       genmat_dep, genmat_indep, indexvector, polynomial_order, translate,
        linear_map, quadratic_map, remove_redundant_generators, reduce_order, ρ
 
 """
@@ -233,6 +233,27 @@ The index vector contains positive integers for the dependent parameters.
 """
 indexvector(P::SPZ) = P.idx
 
+"""
+    polynomial_order(P::SPZ)
+
+Return the polynomial order of a sparse polynomial zonotope.
+
+### Input
+
+- `P` -- sparse polynomial zonotope
+
+### Output
+
+The polynomial order.
+
+### Notes
+
+The polynomial order is the maximum sum of all monomials' parameter exponents.
+"""
+function polynomial_order(P::SPZ)
+    return maximum(sum, eachcol(expmat(P)))
+end
+
 """
     uniqueID(n::Int)
 

src/convert.jl

@@ -1162,3 +1162,117 @@ 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, set_variables, HomogeneousPolynomial,
+                                          in_base, index_table, pos_table, TaylorN
+
+        # 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
+
+        # implements Proposition 3.1.13 in thesis
+        function convert(::Type{Vector{<:TaylorModelN}}, P::SparsePolynomialZonotope{N}) where {N}
+            n = dim(P)
+            p = nparams(P)  # number of parameters
+            q = ngens_indep(P)  # independent/linear parameters
+            h = ngens_dep(P)  # dependent/nonlinear parameters
+            r = p + q
+            c = center(P)
+            G = genmat_dep(P)
+            GI = genmat_indep(P)
+            E = expmat(P)
+            poly_order = polynomial_order(P)
+            z = zeros(Int, q)
+            # we need to rewrite the global variables
+            set_variables("x"; order=poly_order, numvars=r)
+            # the following vectors are shared for each polynomial
+            rem = zero_itv(N)
+            dom = sym_box(r, N)
+            x0 = zero_box(r, N)
+
+            vTM = Vector{TaylorModelN{r,N,N}}(undef, n)
+            @inbounds for i in 1:n
+                coeffs = zeros(HomogeneousPolynomial{N}, poly_order)
+
+                # constant term
+                coeffs[1] = c[i]
+
+                # linear terms
+                for j in 1:q
+                    v = zeros(N, r)
+                    v[p + j] = GI[i, j]
+                    coeffs[2] += HomogeneousPolynomial(v, 1)
+                end
+
+                # nonlinear terms
+                # TODO can be faster if loop is restructured so that index j is
+                # created for each dimension
+                for j in 1:h
+                    Ej = E[:, j]
+                    ord = sum(Ej)
+                    G[i, j]
+                    idx = pos_table[ord + 1][in_base(poly_order, Ej)]
+                    l = length(coeff_table[ord + 1])
+                    v = Vector(SingleEntryVector(idx, l, G[i, j]))
+                    pj = HomogeneousPolynomial(v, ord)
+                    coeffs[ord + 1] += pj
+                end
+
+                pol = TaylorN(coeffs, poly_order)
+                vTM[i] = TaylorModelN(pol, rem, x0, dom)
+            end
+            return vTM
+        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

@@ -105,13 +105,42 @@ for N in [Float64, Float32, Rational{Int}]
     end
 end
 
-SSPZ = SimpleSparsePolynomialZonotope([0.2, -0.6], [1 0; 0 0.4], [1 0; 0 1])
-SPZ = convert(SparsePolynomialZonotope, SSPZ)
-@test center(SPZ) == center(SSPZ)
-@test genmat_dep(SPZ) == genmat(SSPZ)
-@test expmat(SPZ) == expmat(SSPZ)
-@test isempty(genmat_indep(SPZ))
-@test indexvector(SPZ) == 1:2
+for N in [Float64]
+    SSPZ = SimpleSparsePolynomialZonotope(N[0.2, -0.6], N[1 0; 0 0.4], [1 0; 0 1])
+    SPZ = convert(SparsePolynomialZonotope, SSPZ)
+    @test center(SPZ) == center(SSPZ)
+    @test genmat_dep(SPZ) == genmat(SSPZ)
+    @test expmat(SPZ) == expmat(SSPZ)
+    @test isempty(genmat_indep(SPZ))
+    @test indexvector(SPZ) == 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)
+
+    # conversion back to Taylor model
+    vTM2 = convert(Vector{<:TaylorModelN}, PZ)
+    @test vTM == vTM2
+end
 
 for Z in [rand(Zonotope), rand(Hyperrectangle)]
     ZS = convert(SparsePolynomialZonotope, Z)