diff --git a/.travis.yml b/.travis.yml
index caf2f5c..ac31ec4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,8 +4,8 @@ os:
   - linux
   - osx
 julia:
-  - 0.5
-  - 0.6
+  - 0.7
+  - 1.0
   - nightly
 notifications:
   email: false
diff --git a/REQUIRE b/REQUIRE
index 4f32cce..bfd619e 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,4 +1,7 @@
-julia 0.5
+julia 0.7
 
-IntervalArithmetic 0.9.1
-IntervalRootFinding 0.1
+IntervalArithmetic 0.15
+IntervalRootFinding 0.2
+IntervalConstraintProgramming
+DataStructures
+ForwardDiff 0.8
diff --git a/appveyor.yml b/appveyor.yml
index 64262ea..c2588f1 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -1,26 +1,43 @@
 environment:
   matrix:
-  - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
-  - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe"
+  - julia_version: 0.7
+  - julia_version: 1
+  - julia_version: nightly
+
+platform:
+  - x86 # 32-bit
+  - x64 # 64-bit
+
+# # Uncomment the following lines to allow failures on nightly julia
+# # (tests will run but not make your overall status red)
+# matrix:
+#   allow_failures:
+#   - julia_version: nightly
+
 branches:
   only:
     - master
     - /release-.*/
+
 notifications:
   - provider: Email
     on_build_success: false
     on_build_failure: false
     on_build_status_changed: false
+
 install:
-  - ps: (new-object net.webclient).DownloadFile(
-        $("http://s3.amazonaws.com/"+$env:JULIAVERSION),
-        "C:\projects\julia-binary.exe")
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
+  - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))
+
 build_script:
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"IntervalOptimisation\"); Pkg.build(\"IntervalOptimisation\")"
+  - echo "%JL_BUILD_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"
+
 test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"IntervalOptimisation\")"
+  - echo "%JL_TEST_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
+
+# # Uncomment to support code coverage upload. Should only be enabled for packages
+# # which would have coverage gaps without running on Windows
+# on_success:
+#   - echo "%JL_CODECOV_SCRIPT%"
+#   - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%"
diff --git a/examples/griewank.jl b/examples/griewank.jl
new file mode 100644
index 0000000..5873274
--- /dev/null
+++ b/examples/griewank.jl
@@ -0,0 +1,12 @@
+using IntervalArithmetic, IntervalOptimisation
+
+G(X) = 1 + sum(abs2, X) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) )
+
+@time minimise(G, -600..600)
+@time minimise(G, -600..600)
+
+for N in (10, 30, 50)
+    @show N
+    @time minimise(G, IntervalBox(-600..600, N))
+    @time minimise(G, IntervalBox(-600..600, N))
+end
diff --git a/examples/optimise1.jl b/examples/optimise1.jl
new file mode 100644
index 0000000..13892cb
--- /dev/null
+++ b/examples/optimise1.jl
@@ -0,0 +1,56 @@
+using IntervalOptimisation, IntervalArithmetic, IntervalConstraintProgramming
+
+
+# Unconstrained Optimisation
+# Example from Eldon Hansen - Global Optimisation using Interval Analysis, Chapter 11
+f(x, y) = (1.5 - x * (1 - y))^2 + (2.25 - x * (1 - y^2))^2 + (2.625 - x * (1 - y^3))^2
+# f (generic function with 2 methods)
+
+f(X) = f(X...)
+# f (generic function with 2 methods)
+
+X = (-1e6..1e6) × (-1e6..1e6)
+# [-1e+06, 1e+06] × [-1e+06, 1e+06]
+
+minimise_icp(f, X)
+# ([0, 2.39527e-07], IntervalArithmetic.IntervalBox{2,Float64}[[2.99659, 2.99748] × [0.498906, 0.499523], [2.99747, 2.99834] × [0.499198, 0.500185], [2.99833, 2.99922] × [0.499198, 0.500185], [2.99921, 3.00011] × [0.499621, 0.500415], [3.0001, 3.001] × [0.499621, 0.500415], [3.00099, 3.0017] × [0.500169, 0.500566], [3.00169, 3.00242] × [0.500169, 0.500566]])
+
+
+
+f(X) = X[1]^2 + X[2]^2
+# f (generic function with 2 methods)
+
+X = (-∞..∞) × (-∞..∞)
+# [-∞, ∞] × [-∞, ∞]
+
+minimise_icp(f, X)
+# ([0, 0], IntervalArithmetic.IntervalBox{2,Float64}[[-0, 0] × [-0, 0], [0, 0] × [-0, 0]])
+
+
+
+
+# Constrained Optimisation
+# Example 1 from http://archimedes.cheme.cmu.edu/?q=baronexamples
+f(X) = -1 * (X[1] + X[2])
+# f (generic function with 1 method)
+
+X = (-∞..∞) × (-∞..∞)
+# [-∞, ∞] × [-∞, ∞]
+
+constraints = [IntervalOptimisation.constraint(x->(x[1]), -∞..6), IntervalOptimisation.constraint(x->x[2], -∞..4), IntervalOptimisation.constraint(x->(x[1]*x[2]), -∞..4)]
+# 3-element Array{IntervalOptimisation.constraint{Float64},1}:
+#  IntervalOptimisation.constraint{Float64}(#3, [-∞, 6])
+#  IntervalOptimisation.constraint{Float64}(#4, [-∞, 4])
+#  IntervalOptimisation.constraint{Float64}(#5, [-∞, 4])
+
+minimise_icp_constrained(f, X, constraints)
+# ([-6.66676, -6.66541], IntervalArithmetic.IntervalBox{2,Float64}[[5.99918, 6] × [0.666233, 0.666758], [5.99918, 6] × [0.665717, 0.666234], [5.99887, 5.99919] × [0.666233, 0.666826], [5.99984, 6] × [0.665415, 0.665718], [5.99856, 5.99888] × [0.666233, 0.666826], [5.99969, 5.99985] × [0.665415, 0.665718]])
+
+
+
+# One-dimensional case
+minimise1d(x -> (x^2 - 2)^2, -10..11)
+# ([0, 1.33476e-08], IntervalArithmetic.Interval{Float64}[[-1.41426, -1.41358], [1.41364, 1.41429]])
+
+minimise1d_deriv(x -> (x^2 - 2)^2, -10..11)
+# ([0, 8.76812e-08], IntervalArithmetic.Interval{Float64}[[-1.41471, -1.41393], [1.41367, 1.41444]])
diff --git a/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl
index 7845b75..654aa61 100644
--- a/src/IntervalOptimisation.jl
+++ b/src/IntervalOptimisation.jl
@@ -2,25 +2,23 @@
 
 module IntervalOptimisation
 
-export minimise, maximise,
-       minimize, maximize
+using DataStructures
 
-include("SortedVectors.jl")
-using .SortedVectors
+export minimise, maximise,
+       minimize, maximize,
+       minimise1d, minimise1d_deriv,
+       minimise_icp, minimise_icp_constrained
 
-using IntervalArithmetic, IntervalRootFinding
+# include("SortedVectors.jl")
+# using .SortedVectors
 
-if !isdefined(:sup)
-    const sup = supremum
-end
+using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, ForwardDiff
 
-if !isdefined(:inf)
-    const inf = infimum
-end
 
 include("optimise.jl")
+include("optimise1.jl")
 
 const minimize = minimise
-const maximize = maximise 
+const maximize = maximise
 
 end
diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl
index 37f647c..e950342 100644
--- a/src/SortedVectors.jl
+++ b/src/SortedVectors.jl
@@ -1,25 +1,25 @@
 module SortedVectors
 
 import Base: getindex, length, push!, isempty,
-        pop!, shift!, resize!
+        pop!, resize!, popfirst!
 
 export SortedVector
 
 """
 A `SortedVector` behaves like a standard Julia `Vector`, except that its elements are stored in sorted order, with an optional function `by` that determines the sorting order in the same way as `Base.searchsorted`.
 """
-immutable SortedVector{T, F<:Function}
+struct SortedVector{T, F<:Function}
     data::Vector{T}
     by::F
 
-    function SortedVector(data::Vector{T}, by::F)
-        new(sort(data), by)
+    function SortedVector(data::Vector{T}, by::F) where {T,F}
+        new{T,F}(sort(data, by=by), by)
     end
 end
 
 
-SortedVector{T,F}(data::Vector{T}, by::F) = SortedVector{T,F}(data, by)
-SortedVector{T}(data::Vector{T}) = SortedVector{T,typeof(identity)}(data, identity)
+#SortedVector(data::Vector{T}, by::F) where {T,F} = SortedVector(data, by)
+SortedVector(data::Vector{T}) where {T} = SortedVector{T,typeof(identity)}(data, identity)
 
 function show(io::IO, v::SortedVector)
     print(io, "SortedVector($(v.data))")
@@ -30,7 +30,7 @@ end
 getindex(v::SortedVector, i::Int) = v.data[i]
 length(v::SortedVector) = length(v.data)
 
-function push!{T}(v::SortedVector{T}, x::T)
+function push!(v::SortedVector{T}, x::T) where {T}
     i = searchsortedfirst(v.data, x, by=v.by)
     insert!(v.data, i, x)
     return v
@@ -40,7 +40,7 @@ isempty(v::SortedVector) = isempty(v.data)
 
 pop!(v::SortedVector) = pop!(v.data)
 
-shift!(v::SortedVector) = shift!(v.data)
+popfirst!(v::SortedVector) = popfirst!(v.data)
 
 
 function resize!(v::SortedVector, n::Int)
diff --git a/src/optimise.jl b/src/optimise.jl
index 337b8f4..9519c8f 100644
--- a/src/optimise.jl
+++ b/src/optimise.jl
@@ -4,59 +4,88 @@
 
 Find the global minimum of the function `f` over the `Interval` or `IntervalBox` `X` using the Moore-Skelboe algorithm.
 
-For higher-dimensional functions ``f:\mathbb{R}^n \to \mathbb{R}``, `f` must take a single vector argument.
+For higher-dimensional functions ``f:\\mathbb{R}^n \\to \\mathbb{R}``, `f` must take a single vector argument.
 
 Returns an interval containing the global minimum, and a list of boxes that contain the minimisers.
 """
-function minimise{T}(f, X::T, tol=1e-3)
+function minimise(f, X::T, tol=1e-3) where {T}
 
     # list of boxes with corresponding lower bound, ordered by increasing lower bound:
-    working = SortedVector([(X, ∞)], x->x[2])
+    working = PriorityQueue(X => ∞)
 
-    minimizers = T[]
-    global_min = ∞  # upper bound
+    return _minimise(f, working, tol)
+end
+
+function minimise(f, Xs::Vector{T}, tol=1e-3) where {T}
+
+    # list of boxes with corresponding lower bound, ordered by increasing lower bound:
+    working = PriorityQueue(Dict(X => f(X).lo for X in Xs))
+
+    return _minimise(f, working, tol)
+end
+
+function _minimise(f, working::PriorityQueue{K,V}, tol) where {K,V}
+
+    minimizers = K[]
+    global_min = ∞   # upper bound for the global minimum value
 
     num_bisections = 0
 
     while !isempty(working)
 
-        (X, X_min) = shift!(working)
+        #@show working, minimizers
+
+        (XX, X_min) = dequeue_pair!(working)
 
         if X_min > global_min    # X_min == inf(f(X))
             continue
         end
 
         # find candidate for upper bound of global minimum by just evaluating a point in the interval:
-        m = sup(f(Interval.(mid.(X))))   # evaluate at midpoint of current interval
 
-        if m < global_min
+        #@show f(interval.(mid(XX)))
+
+        m = sup(f(interval.(mid(XX))))   # evaluate at midpoint of current interval
+
+        if m < global_min && m > -∞
             global_min = m
         end
 
         # Remove all boxes whose lower bound is greater than the current one:
         # Since they are ordered, just find the first one that is too big
 
-        cutoff = searchsortedfirst(working.data, (X, global_min), by=x->x[2])
-        resize!(working, cutoff-1)
+        # cutoff = searchsortedfirst(working.data, (X, global_min), by=x->x[2])
+        # resize!(working, cutoff)
 
-        if diam(X) < tol
-            push!(minimizers, X)
+        if diam(XX) <= tol
+            push!(minimizers, XX)
 
         else
-            X1, X2 = bisect(X)
-            push!( working, (X1, inf(f(X1))), (X2, inf(f(X2))) )
+            X1, X2 = bisect(XX)
+
+            inf1 = inf(f(X1))
+            inf1 <= global_min && enqueue!(working, X1 => inf1 )
+
+            inf2 = inf(f(X2))
+            inf2 <= global_min && enqueue!(working, X2 => inf2 )
+
             num_bisections += 1
         end
 
     end
 
-    lower_bound = minimum(inf.(f.(minimizers)))
+    # @show global_min
+    # @show minimizers
+
+    new_minimizers = [minimizer for minimizer in minimizers if f(minimizer).lo <= global_min]
+
+    lower_bound = minimum(inf.(f.(new_minimizers)))
 
-    return Interval(lower_bound, global_min), minimizers
+    return Interval(lower_bound, global_min), new_minimizers
 end
 
 
-function maximise{T}(f, X::T, tol=1e-3)
+function maximise(f, X::T, tol=1e-3) where {T}
     bound, minimizers = minimise(x -> -f(x), X, tol)
     return -bound, minimizers
 end
diff --git a/src/optimise1.jl b/src/optimise1.jl
new file mode 100644
index 0000000..8201eb3
--- /dev/null
+++ b/src/optimise1.jl
@@ -0,0 +1,196 @@
+struct IntervalMinimum{T<:Real}
+    interval::Interval{T}
+    minimum::T
+end
+
+Base.isless(a::IntervalMinimum{T}, b::IntervalMinimum{T}) where {T<:Real} = isless(a.minimum, b.minimum)
+
+function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3, use_deriv=false, use_second_deriv=false) where {T<:Real}
+
+    Q = binary_minheap(IntervalMinimum{T})
+
+    global_minimum = f(interval(mid(x))).hi
+
+    arg_minima = Interval{T}[]
+
+    push!(Q, IntervalMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        if use_deriv
+            deriv = ForwardDiff.derivative(f, p.interval)
+            if 0 ∉ deriv
+                continue
+            end
+        end
+
+        # Second derivative contractor
+        if use_second_deriv
+            doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.interval)
+            if doublederiv < 0
+                continue
+            end
+        end
+
+        m = mid(p.interval)
+        current_minimum = f(interval(m)).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+
+        # Contractor 1
+        if use_deriv
+            x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv)
+
+            x = x .∩ p.interval
+        end
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+        else
+
+            if use_deriv && isempty(x[2])
+                x1, x2 = bisect(x[1])
+                push!(Q, IntervalMinimum(x1, f(x1).lo), IntervalMinimum(x2, f(x2).lo))
+                continue
+            end
+
+            x1, x2 = bisect(p.interval)
+            push!(Q, IntervalMinimum(x1, f(x1).lo), IntervalMinimum(x2, f(x2).lo))
+        end
+    end
+    lb = minimum(inf.(f.(arg_minima)))
+
+    return lb..global_minimum, arg_minima
+end
+
+
+struct IntervalBoxMinimum{N, T<:Real}
+    interval::IntervalBox{N, T}
+    minimum::T
+end
+
+"""
+Datatype to provide constraints to Global Optimisation such as:
+```
+Constraint(x->(x^2 - 10), -∞..1)
+```
+"""
+struct Constraint{T<:Real}
+    f::Function
+    c::Interval{T}
+end
+
+Base.isless(a::IntervalBoxMinimum{N, T}, b::IntervalBoxMinimum{N, T}) where {N, T<:Real} = isless(a.minimum, b.minimum)
+
+function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real}
+
+    Q = binary_minheap(IntervalBoxMinimum{N, T})
+
+    global_minimum = ∞
+
+    x = icp(f, x, -∞..global_minimum)
+
+    arg_minima = IntervalBox{N, T}[]
+
+    push!(Q, IntervalBoxMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        current_minimum = sup(f(interval.(mid(p.interval))))
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+        X = icp(f, p.interval, -∞..global_minimum)
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+
+        else
+            x1, x2 = bisect(X)
+            push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo))
+        end
+    end
+
+    #@show arg_minima
+
+    lb = minimum(inf.(f.(arg_minima)))
+
+    return lb..global_minimum, arg_minima
+end
+
+function minimise_icp_constrained(f::Function, x::IntervalBox{N, T}, constraints::Vector{Constraint{T}} = Vector{Constraint{T}}(); reltol=1e-3, abstol=1e-3) where {N, T<:Real}
+
+    Q = binary_minheap(IntervalBoxMinimum{N, T})
+
+    global_minimum = ∞
+
+    for t in constraints
+        x = icp(t.f, x, t.c)
+    end
+
+    x = icp(f, x, -∞..global_minimum)
+
+    arg_minima = IntervalBox{N, T}[]
+
+    push!(Q, IntervalBoxMinimum(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.interval)
+            continue
+        end
+
+        if p.minimum > global_minimum
+            continue
+        end
+
+        # current_minimum = f(interval.(mid(p.interval))).hi
+        current_minimum = f(p.interval).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+        X = icp(f, p.interval, -∞..global_minimum)
+
+        for t in constraints
+            X = icp(t.f, X, t.c)
+        end
+
+        if diam(p.interval) < abstol
+            push!(arg_minima, p.interval)
+        else
+            x1, x2 = bisect(X)
+            push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo))
+        end
+    end
+    lb = minimum(inf.(f.(arg_minima)))
+    return lb..global_minimum, arg_minima
+end