From 2499f85bdb53d950ed69886b1bf89382a6e39f3e Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Mon, 18 Jun 2018 23:19:16 +0530 Subject: [PATCH 01/12] Fix deprecated syntax for inner constructor of SortedVector --- src/SortedVectors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl index 37f647c..c775811 100644 --- a/src/SortedVectors.jl +++ b/src/SortedVectors.jl @@ -12,7 +12,7 @@ immutable SortedVector{T, F<:Function} data::Vector{T} by::F - function SortedVector(data::Vector{T}, by::F) + function SortedVector{T,F}(data::Vector{T}, by::F) where {T,F} new(sort(data), by) end end From 40997354046d9aebb0b155e35cbd5162b763e014 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 19 Jun 2018 22:01:09 +0530 Subject: [PATCH 02/12] Initial implementation of global optimisation --- src/optimise1.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 src/optimise1.jl diff --git a/src/optimise1.jl b/src/optimise1.jl new file mode 100644 index 0000000..c042398 --- /dev/null +++ b/src/optimise1.jl @@ -0,0 +1,52 @@ +using IntervalArithmetic, DataStructures + +import Base.isless + +struct IntervalMinima{T<:Real} + intval::Interval{T} + minima::T +end + +function isless(a::IntervalMinima{T}, b::IntervalMinima{T}) where {T<:Real} + return isless(a.minima, b.minima) +end + +function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} + + Q = binary_minheap(IntervalMinima{T}) + + minima = f(interval(mid(x))).hi + arg_minima = Interval{T}[] + + push!(Q, IntervalMinima(x, minima)) + + while !isempty(Q) + + p = pop!(Q) + + if p.minima > minima + continue + end + # + # if 0 ∉ ForwardDiff.derivative(f, p.intval) + # continue + # end + + # current_minima = f(p.intval).lo + + current_minima = f(interval(mid(p.intval))).hi + + if current_minima < minima + minima = current_minima + end + + if diam(p.intval) < abstol + push!(arg_minima, p.intval) + else + x1, x2 = bisect(p.intval) + push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) + end + end + + return minima, arg_minima +end From a9c1eab527d35dc1333ccb3dc5dff1b7f0db2f6c Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 19 Jul 2018 19:31:23 +0530 Subject: [PATCH 03/12] Add derivative contractors --- src/optimise1.jl | 126 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 125 insertions(+), 1 deletion(-) diff --git a/src/optimise1.jl b/src/optimise1.jl index c042398..96c7b71 100644 --- a/src/optimise1.jl +++ b/src/optimise1.jl @@ -1,4 +1,4 @@ -using IntervalArithmetic, DataStructures +using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures import Base.isless @@ -15,6 +15,83 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where Q = binary_minheap(IntervalMinima{T}) + minima = f(interval(mid(x))).hi + + arg_minima = Interval{T}[] + + push!(Q, IntervalMinima(x, minima)) + + while !isempty(Q) + + p = pop!(Q) + + if isempty(p.intval) + continue + end + + if p.minima > minima + continue + end + + # deriv = ForwardDiff.derivative(f, p.intval) + # if 0 ∉ deriv + # continue + # end + + doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval) + if doublederiv < 0 + continue + end + + m = mid(p.intval) + + current_minima = f(interval(m)).hi + + if current_minima < minima + minima = current_minima + end + + + # # Contractor 1 + # x = m .+ extended_div((interval(-∞, minima) - f(m)), deriv) + # + # x = x .∩ p.intval + + # Contractor 2 + @show p.intval + @show (doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima)) + @show (m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima))) + x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima)) + + x = x .∩ p.intval + + if diam(p.intval) < abstol + push!(arg_minima, p.intval) + else + if length(x) == 1 + x1, x2 = bisect(x[1]) + push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) + else + push!(Q, IntervalMinima.(x, inf.(f.(x)))...) + end + # if isempty(x[2]) + # x1, x2 = bisect(x[1]) + # push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) + # else + # x1, x2, x3 = x + # push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo), IntervalMinima(x3, f(x3).lo)) + # end + end + end + + return minima, arg_minima +end + + +function minimise1d_noderiv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} + + Q = binary_minheap(IntervalMinima{T}) + minima = f(interval(mid(x))).hi arg_minima = Interval{T}[] @@ -50,3 +127,50 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where return minima, arg_minima end + + + +struct IntervalBoxMinima{N, T<:Real} + intval::IntervalBox{N, T} + minima::T +end + +function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real} + return isless(a.minima, b.minima) +end + +function minimisend(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real} + + Q = binary_minheap(IntervalBoxMinima{N, T}) + + minima = f(x).lo + arg_minima = IntervalBox{N, T}[] + + push!(Q, IntervalBoxMinima(x, minima)) + + while !isempty(Q) + + p = pop!(Q) + + if p.minima > minima + continue + end + + if 0 .∉ ForwardDiff.gradient(f, p.intval) + continue + end + + if p.minima < minima + minima = p.minima + end + + if diam(p.intval) < abstol + push!(arg_minima, p.intval) + else + x1, x2 = bisect(p.intval) + push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo)) + end + end + + return minima, arg_minima +end From ae5554199709243e30b7a07efd9690b2e3b281ea Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Mon, 23 Jul 2018 20:39:30 +0530 Subject: [PATCH 04/12] Add icp contractors --- examples/optimise1.jl | 54 ++++++++++ src/IntervalOptimisation.jl | 9 +- src/optimise1.jl | 198 ++++++++++++++++++++++++------------ 3 files changed, 195 insertions(+), 66 deletions(-) create mode 100644 examples/optimise1.jl diff --git a/examples/optimise1.jl b/examples/optimise1.jl new file mode 100644 index 0000000..edd8891 --- /dev/null +++ b/examples/optimise1.jl @@ -0,0 +1,54 @@ +using IntervalOptimisation, IntervalArithmetic, IntervalConstraintProgramming + + +# Unconstrained Optimisation +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 +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..7cee90e 100644 --- a/src/IntervalOptimisation.jl +++ b/src/IntervalOptimisation.jl @@ -3,12 +3,14 @@ module IntervalOptimisation export minimise, maximise, - minimize, maximize + minimize, maximize, + minimise1d, minimise1d_deriv, + minimise_icp, minimise_icp_constrained include("SortedVectors.jl") using .SortedVectors -using IntervalArithmetic, IntervalRootFinding +using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, StaticArrays, ForwardDiff if !isdefined(:sup) const sup = supremum @@ -19,8 +21,9 @@ if !isdefined(:inf) end include("optimise.jl") +include("optimise1.jl") const minimize = minimise -const maximize = maximise +const maximize = maximise end diff --git a/src/optimise1.jl b/src/optimise1.jl index 96c7b71..306f8a1 100644 --- a/src/optimise1.jl +++ b/src/optimise1.jl @@ -1,25 +1,25 @@ -using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures +using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures, IntervalConstraintProgramming import Base.isless struct IntervalMinima{T<:Real} intval::Interval{T} - minima::T + global_minimum::T end function isless(a::IntervalMinima{T}, b::IntervalMinima{T}) where {T<:Real} - return isless(a.minima, b.minima) + return isless(a.global_minimum, b.global_minimum) end function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} Q = binary_minheap(IntervalMinima{T}) - minima = f(interval(mid(x))).hi + global_minimum = f(interval(mid(x))).hi arg_minima = Interval{T}[] - push!(Q, IntervalMinima(x, minima)) + push!(Q, IntervalMinima(x, global_minimum)) while !isempty(Q) @@ -29,51 +29,90 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where continue end - if p.minima > minima + if p.global_minimum > global_minimum continue end - # deriv = ForwardDiff.derivative(f, p.intval) - # if 0 ∉ deriv - # continue - # end + current_minimum = f(interval(mid(p.intval))).hi + + if current_minimum < global_minimum + global_minimum = current_minimum + end + + if diam(p.intval) < abstol + push!(arg_minima, p.intval) + else + x1, x2 = bisect(p.intval) + push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) + end + end + lb = minimum(inf.(f.(arg_minima))) + + return lb..global_minimum, arg_minima +end + +function minimise1d_deriv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} + + Q = binary_minheap(IntervalMinima{T}) + + global_minimum = f(interval(mid(x))).hi + + arg_minima = Interval{T}[] + + push!(Q, IntervalMinima(x, global_minimum)) + + while !isempty(Q) - doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval) - if doublederiv < 0 + p = pop!(Q) + + if isempty(p.intval) + continue + end + + if p.global_minimum > global_minimum continue end + deriv = ForwardDiff.derivative(f, p.intval) + if 0 ∉ deriv + continue + end + # Second derivative contractor + # doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval) + # if doublederiv < 0 + # continue + # end + m = mid(p.intval) - current_minima = f(interval(m)).hi + current_minimum = f(interval(m)).hi - if current_minima < minima - minima = current_minima + if current_minimum < global_minimum + global_minimum = current_minimum end - # # Contractor 1 - # x = m .+ extended_div((interval(-∞, minima) - f(m)), deriv) - # - # x = x .∩ p.intval - - # Contractor 2 - @show p.intval - @show (doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima)) - @show (m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima))) - x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima)) + # Contractor 1 + x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv) x = x .∩ p.intval + # # Contractor 2 (Second derivative, expanding more on the Taylor Series, not beneficial in practice) + # x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, global_minimum)) + # + # x = x .∩ p.intval + if diam(p.intval) < abstol push!(arg_minima, p.intval) else - if length(x) == 1 + if isempty(x[2]) x1, x2 = bisect(x[1]) push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) else push!(Q, IntervalMinima.(x, inf.(f.(x)))...) end + + # Second Deriv contractor # if isempty(x[2]) # x1, x2 = bisect(x[1]) # push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) @@ -83,94 +122,127 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where # end end end + lb = minimum(inf.(f.(arg_minima))) - return minima, arg_minima + return lb..global_minimum, arg_minima end -function minimise1d_noderiv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} - Q = binary_minheap(IntervalMinima{T}) - minima = f(interval(mid(x))).hi - arg_minima = Interval{T}[] +struct IntervalBoxMinima{N, T<:Real} + intval::IntervalBox{N, T} + global_minimum::T +end + +struct constraint{T<:Real} + f::Function + c::Interval{T} +end + +function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real} + return isless(a.global_minimum, b.global_minimum) +end + +function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real} + + Q = binary_minheap(IntervalBoxMinima{N, T}) + + global_minimum = ∞ + + x = icp(f, x, -∞..global_minimum) - push!(Q, IntervalMinima(x, minima)) + arg_minima = IntervalBox{N, T}[] + + push!(Q, IntervalBoxMinima(x, global_minimum)) while !isempty(Q) p = pop!(Q) - if p.minima > minima + if isempty(p.intval) continue end - # - # if 0 ∉ ForwardDiff.derivative(f, p.intval) - # continue - # end - # current_minima = f(p.intval).lo + if p.global_minimum > global_minimum + continue + end - current_minima = f(interval(mid(p.intval))).hi + current_minimum = f(interval.(mid(p.intval))).hi - if current_minima < minima - minima = current_minima + if current_minimum < global_minimum + global_minimum = current_minimum end + # if all(0 .∉ ForwardDiff.gradient(f, p.intval.v)) + # continue + # end + + X = icp(f, p.intval, -∞..global_minimum) if diam(p.intval) < abstol push!(arg_minima, p.intval) + else - x1, x2 = bisect(p.intval) - push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) + x1, x2 = bisect(X) + push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo)) end end - return minima, arg_minima -end + 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} -struct IntervalBoxMinima{N, T<:Real} - intval::IntervalBox{N, T} - minima::T -end + Q = binary_minheap(IntervalBoxMinima{N, T}) -function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real} - return isless(a.minima, b.minima) -end + global_minimum = ∞ -function minimisend(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real} + for t in constraints + x = icp(t.f, x, t.c) + end - Q = binary_minheap(IntervalBoxMinima{N, T}) + x = icp(f, x, -∞..global_minimum) - minima = f(x).lo arg_minima = IntervalBox{N, T}[] - push!(Q, IntervalBoxMinima(x, minima)) + push!(Q, IntervalBoxMinima(x, global_minimum)) while !isempty(Q) p = pop!(Q) - if p.minima > minima + if isempty(p.intval) continue end - if 0 .∉ ForwardDiff.gradient(f, p.intval) + if p.global_minimum > global_minimum continue end - if p.minima < minima - minima = p.minima + # current_minimum = f(interval.(mid(p.intval))).hi + current_minimum = f(p.intval).hi + + if current_minimum < global_minimum + global_minimum = current_minimum + end + # if 0 .∉ ForwardDiff.gradient(f, p.intval.v) + # continue + # end + X = icp(f, p.intval, -∞..global_minimum) + + for t in constraints + X = icp(t.f, X, t.c) end if diam(p.intval) < abstol push!(arg_minima, p.intval) else - x1, x2 = bisect(p.intval) + x1, x2 = bisect(X) push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo)) end end - - return minima, arg_minima + lb = minimum(inf.(f.(arg_minima))) + return lb..global_minimum, arg_minima end From b3ef8117ead5442bbee79b2c254755b45a5dbb35 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Thu, 2 Aug 2018 06:12:39 +0530 Subject: [PATCH 05/12] Incorporate changes from review comments --- REQUIRE | 2 + examples/optimise1.jl | 2 + src/optimise1.jl | 186 +++++++++++++++--------------------------- 3 files changed, 70 insertions(+), 120 deletions(-) diff --git a/REQUIRE b/REQUIRE index 4f32cce..2f74c15 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,3 +2,5 @@ julia 0.5 IntervalArithmetic 0.9.1 IntervalRootFinding 0.1 +IntervalConstraintProgramming +DataStructures diff --git a/examples/optimise1.jl b/examples/optimise1.jl index edd8891..13892cb 100644 --- a/examples/optimise1.jl +++ b/examples/optimise1.jl @@ -2,6 +2,7 @@ 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) @@ -29,6 +30,7 @@ minimise_icp(f, X) # 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) diff --git a/src/optimise1.jl b/src/optimise1.jl index 306f8a1..3dc9f77 100644 --- a/src/optimise1.jl +++ b/src/optimise1.jl @@ -1,90 +1,48 @@ -using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures, IntervalConstraintProgramming - -import Base.isless - -struct IntervalMinima{T<:Real} - intval::Interval{T} - global_minimum::T +struct IntervalMinimum{T<:Real} + interval::Interval{T} + minimum::T end -function isless(a::IntervalMinima{T}, b::IntervalMinima{T}) where {T<:Real} - return isless(a.global_minimum, b.global_minimum) -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) where {T<:Real} +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(IntervalMinima{T}) + Q = binary_minheap(IntervalMinimum{T}) global_minimum = f(interval(mid(x))).hi arg_minima = Interval{T}[] - push!(Q, IntervalMinima(x, global_minimum)) + push!(Q, IntervalMinimum(x, global_minimum)) while !isempty(Q) p = pop!(Q) - if isempty(p.intval) - continue - end - - if p.global_minimum > global_minimum + if isempty(p.interval) continue end - current_minimum = f(interval(mid(p.intval))).hi - - if current_minimum < global_minimum - global_minimum = current_minimum - end - - if diam(p.intval) < abstol - push!(arg_minima, p.intval) - else - x1, x2 = bisect(p.intval) - push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) - end - end - lb = minimum(inf.(f.(arg_minima))) - - return lb..global_minimum, arg_minima -end - -function minimise1d_deriv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real} - - Q = binary_minheap(IntervalMinima{T}) - - global_minimum = f(interval(mid(x))).hi - - arg_minima = Interval{T}[] - - push!(Q, IntervalMinima(x, global_minimum)) - - while !isempty(Q) - - p = pop!(Q) - - if isempty(p.intval) + if p.minimum > global_minimum continue end - if p.global_minimum > global_minimum - continue + if use_deriv + deriv = ForwardDiff.derivative(f, p.interval) + if 0 ∉ deriv + continue + end end - deriv = ForwardDiff.derivative(f, p.intval) - if 0 ∉ deriv - continue - end # Second derivative contractor - # doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval) - # if doublederiv < 0 - # continue - # end - - m = mid(p.intval) + 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 @@ -93,33 +51,24 @@ function minimise1d_deriv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) # Contractor 1 - x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv) - - x = x .∩ p.intval + if use_deriv + x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv) - # # Contractor 2 (Second derivative, expanding more on the Taylor Series, not beneficial in practice) - # x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, global_minimum)) - # - # x = x .∩ p.intval + x = x .∩ p.interval + end - if diam(p.intval) < abstol - push!(arg_minima, p.intval) + if diam(p.interval) < abstol + push!(arg_minima, p.interval) else - if isempty(x[2]) + + if use_deriv && isempty(x[2]) x1, x2 = bisect(x[1]) - push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) - else - push!(Q, IntervalMinima.(x, inf.(f.(x)))...) + push!(Q, IntervalMinimum(x1, f(x1).lo), IntervalMinimum(x2, f(x2).lo)) + continue end - # Second Deriv contractor - # if isempty(x[2]) - # x1, x2 = bisect(x[1]) - # push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo)) - # else - # x1, x2, x3 = x - # push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo), IntervalMinima(x3, f(x3).lo)) - # 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))) @@ -128,25 +77,27 @@ function minimise1d_deriv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) end - - -struct IntervalBoxMinima{N, T<:Real} - intval::IntervalBox{N, T} - global_minimum::T +struct IntervalBoxMinimum{N, T<:Real} + interval::IntervalBox{N, T} + minimum::T end -struct constraint{T<:Real} - f::Function +""" +Datatype to provide constraints to Global Optimisation such as: +``` +Constraint(x->(x^2 - 10), -∞..1) +``` +""" +struct Constraint{F, T<:Real} + f::F c::Interval{T} end -function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real} - return isless(a.global_minimum, b.global_minimum) -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(IntervalBoxMinima{N, T}) + Q = binary_minheap(IntervalBoxMinimum{N, T}) global_minimum = ∞ @@ -154,37 +105,34 @@ function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e- arg_minima = IntervalBox{N, T}[] - push!(Q, IntervalBoxMinima(x, global_minimum)) + push!(Q, IntervalBoxMinimum(x, global_minimum)) while !isempty(Q) p = pop!(Q) - if isempty(p.intval) + if isempty(p.interval) continue end - if p.global_minimum > global_minimum + if p.minimum > global_minimum continue end - current_minimum = f(interval.(mid(p.intval))).hi + current_minimum = f(interval.(mid(p.interval))).hi if current_minimum < global_minimum global_minimum = current_minimum end - # if all(0 .∉ ForwardDiff.gradient(f, p.intval.v)) - # continue - # end - X = icp(f, p.intval, -∞..global_minimum) + X = icp(f, p.interval, -∞..global_minimum) - if diam(p.intval) < abstol - push!(arg_minima, p.intval) + if diam(p.interval) < abstol + push!(arg_minima, p.interval) else x1, x2 = bisect(X) - push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo)) + push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo)) end end @@ -193,9 +141,9 @@ function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e- 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} +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(IntervalBoxMinima{N, T}) + Q = binary_minheap(IntervalBoxMinimum{N, T}) global_minimum = ∞ @@ -207,40 +155,38 @@ function minimise_icp_constrained(f::Function, x::IntervalBox{N, T}, constraints arg_minima = IntervalBox{N, T}[] - push!(Q, IntervalBoxMinima(x, global_minimum)) + push!(Q, IntervalBoxMinimum(x, global_minimum)) while !isempty(Q) p = pop!(Q) - if isempty(p.intval) + if isempty(p.interval) continue end - if p.global_minimum > global_minimum + if p.minimum > global_minimum continue end - # current_minimum = f(interval.(mid(p.intval))).hi - current_minimum = f(p.intval).hi + # current_minimum = f(interval.(mid(p.interval))).hi + current_minimum = f(p.interval).hi if current_minimum < global_minimum global_minimum = current_minimum end - # if 0 .∉ ForwardDiff.gradient(f, p.intval.v) - # continue - # end - X = icp(f, p.intval, -∞..global_minimum) + + X = icp(f, p.interval, -∞..global_minimum) for t in constraints X = icp(t.f, X, t.c) end - if diam(p.intval) < abstol - push!(arg_minima, p.intval) + if diam(p.interval) < abstol + push!(arg_minima, p.interval) else x1, x2 = bisect(X) - push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo)) + push!(Q, IntervalBoxMinimum(x1, f(x1).lo), IntervalBoxMinimum(x2, f(x2).lo)) end end lb = minimum(inf.(f.(arg_minima))) From a5b6b5b120bfdca6ee4e0511b9c0724043583814 Mon Sep 17 00:00:00 2001 From: Eeshan Gupta Date: Tue, 7 Aug 2018 15:58:16 +0530 Subject: [PATCH 06/12] Remove parameterization of F --- src/optimise1.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/optimise1.jl b/src/optimise1.jl index 3dc9f77..592733f 100644 --- a/src/optimise1.jl +++ b/src/optimise1.jl @@ -88,8 +88,8 @@ Datatype to provide constraints to Global Optimisation such as: Constraint(x->(x^2 - 10), -∞..1) ``` """ -struct Constraint{F, T<:Real} - f::F +struct Constraint{T<:Real} + f::Function c::Interval{T} end From ca21191c448ecdee6760cbbb00c95b868f09ea60 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Wed, 5 Sep 2018 17:12:42 -0500 Subject: [PATCH 07/12] Deprecations for Julia 1.0 --- REQUIRE | 7 ++++--- src/IntervalOptimisation.jl | 9 +-------- src/SortedVectors.jl | 12 ++++++------ src/optimise.jl | 6 +++--- 4 files changed, 14 insertions(+), 20 deletions(-) diff --git a/REQUIRE b/REQUIRE index 2f74c15..bfd619e 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,6 +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/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl index 7cee90e..b28aa60 100644 --- a/src/IntervalOptimisation.jl +++ b/src/IntervalOptimisation.jl @@ -10,15 +10,8 @@ export minimise, maximise, include("SortedVectors.jl") using .SortedVectors -using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, StaticArrays, ForwardDiff +using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, ForwardDiff -if !isdefined(:sup) - const sup = supremum -end - -if !isdefined(:inf) - const inf = infimum -end include("optimise.jl") include("optimise1.jl") diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl index c775811..6ca82ea 100644 --- a/src/SortedVectors.jl +++ b/src/SortedVectors.jl @@ -8,18 +8,18 @@ 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{T,F}(data::Vector{T}, by::F) where {T,F} - new(sort(data), by) + function SortedVector(data::Vector{T}, by::F) where {T,F} + new{T,F}(sort(data), 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 diff --git a/src/optimise.jl b/src/optimise.jl index 337b8f4..fb7eb65 100644 --- a/src/optimise.jl +++ b/src/optimise.jl @@ -4,11 +4,11 @@ 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]) @@ -56,7 +56,7 @@ function minimise{T}(f, X::T, tol=1e-3) 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 From 1534646f02ec4d968df0d116ff1a5d625c8864d5 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Wed, 5 Sep 2018 17:14:53 -0500 Subject: [PATCH 08/12] shift! -> popfirst! --- src/SortedVectors.jl | 4 ++-- src/optimise.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl index 6ca82ea..714a70e 100644 --- a/src/SortedVectors.jl +++ b/src/SortedVectors.jl @@ -1,7 +1,7 @@ module SortedVectors import Base: getindex, length, push!, isempty, - pop!, shift!, resize! + pop!, resize!, popfirst! export SortedVector @@ -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 fb7eb65..4629811 100644 --- a/src/optimise.jl +++ b/src/optimise.jl @@ -20,7 +20,7 @@ function minimise(f, X::T, tol=1e-3) where {T} while !isempty(working) - (X, X_min) = shift!(working) + (X, X_min) = popfirst!(working) if X_min > global_min # X_min == inf(f(X)) continue From 3368a0426995ca4f97751ad424f64c8b8f1b624a Mon Sep 17 00:00:00 2001 From: David Sanders Date: Wed, 5 Sep 2018 19:29:47 -0500 Subject: [PATCH 09/12] Update .travis.yml and appveyor.yml --- .travis.yml | 4 ++-- appveyor.yml | 41 +++++++++++++++++++++++++++++------------ 2 files changed, 31 insertions(+), 14 deletions(-) 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/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%" From 09b8c689ede531762926e2c48b8f6bf4088133cf Mon Sep 17 00:00:00 2001 From: David Sanders Date: Sun, 23 Sep 2018 20:27:20 -0500 Subject: [PATCH 10/12] Add Griewank example --- examples/griewank.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 examples/griewank.jl 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 From d46cb15101d2b27d0e4c162d5c066f57f74bd5d9 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Sun, 7 Oct 2018 22:51:57 -0500 Subject: [PATCH 11/12] Change to using a PriorityQueue; 5x speedup for hard problems --- src/IntervalOptimisation.jl | 6 +++-- src/SortedVectors.jl | 2 +- src/optimise.jl | 48 +++++++++++++++++++++++++++---------- src/optimise1.jl | 4 +++- 4 files changed, 44 insertions(+), 16 deletions(-) diff --git a/src/IntervalOptimisation.jl b/src/IntervalOptimisation.jl index b28aa60..654aa61 100644 --- a/src/IntervalOptimisation.jl +++ b/src/IntervalOptimisation.jl @@ -2,13 +2,15 @@ module IntervalOptimisation +using DataStructures + export minimise, maximise, minimize, maximize, minimise1d, minimise1d_deriv, minimise_icp, minimise_icp_constrained -include("SortedVectors.jl") -using .SortedVectors +# include("SortedVectors.jl") +# using .SortedVectors using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, ForwardDiff diff --git a/src/SortedVectors.jl b/src/SortedVectors.jl index 714a70e..e950342 100644 --- a/src/SortedVectors.jl +++ b/src/SortedVectors.jl @@ -13,7 +13,7 @@ struct SortedVector{T, F<:Function} by::F function SortedVector(data::Vector{T}, by::F) where {T,F} - new{T,F}(sort(data), by) + new{T,F}(sort(data, by=by), by) end end diff --git a/src/optimise.jl b/src/optimise.jl index 4629811..4c3f1f4 100644 --- a/src/optimise.jl +++ b/src/optimise.jl @@ -11,23 +11,41 @@ Returns an interval containing the global minimum, and a list of boxes that cont 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) = popfirst!(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 + + #@show f(interval.(mid(XX))) + + m = sup(f(interval.(mid(XX)))) # evaluate at midpoint of current interval if m < global_min global_min = m @@ -36,21 +54,27 @@ function minimise(f, X::T, tol=1e-3) where {T} # 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) + enqueue!(working, X1 => inf(f(X1)) ) + enqueue!(working, X2 => inf(f(X2)) ) 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 end diff --git a/src/optimise1.jl b/src/optimise1.jl index 592733f..8201eb3 100644 --- a/src/optimise1.jl +++ b/src/optimise1.jl @@ -119,7 +119,7 @@ function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e- continue end - current_minimum = f(interval.(mid(p.interval))).hi + current_minimum = sup(f(interval.(mid(p.interval)))) if current_minimum < global_minimum global_minimum = current_minimum @@ -136,6 +136,8 @@ function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e- end end + #@show arg_minima + lb = minimum(inf.(f.(arg_minima))) return lb..global_minimum, arg_minima From ff220147a2728a64980fd0e5989b91501f512c02 Mon Sep 17 00:00:00 2001 From: David Sanders Date: Sat, 27 Oct 2018 17:29:16 -0400 Subject: [PATCH 12/12] Change < to <= --- src/optimise.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/optimise.jl b/src/optimise.jl index 4c3f1f4..9519c8f 100644 --- a/src/optimise.jl +++ b/src/optimise.jl @@ -47,7 +47,7 @@ function _minimise(f, working::PriorityQueue{K,V}, tol) where {K,V} m = sup(f(interval.(mid(XX)))) # evaluate at midpoint of current interval - if m < global_min + if m < global_min && m > -∞ global_min = m end @@ -62,8 +62,13 @@ function _minimise(f, working::PriorityQueue{K,V}, tol) where {K,V} else X1, X2 = bisect(XX) - enqueue!(working, X1 => inf(f(X1)) ) - enqueue!(working, X2 => inf(f(X2)) ) + + 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 @@ -72,11 +77,11 @@ function _minimise(f, working::PriorityQueue{K,V}, tol) where {K,V} # @show global_min # @show minimizers - new_minimizers = [minimizer for minimizer in minimizers if !(f(minimizer).lo > global_min)] + 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