Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

#781 - Intersection emptiness check for two half-spaces #782

Merged
merged 3 commits into from
Oct 14, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/src/lib/binary_functions.md
Original file line number Diff line number Diff line change
Expand Up @@ -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}
```

Expand Down
61 changes: 61 additions & 0 deletions src/helper_functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
89 changes: 89 additions & 0 deletions src/is_intersection_empty.jl
Original file line number Diff line number Diff line change
Expand Up @@ -750,6 +750,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{HPolytope{N}, AbstractHPolygon{N}}
Expand Down
10 changes: 10 additions & 0 deletions test/unit_HalfSpace.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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