Set operations and setdiff (#143)

* Factor out set operations into separate file

* Add setdiff for IntervalBox in 2D

* Fixup

* Fix versions for REQUIRE and test/REQUIRE

REQUIRE

@@ -1,6 +1,5 @@
 julia 0.4
 CRlibm 0.2.2
 Compat 0.7.11
-FixedSizeArrays 0.1
+FixedSizeArrays 0.1.2
 ForwardDiff 0.1.8
-

src/intervals/arithmetic.jl

@@ -1,17 +1,5 @@
 # This file is part of the ValidatedNumerics.jl package; MIT licensed
 
-"""
-    in(x, a)
-    ∈(x, a)
-
-Checks if the number `x` is a member of the interval `a`, treated as a set.
-Corresponds to `isMember` in the ITF-1788 Standard.
-"""
-function in{T<:Real}(x::T, a::Interval)
-    isinf(x) && return false
-    a.lo <= x <= a.hi
-end
-
 
 ## Comparisons
 
@@ -26,36 +14,13 @@ function ==(a::Interval, b::Interval)
 end
 !=(a::Interval, b::Interval) = !(a==b)
 
-"""
-    issubset(a,b)
-    ⊆(a,b)
-
-Checks if all the points of the interval `a` are within the interval `b`.
-"""
-function ⊆(a::Interval, b::Interval)
-    isempty(a) && return true
-    b.lo ≤ a.lo && a.hi ≤ b.hi
-end
 
 # Auxiliary functions: equivalent to </<=, but Inf <,<= Inf returning true
 function islessprime{T<:Real}(a::T, b::T)
     (isinf(a) || isinf(b)) && a==b && return true
     a < b
 end
 
-# Interior
-function interior(a::Interval, b::Interval)
-    isempty(a) && return true
-    islessprime(b.lo, a.lo) && islessprime(a.hi, b.hi)
-end
-const ⪽ = interior  # \subsetdot
-
-# Disjoint:
-function isdisjoint(a::Interval, b::Interval)
-    (isempty(a) || isempty(b)) && return true
-    islessprime(b.hi, a.lo) || islessprime(a.hi, b.lo)
-end
-
 # Weakly less, \le, <=
 function <=(a::Interval, b::Interval)
     isempty(a) && isempty(b) && return true
@@ -267,43 +232,6 @@ function max(a::Interval, b::Interval)
 end
 
 
-## Set operations
-"""
-    intersect(a, b)
-    ∩(a,b)
-
-Returns the intersection of the intervals `a` and `b`, considered as
-(extended) sets of real numbers. That is, the set that contains
-the points common in `a` and `b`.
-"""
-function intersect{T}(a::Interval{T}, b::Interval{T})
-    isdisjoint(a,b) && return emptyinterval(T)
-
-    Interval(max(a.lo, b.lo), min(a.hi, b.hi))
-end
-# Specific promotion rule for intersect:
-intersect{T,S}(a::Interval{T}, b::Interval{S}) = intersect(promote(a,b)...)
-
-
-## Hull
-"""
-    hull(a, b)
-
-Returns the "convex hull" of the intervals `a` and `b`, considered as
-(extended) sets of real numbers. That is, the minimum set that contains
-all points in `a` and `b`.
-"""
-hull{T}(a::Interval{T}, b::Interval{T}) = Interval(min(a.lo, b.lo), max(a.hi, b.hi))
-
-"""
-    union(a, b)
-    ∪(a,b)
-
-Returns the union (convex hull) of the intervals `a` and `b`; it is equivalent
-to `hull(a,b)`.
-"""
-union(a::Interval, b::Interval) = hull(a, b)
-
 
 dist(a::Interval, b::Interval) = max(abs(a.lo-b.lo), abs(a.hi-b.hi))
 eps(a::Interval) = max(eps(a.lo), eps(a.hi))

src/intervals/intervals.jl

@@ -44,6 +44,7 @@ include("special.jl")
 include("macros.jl")
 include("conversion.jl")
 include("precision.jl")
+include("set_operations.jl")
 include("arithmetic.jl")
 include("functions.jl")
 include("trigonometric.jl")

src/intervals/set_operations.jl

@@ -0,0 +1,109 @@
+# This file is part of the ValidatedNumerics.jl package; MIT licensed
+
+
+"""
+    in(x, a)
+    ∈(x, a)
+
+Checks if the number `x` is a member of the interval `a`, treated as a set.
+Corresponds to `isMember` in the ITF-1788 Standard.
+"""
+function in{T<:Real}(x::T, a::Interval)
+    isinf(x) && return false
+    a.lo <= x <= a.hi
+end
+
+
+
+"""
+    issubset(a,b)
+    ⊆(a,b)
+
+Checks if all the points of the interval `a` are within the interval `b`.
+"""
+function ⊆(a::Interval, b::Interval)
+    isempty(a) && return true
+    b.lo ≤ a.lo && a.hi ≤ b.hi
+end
+
+
+# Interior
+function interior(a::Interval, b::Interval)
+    isempty(a) && return true
+    islessprime(b.lo, a.lo) && islessprime(a.hi, b.hi)
+end
+const ⪽ = interior  # \subsetdot
+
+# Disjoint:
+function isdisjoint(a::Interval, b::Interval)
+    (isempty(a) || isempty(b)) && return true
+    islessprime(b.hi, a.lo) || islessprime(a.hi, b.lo)
+end
+
+
+# Intersection
+"""
+    intersect(a, b)
+    ∩(a,b)
+
+Returns the intersection of the intervals `a` and `b`, considered as
+(extended) sets of real numbers. That is, the set that contains
+the points common in `a` and `b`.
+"""
+function intersect{T}(a::Interval{T}, b::Interval{T})
+    isdisjoint(a,b) && return emptyinterval(T)
+
+    Interval(max(a.lo, b.lo), min(a.hi, b.hi))
+end
+# Specific promotion rule for intersect:
+intersect{T,S}(a::Interval{T}, b::Interval{S}) = intersect(promote(a,b)...)
+
+
+## Hull
+"""
+    hull(a, b)
+
+Returns the "convex hull" of the intervals `a` and `b`, considered as
+(extended) sets of real numbers. That is, the minimum set that contains
+all points in `a` and `b`.
+"""
+hull{T}(a::Interval{T}, b::Interval{T}) = Interval(min(a.lo, b.lo), max(a.hi, b.hi))
+
+"""
+    union(a, b)
+    ∪(a,b)
+
+Returns the union (convex hull) of the intervals `a` and `b`; it is equivalent
+to `hull(a,b)`.
+"""
+union(a::Interval, b::Interval) = hull(a, b)
+
+
+
+
+doc"""
+    setdiff(x::Interval, y::Interval)
+
+Calculate the set difference `x \ y`, i.e. the set of values
+that are inside the interval `x` but not inside `y`.
+
+Returns an array of intervals.
+The array may:
+
+- be empty if `x ⊆ y`;
+- contain a single interval, if `y` overlaps `x`
+- contain two intervals, if `y` is strictly contained within `x`.
+"""
+function setdiff(x::Interval, y::Interval)
+    intersection = x ∩ y
+
+    isempty(intersection) && return [x]
+    intersection == x && return typeof(x)[]  # x is subset of y; setdiff is empty
+
+    x.lo == intersection.lo && return [Interval(intersection.hi, x.hi)]
+    x.hi == intersection.hi && return [Interval(x.lo, intersection.lo)]
+
+    return [Interval(x.lo, y.lo),
+            Interval(y.hi, x.hi)]
+
+end

src/multidim/intervalbox.jl

@@ -16,8 +16,11 @@ mid(X::IntervalBox) = [mid(x) for x in X]
 ⊆(X::IntervalBox, Y::IntervalBox) = all([x ⊆ y for (x,y) in zip(X, Y)])
 
 ∩(X::IntervalBox, Y::IntervalBox) = IntervalBox([x ∩ y for (x,y) in zip(X, Y)]...)
+
 isempty(X::IntervalBox) = any([isempty(x) for x in X])
 
+diam(X::IntervalBox) = maximum([diam(x) for x in X])
+
 
 function show(io::IO, X::IntervalBox)
     for (i, x) in enumerate(X)

src/multidim/multidim.jl

@@ -1,2 +1,3 @@
 include("intervalbox.jl")
 include("intervalbox_macro.jl")
+include("setdiff.jl")

src/multidim/setdiff.jl

@@ -0,0 +1,56 @@
+doc"""
+Returns a list of pairs (interval, label)
+label is 1 if the interval is *excluded* from the setdiff
+label is 0 if the interval is included in the setdiff
+label is -1 if the intersection of the two intervals was empty
+"""
+function labelled_setdiff{T}(x::Interval{T}, y::Interval{T})
+    intersection = x ∩ y
+
+    isempty(intersection) && return [(x, -1)]
+    intersection == x && return [(x, 1)]
+
+    x.lo == intersection.lo && return [(intersection,1), (Interval(intersection.hi, x.hi), 0)]
+    x.hi == intersection.hi && return [(Interval(x.lo, intersection.lo), 0), (intersection, 1)]
+
+    return [(Interval(x.lo, y.lo), 0),
+            (y, 1),
+            (Interval(y.hi, x.hi), 0)]
+
+end
+
+doc"""
+    setdiff(A::IntervalBox{2,T}, B::IntervalBox{2,T})
+
+Returns a vector of `IntervalBox`es that are in the set difference `A \ B`,
+i.e. the set of `x` that are in `A` but not in `B`.
+"""
+function setdiff{T}(A::IntervalBox{2,T}, B::IntervalBox{2,T})
+    X = labelled_setdiff(A[1], B[1])
+    Y = labelled_setdiff(A[2], B[2])
+
+    results_list = typeof(A)[]
+
+    for (x, label) in X
+        label == -1 && return [A]   # intersection in one direction empty, so total intersection empty
+
+        if label == 0
+            push!(results_list, IntervalBox(x, A[2]))
+            continue
+        end
+
+        # label is 1 here, so there is some intersection in the x direction
+        for (y, label) in Y
+            label == -1 && return [A]
+
+            if label == 0
+                push!(results_list, IntervalBox(x, y))
+                continue
+            end
+
+            # label == 1:  exclude this box since all labels are 1
+        end
+    end
+
+    return results_list
+end

test/REQUIRE

@@ -1,2 +1,2 @@
 FactCheck 0.3
-Polynomials
+Polynomials 0.0.5

test/interval_tests/set_operations.jl

@@ -0,0 +1,34 @@
+# This file is part of the ValidatedNumerics.jl package; MIT licensed
+
+using ValidatedNumerics
+using FactCheck
+
+setprecision(Interval, Float64)
+
+facts("setdiff") do
+    x = 2..4
+    y = 3..5
+
+    d = setdiff(x, y)
+
+    @fact typeof(d) --> Vector{Interval{Float64}}
+    @fact length(d) --> 1
+    @fact d --> [2..3]
+    @fact setdiff(y, x) --> [4..5]
+
+    x = 2..4
+    y = 2..5
+
+    @fact typeof(d) --> Vector{Interval{Float64}}
+    @fact length(setdiff(x, y)) --> 0
+    @fact setdiff(y, x) --> [4..5]
+
+    x = 2..5
+    y = 3..4
+    @fact setdiff(x, y) --> [2..3, 4..5]
+
+#    X = IntervalBox(2..4, 3..5)
+#    Y = IntervalBox(3..5, 4..6)
+
+    #@fact setdiff(X, Y) --> IntervalBox(2..3, 3..4)
+end