Decomposed linear map of CartesianProductArray

src/Approximations/overapproximate.jl

@@ -516,3 +516,67 @@ function overapproximate(cap::Intersection{N,
     return overapproximate(swap(cap), dir; kwargs...)
 end
 
+"""
+    overapproximate(lm::LinearMap{N, <:CartesianProductArray},
+                      ::Type{<:CartesianProductArray}, oa=Hyperrectangle)  where {N}
+
+Decompose a lazy linear map of a cartesian product array.
+
+### Input
+
+- `lm`   -- lazy linear map of cartesian product array
+- `CartesianProductArray` -- type for dispatch
+- `oa` -- approximation option for decomposition
+
+### Output
+
+A `CartesianProductArray` representing the decomposed linear map.
+"""
+function overapproximate(lm::LinearMap{N, <:CartesianProductArray},
+                           ::Type{<:CartesianProductArray}, oa=Hyperrectangle) where {N}
+
+    cpa = lm.X
+    M = lm.M
+
+    array = allocate_result(oa, N)
+    sizehint!(array, length(cpa.array))
+
+    row_start_ind, row_end_ind = 1, 0
+    for bi in cpa.array
+        row_end_ind += dim(bi)
+        ms = blocks_linear_map(cpa, M, row_start_ind : row_end_ind)
+        push!(array, overapproximate(ms, oa))
+        row_start_ind += row_end_ind
+    end
+
+    result = CartesianProductArray(array)
+    return result
+end
+
+
+function blocks_linear_map(cpa::CartesianProductArray,
+                           M::AbstractMatrix{N},
+                           row_range::UnitRange{Int}) where {N}
+    col_start_ind, col_end_ind = 1, 0
+    array = Vector{LazySet{N}}()
+    sizehint!(array, length(cpa.array))
+    for bi in cpa.array
+        col_end_ind += dim(bi)
+        push!(array, LinearMap(M[row_range, col_start_ind : col_end_ind], bi))
+        col_start_ind += col_end_ind
+    end
+
+    return MinkowskiSumArray(array)
+end
+
+function allocate_result(oa, N)
+    return Vector{LazySet{N}}()
+end
+
+function allocate_result(oa::Type{<:LazySet}, N)
+    return Vector{oa{N}}()
+end
+
+function allocate_result(oa::Type{<:AbstractDirections}, N)
+    return Vector{HPolytope{N}}()
+end

test/unit_overapproximate.jl

@@ -47,12 +47,36 @@ 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)
+    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
 end
 
 # tests that do not work with Rational{Int}
@@ -129,4 +153,5 @@ for N in [Float64, Float32]
     Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope
     @test Y_polygon ⊆ Y_zonotope
     @test !(Y_zonotope ⊆ Y_polygon)
+
 end