Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
IntervalOptimisation problem benchmark
# 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
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)
----------------------------------------------------------------
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)
----------------------------------------------------------------
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
@ericphanson
Copy link
Author

ericphanson commented Dec 10, 2019

f_vs_fp

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment