Merge pull request #2241 from JuliaReach/schillic/1227

#1227 - Efficient subset check between CartesianProducts

docs/src/lib/binary_functions.md

@@ -122,6 +122,7 @@ pontryagin_difference
 ⊆(::LazySet{N}, ::Universe{N}, ::Bool=false) where {N<:Real}
 ⊆(::Universe{N}, ::LazySet{N}, ::Bool=false) where {N<:Real}
 ⊆(::LazySet{N}, ::Complement{N}, ::Bool=false) where {N<:Real}
+⊆(::CartesianProduct{N}, ::CartesianProduct{N}, ::Bool=false) where {N<:Real}
 ⊆(::CartesianProductArray{N}, ::CartesianProductArray{N}, ::Bool=false) where {N<:Real}
 ```
 

src/ConcreteOperations/issubset.jl

@@ -957,6 +957,77 @@ end
 # --- CartesianProduct ---
 
 
+"""
+    ⊆(X::CartesianProduct{N}, Y::CartesianProduct{N}, witness::Bool=false;
+      check_block_equality::Bool=true) where {N<:Real}
+
+Check whether a Cartesian product of two convex sets is contained in another
+Cartesian product of two convex sets, and otherwise optionally compute a
+witness.
+
+### Input
+
+- `X`       -- Cartesian product of two convex sets
+- `Y`       -- Cartesian product of two convex sets
+- `witness` -- (optional, default: `false`) compute a witness if activated
+- `check_block_equality` -- (optional, default: `true`) flag for checking that
+               the block structure of the two sets is identical
+
+### Output
+
+* If `witness` option is deactivated: `true` iff ``X ⊆ Y``
+* If `witness` option is activated:
+  * `(true, [])` iff ``X ⊆ Y``
+  * `(false, v)` iff ``X \\not\\subseteq Y`` and
+    ``v ∈ X \\setminus Y``
+
+### Notes
+
+This algorithm requires that the two Cartesian products share the same block
+structure.
+If `check_block_equality` is activated, we check this property and, if it does
+not hold, we use a fallback implementation based on conversion to constraint
+representation (assuming that the sets are polyhedral).
+
+### Algorithm
+
+We check for inclusion for each block of the Cartesian products.
+
+For witness production, we obtain a witness in one of the blocks.
+We then construct a high-dimensional witness by obtaining any point in the other
+blocks (using `an_element`) and concatenating these points.
+"""
+function ⊆(X::CartesianProduct{N}, Y::CartesianProduct{N}, witness::Bool=false;
+           check_block_equality::Bool=true) where {N<:Real}
+    n1 = dim(X.X)
+    n2 = dim(X.Y)
+    if check_block_equality && (n1 != dim(Y.X) || n2 != dim(Y.Y))
+        return invoke(⊆, Tuple{LazySet{N}, LazySet{N}, Bool}, X, Y, witness)
+    end
+
+    # check first block
+    result = ⊆(X.X, Y.X, witness)
+    if !witness && !result
+        return false
+    elseif witness && !result[1]
+        # construct a witness
+        w = vcat(result[2], an_element(X.Y))
+        return (false, w)
+    end
+
+    # check second block
+    result = ⊆(X.Y, Y.Y, witness)
+    if !witness && !result
+        return false
+    elseif witness && !result[1]
+        # construct a witness
+        w = vcat(an_element(X.X), result[2])
+        return (false, w)
+    end
+
+    return witness ? (true, N[]) : true
+end
+
 """
     ⊆(X::CartesianProductArray{N}, Y::CartesianProductArray{N},
       witness::Bool=false; check_block_equality::Bool=true
@@ -986,7 +1057,9 @@ compute a witness.
 
 This algorithm requires that the two Cartesian products share the same block
 structure.
-Depending on the value of `check_block_equality`, we check this property.
+If `check_block_equality` is activated, we check this property and, if it does
+not hold, we use a fallback implementation based on conversion to constraint
+representation (assuming that the sets are polyhedral).
 
 ### Algorithm
 
@@ -1004,8 +1077,7 @@ function ⊆(X::CartesianProductArray{N},
     aX = array(X)
     aY = array(Y)
     if check_block_equality && !same_block_structure(aX, aY)
-        throw(ArgumentError("this inclusion check requires Cartesian products" *
-                            "with the same block structure"))
+        return invoke(⊆, Tuple{LazySet{N}, LazySet{N}, Bool}, X, Y, witness)
     end
 
     for i in 1:length(aX)

test/unit_CartesianProduct.jl

@@ -139,6 +139,22 @@ for N in [Float64, Float32, Rational{Int}]
         @test concretize(cp) === cp
     end
 
+    # inclusion
+    # same block structure
+    cp1 = CartesianProduct(Ball1(N[2], N(1)), Ball1(N[0, 0], N(2)))
+    cp2 = CartesianProduct(Interval(N(1), N(3)), BallInf(N[0, 0], N(3)))
+    res, w = ⊆(cp1, cp2, true)
+    @test cp1 ⊆ cp2 && res && w == N[]
+    res, w = ⊆(cp2, cp1, true)
+    @test cp2 ⊈ cp1 && !res && w ∈ cp2 && w ∉ cp1
+    # different block structure
+    cp1 = CartesianProduct(Interval(N(1), N(2)), BallInf(N[1, 1], N(1)))
+    cp2 = CartesianProduct(BallInf(N[1, 1], N(2)), Interval(N(0), N(2)))
+    res, w = ⊆(cp1, cp2, true)
+    @test cp1 ⊆ cp2 && res && w == N[]
+    res, w = ⊆(cp2, cp1, true)
+    @test cp2 ⊈ cp1 && !res && w ∈ cp2 && w ∉ cp1
+
     # ==================================
     # Conversions of Cartesian Products
     # ==================================
@@ -301,6 +317,13 @@ for N in [Float64, Float32, Rational{Int}]
     @test cpa1 ⊆ cpa2 && res && w == N[]
     res, w = ⊆(cpa2, cpa1, true)
     @test cpa2 ⊈ cpa1 && !res && w ∈ cpa2 && w ∉ cpa1
+    # different block structure
+    cpa1 = CartesianProductArray([Interval(N(1), N(2)), BallInf(N[1, 1], N(1))])
+    cpa2 = CartesianProductArray([BallInf(N[1, 1], N(2)), Interval(N(0), N(2))])
+    res, w = ⊆(cpa1, cpa2, true)
+    @test cpa1 ⊆ cpa2 && res && w == N[]
+    res, w = ⊆(cpa2, cpa1, true)
+    @test cpa2 ⊈ cpa1 && !res && w ∈ cpa2 && w ∉ cpa1
 
     # convert a hyperrectangle to the cartesian product array of intervals
     # convert a hyperrectangle to a cartesian product of intervals