IntervalOptimisation problem benchmark

IntervalSpecialFunctions.jl

# This file contains code taken from https://github.com/JuliaIntervals/IntervalSpecialFunctions.jl
# which is available under the following MIT license:

# > Copyright (c) 2018: David Sanders.
# >
# > Permission is hereby granted, free of charge, to any person obtaining a copy
# > of this software and associated documentation files (the "Software"), to deal
# > in the Software without restriction, including without limitation the rights
# > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# > copies of the Software, and to permit persons to whom the Software is
# > furnished to do so, subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be included in all
# > copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.

using IntervalArithmetic
using IntervalArithmetic: @round, big53

import SpecialFunctions: erf, erfc, erfinv, erfcinv

for f in (:erf, :erfc)
    @eval function($f)(x::BigFloat, r::RoundingMode)
        setrounding(BigFloat, r) do
            ($f)(x)
        end
    end
    @eval ($f)(a::Interval{Float64}) = convert(Interval{Float64}, ($f)(big53(a)))
end

function erf(a::Interval{T}) where T
    isempty(a) && return a
    @round( erf(a.lo), erf(a.hi) )
end

function erfc(a::Interval{T}) where T
    isempty(a) && return a
    @round( erfc(a.hi), erfc(a.lo) )
end

lb=0.9_ub=1.5_T=BigFloat.txt

Run with T=BigFloat and tol=1.000000000000000000000000000000000000000000000000000000000000000000000000000007e-05
Run took 13.0 seconds
Run allocated 1492.78 MB
Found maximum:
Interval(0.3468130470973156613776397713904621493388781910433821366603606425257790080787701, 0.3468359173471710235926038024080365516878940622795836636770898231387642913645044)
Found maximizer:
Interval(1.145316457021121284721704055920273340077874567647491198384260011124115290194354, 1.178454475936854016988366656724494955023371118549792529543291444959419755814886)
----------------------------------------------------------------

Run with T=BigFloat and tol=1.000000000000000000000000000000000000000000000000000000000000000000000000000004e-06
Run took 309.92 seconds
Run allocated 18367.57 MB
Found maximum:
Interval(0.3468130470974660080993051747543897684483893939381928306214843688667648775109691, 0.346814489964152696730616727757849799182536517888315716393463469986333908951801)
Found maximizer:
Interval(1.15740957984836702385748742324209963220736599774267720256456516909958175999172, 1.165673402878772514217799088442116636597171320206094700284711092283516209740803)
----------------------------------------------------------------

Run with T=BigFloat and tol=1.000000000000000000000000000000000000000000000000000000000000000000000000000006e-07
Run took 2983.69 seconds
Run allocated 143912.46 MB
Found maximum:
Interval(0.3468130470974666320858031464702992090647364338224003842291257587540849098263324, 0.3468132298368594581764267081993842072687556851927528251870341602452707027458422)
Found maximizer:
Interval(1.160092254271448061513436838888535612967251556384308869749386338427798708567232, 1.16300186634084746632919871169566913121440431548028997289468327571227109704459)
----------------------------------------------------------------

Run with T=BigFloat and tol=1.000000000000000000000000000000000000000000000000000000000000000000000000000008e-08
Run took 21465.83 seconds
Run allocated 1.09770966e6 MB
Found maximum:
Interval(0.3468130470974666515631888744240971365907409439486076174094086281975913392036961, 0.3468130702656272851004971930136917978269813439718454979114632475214503429354717)
Found maximizer:
Interval(1.161018439614974904085783497830374188956854291805951517754182967951945647096709, 1.162036724367762351719120807412913001540914655485962803741830373780133111409785)
----------------------------------------------------------------

lb=0.9_ub=1.5_T=Float64.txt

Run with T=Float64 and tol=1.0e-5
Run took 8.22 seconds
Run allocated 2584.66 MB
Found maximum:
Interval(0.3468130470973151, 0.3468359173471721)
Found maximizer:
Interval(1.1453164570211216, 1.1784544759368543)
----------------------------------------------------------------

Run with T=Float64 and tol=1.0000000000000002e-6
Run took 86.35 seconds
Run allocated 37945.33 MB
Found maximum:
Interval(0.3468130470974652, 0.34681448996415387)
Found maximizer:
Interval(1.1574095798483668, 1.1656734028787725)
----------------------------------------------------------------

Run with T=Float64 and tol=1.0000000000000002e-7
Run took 725.95 seconds
Run allocated 308646.46 MB
Found maximum:
Interval(0.34681304709746597, 0.3468132298368607)
Found maximizer:
Interval(1.160092254271448, 1.1630018663408475)
----------------------------------------------------------------

Run with T=Float64 and tol=1.0000000000000002e-8
Run took 6213.52 seconds
Run allocated 1.95922189e6 MB
Found maximum:
Interval(0.3468130470974662, 0.34681307026562846)
Found maximizer:
Interval(1.1610184396149748, 1.1620367243677623)
----------------------------------------------------------------

optimise.jl

using IntervalArithmetic, IntervalOptimisation, SpecialFunctions
setformat(:full)

include("IntervalSpecialFunctions.jl")

function Phic(x::T) where {T}
    invsqrt2 = inv(sqrt(T(2)))
    erfc(x * invsqrt2)/2
end

function phi(x::T) where {T}
    invsqrt2π = inv(sqrt(2*T(π)))
    exp(-abs2(x)/2) * invsqrt2π
end
M(x) = Phic(x)/phi(x)
f(x) = x - x^2 * M(x)
(m, minimizers), t, bytes, _ = @timed minimise(x -> -f(x), T(.9)..T(1.5); tol=tol)

function run(io, tol, T)
    (m, minimizers), t, bytes, _ = @timed minimise(x -> -f(x), T(.9)..T(1.5); tol=tol)
    X = reduce(∪, minimizers)
    println(io, "Run with T=$T and tol=$tol")
    println(io, "Run took $(round(t, digits=2)) seconds")
    println(io, "Run allocated $(round(bytes*1e-6, digits=2)) MB")
    println(io, "Found maximum:")
    println(io, -m)
    println(io, "Found maximizer:")
    println(io, X)
    println(io, "----------------------------------------------------------------\n")
    flush(io)
    return diam(X)
end

function repeat_run(io, tol, T)
    d = Inf
    while d > 1e-10
        tol /= 10
        d = run(io, tol, T)
    end
end

# T = Float64
T = BigFloat
open("lb=0.9_ub=1.5_T=$T.txt", "w") do io
    tol = T(1.0)/10^4
    repeat_run(io, tol, T)
end