Commit 00f904b0 authored by Antoine Cully's avatar Antoine Cully
Browse files

first version of the benchmark integrated in the wscript

parent 971be432
#define SKIP_TRICKS
#include "testfunctions.hpp"
#include <chrono>
#include <fstream>
#include <string>
template <typename Function>
void benchmark(const bopt_params& par, const std::string& name)
{
auto t1 = std::chrono::steady_clock::now();
Benchmark<Function> benchmark(par);
vectord result(Function::dim_in);
benchmark.optimize(result);
auto time_running = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - t1).count();
std::cout.precision(17);
std::cout << std::endl;
auto best = benchmark.evaluateSample(result);
double accuracy = benchmark.accuracy(best);
std::cout << name << std::endl;
std::cout << "Result: " << std::fixed << result << " -> " << best << std::endl;
std::cout << "Smallest difference: " << accuracy << std::endl;
std::cout << "Time running: " << time_running << "ms" << std::endl << std::endl;
std::ofstream res_file(name + ".dat",std::ios_base::out | std::ios_base::app);
res_file.precision(17);
res_file << std::fixed << accuracy << " "<< time_running << std::endl;
}
int main(int nargs, char *args[])
{
srand(time(NULL));
bopt_params par = initialize_parameters_to_default();
par.n_iterations = 190;
par.n_inner_iterations = 250;
par.n_iter_relearn = 0;
//par.random_seed = 0;
par.verbose_level = 0;
par.noise = 1e-10;
par.sigma_s = 1;
par.sc_type = SC_ML;
par.init_method = 3;
strcpy(par.crit_name,"cLCB");
par.crit_params[0] = 0.125;
par.n_crit_params = 1;
par.force_jump = 0;
strcpy(par.kernel.name,"kMaternISO5");
benchmark<BraninNormalized>(par, "branin");
benchmark<Hartman6>(par, "hartman6");
benchmark<Hartman3>(par, "hartman3");
benchmark<Rastrigin>(par, "rastrigin");
benchmark<Sphere>(par, "sphere");
benchmark<Ellipsoid>(par, "ellipsoid");
benchmark<GoldenPrice>(par, "goldenprice");
benchmark<SixHumpCamel>(par, "sixhumpcamel");
return 0;
}
/*
-------------------------------------------------------------------------
This file is part of BayesOpt, an efficient C++ library for
Bayesian optimization.
Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
BayesOpt is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
BayesOpt is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with BayesOpt. If not, see <http://www.gnu.org/licenses/>.
------------------------------------------------------------------------
*/
#define _USE_MATH_DEFINES
#include <cmath>
#include <algorithm>
#include <boost/numeric/ublas/assignment.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <bayesopt/bayesopt.hpp>
using namespace bayesopt;
// support functions
inline double sign(double x)
{
if (x < 0)
return -1;
if (x > 0)
return 1;
return 0;
}
inline double sqr(double x)
{
return x * x;
};
inline double hat(double x)
{
if (x != 0)
return log(fabs(x));
return 0;
}
inline double c1(double x)
{
if (x > 0)
return 10;
return 5.5;
}
inline double c2(double x)
{
if (x > 0)
return 7.9;
return 3.1;
}
inline vectord t_osz(const vectord& x)
{
vectord r = x;
for (int i = 0; i < x.size(); i++)
r(i) = sign(x(i)) * exp(hat(x(i)) + 0.049 * sin(c1(x(i)) * hat(x(i))) + sin(c2(x(i)) * hat(x(i))));
return r;
}
struct Sphere {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
vectord opt(2);
opt <<= 0.5, 0.5;
return sqr(norm_2(x - opt));
}
matrixd solutions() const
{
matrixd sols(1, 2);
sols <<= 0.5, 0.5;
return sols;
}
};
struct Ellipsoid {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
vectord opt(2);
opt <<= 0.5, 0.5;
vectord z = t_osz(x - opt);
double r = 0;
for (size_t i = 0; i < dim_in; ++i)
r += std::pow(10, ((double)i) / (dim_in - 1.0)) * z(i) * z(i) + 1;
return r;
}
matrixd solutions() const
{
matrixd sols(1, 2);
sols <<= 0.5, 0.5;
return sols;
}
};
struct Rastrigin {
static constexpr size_t dim_in = 4;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
double f = 10 * dim_in;
for (size_t i = 0; i < dim_in; ++i)
f += x(i) * x(i) - 10 * cos(2 * M_PI * x(i));
return f;
}
matrixd solutions() const
{
matrixd sols(1, dim_in);
for (size_t i = 0; i < dim_in; ++i)
sols(0, i) = 0;
return sols;
}
};
// see : http://www.sfu.ca/~ssurjano/hart3.html
struct Hartman3 {
static constexpr size_t dim_in = 3;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
matrixd a(4, 3);
matrixd p(4, 3);
a <<= 3.0, 10, 30, 0.1, 10, 35, 3.0, 10, 30, 0.1, 10, 36;
p <<= 0.3689, 0.1170, 0.2673, 0.4699, 0.4387, 0.7470, 0.1091, 0.8732, 0.5547,
0.0382, 0.5743, 0.8828;
vectord alpha(4);
alpha <<= 1.0, 1.2, 3.0, 3.2;
double res = 0;
for (int i = 0; i < 4; i++) {
double s = 0.0f;
for (size_t j = 0; j < 3; j++) {
s += a(i, j) * (x(j) - p(i, j)) * (x(j) - p(i, j));
}
res += alpha(i) * exp(-s);
}
return -res;
}
matrixd solutions() const
{
matrixd sols(1, 3);
sols <<= 0.114614, 0.555649, 0.852547;
return sols;
}
};
// see : http://www.sfu.ca/~ssurjano/hart6.html
struct Hartman6 {
static constexpr size_t dim_in = 6;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
matrixd a(4, 6);
matrixd p(4, 6);
a <<= 10, 3, 17, 3.5, 1.7, 8, 0.05, 10, 17, 0.1, 8, 14, 3, 3.5, 1.7, 10, 17,
8, 17, 8, 0.05, 10, 0.1, 14;
p <<= 0.1312, 0.1696, 0.5569, 0.0124, 0.8283, 0.5886, 0.2329, 0.4135, 0.8307,
0.3736, 0.1004, 0.9991, 0.2348, 0.1451, 0.3522, 0.2883, 0.3047, 0.665,
0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381;
vectord alpha(4);
alpha <<= 1.0, 1.2, 3.0, 3.2;
double res = 0;
for (int i = 0; i < 4; i++) {
double s = 0.0f;
for (size_t j = 0; j < 6; j++) {
s += a(i, j) * sqr(x(j) - p(i, j));
}
res += alpha(i) * exp(-s);
}
return -res;
}
matrixd solutions() const
{
matrixd sols(1, 6);
sols <<= 0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573;
return sols;
}
};
// see : http://www.sfu.ca/~ssurjano/goldpr.html
// (with ln, as suggested in Jones et al.)
struct GoldenPrice {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const vectord& xx) const
{
vectord x = (4.0 * xx);
x(0) -= 2.0;
x(1) -= 2.0;
double r = (1 + (x(0) + x(1) + 1) * (x(0) + x(1) + 1) * (19 - 14 * x(0) + 3 * x(0) * x(0) - 14 * x(1) + 6 * x(0) * x(1) + 3 * x(1) * x(1))) * (30 + (2 * x(0) - 3 * x(1)) * (2 * x(0) - 3 * x(1)) * (18 - 32 * x(0) + 12 * x(0) * x(0) + 48 * x(1) - 36 * x(0) * x(1) + 27 * x(1) * x(1)));
return log(r) - 5;
}
matrixd solutions() const
{
matrixd sols(1, 2);
sols <<= 0.5, 0.25;
return sols;
}
};
struct BraninNormalized {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
double a = x(0) * 15 - 5;
double b = x(1) * 15;
return sqr(b - (5.1 / (4 * sqr(M_PI))) * sqr(a) + 5 * a / M_PI - 6) + 10 * (1 - 1 / (8 * M_PI)) * cos(a) + 10;
}
matrixd solutions() const
{
matrixd sols(3, 2);
sols <<= 0.1238938, 0.818333,
0.5427728, 0.151667,
0.961652, 0.1650;
return sols;
}
};
struct SixHumpCamel {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const vectord& x) const
{
double x1_2 = x(0)*x(0);
double x2_2 = x(1)*x(1);
double tmp1 = (4 - 2.1 * x1_2 + (x1_2*x1_2)/3) * x1_2;
double tmp2 = x(0)*x(1);
double tmp3 = (-4 + 4 * x2_2) * x2_2;
return tmp1 + tmp2 + tmp3;
}
matrixd solutions() const
{
matrixd sols(2, 2);
sols <<= 0.0898, -0.7126,
-0.0898, 0.7126;
return sols;
}
};
template <typename Function>
class Benchmark : public bayesopt::ContinuousModel
{
public:
Benchmark(bopt_params par):
ContinuousModel(Function::dim_in, par) {}
double evaluateSample(const vectord& xin)
{
return f(xin);
}
bool checkReachability(const vectord &query)
{return true;};
double accuracy(double x)
{
matrixd sols = f.solutions();
double diff = std::abs(x - f(row(sols, 0)));
double min_diff = diff;
for (size_t i = 1; i < sols.size1(); i++) {
diff = std::abs(x - f(row(sols, i)));
if (diff < min_diff)
min_diff = diff;
}
return min_diff;
}
Function f;
};
#include <iostream>
#include <chrono>
#include <limbo/limbo.hpp>
#include "testfunctions.hpp"
using namespace limbo;
struct Params {
struct bayes_opt_bobase {
BO_PARAM(bool, stats_enabled, true);
};
struct bayes_opt_boptimizer {
BO_PARAM(double, noise, 1e-10);
BO_PARAM(int, dump_period, -1);
};
struct stop_maxiterations {
BO_PARAM(int, iterations, 190);
};
struct kernel_maternfivehalfs {
BO_PARAM(double, sigma, 1);
BO_PARAM(double, l, 1);
};
struct acqui_ucb {
BO_PARAM(double, alpha, 0.125);
};
struct init_randomsampling {
BO_PARAM(int, samples, 10);
};
struct opt_nloptnograd {
BO_DYN_PARAM(int, iterations);
BO_PARAM(double, epsilon, 1e-8);
};
struct mean_constant {
BO_PARAM(double, constant, 1);
};
};
BO_DECLARE_DYN_PARAM(int, Params::opt_nloptnograd, iterations);
template <typename Params>
struct StaticInit {
template <typename StateFunction, typename AggregatorFunction, typename Opt>
void operator()(const StateFunction& seval, const AggregatorFunction&, Opt& opt) const
{
Eigen::Vector2d x;
x(0) = -0.055497063038156969690135505676;
x(1) = -0.351795677434771201331992637257;
opt.add_new_sample(x, seval(x));
x(0) = 0.240911582092300627645879219810;
x(1) = -0.351795677434771201331992637257;
opt.add_new_sample(x, seval(x));
x(0) = -0.351795677434771201331992637257;
x(1) = -0.343720034134076580133895794673;
opt.add_new_sample(x, seval(x));
x(0) = 0.055387032304313603995977911771;
x(1) = -0.076908787408790806346486820340;
opt.add_new_sample(x, seval(x));
x(0) = 0.351795677434771201331992637257;
x(1) = -0.076908787408790806346486820339;
opt.add_new_sample(x, seval(x));
x(0) = -0.351795677434771201331992637257;
x(1) = -0.047311389003618982797881069187;
opt.add_new_sample(x, seval(x));
x(0) = -0.132623434010086730621834242216;
x(1) = 0.152242144215576109267055784036;
opt.add_new_sample(x, seval(x));
x(0) = 0.351795677434771201331992637257;
x(1) = 0.219499857721666790989527905147;
opt.add_new_sample(x, seval(x));
x(0) = -0.351795677434771201331992637257;
x(1) = 0.351795677434771201331992637257;
opt.add_new_sample(x, seval(x));
x(0) = 0.086548809414597740088324152827;
x(1) = 0.351795677434771201331992637257;
opt.add_new_sample(x, seval(x));
}
};
template <typename Optimizer, typename Function>
void benchmark(const std::string& name)
{
int iters_base = 250;
Params::opt_nloptnograd::set_iterations(iters_base * Function::dim_in);
auto t1 = std::chrono::steady_clock::now();
Optimizer opt;
Benchmark<Function> target;
opt.optimize(target);
auto time_running = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - t1).count();
std::cout.precision(17);
std::cout << std::endl;
auto best = opt.best_observation();
double accuracy = target.accuracy(best);
std::cout << name << std::endl;
std::cout << "Result: " << std::fixed << opt.best_sample().transpose() << " -> " << best << std::endl;
std::cout << "Smallest difference: " << accuracy << std::endl;
std::cout << "Time running: " << time_running << "ms" << std::endl
<< std::endl;
std::ofstream res_file(name + ".dat",std::ios_base::out | std::ios_base::app);
res_file.precision(17);
res_file << std::fixed << accuracy << " " << time_running << std::endl;
}
int main()
{
srand(time(NULL));
typedef kernel::MaternFiveHalfs<Params> Kernel_t;
typedef opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L> AcquiOpt_t;
typedef boost::fusion::vector<stop::MaxIterations<Params>> Stop_t;
// typedef mean_functions::MeanFunctionARD<Params, mean_functions::MeanData<Params>> Mean_t;
typedef mean::Constant<Params> Mean_t;
typedef boost::fusion::vector<> Stat_t;
typedef init::RandomSampling<Params> Init_t;
typedef model::GP<Params, Kernel_t, Mean_t> GP_t;
typedef acqui::UCB<Params, GP_t> Acqui_t;
typedef bayes_opt::BOptimizer<Params, modelfun<GP_t>, initfun<Init_t>, acquifun<Acqui_t>, acquiopt<AcquiOpt_t>, statsfun<Stat_t>, stopcrit<Stop_t>> Opt_t;
benchmark<Opt_t, BraninNormalized>("branin");
benchmark<Opt_t, Hartman6>("hartman6");
benchmark<Opt_t, Hartman3>("hartman3");
benchmark<Opt_t, Rastrigin>("rastrigin");
benchmark<Opt_t, Sphere>("sphere");
benchmark<Opt_t, Ellipsoid>("ellipsoid");
benchmark<Opt_t, GoldenPrice>("goldenprice");
benchmark<Opt_t, SixHumpCamel>("sixhumpcamel");
return 0;
}
#define _USE_MATH_DEFINES
#include <Eigen/Core>
// support functions
inline double sign(double x)
{
if (x < 0)
return -1;
if (x > 0)
return 1;
return 0;
}
inline double sqr(double x)
{
return x * x;
};
inline double hat(double x)
{
if (x != 0)
return log(fabs(x));
return 0;
}
inline double c1(double x)
{
if (x > 0)
return 10;
return 5.5;
}
inline double c2(double x)
{
if (x > 0)
return 7.9;
return 3.1;
}
inline Eigen::VectorXd t_osz(const Eigen::VectorXd& x)
{
Eigen::VectorXd r = x;
for (int i = 0; i < x.size(); i++)
r(i) = sign(x(i)) * exp(hat(x(i)) + 0.049 * sin(c1(x(i)) * hat(x(i))) + sin(c2(x(i)) * hat(x(i))));
return r;
}
struct Sphere {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd opt(2);
opt << 0.5, 0.5;
return (x - opt).squaredNorm();
}
Eigen::MatrixXd solutions() const
{
Eigen::MatrixXd sols(1, 2);
sols << 0.5, 0.5;
return sols;
}
};
struct Ellipsoid {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 1;
double operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd opt(2);
opt << 0.5, 0.5;
Eigen::VectorXd z = t_osz(x - opt);
double r = 0;
for (size_t i = 0; i < dim_in; ++i)
r += std::pow(10, ((double)i) / (dim_in - 1.0)) * z(i) * z(i) + 1;
return r;
}
Eigen::MatrixXd solutions() const
{