Decomposed linear map of CartesianProductArray

src/Approximations/overapproximate.jl

@@ -334,7 +334,7 @@ function overapproximate(X::LazySet{N},
 end
 
 """
-    overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real}
+    overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real}
 
 Return the overapproximation of a real unidimensional set with an interval.
 
@@ -352,11 +352,10 @@ An interval.
 The method relies on the exact conversion to `Interval`. Two support
 function evaluations are needed in general.
 """
-function overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real}
+function overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real}
     @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`"
     return convert(Interval, S)
 end
-
 function overapproximate_cap_helper(X::LazySet{N},             # compact set
                                     P::AbstractPolyhedron{N},  # polyhedron
                                     dir::AbstractDirections{N};
@@ -516,3 +515,85 @@ function overapproximate(cap::Intersection{N,
     return overapproximate(swap(cap), dir; kwargs...)
 end
 
+"""
+    overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}},
+                      ::Type{CartesianProductArray{N, O}}) where {N, T<:LazySet{N}, O<:LazySet{N}}
+
+Decompose a lazy linear map of a cartesian product array keeping original block structure.
+
+### Input
+
+- `lm`   -- lazy linear map of cartesian product array
+- `CartesianProductArray` -- type for dispatch
+
+### Output
+
+A `CartesianProductArray` representing the decomposed linear map.
+"""
+function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}},
+                           ::Type{CartesianProductArray{N, O}}) where {N, T<:LazySet{N}, O<:LazySet{N}}
+
+    M, cpa = lm.M, lm.X
+    array = Vector{O}(undef, length(cpa.array))
+    return _overapproximate_lm_cpa!(M, cpa, array, O)
+end
+
+"""
+    overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}},
+                      ::Type{CartesianProductArray{N, HPolytope{N}}},
+                   dir::AbstractDirections{N}) where {N, T<:LazySet{N}}
+
+Decompose a lazy linear map of a cartesian product array.
+
+### Input
+
+- `lm`   -- lazy linear map of cartesian product array
+- `CartesianProductArray` -- type for dispatch
+- `dir` -- approximation option for decomposition
+
+### Output
+
+A `CartesianProductArray` representing the decomposed linear map.
+"""
+function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}},
+                           ::Type{<:CartesianProductArray},
+                         dir::Type{<:AbstractDirections}) where {N, T<:LazySet{N}}
+    M, cpa = lm.M, lm.X
+    array = Vector{HPolytope{N}}(undef, length(cpa.array))
+    return _overapproximate_lm_cpa!(M, cpa, array, dir)
+end
+
+#same for <:AbstractHyperrectangle
+function overapproximate(lm::LinearMap{N, CartesianProductArray{N, T}},
+                           ::Type{<:CartesianProductArray},
+                         oa::Type{<:AbstractHyperrectangle}) where {N, T<:LazySet{N}}
+    M, cpa = lm.M, lm.X
+    array = Vector{oa{N}}(undef, length(cpa.array))
+    return _overapproximate_lm_cpa!(M, cpa, array, oa)
+end
+
+function _overapproximate_lm_cpa!(M, cpa, array, O)
+    row_start_ind, row_end_ind = 1, 0
+    @inbounds for (i, bi) in enumerate(cpa.array)
+        row_end_ind += dim(bi)
+        ms = blocks_linear_map(cpa, M, row_start_ind : row_end_ind)
+        array[i] = overapproximate(ms, O)
+        row_start_ind = row_end_ind + 1
+    end
+    return CartesianProductArray(array)
+end
+
+function blocks_linear_map(cpa::CartesianProductArray{N, T},
+                           M::AbstractMatrix{N},
+                           row_range::UnitRange{Int}) where {N, T<:LazySet{N}}
+    array = Vector{LinearMap{N, <:T}}(undef, length(cpa.array))
+
+    col_start_ind, col_end_ind = 1, 0
+    @inbounds for (i, bi) in enumerate(cpa.array)
+        col_end_ind += dim(bi)
+        array[i] = LinearMap(M[row_range, col_start_ind : col_end_ind], bi)
+        col_start_ind = col_end_ind + 1
+    end
+
+    return MinkowskiSumArray(array)
+end

test/unit_overapproximate.jl

@@ -47,12 +47,48 @@ for N in [Float64, Rational{Int}, Float32]
     h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
     H = overapproximate(h1 × h2, Hyperrectangle) # defaults to convert method
     @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]
-    
+
     # overapproximation of the lazy linear map of a hyperrectangular set
     H = Hyperrectangle(N[0, 0], N[1/2, 1])
     M = Diagonal(N[2, 2])
     OA = overapproximate(M*H, Hyperrectangle)
     @test OA isa Hyperrectangle && OA.center == N[0, 0] && OA.radius == N[1, 2]
+
+    #overapproximation of Minkowski sum of linear maps for each block in the row block
+    i1 = Interval(N[0, 1])
+    h = Hyperrectangle(low=N[3, 4], high=N[5, 7])
+    M = N[1 2 3; 4 5 6; 7 8 9]
+    cpa = CartesianProductArray([i1, h])
+    lm = M * cpa
+
+    oa = overapproximate(lm, Hyperrectangle)
+    oa_box = overapproximate(lm, Approximations.BoxDirections)
+    d_oa_d_hp = overapproximate(lm, CartesianProductArray{N, Hyperrectangle{N}})
+    d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections)
+    oa_d_hp = overapproximate(d_oa_d_hp)
+    oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections)
+
+    @test oa == oa_d_hp
+    @test oa_box == oa_d_box
+
+    for (oax, set_type) in [(d_oa_d_hp, Hyperrectangle), (d_oa_d_box, HPolytope)]
+        @test oax isa CartesianProductArray
+        arr = oax.array
+        @test length(arr) == 2 && dim(arr[1]) == 1 && dim(arr[2]) == 2
+        @test all(X -> X isa set_type, arr)
+    end
+
+    i1 = Interval(N[0, 1])
+    i2 = Interval(N[2, 3])
+    i3 = Interval(N[1, 4])
+    cpa = CartesianProductArray([i1, i2, i3])
+    M = N[1 2 0; 0 1 0; 0 1 1]
+    lm = M * cpa
+    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
+    oa = overapproximate(lm)
+    @test overapproximate(d_oa) == oa
+    @test typeof(d_oa) == CartesianProductArray{N, Interval{N}}
+
 end
 
 # tests that do not work with Rational{Int}
@@ -130,3 +166,15 @@ for N in [Float64, Float32]
     @test Y_polygon ⊆ Y_zonotope
     @test !(Y_zonotope ⊆ Y_polygon)
 end
+
+for N in [Float64] # due to sparse vectors: a = sparse(Float32[1 -1; 1 1];); a \ Float32[4, 10]
+    #decomposed linear map approximation
+    i1 = Interval(N[0, 1])
+    i2 = Interval(N[2, 3])
+    M = N[1 2; 0 1]
+    cpa = CartesianProductArray([i1, i2])
+    lm = M * cpa
+    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
+    oa = overapproximate(lm, OctDirections)
+    @test oa ⊆ d_oa
+end