Adding HeapedVector.jl (#28)

* Adding HeapedVector.jl and adding benchmarks

* trying to import Strategy.jl but cant do it

* Error:In Extending the function of StrategyBase.jl

* Tests and some bugs are fixed

* Benchmarks is shift to another branch

* .DS_store is removedand made few changes in Heaped and SortedVector

* All .DS_Store files are removed

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

src/IntervalOptimisation.jl

@@ -1,18 +1,23 @@
-
-
 module IntervalOptimisation
 
 export minimise, maximise,
        minimize, maximize
+export HeapedVector, SortedVector
+
+include("StrategyBase.jl")
+using .StrategyBase
 
 include("SortedVectors.jl")
 using .SortedVectors
 
-using IntervalArithmetic, IntervalRootFinding
+include("HeapedVectors.jl")
+using .HeapedVectors
 
+using IntervalArithmetic, IntervalRootFinding
 
 include("optimise.jl")
 
+
 const minimize = minimise
 const maximize = maximise
 

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
 

src/StrategyBase.jl

@@ -0,0 +1,7 @@
+module StrategyBase
+
+  export filter_elements!, Strategy
+  abstract type Strategy end
+  function filter_elements!(s::Strategy)end
+
+end

src/optimise.jl

@@ -1,18 +1,22 @@
+"""
+    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"
 
-"""
-    minimise(f, X, tol=1e-3)
-
-Find the global minimum of the function `f` over the `Interval` or `IntervalBox` `X` using the Moore-Skelboe algorithm.
+Find the global minimum of the function `f` over the `Interval` or `IntervalBox` `X` using the Moore-Skelboe algorithm. By specifing the way in which vector element are kept arranged which is in heaped array or in 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}
+
+    # list of boxes with corresponding lower bound, arranged according to selected structure :
+    working = structure([(X, ∞)], x->x[2])
     minimizers = T[]
     global_min = ∞  # upper bound
 
@@ -34,17 +38,14 @@ function minimise(f, X::T, tol=1e-3) where {T}
         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
 
@@ -56,7 +57,7 @@ function minimise(f, X::T, tol=1e-3) where {T}
 end
 
 
-function maximise(f, X::T, tol=1e-3) where {T}
-    bound, minimizers = minimise(x -> -f(x), X, tol)
+function maximise(f, X::T; structure = HeapedVector, tol=1e-3 ) where {T}
+    bound, minimizers = minimise(x -> -f(x), X, structure, tol)
     return -bound, minimizers
 end

test/heaped_vector.jl

@@ -0,0 +1,36 @@
+using IntervalOptimisation
+using Test
+
+const HeapedVector = 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

test/optimise.jl

@@ -1,42 +1,62 @@
-using IntervalArithmetic, IntervalOptimisation
+using IntervalArithmetic, IntervalOptimisation
 using Test
 
 @testset "IntervalOptimisation tests" begin
 
-    @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)
+        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
 
         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 SoretedVector" 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
+            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 -> 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
+            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 "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)
+        @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)
+        end
+
     end
 
-end
+end

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")

test/sorted_vector.jl

@@ -47,4 +47,4 @@ const SortedVector = IntervalOptimisation.SortedVector
         @test v.data == [(1, "a"), (4, "c"), (3, "x"), (5, "y"), (2, "z")]
     end
 
-end
+end