added new tests and refactoring

added functionality to keep the original structure and preallocate corresponding  set

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}, ::Union{Type{Interval}, Type{Interval{N}}}) 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}, ::Union{Type{Interval}, Type{Interval{N}}}) 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};
@@ -517,66 +516,84 @@ function overapproximate(cap::Intersection{N,
 end
 
 """
-    overapproximate(lm::LinearMap{N, <:CartesianProductArray},
-                      ::Type{<:CartesianProductArray}, oa=Hyperrectangle)  where {N}
+    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.
+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
-- `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}
+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
 
-    cpa = lm.X
-    M = lm.M
+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
 
-    array = allocate_result(oa, N)
-    sizehint!(array, length(cpa.array))
+#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
-    for bi in cpa.array
+    @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)
-        push!(array, overapproximate(ms, oa))
-        row_start_ind += row_end_ind
+        array[i] = overapproximate(ms, O)
+        row_start_ind = row_end_ind + 1
     end
-
-    result = CartesianProductArray(array)
-    return result
+    return CartesianProductArray(array)
 end
 
-
-function blocks_linear_map(cpa::CartesianProductArray,
+function blocks_linear_map(cpa::CartesianProductArray{N, T},
                            M::AbstractMatrix{N},
-                           row_range::UnitRange{Int}) where {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
-    array = Vector{LazySet{N}}()
-    sizehint!(array, length(cpa.array))
-    for bi in cpa.array
+    @inbounds for (i, bi) in enumerate(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
+        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
-
-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

@@ -63,7 +63,7 @@ for N in [Float64, Rational{Int}, Float32]
 
     oa = overapproximate(lm, Hyperrectangle)
     oa_box = overapproximate(lm, Approximations.BoxDirections)
-    d_oa_d_hp = overapproximate(lm, CartesianProductArray)
+    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)
@@ -77,6 +77,26 @@ for N in [Float64, Rational{Int}, Float32]
         @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])
+    M = N[1 2; 0 1]
+    cpa = CartesianProductArray([i1, i2])
+    lm = M * cpa
+    d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}})
+    oa = overapproximate(lm, 1e-8)
+    @test oa ⊆ d_oa
+
+    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}