optimization test functions including from Kandasamy et al 2015, "High Dimensional Bayesian Optimization and Bandits"

gp_opt_testfun.jl

using Distributions
using Plots
using StatsFuns

# function from Kandasamy et al

# each Fd is trimodal
struct Fd <: Function
        d::Int64
        v1::Vector{Float64}
        v2::Vector{Float64}
        v3::Vector{Float64}
end


# note: paper does not have negative signs on exponential terms,
# but this blows up otherwise (and anyway, isn't it effectively sampling from a weird GP?)
function (f::Fd)(x)
        h = 0.01 * f.d^(0.1)
        components_before_exp= [(  -(norm(x - f.v1)^2) / (2*h^2) + log(0.1 * (1.0 / h^f.d)))
        ,
        (  -(norm(x - f.v2)^2) / (2*h^2) + log(0.1 * (1.0 / h^f.d)))
        ,
        (  -(norm(x - f.v3)^2) / (2*h^2) + log(0.8 * (1.0 / h^f.d)))]
        return logsumexp(components_before_exp)
end

function sample_fd(d)
        v1 = rand(d)
        v2 = rand(d)
        v3 = rand(d)
        Fd(d, v1, v2, v3)
end

# sum many trimodal functions over random triples of input dimensions
struct TestFunction <: Function
        component_functions::Vector{Fd}
        dimension_choices::Vector{Vector{Int64}}
end

function (f::TestFunction)(x)
        total = 0.0
        dim = length(f.dimension_choices[1])
        for i=1:length(f.component_functions)
                dims = f.dimension_choices[i]
                subset = [x[j] for j in dims]
                total += f.component_functions[i](subset)
        end
        total
end

function sample_full_f(all_d, group_d, M)
        funcs = [sample_fd(group_d) for _=1:M]
        dim_choices = [rand(1:all_d, group_d) for _=1:M]
        TestFunction(funcs, dim_choices)
end

full_f = sample_full_f(2, 2, 5)
contour(0:0.01:1, 0:0.01:1, (x, y) -> full_f([x, y]), fill=true)



struct Shekel <: Function
        a::Array{Float64, 2}
        c::Vector{Float64}
end

function (f::Shekel)(x)
        total = 0.0
        current = 0.0
        for i=1:length(f.c)
                current = f.c[i]
                for j=1:length(x)
                        current += (x[j] - f.a[i, j])^2
                end
                total += 1.0 / current
        end
        total
end

sample_random_shekel(num_maxima, num_dim) = Shekel(rand(num_maxima, num_dim), rand(num_dim))

# sort of from https://www.sfu.ca/~ssurjano/shekel.html
example_shekel = Shekel([4.0 1.0 8.0 6.0 3.0 2.0 5.0 8.0 6.0 7.0;4.0 1.0 8.0 6.0 7.0 9.0 3.0 1.0 2.0 3.6]', [.1,.2,.2,.4,.4,.6,.3,.7,.5,.5])

surface(collect(0:0.1:10), collect(0:0.1:10), (x, y) -> example_shekel([x, y]))
# test out simple optimization techniques
# ReverseDiff and gradient-based stuff from Optim might be more efficient but defeats the purpose

using BlackBoxOptim

opt_f(x) = -full_f(x)
bboptimize(opt_f; SearchRange=(0.0,1.0), NumDimensions = 2)

opt_shekel(x) = -example_shekel(x)
bboptimize(opt_shekel; SearchRange=(0.0, 10.0), NumDimensions = 2)