Global Optimisation

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

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

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

src/optimise1.jl

@@ -0,0 +1,248 @@
+using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures, IntervalConstraintProgramming
+
+import Base.isless
+
+struct IntervalMinima{T<:Real}
+    intval::Interval{T}
+    global_minimum::T
+end
+
+function isless(a::IntervalMinima{T}, b::IntervalMinima{T}) where {T<:Real}
+    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})
+
+    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)
+            continue
+        end
+
+        if p.global_minimum > global_minimum
+            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)
+            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_minimum = f(interval(m)).hi
+
+        if current_minimum < global_minimum
+            global_minimum = current_minimum
+        end
+
+
+        # 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 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))
+            # 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
+    lb = minimum(inf.(f.(arg_minima)))
+
+    return lb..global_minimum, arg_minima
+end
+
+
+
+
+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)
+
+    arg_minima = IntervalBox{N, T}[]
+
+    push!(Q, IntervalBoxMinima(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.intval)
+            continue
+        end
+
+        if p.global_minimum > global_minimum
+            continue
+        end
+
+        current_minimum = f(interval.(mid(p.intval))).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)
+
+        if diam(p.intval) < abstol
+            push!(arg_minima, p.intval)
+
+        else
+            x1, x2 = bisect(X)
+            push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo))
+        end
+    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}
+
+    Q = binary_minheap(IntervalBoxMinima{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, IntervalBoxMinima(x, global_minimum))
+
+    while !isempty(Q)
+
+        p = pop!(Q)
+
+        if isempty(p.intval)
+            continue
+        end
+
+        if p.global_minimum > global_minimum
+            continue
+        end
+
+        # 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(X)
+            push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo))
+        end
+    end
+    lb = minimum(inf.(f.(arg_minima)))
+    return lb..global_minimum, arg_minima
+end