Merge pull request #782 from JuliaReach/schillic/781

#781 - Intersection emptiness check for two half-spaces

docs/src/lib/binary_functions.md

@@ -49,6 +49,7 @@ is_intersection_empty(::AbstractPolytope{N}, ::AbstractPolytope{N}, ::Bool=false
 is_intersection_empty(::AbstractSingleton{N}, ::AbstractPolytope{N}, ::Bool=false) where {N<:Real}
 is_intersection_empty(::LazySet{N}, ::Union{Hyperplane{N}, Line{N}}, ::Bool=false) where {N<:Real}
 is_intersection_empty(::LazySet{N}, ::HalfSpace{N}, ::Bool=false) where {N<:Real}
+is_intersection_empty(::HalfSpace{N}, ::HalfSpace{N}, ::Bool=false) where {N<:Real}
 is_intersection_empty(::LazySet{N}, ::Union{HPolytope{N}, AbstractHPolygon{N}}) where {N<:Real}
 ```
 

src/helper_functions.jl

@@ -86,6 +86,67 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T
     return true
 end
 
+"""
+    samedir(u::AbstractVector{N},
+            v::AbstractVector{N})::Tuple{Bool, Real} where N<:Real
+
+Check whether two vectors point in the same direction.
+
+### Input
+
+- `u` -- first vector
+- `v` -- second vector
+
+### Output
+
+`true` iff the vectors are identical up to a positive scaling factor.
+
+### Examples
+
+```jldoctest
+julia> LazySets.samedir([1, 2, 3], [2, 4, 6])
+true, 0.5
+
+julia> LazySets.samedir([1, 2, 3], [3, 2, 1])
+false
+
+julia> LazySets.samedir([1, 2, 3], [-1, -2, -3])
+false
+
+```
+"""
+function samedir(u::AbstractVector{N},
+                 v::AbstractVector{N}
+                )::Tuple{Bool, Real} where N<:Real
+    @assert length(u) == length(v) "wrong dimension"
+    no_factor = true
+    factor = 0
+    @inbounds for i in 1:length(u)
+        if u[i] == 0
+            if v[i] != 0
+                return (false, 0)
+            end
+            continue
+        elseif v[i] == 0
+            return (false, 0)
+        end
+        if no_factor
+            no_factor = false
+            factor = u[i] / v[i]
+            if factor < 0
+                return (false, 0)
+            end
+        elseif factor != u[i] / v[i]
+            return (false, 0)
+        end
+    end
+    if no_factor
+        # both vectors are zero
+        return (true, 0)
+    end
+    return (true, factor)
+end
+
 """
     nonzero_indices(v::AbstractVector{N})::Vector{Int} where N<:Real
 

src/is_intersection_empty.jl

@@ -748,6 +748,95 @@ function is_intersection_empty(hs::HalfSpace{N},
     return is_intersection_empty(X, hs, witness; kwargs...)
 end
 
+"""
+    is_intersection_empty(hs1::HalfSpace{N},
+                          hs2::HalfSpace{N},
+                          [witness]::Bool=false
+                         )::Union{Bool, Tuple{Bool,Vector{N}}} where N<:Real
+
+Check whether two half-spaces do not intersect, and otherwise optionally compute
+a witness.
+
+### Input
+
+- `hs1`     -- half-space
+- `hs2`     -- half-space
+- `witness` -- (optional, default: `false`) compute a witness if activated
+
+### Output
+
+* If `witness` option is deactivated: `true` iff ``hs1 ∩ hs2 = ∅``
+* If `witness` option is activated:
+  * `(true, [])` iff ``hs1 ∩ hs2 = ∅``
+  * `(false, v)` iff ``hs1 ∩ hs2 ≠ ∅`` and ``v ∈ hs1 ∩ hs2``
+
+### Algorithm
+
+Two half-spaces do not intersect if and only if their normal vectors point in
+the opposite direction and there is a gap between the two defining hyperplanes.
+
+The latter can be checked as follows:
+Let ``hs_1 : a_1⋅x = b_1`` and ``hs2 : a_2⋅x = b_2``.
+Then we already know that ``a_2 = -k⋅a_1`` for some positive scaling factor
+``k``.
+Let ``x_1`` be a point on the defining hyperplane of ``hs_1``.
+We construct a line segment from ``x_1`` to the point ``x_2`` on the defining
+hyperplane of ``hs_2`` by shooting a ray from ``x_1`` with direction ``a_1``.
+Thus we look for a factor ``s`` such that ``(x_1 + s⋅a_1)⋅a_2 = b_2``.
+This gives us ``s = (b_2 - x_1⋅a_2) / (-k a_1⋅a_1)``.
+The gap exists if and only if ``s`` is positive.
+
+If the normal vectors do not point in opposite directions, then the defining
+hyperplanes intersect and we can produce a witness as follows.
+All points ``x`` in this intersection satisfy ``a_1⋅x = b_1`` and ``a_2⋅x = b_2``.
+Thus we have ``(a_1 + a_2)⋅x = b_1+b_2``.
+We now find a dimension where ``a_1 + a_2`` is non-zero, say, ``i``.
+Then the result is a vector with one non-zero entry in dimension ``i``, defined
+as ``[0, …, 0, (b_1 + b_2)/(a_1[i] + a_2[i]), 0, …, 0]``.
+Such a dimension ``i`` always exists.
+"""
+function is_intersection_empty(hs1::HalfSpace{N},
+                               hs2::HalfSpace{N},
+                               witness::Bool=false;
+                               kwargs...
+                              )::Union{Bool, Tuple{Bool,Vector{N}}} where N<:Real
+    a1 = hs1.a
+    a2 = hs2.a
+    issamedir, k = samedir(a1, -a2)
+    if issamedir
+        x1 = an_element(Hyperplane(a1, hs1.b))
+        b2 = hs2.b
+        s = (b2 - dot(x1, a2)) / (-k * dot(a1, a1))
+        empty_intersection = s > 0
+        if witness
+            if empty_intersection
+                v = N[]
+            else
+                # both defining hyperplanes are contained in each half-space
+                v = x1
+            end
+        end
+    else
+        empty_intersection = false
+        if witness
+            v = zeros(N, length(a1))
+            for i in 1:length(a1)
+                a_sum_i = a1[i] + a2[i]
+                if a_sum[i] != 0
+                    v[i] = (hs1.b + hs2.b) / a_sum_i
+                    break
+                end
+            end
+        end
+    end
+
+    if !witness
+        return empty_intersection
+    else
+        return (empty_intersection, v)
+    end
+end
+
 """
     is_intersection_empty(X::LazySet{N},
                           P::Union{HPolyhedron{N}, HPolytope{N}, AbstractHPolygon{N}}

test/unit_HalfSpace.jl

@@ -48,4 +48,14 @@ for N in [Float64, Rational{Int}, Float32]
     b = BallInf(N[1, 1, 1], N(1))
     empty_intersection, v = is_intersection_empty(b, hs, true)
     @test !is_intersection_empty(b, hs) && !empty_intersection && v ∈ hs
+    hs1 = HalfSpace(N[1, 0], N(1)) # x <= 1
+    hs2 = HalfSpace(N[-1, 0], N(-2)) # x >= 2
+    empty_intersection, v = is_intersection_empty(hs1, hs2, true)
+    @test is_intersection_empty(hs1, hs2) && empty_intersection && v == N[]
+    hs3 = HalfSpace(N[-1, 0], N(-1)) # x >= 1
+    empty_intersection, v = is_intersection_empty(hs1, hs3, true)
+    @test !is_intersection_empty(hs1, hs3) && !empty_intersection && v == N[1, 0]
+    hs4 = HalfSpace(N[-1, 0], N(0)) # x >= 0
+    empty_intersection, v = is_intersection_empty(hs1, hs4, true)
+    @test !is_intersection_empty(hs1, hs4) && !empty_intersection && v ∈ hs1 && v ∈ hs4
 end