#2522 - Concrete projection of interval x zonotope (#2525)

* add projection to cp of interval and zonotope * Update src/LazyOperations/CartesianProduct.jl Co-authored-by: Christian Schilling <[email protected]> Co-authored-by: Christian Schilling <[email protected]>
JuliaReach · Jan 29, 2021 · ede6c23 · ede6c23
1 parent 00d75fd
commit ede6c23
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 0 deletions.
diff --git a/src/LazyOperations/CartesianProduct.jl b/src/LazyOperations/CartesianProduct.jl
@@ -325,3 +325,30 @@ function project(cp::CartesianProduct{N, <:Interval, <:AbstractHyperrectangle{N}
     end
     return Hyperrectangle(cH[block_vec], rH[block_vec], check_bounds=false)
 end
+
+"""
+    project(cp::CartesianProduct{N, <:Interval, <:AbstractZonotope{N}}, block::AbstractVector{Int}) where {N}
+
+Concrete projection of a Cartesian product between an interval and a zonotopic set.
+
+### Input
+
+- `cp`       -- Cartesian product between an interval and a zonotopic set
+- `block`    -- block structure, a vector with the dimensions of interest
+
+### Output
+
+A set representing the projection of the cartesian product `cp` on the
+dimensions specified by `block`.
+"""
+function project(cp::CartesianProduct{N, <:Interval, <:AbstractZonotope{N}}, block::AbstractVector{Int}) where {N}
+    block_vec = collect(block)
+    Z = cp.Y
+    if 1 ∉ block_vec
+        block_vec .-= 1
+    else
+        Z = convert(Zonotope, cp)
+    end
+    M = projection_matrix(block_vec, dim(Z), N)
+    return linear_map(M, Z)
+end
diff --git a/test/unit_CartesianProduct.jl b/test/unit_CartesianProduct.jl
@@ -387,6 +387,10 @@ for N in [Float64, Float32]
     cp = Interval(N(0), N(2)) × Hyperrectangle(N[2, 3], N[1, 1])
     @test project(cp, 1:2) == Hyperrectangle(N[1, 2], N[1, 1])
     @test project(cp, 2:3) == Hyperrectangle(N[2, 3], N[1, 1])
+
+    cp = Interval(N(0), N(2)) × Zonotope(N[2, 3], N[1 0; 0 1])
+    @test isequivalent(project(cp, 1:2), Zonotope(N[1, 2], N[1 0; 0 1]))
+    @test isequivalent(project(cp, 2:3), Zonotope(N[2, 3], N[1 0; 0 1]))
 end
 
 for N in [Float64]