diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
new file mode 100644
index 0000000..d77d3a0
--- /dev/null
+++ b/.github/workflows/TagBot.yml
@@ -0,0 +1,11 @@
+name: TagBot
+on:
+  schedule:
+    - cron: 0 * * * *
+jobs:
+  TagBot:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
diff --git a/.travis.yml b/.travis.yml
index a6fbab2..01b2f71 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -6,8 +6,8 @@ os:
   - osx
 
 julia:
-  - 0.7
   - 1.0
+  - 1.5
   - nightly
 
 notifications:
@@ -15,7 +15,6 @@ notifications:
 
 matrix:
   allow_failures:
-    - julia: 0.7
     - julia: nightly
 
 after_success:
diff --git a/NEWS.md b/NEWS.md
index 97b117c..3dc4f92 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,8 @@
 # NEWS for IntervalOptimisation.jl
 
+## v0.4
+- Drop support for Julia 0.7  
+
 ## v0.3
 - Drop support for Julia 0.6. The package is now fully compatible with Julia 1.0.
 
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..d59c4fc
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,19 @@
+name = "IntervalOptimisation"
+uuid = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2"
+version = "0.4.2"
+
+[deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[compat]
+IntervalArithmetic = "0.15, 0.16, 0.17"
+julia = "1.0"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/README.md b/README.md
index 55f9dd0..ab15ffc 100644
--- a/README.md
+++ b/README.md
@@ -4,8 +4,6 @@
 [![appveyor badge][appveyor_badge]][appveyor_url]
 [![codecov badge][codecov_badge]][codecov_url]
 
-## Documentation [here][documenter_latest]
-
 [travis_badge]: https://travis-ci.org/JuliaIntervals/IntervalOptimisation.jl.svg?branch=master
 [travis_url]: https://travis-ci.org/JuliaIntervals/IntervalOptimisation.jl
 
@@ -35,7 +33,7 @@ They return an `Interval` that is guaranteed to contain the global minimum (maxi
 
 
 #### 1D
-```
+```julia
 using IntervalArithmetic, IntervalOptimisation
 
 julia> @time global_min, minimisers = minimise(x -> (x^2 - 2)^2, -10..11);
@@ -52,7 +50,7 @@ julia> minimisers
 
 #### 2D
 
-```
+```julia
 julia> @time global_min, minimisers = minimise(  X -> ( (x,y) = X; x^2 + y^2 ),
                                                         (-10000..10001) × (-10000..10001) );
   0.051122 seconds (46.80 k allocations: 2.027 MiB)
diff --git a/REQUIRE b/REQUIRE
deleted file mode 100644
index 3161ea8..0000000
--- a/REQUIRE
+++ /dev/null
@@ -1,7 +0,0 @@
-julia 0.7
-
-IntervalArithmetic 0.15
-IntervalRootFinding 0.4
-IntervalConstraintProgramming 0.9
-DataStructures 0.9
-ForwardDiff 0.8
diff --git a/appveyor.yml b/appveyor.yml
index c2588f1..c799c17 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -1,7 +1,7 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1
+  - julia_version: 1.5
   - julia_version: nightly
 
 platform:
diff --git a/src/HeapedVectors.jl b/src/HeapedVectors.jl
new file mode 100644
index 0000000..98f2092
--- /dev/null
+++ b/src/HeapedVectors.jl
@@ -0,0 +1,107 @@
+__precompile__()
+
+module HeapedVectors
+
+import Base: getindex, length, push!, isempty,
+        pop!, filter!, popfirst!
+
+export HeapedVector
+
+import ..StrategyBase:filter_elements!
+using ..StrategyBase
+
+struct HeapedVector{T, F<:Function} <: Strategy
+    data::Vector{T}
+    by::F
+    function HeapedVector(v::Vector{T}, by::F) where {T, F}
+        new{T, F}(heaping(v, by), by)
+    end
+end
+
+HeapedVector(data::Vector{T}) where {T} = HeapedVector(data, identity)
+
+
+function heaping(v, by)
+    ar = typeof(v[1])[]
+    for i = 1:length(v)
+        insert!(ar, i, v[i])
+        floatup!(ar, length(ar), by)
+    end
+    return ar
+end
+
+function floatup!(ar, index, by)
+    par = convert(Int, floor(index/2))
+    if index <= 1
+        return ar
+    end
+    if by(ar[index]) < by(ar[par])
+       ar[par], ar[index] = ar[index], ar[par]
+    end
+    floatup!(ar, par, by)
+end
+
+
+function push!(v::HeapedVector{T}, x::T) where {T}
+    insert!(v.data, length(v.data)+1, x)
+    floatup!(v.data, length(v.data), v.by)
+    return v
+end
+
+
+isempty(v::HeapedVector) = isempty(v.data)
+
+
+function popfirst!(v::HeapedVector{T}) where {T}
+    if length(v.data) == 0
+         return
+    end
+
+    if length(v.data) > 2
+        v.data[length(v.data)], v.data[1] = v.data[1], v.data[length(v.data)]
+    	minm = pop!(v.data)
+    	bubbledown!(v::HeapedVector{T}, 1)
+
+    elseif length(v.data) == 2
+        v.data[length(v.data)], v.data[1] = v.data[1], v.data[length(v.data)]
+    	minm = pop!(v.data)
+    else
+    	minm = pop!(v.data)
+    end
+    return minm
+end
+
+
+function bubbledown!(v::HeapedVector{T}, index) where{T}
+    left = index*2
+    right = index*2+1
+    smallest = index
+
+    if length(v.data)+1 > left && v.by(v.data[smallest]) > v.by(v.data[left])
+        smallest = left
+    end
+
+    if length(v.data)+1 > right && v.by(v.data[smallest]) > v.by(v.data[right])
+	   smallest = right
+    end
+
+    if smallest != index
+       v.data[index], v.data[smallest] = v.data[smallest], v.data[index]
+	   bubbledown!(v, smallest)
+    end
+end
+
+function filter_elements!(A::HeapedVector{T}, x::T) where{T}
+    func(y) = A.by(y) < A.by(x)
+    filter!(func, A.data)
+
+    if length(A.data) == 0
+        return A
+    end
+
+    heaping(A.data, A.by)
+    return A
+
+end
+
+end
diff --git a/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl
index 65cd52c..aaf3570 100644
--- a/src/IntervalOptimisation.jl
+++ b/src/IntervalOptimisation.jl
@@ -1,17 +1,24 @@
-
-
 module IntervalOptimisation
 
 export minimise, maximise,
        minimize, maximize
+export HeapedVector, SortedVector
+export mean_value_form_scalar, third_order_taylor_form_scalar
+
+include("StrategyBase.jl")
+using .StrategyBase
 
 include("SortedVectors.jl")
 using .SortedVectors
 
-using IntervalArithmetic, IntervalRootFinding
+include("HeapedVectors.jl")
+using .HeapedVectors
 
+using IntervalArithmetic, IntervalRootFinding
+using LinearAlgebra
 
 include("optimise.jl")
+include("centered_forms.jl") 
 
 const minimize = minimise
 const maximize = maximise
diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl
index 93b8dfb..c42dea7 100644
--- a/src/SortedVectors.jl
+++ b/src/SortedVectors.jl
@@ -3,31 +3,29 @@ __precompile__()
 module SortedVectors
 
 import Base: getindex, length, push!, isempty,
-        pop!, resize!, popfirst!
+        pop!, 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`.
-"""
-struct SortedVector{T, F<:Function}
+import ..StrategyBase:filter_elements!
+using ..StrategyBase
+
+struct SortedVector{T, F<:Function} <: Strategy
+
     data::Vector{T}
     by::F
 
     function SortedVector(data::Vector{T}, by::F) where {T,F}
-        new{T,F}(sort(data,by=by), by)
+        new{T,F}(sort(data, by = by), by)
     end
 end
 
-
 SortedVector(data::Vector{T}) where {T} = SortedVector(data, identity)
 
 function show(io::IO, v::SortedVector)
     print(io, "SortedVector($(v.data))")
 end
 
-
-
 getindex(v::SortedVector, i::Int) = v.data[i]
 length(v::SortedVector) = length(v.data)
 
@@ -43,9 +41,9 @@ pop!(v::SortedVector) = pop!(v.data)
 
 popfirst!(v::SortedVector) = popfirst!(v.data)
 
-
-function resize!(v::SortedVector, n::Int)
-    resize!(v.data, n)
+function filter_elements!(v::SortedVector{T}, x::T) where {T}
+    cutoff = searchsortedfirst(v.data, x, by=v.by)
+    resize!(v.data, cutoff-1)
     return v
 end
 
diff --git a/src/StrategyBase.jl b/src/StrategyBase.jl
new file mode 100644
index 0000000..b88e833
--- /dev/null
+++ b/src/StrategyBase.jl
@@ -0,0 +1,7 @@
+module StrategyBase
+
+  export filter_elements!, Strategy
+  abstract type Strategy end
+  function filter_elements!(s::Strategy)end
+
+end
diff --git a/src/centered_forms.jl b/src/centered_forms.jl
new file mode 100644
index 0000000..42278a7
--- /dev/null
+++ b/src/centered_forms.jl
@@ -0,0 +1,33 @@
+import ForwardDiff: gradient, jacobian, hessian
+
+gradient(f, X::IntervalBox) = gradient(f, X.v)
+# jacobian(f, X::IntervalBox) = jacobian(f, X.v)
+hessian(f, X::IntervalBox) = hessian(f, X.v)
+
+"""
+Calculate the mean-value form of a function ``f:\\mathbb{R}^n \\to \\mathbb{R}``
+using the gradient ``\nabla f``;
+this gives a second-order approximation.
+"""
+function mean_value_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    return f(m) + gradient(f, X.v) ⋅ (X - m)
+end
+
+mean_value_form_scalar(f) = X -> mean_value_form_scalar(f, X)
+
+
+"""
+Calculate a third-order Taylor form of ``f:\\mathbb{R}^n \\to \\mathbb{R}`` using the Hessian.
+"""
+function third_order_taylor_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    H = hessian(f, X)
+    δ = X - m
+
+    return f(m) + gradient(f, m) ⋅ δ + 0.5 * sum(δ[i]*H[i,j]*δ[j] for i in 1:length(X) for j in 1:length(X))
+end
+
+third_order_taylor_form_scalar(f) = X -> third_order_taylor_form_scalar(f, X)
diff --git a/src/optimise.jl b/src/optimise.jl
index 4629811..e244162 100644
--- a/src/optimise.jl
+++ b/src/optimise.jl
@@ -1,20 +1,33 @@
+numeric_type(::Interval{T}) where {T} = T
+numeric_type(::IntervalBox{N, T}) where {N, T} = T
 
-"""
-    minimise(f, X, tol=1e-3)
+interval_mid(X::Interval) = Interval(mid(X))
+interval_mid(X::IntervalBox) = IntervalBox(mid(X))
 
-Find the global minimum of the function `f` over the `Interval` or `IntervalBox` `X` using the Moore-Skelboe algorithm.
+"""
+    minimise(f, X, structure = SortedVector, tol=1e-3)
+    or
+    minimise(f, X, structure = HeapedVector, tol=1e-3)
+    or
+    minimise(f, X, tol=1e-3) in this case the default value of "structure" is "HeapedVector"
+
+Find the global minimum of the function `f` over the `Interval` or `IntervalBox` `X`
+using the Moore-Skelboe algorithm. The way in which vector elements are
+kept arranged can be either a heaped array or a sorted array.
+If you not specify any particular strategy to keep vector elements arranged then
+by default heaped array is used.
 
 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(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])
-
+function minimise(f, X::T; structure = HeapedVector, tol=1e-3) where {T}
+    nT = numeric_type(X)
+    
+    # list of boxes with corresponding lower bound, arranged according to selected structure :
+    working = structure([(X, nT(∞))], x->x[2])
     minimizers = T[]
-    global_min = ∞  # upper bound
+    global_min = nT(∞)  # upper bound
 
     num_bisections = 0
 
@@ -27,24 +40,21 @@ function minimise(f, X::T, tol=1e-3) where {T}
         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
+        m = sup(f(interval_mid(X)))   # evaluate at midpoint of current interval
 
         if m < global_min
             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)
+        filter_elements!(working , (X, global_min) )
 
         if diam(X) < tol
             push!(minimizers, X)
-
         else
             X1, X2 = bisect(X)
-            push!( working, (X1, inf(f(X1))), (X2, inf(f(X2))) )
+            push!( working, (X1, inf(f(X1))) )
+            push!( working, (X2, inf(f(X2))) )
             num_bisections += 1
         end
 
@@ -55,8 +65,18 @@ function minimise(f, X::T, tol=1e-3) where {T}
     return Interval(lower_bound, global_min), minimizers
 end
 
-
-function maximise(f, X::T, tol=1e-3) where {T}
-    bound, minimizers = minimise(x -> -f(x), X, tol)
+"""
+    maximise(f, X, structure = SortedVector, tol=1e-3)
+    or
+    maximise(f, X, structure = HeapedVector, tol=1e-3)
+    or
+    maximise(f, X, tol=1e-3) in this case the default value of "structure" is "HeapedVector"
+
+Find the global maximum of the function `f` over the `Interval` or `IntervalBox` `X`
+using the Moore-Skelboe algorithm. See [`minimise`](@ref) for a description
+of the available options.
+"""
+function maximise(f, X::T; structure=HeapedVector, tol=1e-3) where {T}
+    bound, minimizers = minimise(x -> -f(x), X, structure=structure, tol=tol)
     return -bound, minimizers
 end
diff --git a/test/heaped_vector.jl b/test/heaped_vector.jl
new file mode 100644
index 0000000..02e0a38
--- /dev/null
+++ b/test/heaped_vector.jl
@@ -0,0 +1,36 @@
+using IntervalOptimisation
+using Test
+
+using IntervalOptimisation: HeapedVector
+
+@testset "heap" begin
+	
+	@testset "with standard heaping " begin
+
+    h=HeapedVector([9, 7, 5, 2, 3]) 
+
+    @test h.data == [2, 3, 7, 9, 5]
+
+    push!(h, 1)
+    @test h.data == [1, 3, 2, 9, 5, 7]
+
+    popfirst!(h)
+    @test h.data == [2, 3, 7, 9, 5]
+
+  end 
+
+  @testset "with ordering function" begin
+
+    h=HeapedVector( [(9,"f"), (7,"c"), (5,"e"), (2,"d"), (3,"b")] , x->x[2])
+        
+    @test h.data == [(3, "b"), (7, "c"), (5, "e"), (9, "f"), (2, "d")]
+
+    push!(h, (1,"a"))
+    @test h.data == [(1, "a"), (7, "c"), (3, "b"), (9, "f"), (2, "d"), (5, "e")]
+        
+    popfirst!(h)
+    @test h.data == [(3, "b"), (7, "c"), (5, "e"), (9, "f"), (2, "d")]
+
+  end 
+
+end
diff --git a/test/optimise.jl b/test/optimise.jl
index 641c1f2..3377040 100644
--- a/test/optimise.jl
+++ b/test/optimise.jl
@@ -1,15 +1,36 @@
 using IntervalArithmetic, IntervalOptimisation
 using Test
+using IntervalOptimisation: numeric_type
 
 @testset "IntervalOptimisation tests" begin
+    @testset "numeric_type" begin
+        x = -10..10
+        big_x = big(x)
+        @test numeric_type(x) == Float64
+        @test numeric_type(big_x) == BigFloat
+        @test numeric_type(IntervalBox(x, x)) == Float64
+        @test numeric_type(IntervalBox(big_x, big_x)) == BigFloat
+    end
 
-    @testset "Minimise in 1D" begin
+    @testset "Minimise in 1D using default data structure i.e HeapedVector" begin
         global_min, minimisers = minimise(x->x, -10..10)
         @test global_min ⊆ -10 .. -9.999
         @test length(minimisers) == 1
         @test minimisers[1] ⊆ -10 .. -9.999
 
-        global_min, minimisers = minimise(x->x^2, -10..11, 1e-10)
+        # same but with maximise
+        global_max, maximisers = maximise(x->x, -10..10)
+        @test global_max ⊆ 9.999 .. 10
+        @test length(maximisers) == 1
+        @test maximisers[1] ⊆ 9.999 .. 10
+
+        # same but with BigFloats
+        global_min, minimisers = minimise(x->x, -big(10.0)..big(10.0))
+        @test global_min ⊆ -10 .. -9.999
+        @test length(minimisers) == 1
+        @test minimisers[1] ⊆ -10 .. -9.999
+
+        global_min, minimisers = minimise(x->x^2, -10..11, tol = 1e-10)
         @test global_min ⊆ 0..1e-20
         @test length(minimisers) == 1
         @test minimisers[1] ⊆ -0.1..0.1
@@ -17,26 +38,70 @@ using Test
         global_min, minimisers = minimise(x->(x^2-2)^2, -10..11)
         @test global_min ⊆ 0..1e-7
         @test length(minimisers) == 2
-        @test sqrt(2) ∈ minimisers[1]
+        @test sqrt(2) ∈ minimisers[2]
     end
 
+    for Structure in (SortedVector, HeapedVector)
 
-    @testset "Discontinuous function in 1D" begin
+        @testset "Minimise in 1D using SortedVector" begin
+            global_min, minimisers = minimise(x->x, -10..10, structure = Structure)
+            @test global_min ⊆ -10 .. -9.999
+            @test length(minimisers) == 1
+            @test minimisers[1] ⊆ -10 .. -9.999
 
-        H(x) = (sign(x) + 1) / 2   # Heaviside function except at 0, where H(0) = 0.5
+            # same but with maximise
+            global_max, maximisers = maximise(x->x, -10..10, structure = Structure)
+            @test global_max ⊆ 9.999 .. 10
+            @test length(maximisers) == 1
+            @test maximisers[1] ⊆ 9.999 .. 10
 
-        global_min, minimisers = minimise(x -> abs(x) + H(x) - 1, -10..11, 1e-5)
-        @test global_min ⊆ -1 .. -0.9999
-        @test length(minimisers) == 1
-        @test 0 ∈ minimisers[1]
-        @test diam(minimisers[1]) <= 1e-5
-    end
+            # same but with BigFloats
+            global_min, minimisers = minimise(x->x, -big(10.0)..big(10.0), structure = Structure)
+            @test global_min ⊆ -10 .. -9.999
+            @test length(minimisers) == 1
+            @test minimisers[1] ⊆ -10 .. -9.999
 
+            global_min, minimisers = minimise(x->x^2, -10..11, tol=1e-10, structure = Structure)
+            @test global_min ⊆ 0..1e-20
+            @test length(minimisers) == 1
+            @test minimisers[1] ⊆ -0.1..0.1
+
+            global_min, minimisers = minimise(x->(x^2-2)^2, -10..11, structure = Structure)
+            @test global_min ⊆ 0..1e-7
+            @test length(minimisers) == 2
+            @test sqrt(2) ∈ max(minimisers[1], minimisers[2])
+        end
+
+
+        @testset "Discontinuous function in 1D" begin
+
+            H(x) = (sign(x) + 1) / 2   # Heaviside function except at 0, where H(0) = 0.5
+            global_min, minimisers = minimise(x -> abs(x) + H(x) - 1, -10..11, tol=1e-5, structure = Structure)
+            @test global_min ⊆ -1 .. -0.9999
+            @test length(minimisers) == 1
+            @test 0 ∈ minimisers[1]
+            @test diam(minimisers[1]) <= 1e-5
+        end
+
+
+        @testset "Smooth function in 2D" begin
+            global_min, minimisers = minimise( X -> ( (x,y) = X; x^2 + y^2 ), (-10..10) × (-10..10), structure = Structure )
+            @test global_min ⊆ 0..1e-7
+            @test all(X ⊆ (-1e-3..1e3) × (-1e-3..1e-3) for X in minimisers)
+
+            # same but with maximise
+            global_max, maximisers = maximise( X -> ( (x,y) = X; x^2 + y^2 ), (-10..10) × (-10..10), structure = Structure )
+            @test global_max ⊆ 199.9..200
+            m = (9.99..10)
+            @test all(X ⊆ m × m || X ⊆ -m × m || X ⊆ m × -m || X ⊆ -m × -m for X in maximisers)
+
+            # same but with BigFloats
+            global_min, minimisers = minimise( X -> ( (x,y) = X; x^2 + y^2 ), (-big(10.0)..big(10.0)) × (-big(10.0)..big(10.0)), structure = Structure )
+            @test global_min ⊆ 0..1e-7
+            @test all(X ⊆ big(-1e-3..1e3) × big(-1e-3..1e-3) for X in minimisers)
+
+        end
 
-    @testset "Smooth function in 2D" begin
-        global_min, minimisers = minimise( X -> ( (x,y) = X; x^2 + y^2 ), (-10..10) × (-10..10) )
-        @test global_min ⊆ 0..1e-7
-        @test all(X ⊆ (-1e-3..1e3) × (-1e-3..1e-3) for X in minimisers)
     end
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index cadea80..e50374e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,4 +2,5 @@ using IntervalOptimisation
 using Test
 
 include("sorted_vector.jl")
-include("optimise.jl")
+include("heaped_vector.jl")
+include("optimise.jl")
\ No newline at end of file
diff --git a/test/sorted_vector.jl b/test/sorted_vector.jl
index 70413f2..1358455 100644
--- a/test/sorted_vector.jl
+++ b/test/sorted_vector.jl
@@ -1,7 +1,7 @@
 using IntervalOptimisation
 using Test
 
-const SortedVector = IntervalOptimisation.SortedVector
+using IntervalOptimisation: SortedVector
 
 @testset "SortedVector" begin