refactoring; use TimerOutputs; repeated function eval option

Author: jbrea

README.md

@@ -30,12 +30,12 @@ modeloptimizer = MLGPOptimizer(every = 50, noisebounds = [-4, 3],       # bounds
                                maxeval = 40)
 opt = BOpt(f,
            model,
-           ExpectedImprovement(),                 # type of acquisition
+           UpperConfidenceBound(),                # type of acquisition
            modeloptimizer,                        
            [-5., -5.], [5., 5.],                  # lowerbounds, upperbounds         
-           maxiterations = 500,
-           sense = Min,
-           gradientfree = false,                  # use gradient information
+           repetitions = 5,                       # evaluate the function for each input 5 times
+           maxiterations = 100,                   # evaluate at 100 input positions
+           sense = Min,                           # minimize the function
            verbosity = Progress)
 
 result = boptimize!(opt)

REQUIRE

@@ -2,6 +2,7 @@ julia 1.0
 GaussianProcesses 0.9.0
 NLopt
 SpecialFunctions
+TimerOutputs
 ElasticPDMats 0.2.1
 ForwardDiff
 DiffResults

examples/1dplot.jl

@@ -0,0 +1,44 @@
+using BayesianOptimization, GaussianProcesses, PGFPlotsX, Random
+push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{fillbetween}")
+
+Random.seed!(13)
+f(x, noisevariance = 1) = .1*sum((x .- 2).^2) + cos(sum(π/2 * x)) + noisevariance * randn()
+model = ElasticGPE(1, mean = MeanConst(0.),
+                   kernel = SEArd([0.], 5.), logNoise = 0.)
+modeloptimizer = MLGPOptimizer(every = 50, noisebounds = [-2., 3],
+                               kernbounds = [[-1, 0], [4, 10]], maxeval = 40)
+opt = BOpt(f, model, ExpectedImprovement(),
+           modeloptimizer, [-5.], [5.],
+           maxiterations = 5, sense = Min, repetitions = 5,
+           acquisitionoptions = (maxeval = 4000, restarts = 50),
+           verbosity = Progress)
+result = boptimize!(opt)
+
+acqfunc = BayesianOptimization.acquisitionfunction(opt.acquisition, model)
+xs = -5:.02:5
+ms, var = predict_f(model, xs)
+sig = sqrt.(var)
+fmax, xmax = BayesianOptimization.acquire_max(opt.opt,
+                        opt.lowerbounds, opt.upperbounds,
+                        opt.acquisitionoptions.restarts)
+@pgf GroupPlot({group_style = {group_size = "1 by 2", vertical_sep = "4mm"},
+                height = "4cm", width = "8cm", legend_pos = "outer north east",
+                legend_style = {draw = "none"}},
+               {legend_columns = 1, xticklabels = ""},
+             Plot({only_marks}, Coordinates(model.x[:], -model.y[:])),
+             "\\addlegendentry{observations}",
+             Plot({no_marks, "red!20", name_path = "A", forget_plot}, Coordinates(xs, -ms .+ sig)),
+             Plot({no_marks, "red!20", name_path = "B", forget_plot}, Coordinates(xs, -ms .- sig)),
+             "\\addplot[red!20] fill between [of=A and B];",
+             "\\addlegendentry{model std}",
+             Plot({no_marks, blue, very_thick}, Coordinates(xs, -ms)),
+             "\\addlegendentry{model mean}",
+             Plot({no_marks, "green!80!black", very_thick}, Coordinates(xs, f.(xs, 0))),
+             "\\addlegendentry{noisefree target}",
+             {height = "3cm", ytick = [0, 4e-2, 8e-2]},
+             Plot({only_marks, red, mark="triangle*",
+                   mark_size = "5pt", mark_options = {rotate = "0"}},
+                  Coordinates(xmax, [fmax])),
+             Plot({no_marks, very_thick}, Coordinates(xs, (x -> acqfunc([x])).(xs))),
+             "\\addlegendentry{next acquisition}",
+             "\\addlegendentry{acquisition function}")

examples/branin_hartmann.jl

@@ -0,0 +1,95 @@
+using BayesianOptimization, GaussianProcesses, Random
+
+# function and regret definitions
+branin(x::Vector; kwargs...) = branin(x[1], x[2]; kwargs...)
+branin(x1, x2; a = 1, b = 5.1/(4π^2), c = 5/π, r = 6, s = 10, t = 1/(8π),
+       noiselevel = 0) =
+    a * (x2 - b*x1^2 + c*x1 - r)^2 + s*(1 - t)*cos(x1) + s + noiselevel * randn()
+
+minima(::typeof(branin)) = [[-π, 12.275], [π, 2.275], [9.42478, 2.475]], 0.397887
+
+hartmann(x; α = [1., 1.2, 3., 3.2], A = [10 3 17 3.5 1.7 8;
+                                         .05 10 17 .1 8 14;
+                                         3 3.5 1.7 10 17 8;
+                                         17 8 .05 10 .1 14],
+         P = 1e-4 * [1312 1696 5569 124 8283 5886;
+                     2329 4135 8307 3736 1004 9991;
+                     2348 1451 3522 2883 3047 6650;
+                     4047 8828 8732 5743 1091 381]) =
+ -sum([α[i] * exp(-sum([A[i, j] * (x[j] - P[i, j])^2 for j in 1:6])) for i in 1:4])
+
+minima(::typeof(hartmann)) = [[.20169, .150011, .476874, .275332, .311652, .6573]], -3.32237
+
+euclidean(x, y) = √sum((x .- y).^2)
+function regret(opt, func)
+    mins, fmin = minima(func)
+    (observed_dist = minimum(map(m -> euclidean(m, opt.observed_optimizer), mins)),
+     observed_regret = abs(opt.observed_optimum - fmin),
+     model_dist = minimum(map(m -> euclidean(m, opt.model_optimizer), mins)),
+     model_regret = abs(Int(opt.sense) * opt.model_optimum - fmin))
+end
+
+# optimize noise-free branin
+
+opt = BOpt(x -> branin(x, noiselevel = 0),
+           ElasticGPE(2, mean = MeanConst(-10.), kernel = SEArd([0., 0.], 5.),
+                      logNoise = -2., capacity = 3000),
+           ExpectedImprovement(),
+           MLGPOptimizer(every = 50, noisebounds = [-4, 3],
+                         kernbounds = [[-1, -1, 0], [4, 4, 10]],
+                         maxeval = 40),
+           [-5., 0.], [10., 15.], maxiterations = 200,
+           sense = Min)
+@time boptimize!(opt)
+regret(opt, branin)
+
+# optimize hartman
+
+opt = BOpt(hartmann,
+           ElasticGPE(6, mean = MeanConst(0.), kernel = Mat52Ard(zeros(6), 0.),
+                      logNoise = -2., capacity = 3000),
+           ExpectedImprovement(),
+           MLGPOptimizer(every = 20, noisebounds = [-4, 3],
+                         kernbounds = [[-3*ones(6); -3], [4*ones(6); 3]],
+                         maxeval = 100),
+           zeros(6), ones(6), maxiterations = 300,
+           sense = Min)
+@time boptimize!(opt)
+regret(opt, hartmann)
+
+# compare model and observation optimizer on noisy branin
+
+all_obs_regs = []
+all_mod_regs = []
+for _ in 1:10
+    opt = BOpt(x -> branin(x, noiselevel = 1),
+               ElasticGPE(2, mean = MeanConst(-50.), kernel = SEArd(zeros(2), 4.),
+                          logNoise = 2., capacity = 1000),
+               UpperConfidenceBound(),
+               MLGPOptimizer(every = 100, noisebounds = [-4, 3],
+                             kernbounds = [[-1, -1, 0], [4, 4, 10]],
+                             f_calls_limit = 40), repetitions = 5,
+               [-5., 0.], [10., 15.], maxiterations = 20,
+               sense = Min, verbosity = Silent)
+    obs_regs = []
+    mod_regs = []
+    for i in 1:10
+        res = boptimize!(opt);
+        _, obs_reg, _, mod_reg = regret(opt, branin)
+        push!(obs_regs, obs_reg)
+        push!(mod_regs, mod_reg)
+    end
+    push!(all_obs_regs, obs_regs)
+    push!(all_mod_regs, mod_regs)
+end
+using PGFPlotsX
+x = (1:length(all_obs_regs[1])) * 20 * 5
+@pgf Axis({ymode = "log", legend_entries = ["average observation regret",
+                                            "average model regret"],
+           legend_columns = 1, legend_pos = "north east",
+           ylabel = "number of observations",
+           title = "model optimizers become more accurate for noisy objectives"},
+          Plot({red, very_thick, no_marks}, Coordinates(x, mean(all_obs_regs))),
+          Plot({blue, very_thick, no_marks}, Coordinates(x, mean(all_mod_regs))),
+          [Plot({red, thin, forget_plot}, Coordinates(x, y)) for y in all_obs_regs]...,
+          [Plot({blue, thin, forget_plot}, Coordinates(x, y)) for y in all_mod_regs]...)

src/BayesianOptimization.jl

@@ -2,65 +2,34 @@ module BayesianOptimization
 import NLopt, GaussianProcesses
 import GaussianProcesses: GPBase, GPE
 import ElasticPDMats: ElasticPDMat
-using ForwardDiff, DiffResults, Random, Dates, SpecialFunctions
+using ForwardDiff, DiffResults, Random, Dates, SpecialFunctions, TimerOutputs
 export BOpt, ExpectedImprovement, ProbabilityOfImprovement, UpperConfidenceBound,
 ThompsonSamplingSimple, MutualInformation, boptimize!, MLGPOptimizer, NoOptimizer,
 Min, Max, BrochuBetaScaling, NoBetaScaling, Silent, Timings, Progress
 
+ENABLE_TIMINGS = true
 
-mutable struct IterationCounter
-    c::Int
-    i::Int
-    N::Int
-end
-isdone(s::IterationCounter) = s.c == s.N
-step!(s::IterationCounter) = (s.c += 1; s.i += 1)
-init!(s::IterationCounter) = s.c = 0
-mutable struct DurationCounter
-    starttime::Float64
-    duration::Float64
-    now::Float64
-    endtime::Float64
-end
-function init!(s::DurationCounter)
-    s.starttime = time()
-    s.endtime = s.starttime + s.duration
-end
-isdone(s::DurationCounter) = (s.now = time()) >= s.endtime
 abstract type ModelOptimizer end
-mutable struct MLGPOptimizer{NT} <: ModelOptimizer
-    i::Int
-    every::Int
-    options::NT
-end
-"""
-    MLGPOptimizer(; every = 10, kwargs...)
-
-Set the GP hyperparameters to the maximum likelihood estimate `every` number of steps.
-"""
-MLGPOptimizer(; every = 10, kwargs...) = MLGPOptimizer(0, every, kwargs.data)
-function optimizemodel!(o::MLGPOptimizer, model::GPBase)
-    if o.i % o.every == 0
-        optimizemodel!(model; o.options...)
-    end
-    o.i += 1
-end
-
 """
 Don't optimize the model ever.
 """
 struct NoModelOptimizer <: ModelOptimizer end
 optimizemodel!(o::NoModelOptimizer, model) = Nothing
 
+include("utils.jl")
+include("acquisitionfunctions.jl")
+include("acquisition.jl")
+include("models/gp.jl")
+
 @enum Sense Min=-1 Max=1
 @enum Verbosity Silent=0 Timings=1 Progress=2
 
-mutable struct BOpt{F,M,A,MO}
+mutable struct BOpt{F,M,A,AO,MO}
     func::F
     sense::Sense
     model::M
     acquisition::A
-    acquisitionoptions::NamedTuple
+    acquisitionoptions::AO
     modeloptimizer::MO
     lowerbounds::Array{Float64, 1}
     upperbounds::Array{Float64, 1}
@@ -73,34 +42,35 @@ mutable struct BOpt{F,M,A,MO}
     opt::NLopt.Opt
     verbosity::Verbosity
     lhs_iterations::Int
+    repetitions::Int
+    timeroutput::TimerOutput
 end
-isdone(o::BOpt) = isdone(o.iterations) || isdone(o.duration)
 """
     BOpt(func, model, acquisition, modeloptimizer, lowerbounds, upperbounds;
               sense = Max, maxiterations = 10^4, maxduration = Inf,
-              acquisitionoptions = NamedTuple(), gradientfree = false,
+              acquisitionoptions = NamedTuple(), repetitions = 1,
               verbosity = Progress, lhs_iterations = 5*length(lowerbounds))
 """
 function BOpt(func, model, acquisition, modeloptimizer, lowerbounds, upperbounds;
               sense = Max, maxiterations = 10^4, maxduration = Inf,
-              acquisitionoptions = NamedTuple(), gradientfree = false,
-              verbosity = Progress, lhs_iterations = 5*length(lowerbounds))
-    if gradientfree
-        default_acquisitionoptions = (method = :GN_DIRECT_L, restarts = 1, maxeval = 500)
-    else
-        default_acquisitionoptions = (method = :LD_LBFGS, restarts = 10, maxeval = 500)
-    end
-    acquisitionoptions = merge(default_acquisitionoptions, acquisitionoptions)
+              acquisitionoptions = NamedTuple(),
+              repetitions = 1, verbosity = Progress,
+              lhs_iterations = 5*length(lowerbounds))
     now = time()
-    BOpt(func, sense, model, acquisition, acquisitionoptions,
+    acquisitionoptions = merge(defaultoptions(typeof(model), typeof(acquisition)),
+               acquisitionoptions)
+    maxiterations < lhs_iterations && @error("maxiterations = $maxiterations < lhs_iterations = $lhs_iterations")
+    BOpt(func, sense, model, acquisition,
+         acquisitionoptions,
          modeloptimizer, lowerbounds, upperbounds,
          -Inf64*Int(sense), Array{Float64}(undef, length(lowerbounds)),
          -Inf64*Int(sense), Array{Float64}(undef, length(lowerbounds)),
          IterationCounter(0, 0, maxiterations),
          DurationCounter(now, maxduration, now, now + maxduration),
          NLopt.Opt(acquisitionoptions.method, length(lowerbounds)),
-         verbosity, lhs_iterations)
+         verbosity, lhs_iterations, repetitions, TimerOutput())
 end
+isdone(o::BOpt) = isdone(o.iterations) || isdone(o.duration)
 import Base: show
 function show(io::IO, mime::MIME"text/plain", o::BOpt)
     println(io, "Bayesian Optimization object\n\nmodel:")
@@ -119,54 +89,52 @@ function show(io::IO, mime::MIME"text/plain", o::BOpt)
     end
 end
 
-sample(lowerbounds, upperbounds) =
-    rand(length(lowerbounds)) .* (upperbounds .- lowerbounds) .+ lowerbounds
-
 function initialise_model!(o)
-    dac = @elapsed x = latin_hypercube_sampling(o.lowerbounds, o.upperbounds,
-                                                o.lhs_iterations)
-    dfunc = @elapsed y = Int(o.sense) .* o.func.([x[:, i] for i in 1:size(x, 2)])
-    o.iterations.i = o.iterations.c = length(y)
-    dmu = @elapsed update!(o.model, x, y)
-    dom = @elapsed optimizemodel!(o.modeloptimizer, o.model)
-    o.opt = nlopt_setup(o.acquisition, o.model, o.lowerbounds, o.upperbounds;
-                        o.acquisitionoptions...)
-    dac, dfunc, dmu, dom
+    @mytimeit o.timeroutput "acquisition" x = latin_hypercube_sampling(o.lowerbounds, o.upperbounds, o.lhs_iterations)
+    y = Float64[]
+    for i in 1:size(x, 2)
+        for j in 1:o.repetitions
+            @mytimeit o.timeroutput "function evaluation" push!(y, Int(o.sense) * o.func(x[:, i]))
+        end
+    end
+    o.iterations.i = o.iterations.c = length(y)/o.repetitions
+    @mytimeit o.timeroutput "model update" update!(o.model,
+                           hcat(hcat([fill(x[:, i], o.repetitions) for i in 1:size(x, 2)]...)...),
+                           y)
+    @mytimeit o.timeroutput "model hyperparameter optimization" optimizemodel!(o.modeloptimizer, o.model)
+    o.opt = nlopt_setup(o.acquisition, o.model, o.lowerbounds, o.upperbounds,
+                        o.acquisitionoptions)
 end
 """
     boptimize!(o::BOpt)
 """
 function boptimize!(o::BOpt)
     init!(o.duration)
     init!(o.iterations)
-    dfunc = dom = dac = dmu = 0.
-    if o.iterations.i == 0 dac, dfunc, dmu, dom = initialise_model!(o) end
+    reset_timer!(o.timeroutput)
+    o.iterations.i == 0 && initialise_model!(o)
     while !isdone(o)
         o.verbosity >= Progress && @info("$(now())\titeration: $(o.iterations.i)\tcurrent optimum: $(o.observed_optimum)")
         setparams!(o.acquisition, o.model)
-        dac += @elapsed f, x = acquire_max(o.opt, o.lowerbounds, o.upperbounds,
-                                           o.acquisitionoptions.restarts)
-        dfunc += @elapsed y = Int(o.sense) * o.func(x)
+        @mytimeit o.timeroutput "acquisition" f, x = acquire_max(o.opt, o.lowerbounds, o.upperbounds, o.acquisitionoptions.restarts)
+        ys = Float64[]
         step!(o.iterations)
-        if y > Int(o.sense) * o.observed_optimum
-            o.observed_optimum = Int(o.sense) * y
-            o.observed_optimizer = x
+        for _ in 1:o.repetitions
+            @mytimeit o.timeroutput "function evaluation" y = Int(o.sense) * o.func(x)
+            push!(ys, y)
+            if y > Int(o.sense) * o.observed_optimum
+                o.observed_optimum = Int(o.sense) * y
+                o.observed_optimizer = x
+            end
         end
-        dmu += @elapsed update!(o.model, x, y)
-        dom += @elapsed optimizemodel!(o.modeloptimizer, o.model)
+        @mytimeit o.timeroutput "model update" update!(o.model, hcat(fill(x, o.repetitions)...), ys)
+        @mytimeit o.timeroutput "model hyperparameter optimization" optimizemodel!(o.modeloptimizer, o.model)
     end
+     @mytimeit o.timeroutput "acquisition" o.model_optimum, o.model_optimizer = acquire_model_max(o)
     o.duration.now = time()
-    o.verbosity >= Timings && @info("time spent for:
-        function evaluation \t $dfunc s
-        model update \t\t $dmu s
-        model optimization \t $dom s
-        acquisition \t\t $dac s")
-    dom += @elapsed o.model_optimum, o.model_optimizer = acquire_model_max(o, restarts = 10, maxeval = 2000)
+    o.verbosity >= Timings && @info(o.timeroutput)
     (observerd_optimum = o.observed_optimum, observed_optimizer = o.observed_optimizer,
      model_optimum = Int(o.sense) * o.model_optimum, model_optimizer = o.model_optimizer)
 end
 
-include("acquisition.jl")
-include("models/gp.jl")
-
 end # module