Merge pull request #3269 from JuliaReach/schillic/1721

#1721 - Concrete minkowski difference for hyperrectangular sets

docs/src/lib/binary_functions.md

@@ -120,6 +120,8 @@ minkowski_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope)
 pontryagin_difference
 minkowski_difference(::LazySet, ::LazySet)
 minkowski_difference(::AbstractZonotope, ::AbstractZonotope)
+minkowski_difference(::AbstractHyperrectangle, ::AbstractHyperrectangle)
+minkowski_difference(::Interval, ::Interval)
 ```
 
 ## Subset check

src/ConcreteOperations/minkowski_difference.jl

@@ -76,7 +76,7 @@ Mathematically:
 """
 const pontryagin_difference = minkowski_difference
 
-for ST in [:LazySet, :AbstractZonotope]
+for ST in [:LazySet, :AbstractZonotope, :AbstractHyperrectangle]
     # Minkowski difference with singleton is a translation
     @eval minkowski_difference(X::($ST), S::AbstractSingleton) =
         translate(X, -element(S))
@@ -85,6 +85,62 @@ for ST in [:LazySet, :AbstractZonotope]
     @eval minkowski_difference(X::($ST), ::ZeroSet) = X
 end
 
+"""
+    minkowski_difference(I1::Interval, I2::Interval)
+
+Compute the Minkowski difference of two intervals.
+
+### Input
+
+- `I1` -- interval
+- `I2` -- interval
+
+### Output
+
+An `Interval` that corresponds to the Minkowski difference of `I1` minus `I2`,
+or an `EmptySet` if the difference is empty.
+"""
+function minkowski_difference(I1::Interval, I2::Interval)
+    l = min(I1) - min(I2)
+    h = max(I1) - max(I2)
+    if h < l
+        N = promote_type(eltype(I1), eltype(I2))
+        return EmptySet{N}(1)
+    end
+    return Interval(l, h)
+end
+
+"""
+    minkowski_difference(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle)
+
+Compute the Minkowski difference of two hyperrectangular sets.
+
+### Input
+
+- `H1` -- hyperrectangular set
+- `H2` -- hyperrectangular set
+
+### Output
+
+A `Hyperrectangle` that corresponds to the Minkowski difference of `H1` minus
+`H2`, or an `EmptySet` if the difference is empty.
+"""
+function minkowski_difference(H1::AbstractHyperrectangle,
+                              H2::AbstractHyperrectangle)
+    n = dim(H1)
+    @assert n == dim(H2) "incompatible dimensions"
+
+    N = promote_type(eltype(H1), eltype(H2))
+    r = Vector{N}(undef, n)
+    for i in 1:n
+        r[i] = radius_hyperrectangle(H1, i) - radius_hyperrectangle(H2, i)
+        if r[i] < zero(N)
+            return EmptySet{N}(n)
+        end
+    end
+    return Hyperrectangle(center(H1) - center(H2), r)
+end
+
 """
     minkowski_difference(Z1::AbstractZonotope, Z2::AbstractZonotope)
 

test/ConcreteOperations/minkowski_difference.jl

@@ -0,0 +1,23 @@
+for N in [Float64, Float32, Rational{Int}]
+    # intervals, nonempty difference
+    X = Interval(N(1), N(3))
+    Y = Interval(N(2), N(3))
+    @test minkowski_difference(X, Y) == Interval(N(-1), N(0))
+
+    # intervals, empty difference
+    X = Interval(N(1), N(3))
+    Y = Interval(N(2), N(5))
+    D = minkowski_difference(X, Y)
+    @test D isa EmptySet{N} && dim(D) == 1
+
+    # hyperrectangles, nonempty difference
+    H1 = Hyperrectangle(N[2, 3], N[1, 3])
+    H2 = Hyperrectangle(N[5/2, -1], N[1/2, 2])
+    @test minkowski_difference(H1, H2) == Hyperrectangle(N[-1/2, 4], N[1/2, 1])
+
+    # hyperrectangles, empty difference
+    H1 = Hyperrectangle(N[2, 3], N[1, 3])
+    H2 = Hyperrectangle(N[5, -1], N[2, 2])
+    D = minkowski_difference(H1, H2)
+    @test D isa EmptySet{N} && dim(D) == 2
+end

test/runtests.jl

@@ -138,6 +138,7 @@ if test_suite_basic
     @time @testset "LazySets.isdisjoint" begin include("ConcreteOperations/isdisjoint.jl") end
     @time @testset "LazySets.distance" begin include("ConcreteOperations/distance.jl") end
     @time @testset "LazySets.isstrictsubset" begin include("ConcreteOperations/isstrictsubset.jl") end
+    @time @testset "LazySets.minkowski_difference" begin include("ConcreteOperations/minkowski_difference.jl") end
     @time @testset "LazySets.minkowski_sum" begin include("ConcreteOperations/minkowski_sum.jl") end
     @time @testset "LazySets.samples" begin include("ConcreteOperations/samples.jl") end