Working n-dimensional setdiff

src/ValidatedNumerics.jl

@@ -25,7 +25,7 @@ import Base:
     precision,
     isfinite, isnan,
     show, showall,
-    isinteger
+    isinteger, setdiff
 
 export
     Interval, AbstractInterval,

src/intervals/functions.jl

@@ -155,9 +155,6 @@ function ^(a::Interval{BigFloat}, x::Interval)
 end
 
 
-Base.inf(x::Rational) = 1//0  # to allow sqrt()
-
-
 function sqrt{T}(a::Interval{T})
     domain = Interval(zero(T), convert(T, Inf))
     a = a ∩ domain

src/multidim/setdiff.jl

@@ -10,47 +10,98 @@ function labelled_setdiff{T}(x::Interval{T}, y::Interval{T})
     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)]
+    x.lo == intersection.lo && return [(intersection, 1), (Interval(intersection.hi, x.hi), 0)]
+    x.hi == intersection.hi && return [(intersection, 1), (Interval(x.lo, intersection.lo), 0)]
 
-    return [(Interval(x.lo, y.lo), 0),
-            (y, 1),
+    return [(y, 1),
+            (Interval(x.lo, y.lo), 0),
             (Interval(y.hi, x.hi), 0)]
 
 end
 
+# function setdiff{N,T}(A::IntervalBox{N,T}, B::IntervalBox{N,T})
+#     X = [labelled_setdiff(a,b) for (a, b) in zip(A, B)]
+#     lengths = map(length, X)
+#     index = ones(N)
+#
+#     # index[j] represents which set we are looking at in direction j
+#``
+#
+#     while index[1] <= N
+#         current_direction = 1
+#         current_piece = [ X[1][index[1]] ]
+#
+#     end
+# end
 doc"""
-    setdiff(A::IntervalBox{2,T}, B::IntervalBox{2,T})
+    setdiff(A::IntervalBox{N,T}, B::IntervalBox{N,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`.
+
+Algorithm: Start from the total overlap (in all directions);
+expand each direction in turn.
 """
-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
+function setdiff{N,T}(A::IntervalBox{N,T}, B::IntervalBox{N,T})
+    X = [labelled_setdiff(a,b) for (a, b) in zip(A, B)]
+    # ordered such that the first in each is the excluded interval
+
+    first = [ i[1] for i in X ]
+    labels = [i[2] for i in first]
+
+    any(labels .== -1) && return [A]  # no overlap
 
-        # label is 1 here, so there is some intersection in the x direction
-        for (y, label) in Y
-            label == -1 && return [A]
+    # @assert all(labels .== 1)
 
-            if label == 0
-                push!(results_list, IntervalBox(x, y))
-                continue
-            end
+    excluded = [i[1] for i in first]
 
-            # label == 1:  exclude this box since all labels are 1
+    result_list = IntervalBox{N,T}[]
+
+    for dimension in N:-1:1
+        for which in X[dimension][2:end]
+            excluded[dimension] = which[1]
+            push!(result_list,
+                IntervalBox(excluded[1:dimension]..., A[dimension+1:end]...))
         end
     end
 
-    return results_list
+    result_list
+
 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/interval_tests/consistency.jl

@@ -288,7 +288,3 @@ facts("Interval power of an interval") do
     @fact b^@interval(0.3) == Interval(1.3903891703159093, 1.5157165665103982) --> true
 
 end
-
-facts("Rational infinity") do
-    @fact inf(3//4) == 1//0 --> true
-end

test/multidim_tests/multidim.jl

@@ -35,6 +35,54 @@ facts("@intervalbox tests") do
 
     @intervalbox g(x, y) = x - y
     @fact isa(g(X), IntervalBox) --> true
+end
+
+facts("setdiff for IntervalBox") do
+    X = IntervalBox(2..4, 3..5)
+    Y = IntervalBox(3..5, 4..6)
+    @fact setdiff(X, Y) --> [ IntervalBox(3..4, 3..4),
+                              IntervalBox(2..3, 3..5) ]
+
+
+    X = IntervalBox(2..5, 3..6)
+    Y = IntervalBox(-10..10, 4..5)
+    @fact setdiff(X, Y) --> [ IntervalBox(2..5, 3..4),
+                              IntervalBox(2..5, 5..6) ]
+
+
+    X = IntervalBox(2..5, 3..6)
+    Y = IntervalBox(4..6, 4..5)
+    @fact setdiff(X, Y) --> [ IntervalBox(4..5, 3..4),
+                              IntervalBox(4..5, 5..6),
+                              IntervalBox(2..4, 3..6) ]
+
+
+    X = IntervalBox(2..5, 3..6)
+    Y = IntervalBox(3..4, 4..5)
+    @fact setdiff(X, Y) --> [ IntervalBox(3..4, 3..4),
+                              IntervalBox(3..4, 5..6),
+                              IntervalBox(2..3, 3..6),
+                              IntervalBox(4..5, 3..6) ]
+
+
+    X = IntervalBox(2..5, 3..6)
+    Y = IntervalBox(2..4, 10..20)
+    @fact setdiff(X, Y) --> typeof(X)[X]
+
+
+    X = IntervalBox(2..5, 3..6)
+    Y = IntervalBox(-10..10, -10..10)
+    @fact setdiff(X, Y) --> typeof(X)[]
+
+
+    X = IntervalBox(1..4, 3..6, 7..10)
+    Y = IntervalBox(2..3, 4..5, 8..9)
+    @fact setdiff(X, Y) --> [ IntervalBox(2..3, 4..5, 7..8),
+                              IntervalBox(2..3, 4..5, 9..10),
+                              IntervalBox(2..3, 3..4, 7..10),
+                              IntervalBox(2..3, 5..6, 7..10),
+                              IntervalBox(1..2, 3..6, 7..10),
+                              IntervalBox(3..4, 3..6, 7..10) ]
 
 
 end