Basic support for multidimensional IntervalBox's

REQUIRE

@@ -1,5 +1,4 @@
 julia 0.4
-FactCheck 0.3
-Polynomials
 CRlibm 0.2.2
 Compat 0.7.11
+FixedSizeArrays

examples/draw_function_image.jl

@@ -0,0 +1,30 @@
+using ValidatedNumerics
+
+using PyCall
+using PyPlot
+
+@pyimport matplotlib.patches as patches
+@pyimport matplotlib.collections as collections
+
+
+const rectangle = patches.Rectangle
+
+function make_rectangle(x, y, xwidth, ywidth, color="grey", alpha=0.5)
+    patches.Rectangle(
+        (x, y), xwidth, ywidth, facecolor=color, alpha=alpha
+    )
+end
+
+import PyPlot.draw
+function draw{T<:IntervalBox}(box_list::Vector{T}, color="grey", alpha=0.5)
+    patch_list = []
+    for box in box_list
+        x, y = box
+        push!(patch_list,
+            make_rectangle(x.lo, y.lo, x.hi-x.lo, y.hi-y.lo, color, alpha))
+    end
+
+    ax = gca()
+    ax[:add_collection](collections.PatchCollection(patch_list, color=color, alpha=alpha,
+    edgecolor="black"))
+end

examples/standard_map.jl

@@ -0,0 +1,69 @@
+using ValidatedNumerics
+
+import Base.mod2pi
+
+function standard_map(X::IntervalBox, k = 1.0)
+    p, θ = X
+
+    p′ = p + k*sin(θ)
+    @show p′
+    θ′ = mod2pi( θ + p′ )
+    p′ = mod2pi(p′)
+
+    @show p′, θ′
+    @show typeof(p′), typeof(θ′)
+
+    IntervalBox(p′, θ′)
+end
+
+function IntervalBox{T}(X::Vector{Interval{T}}, Y::Vector{Interval{T}})
+    vec([IntervalBox(x, y) for x in X, y in Y])
+end
+
+function iterate(f, n, x)
+    for i in 1:n
+        x = f(x)
+    end
+    x
+end
+
+
+import Base.mod
+function mod(X::Interval, width::Real)
+
+
+    @show X, width
+
+    X /= width
+
+    if diam(X) >= 1.
+
+        return [Interval(0, 1) * width]
+    end
+
+    a = X.lo - floor(X.lo)
+    b = X.hi - floor(X.hi)
+
+    if a < b
+        return [Interval(a, b)*width]
+
+    end
+
+    return [Interval(0, b)*width, Interval(a, 1)*width]
+
+end
+
+
+function mod(X::IntervalBox, width::Real)
+    x, y = X
+
+    xx = mod(x, width)
+    yy = mod(y, width)
+
+    vec([IntervalBox(x, y) for x in xx, y in yy])
+end
+
+
+mod2pi{T}(x::Interval{T}) = mod(x, ValidatedNumerics.two_pi(T))
+
+mod2pi{T}(X::Vector{Interval{T}}) = vcat(map(mod2pi, X)...)

src/ValidatedNumerics.jl

@@ -4,10 +4,10 @@ __precompile__(true)
 
 module ValidatedNumerics
 
+using CRlibm
 using Compat
-#using FactCheck
+using FixedSizeArrays
 
-using CRlibm
 
 setrounding(BigFloat, RoundNearest)
 setrounding(Float64, RoundNearest)
@@ -25,7 +25,8 @@ import Base:
     BigFloat, float, widen,
     ⊆, eps,
     floor, ceil, trunc, sign, round,
-    expm1, log1p
+    expm1, log1p,
+    isfinite, isnan
 
 
 export
@@ -43,7 +44,13 @@ export
     pi_interval,
     midpoint_radius, interval_from_midpoint_radius,
     RoundTiesToEven, RoundTiesToAway,
-    cancelminus, cancelplus, isunbounded
+    cancelminus, cancelplus, isunbounded,
+    .., @I_str
+
+
+## Multidimensional
+export
+    IntervalBox, @intervalbox
 
 ## Root finding
 export
@@ -63,6 +70,8 @@ end
 
 include("misc.jl")
 include("intervals/intervals.jl")
+include("multidim/multidim.jl")
+
 include("root_finding/root_finding.jl")
 
 

src/intervals/intervals.jl

@@ -50,3 +50,4 @@ include("precision.jl")
 include("functions.jl")
 include("trigonometric_functions.jl")
 include("hyperbolic_functions.jl")
+include("syntax.jl")

src/intervals/special_intervals.jl

@@ -33,6 +33,8 @@ nai() = nai(get_interval_precision()[1])
 
 isnai(x::Interval) = isnan(x.lo) || isnan(x.hi)
 
+isfinite(x::Interval) = isfinite(x.lo) && isfinite(x.hi)
+isnan(x::Interval) = isnai(x)
 
 doc"""`isthin(x)` corresponds to `isSingleton`, i.e. it checks if `x` is the set consisting of a single exactly representable float. Thus any float which is not exactly representable does *not* yield a thin interval."""
 function isthin(x::Interval)

src/intervals/syntax.jl

@@ -0,0 +1,7 @@
+# Syntax for intervals
+
+a..b = @interval(a, b)
+
+macro I_str(ex)  # I"[3,4]"
+    @interval(ex)
+end

src/multidim/box.jl

@@ -0,0 +1,32 @@
+# Multidimensional intervals, called IntervalBox
+
+#using ValidatedNumerics
+#using FixedSizeArrays
+
+import Base:
+    ⊆, ∩, isempty
+
+
+doc"""An `IntervalBox` is a Cartesian product of an arbitrary number of `Interval`s,
+representing an $N$-dimensional rectangular IntervalBox."""
+
+immutable IntervalBox{N, T} <: FixedVector{N, Interval{T}} # uses FixedSizeArrays package
+    intervals :: NTuple{N, Interval{T}}
+end
+
+mid(X::IntervalBox) = [mid(x) for x in X.intervals]
+
+⊆(X::IntervalBox, Y::IntervalBox) = all([x ⊆ y for (x,y) in zip(X.intervals, Y.intervals)])
+
+∩(X::IntervalBox, Y::IntervalBox) = IntervalBox([x ∩ y for (x,y) in zip(X.intervals, Y.intervals)]...)
+isempty(X::IntervalBox) = any(map(isempty, X.intervals))
+
+
+function Base.show(io::IO, X::IntervalBox)
+    for (i, x) in enumerate(X)
+        print(io, x)
+        if i != length(X)
+            print(io, " × ")
+        end
+    end
+end

src/multidim/box_macro.jl

@@ -0,0 +1,42 @@
+using ValidatedNumerics
+
+
+"""Macro to make a version of a multidimensional function that acts on
+an `IntervalBox` and returns an `IntervalBox`.
+Both the original n-argument function and the `IntervalBox` version are defined.
+
+Example:
+
+    @intervalbox f(x, y) = (x + y, x - y)
+
+    X = IntervalBox(1..1, 2..2)
+    f(X)
+
+(No significant error checking is performed!)
+"""
+
+macro intervalbox(ex)
+    # @show ex
+    # @show ex.head
+    if ex.head != :(=)
+        throw(ArgumentError("Not a function definition."))
+    end
+
+    if ex.args[1].head != :call
+        throw(ArgumentError("Not a function definition."))
+    end
+
+    f = ex.args[1].args[1]  # function name
+    f = esc(f)
+    ex.args[1].args[1] = f  # escape the function name in the original code
+
+
+    new_ex = Expr(:block)
+
+    push!(new_ex.args, ex)  # the function call
+
+    push!(new_ex.args,  :( ($f)(X::IntervalBox) = IntervalBox( $f(X...)...) ) )
+
+    #@show new_ex
+    new_ex
+end

src/multidim/multidim.jl

@@ -0,0 +1,2 @@
+include("box.jl")
+include("box_macro.jl")

test/REQUIRE

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

test/multidim_tests/multidim_tests.jl

@@ -0,0 +1,31 @@
+using FactCheck
+
+facts("Operations on boxes") do
+    A = IntervalBox(1..2, 3..4)
+    B = IntervalBox(0..2, 3..6)
+
+    @fact 2*A --> IntervalBox(2..4, 6..8)
+    @fact A + B --> IntervalBox(1..4, 6..10)
+    @fact dot(A, B) --> @interval(9, 28)
+
+    @fact A ⊆ B --> true
+
+    X = IntervalBox(1..2, 3..4)
+    Y = IntervalBox(3..4, 3..4)
+
+    @fact isempty(X ∩ Y) --> true
+
+    v = [@interval(i, i+1) for i in 1:10]
+    V = IntervalBox(v...)
+
+    @fact length(V) --> 10
+
+end
+
+facts("@box tests") do
+    @intervalbox f(x, y) = (x + y, x - y)
+
+    X = IntervalBox(1..1, 2..2)
+    @fact f(X) --> IntervalBox(3..3, -1 .. -1)
+
+end

test/runtests.jl

@@ -14,13 +14,16 @@ include("misc_tests.jl")
 # Interval tests:
 
 include("interval_tests/interval_tests.jl")
+include("multidim_tests/multidim_tests.jl")
 
 
 # Root-finding tests:
 
 include("root_finding_tests/auto_diff_tests.jl")
 include("root_finding_tests/root_finding_tests.jl")
 
+# Multidimensional boxes tests:
+
 # ITF1788 tests
 
 include("ITF1788_tests/ITF1788_tests.jl")