#2521 - Concrete projection a hyperrectangular cartesian product (#2524)

* add projection to cp of interval and box * Update test/unit_CartesianProduct.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/LazyOperations/CartesianProduct.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/LazyOperations/CartesianProduct.jl Co-authored-by: Christian Schilling <[email protected]> * update * fix error Co-authored-by: Christian Schilling <[email protected]>
JuliaReach · Jan 28, 2021 · 00d75fd · 00d75fd
1 parent 8ab5b11
commit 00d75fd
Show file tree

Hide file tree

Showing 2 changed files with 35 additions and 0 deletions.
diff --git a/src/LazyOperations/CartesianProduct.jl b/src/LazyOperations/CartesianProduct.jl
@@ -294,3 +294,34 @@ end
 function concretize(cp::CartesianProduct)
     return cartesian_product(concretize(cp.X), concretize(cp.Y))
 end
+
+"""
+    project(cp::CartesianProduct{N, <:Interval, <:AbstractHyperrectangle{N}}, block::AbstractVector{Int}) where {N}
+
+Concrete projection of a Cartesian product between an interval and a
+hyperrectangle.
+
+### Input
+
+- `cp`       -- Cartesian product between an interval and a hyperrectangle
+- `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, <:AbstractHyperrectangle{N}}, block::AbstractVector{Int}) where {N}
+    I = cp.X
+    H = cp.Y
+    block_vec = collect(block)
+    if 1 ∉ block_vec
+        block_vec .-= 1
+        cH = center(H)
+        rH = radius_hyperrectangle(H)
+    else
+        cH = vcat(center(I), center(H))
+        rH = vcat(radius_hyperrectangle(I), radius_hyperrectangle(H))
+    end
+    return Hyperrectangle(cH[block_vec], rH[block_vec], check_bounds=false)
+end
diff --git a/test/unit_CartesianProduct.jl b/test/unit_CartesianProduct.jl
@@ -383,6 +383,10 @@ for N in [Float64, Float32]
     @test_throws AssertionError is_intersection_empty(cpa1, Universe{N}(3))
     @test_throws AssertionError is_intersection_empty(Universe{N}(5), cpa2)
 
+    # projection
+    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])
 end
 
 for N in [Float64]