Merge remote-tracking branch 'origin' into ehvi_multi_experimental

Conflicts:
	src/limbo/experimental/bayes_opt/bo_multi.hpp

.gitignore

@@ -17,6 +17,7 @@ waf_xcode.sh
 exp
 src/tests/combinations
 params_*
+benchmark_results
 
 # Ignored folders for the documentation
 _build

docs/update_limbo_docs.sh

+set -x
+set -e
+
+
+export SPHINX_RESIBOTS_THEME="$HOME/git/sphinx_resibots_theme"
+rm -rf /tmp/doc_limbo
+
+mkdir /tmp/doc_limbo/
+cd /tmp/doc_limbo
+git clone git@github.com:resibots/limbo.git
+cd limbo/docs
+make html
+mv doxygen_doc /tmp
+git checkout gh-pages
+rm -rf doxygen_doc
+cp -r /tmp/doxygen_doc .
+cd ..
+cp -r docs/_build/html/* .
+git add .
+git commit -m 'automatic update of the doc [ci skip]'
+git push origin gh-pages
+
+rm -rf /tmp/doxygen_doc
+rm -rf /tmp/doc_limbo
+
+

src/benchmarks/bayesopt/hp_opt.cpp

+#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 = 50;
+    //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, "kSEARD");
+
+    benchmark<BraninNormalized>(par, "branin");
+    benchmark<Hartmann6>(par, "hartmann6");
+    benchmark<Hartmann3>(par, "hartmann3");
+    benchmark<Rastrigin>(par, "rastrigin");
+    benchmark<Sphere>(par, "sphere");
+    benchmark<Ellipsoid>(par, "ellipsoid");
+    benchmark<GoldsteinPrice>(par, "goldsteinprice");
+    benchmark<SixHumpCamel>(par, "sixhumpcamel");
+
+    return 0;
+}

src/benchmarks/bayesopt/simple.cpp

+#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<Hartmann6>(par, "hartmann6");
+    benchmark<Hartmann3>(par, "hartmann3");
+    benchmark<Rastrigin>(par, "rastrigin");
+    benchmark<Sphere>(par, "sphere");
+    benchmark<Ellipsoid>(par, "ellipsoid");
+    benchmark<GoldsteinPrice>(par, "goldsteinprice");
+    benchmark<SixHumpCamel>(par, "sixhumpcamel");
+
+    return 0;
+}

src/benchmarks/bayesopt/testfunctions.hpp

+/*
+-------------------------------------------------------------------------
+   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 Hartmann3 {
+    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 Hartmann6 {
+    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 GoldsteinPrice {
+    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 = -3 + 6 * x(0);
+        double x2 = -2 + 4 * x(1);
+        double x1_2 = x1 * x1;
+        double x2_2 = x2 * x2;
+
+        double tmp1 = (4 - 2.1 * x1_2 + (x1_2 * x1_2) / 3) * x1_2;
+        double tmp2 = x1 * x2;
+        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;
+        sols(0, 0) = (sols(0, 0) + 3.0) / 6.0;
+        sols(1, 0) = (sols(1, 0) + 3.0) / 6.0;
+        sols(0, 1) = (sols(0, 1) + 2.0) / 4.0;
+        sols(1, 1) = (sols(1, 1) + 2.0) / 4.0;
+        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;
+};

src/benchmarks/limbo/hp_opt.cpp

+#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, false);
+    };
+    struct bayes_opt_boptimizer {
+        BO_PARAM(double, noise, 1e-10);
+        BO_PARAM(int, hp_period, 50);
+    };
+    struct stop_maxiterations {
+        BO_PARAM(int, iterations, 190);
+    };
+    struct kernel_maternfivehalfs {
+        BO_PARAM(double, sigma_sq, 1);
+        BO_PARAM(double, l, 1);
+    };
+    struct acqui_ucb {
+        BO_PARAM(double, alpha, 0.125);
+    };
+    struct init_randomsampling {
+        BO_PARAM(int, samples, 10);
+    };
+    struct kernel_squared_exp_ard : public defaults::kernel_squared_exp_ard {
+    };
+    struct mean_constant {
+        BO_PARAM(double, constant, 1);
+    };
+    struct opt_rprop : public defaults::opt_rprop {
+    };
+};
+
+struct DirectParams {
+    struct opt_nloptnograd {
+        BO_DYN_PARAM(int, iterations);
+    };
+};
+
+struct BobyqaParams {
+    struct opt_nloptnograd {
+        BO_DYN_PARAM(int, iterations);
+    };
+};
+struct BobyqaParams_HP {
+    struct opt_nloptnograd {
+        BO_DYN_PARAM(int, iterations);
+    };
+};
+
+BO_DECLARE_DYN_PARAM(int, DirectParams::opt_nloptnograd, iterations);
+BO_DECLARE_DYN_PARAM(int, BobyqaParams::opt_nloptnograd, iterations);
+BO_DECLARE_DYN_PARAM(int, BobyqaParams_HP::opt_nloptnograd, iterations);
+
+template <typename Optimizer, typename Function>
+void benchmark(const std::string& name)
+{
+    int iters_base = 250;
+    DirectParams::opt_nloptnograd::set_iterations(static_cast<int>(iters_base * Function::dim_in * 0.9));
+    BobyqaParams::opt_nloptnograd::set_iterations(iters_base * Function::dim_in - DirectParams::opt_nloptnograd::iterations());
+
+    BobyqaParams_HP::opt_nloptnograd::set_iterations(10 * Function::dim_in * 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::SquaredExpARD<Params> Kernel_t;
+    typedef opt::Chained<Params, opt::NLOptNoGrad<DirectParams, nlopt::GN_DIRECT_L>, opt::NLOptNoGrad<BobyqaParams, nlopt::LN_BOBYQA>> 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, model::gp::KernelLFOpt<Params, opt::NLOptNoGrad<BobyqaParams_HP, nlopt::LN_BOBYQA>>> 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, Hartmann6>("hartmann6");
+    benchmark<Opt_t, Hartmann3>("hartmann3");
+    benchmark<Opt_t, Rastrigin>("rastrigin");
+    benchmark<Opt_t, Sphere>("sphere");
+    benchmark<Opt_t, Ellipsoid>("ellipsoid");
+    benchmark<Opt_t, GoldsteinPrice>("goldsteinprice");
+    benchmark<Opt_t, SixHumpCamel>("sixhumpcamel");