Commit 0a15d627 authored by Vaios Papaspyros's avatar Vaios Papaspyros
Browse files

- First attempt for constrained bopt using modified expected improvement and a modified boptimizer

- Adding a new example
parent ca836aeb
#include <limbo/limbo.hpp>
#include <limbo/experimental/bayes_opt/cboptimizer.hpp>
#include <limbo/experimental/acqui/cei.hpp>
using namespace limbo;
struct Params {
struct cbayes_opt_boptimizer : public defaults::cbayes_opt_boptimizer {
BO_PARAM(double, noise, 0.01);
};
struct init_randomsampling {
BO_PARAM(int, samples, 1);
};
struct kernel_exp : public defaults::kernel_exp {
};
struct bayes_opt_bobase {
BO_PARAM(bool, stats_enabled, true);
};
struct stop_maxiterations {
BO_PARAM(int, iterations, 100);
};
struct acqui_cei : public defaults::acqui_cei {
};
struct mean_constant {
BO_PARAM(double, constant, 0.0);
};
#ifdef USE_NLOPT
struct opt_nloptnograd : public defaults::opt_nloptnograd {
};
#elif defined(USE_LIBCMAES)
struct opt_cmaes : public defaults::opt_cmaes {
};
#else
struct opt_gridsearch : public defaults::opt_gridsearch {
};
#endif
};
#ifdef DIM6
#define ZDT_DIM 6
#elif defined(DIM2)
#define ZDT_DIM 2
#else
#define ZDT_DIM 30
#endif
struct zdt1 {
static constexpr size_t dim_in = ZDT_DIM;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd res(2);
double f1 = x(0);
double g = 1.0;
for (int i = 1; i < x.size(); ++i)
g += 9.0 / (x.size() - 1) * x(i);
double h = 1.0f - sqrtf(f1 / g);
double f2 = g * h;
res(0) = 1.0 - f1;
res(1) = 1.0 - f2;
return res;
}
};
struct zdt2 {
static constexpr size_t dim_in = ZDT_DIM;
static constexpr size_t dim_out = 2;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd res(2);
double f1 = x(0);
double g = 1.0;
for (int i = 1; i < x.size(); ++i)
g += 9.0 / (x.size() - 1) * x(i);
double h = 1.0f - pow((f1 / g), 2.0);
double f2 = g * h;
res(0) = 1.0 - f1;
res(1) = 1.0 - f2;
return res;
}
};
struct zdt3 {
static constexpr size_t dim_in = ZDT_DIM;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd res(2);
double f1 = x(0);
double g = 1.0;
for (int i = 1; i < x.size(); ++i)
g += 9.0 / (x.size() - 1) * x(i);
double h = 1.0f - sqrtf(f1 / g) - f1 / g * sin(10 * M_PI * f1);
double f2 = g * h;
res(0) = 1.0 - f1;
res(1) = 1.0 - f2;
return res;
}
};
struct mop2 {
static constexpr size_t dim_in = 2;
static constexpr size_t dim_out = 2;
static constexpr size_t nb_constraints = 2;
Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
{
Eigen::VectorXd res(4);
// scale to [-2, 2]
Eigen::VectorXd xx = (x * 4.0).array() - 2.0;
// f1, f2
Eigen::VectorXd v1 = (xx.array() - 1.0 / sqrt(xx.size())).array().square();
Eigen::VectorXd v2 = (xx.array() + 1.0 / sqrt(xx.size())).array().square();
double f1 = 1.0 - exp(-v1.sum());
double f2 = 1.0 - exp(-v2.sum());
// we _maximize in [0:1]
res(0) = 1 - f1;
res(1) = 1 - f2;
// testing the constraints
(res(0) > 0.4) ? res(2) = 1 : res(2) = 0;
(res(1) > 0.4) ? res(3) = 1 : res(3) = 0;
std::cout << res(0) << " " << res(1) << std::endl;
return res;
}
};
int main()
{
tools::par::init();
// typedef zdt1 func_t;
// typedef zdt2 func_t;
// typedef zdt3 func_t;
typedef mop2 func_t;
using Stop_t = boost::fusion::vector<stop::MaxIterations<Params>>;
using Stat_t = boost::fusion::vector<stat::Samples<Params>,
stat::BlSamples<Params>,
stat::BestObservations<Params>,
stat::AggregatedObservations<Params>>;
using Mean_t = mean::Constant<Params>;
using Kernel_t = kernel::Exp<Params>;
using GP_t = model::GP<Params, Kernel_t, Mean_t>;
using Acqui_t = experimental::acqui::CEI<Params, GP_t>;
using Init_t = init::RandomSampling<Params>;
experimental::bayes_opt::CBOptimizer<Params,
modelfun<GP_t>,
acquifun<Acqui_t>,
statsfun<Stat_t>,
initfun<Init_t>,
stopcrit<Stop_t>>
opt;
opt.optimize(func_t());
return 0;
}
...@@ -4,32 +4,32 @@ ...@@ -4,32 +4,32 @@
#| This project has received funding from the European Research Council (ERC) under #| This project has received funding from the European Research Council (ERC) under
#| the European Union's Horizon 2020 research and innovation programme (grant #| the European Union's Horizon 2020 research and innovation programme (grant
#| agreement No 637972) - see http://www.resibots.eu #| agreement No 637972) - see http://www.resibots.eu
#| #|
#| Contributor(s): #| Contributor(s):
#| - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr) #| - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
#| - Antoine Cully (antoinecully@gmail.com) #| - Antoine Cully (antoinecully@gmail.com)
#| - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr) #| - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
#| - Federico Allocati (fede.allocati@gmail.com) #| - Federico Allocati (fede.allocati@gmail.com)
#| - Vaios Papaspyros (b.papaspyros@gmail.com) #| - Vaios Papaspyros (b.papaspyros@gmail.com)
#| #|
#| This software is a computer library whose purpose is to optimize continuous, #| This software is a computer library whose purpose is to optimize continuous,
#| black-box functions. It mainly implements Gaussian processes and Bayesian #| black-box functions. It mainly implements Gaussian processes and Bayesian
#| optimization. #| optimization.
#| Main repository: http://github.com/resibots/limbo #| Main repository: http://github.com/resibots/limbo
#| Documentation: http://www.resibots.eu/limbo #| Documentation: http://www.resibots.eu/limbo
#| #|
#| This software is governed by the CeCILL-C license under French law and #| This software is governed by the CeCILL-C license under French law and
#| abiding by the rules of distribution of free software. You can use, #| abiding by the rules of distribution of free software. You can use,
#| modify and/ or redistribute the software under the terms of the CeCILL-C #| modify and/ or redistribute the software under the terms of the CeCILL-C
#| license as circulated by CEA, CNRS and INRIA at the following URL #| license as circulated by CEA, CNRS and INRIA at the following URL
#| "http://www.cecill.info". #| "http://www.cecill.info".
#| #|
#| As a counterpart to the access to the source code and rights to copy, #| As a counterpart to the access to the source code and rights to copy,
#| modify and redistribute granted by the license, users are provided only #| modify and redistribute granted by the license, users are provided only
#| with a limited warranty and the software's author, the holder of the #| with a limited warranty and the software's author, the holder of the
#| economic rights, and the successive licensors have only limited #| economic rights, and the successive licensors have only limited
#| liability. #| liability.
#| #|
#| In this respect, the user's attention is drawn to the risks associated #| In this respect, the user's attention is drawn to the risks associated
#| with loading, using, modifying and/or developing or reproducing the #| with loading, using, modifying and/or developing or reproducing the
#| software by the user in light of its specific status of free software, #| software by the user in light of its specific status of free software,
...@@ -40,10 +40,10 @@ ...@@ -40,10 +40,10 @@
#| requirements in conditions enabling the security of their systems and/or #| requirements in conditions enabling the security of their systems and/or
#| data to be ensured and, more generally, to use and operate it in the #| data to be ensured and, more generally, to use and operate it in the
#| same conditions as regards security. #| same conditions as regards security.
#| #|
#| The fact that you are presently reading this means that you have had #| The fact that you are presently reading this means that you have had
#| knowledge of the CeCILL-C license and that you accept its terms. #| knowledge of the CeCILL-C license and that you accept its terms.
#| #|
def build(bld): def build(bld):
...@@ -59,10 +59,15 @@ def build(bld): ...@@ -59,10 +59,15 @@ def build(bld):
target = 'parego', target = 'parego',
uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT', uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
use = 'limbo') use = 'limbo')
obj = bld.program(features = 'cxx', obj = bld.program(features = 'cxx',
source = 'multi.cpp', source = 'multi.cpp',
includes = '.. ../.. ../../../', includes = '.. ../.. ../../../',
target = 'multi', target = 'multi',
uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT', uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
use = 'limbo') use = 'limbo')
obj = bld.program(features = 'cxx',
source = 'constrained_ei.cpp',
includes = '.. ../.. ../../../',
target = 'constrained_ei',
uselib = 'BOOST EIGEN TBB SFERES LIBCMAES NLOPT',
use = 'limbo')
//| Copyright Inria May 2015
//| This project has received funding from the European Research Council (ERC) under
//| the European Union's Horizon 2020 research and innovation programme (grant
//| agreement No 637972) - see http://www.resibots.eu
//|
//| Contributor(s):
//| - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
//| - Antoine Cully (antoinecully@gmail.com)
//| - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
//| - Federico Allocati (fede.allocati@gmail.com)
//| - Vaios Papaspyros (b.papaspyros@gmail.com)
//|
//| This software is a computer library whose purpose is to optimize continuous,
//| black-box functions. It mainly implements Gaussian processes and Bayesian
//| optimization.
//| Main repository: http://github.com/resibots/limbo
//| Documentation: http://www.resibots.eu/limbo
//|
//| This software is governed by the CeCILL-C license under French law and
//| abiding by the rules of distribution of free software. You can use,
//| modify and/ or redistribute the software under the terms of the CeCILL-C
//| license as circulated by CEA, CNRS and INRIA at the following URL
//| "http://www.cecill.info".
//|
//| As a counterpart to the access to the source code and rights to copy,
//| modify and redistribute granted by the license, users are provided only
//| with a limited warranty and the software's author, the holder of the
//| economic rights, and the successive licensors have only limited
//| liability.
//|
//| In this respect, the user's attention is drawn to the risks associated
//| with loading, using, modifying and/or developing or reproducing the
//| software by the user in light of its specific status of free software,
//| that may mean that it is complicated to manipulate, and that also
//| therefore means that it is reserved for developers and experienced
//| professionals having in-depth computer knowledge. Users are therefore
//| encouraged to load and test the software's suitability as regards their
//| requirements in conditions enabling the security of their systems and/or
//| data to be ensured and, more generally, to use and operate it in the
//| same conditions as regards security.
//|
//| The fact that you are presently reading this means that you have had
//| knowledge of the CeCILL-C license and that you accept its terms.
//|
#ifndef LIMBO_ACQUI_CEI_HPP
#define LIMBO_ACQUI_CEI_HPP
#include <cmath>
#include <vector>
#include <Eigen/Core>
#include <limbo/tools/macros.hpp>
namespace limbo {
namespace defaults {
struct acqui_cei {
/// @ingroup acqui_defaults
BO_PARAM(double, jitter, 0.0);
};
}
namespace experimental {
namespace acqui {
/** @ingroup acqui
\rst
Constrained EI (Expected Improvement). See :cite:`TODO`, p. TODO
.. math::
TODO
Parameters:
- ``double jitter`` - :math:`\xi`
\endrst
*/
template <typename Params, typename Model>
class CEI {
public:
CEI(const std::vector<Model>& models, int iteration = 0) : _models(models) {}
size_t dim_in() const { return _models[0].dim_in(); }
size_t dim_out() const { return _models[0].dim_out(); }
template <typename AggregatorFunction>
double operator()(const Eigen::VectorXd& v, const AggregatorFunction& afun) const
{
Eigen::VectorXd mu;
double sigma_sq;
std::tie(mu, sigma_sq) = _models[0].query(v);
double sigma = std::sqrt(sigma_sq);
// If \sigma(x) = 0 or we do not have any observation yet we return 0
if (sigma < 1e-10 || _models[0].samples().size() < 1)
return 0.0;
// Compute constrained EI(x)
// First find the best so far (predicted) observation
// (We are zeroing infeasible samples subject to the constraint value)
std::vector<double> rewards;
for (auto s : _models[0].samples())
rewards.push_back(afun(_models[0].mu(s)));
double f_max = *std::max_element(rewards.begin(), rewards.end());
// Calculate Z and \Phi(Z) and \phi(Z)
double X = afun(mu) - f_max - Params::acqui_cei::jitter();
double Z = X / sigma;
double phi = std::exp(-0.5 * std::pow(Z, 2.0)) / std::sqrt(2.0 * M_PI);
double Phi = 0.5 * std::erfc(-Z / std::sqrt(2)); //0.5 * (1.0 + std::erf(Z / std::sqrt(2)));
return _pf(v, afun) * (X * Phi + sigma * phi);
}
protected:
const std::vector<Model>& _models;
template <typename AggregatorFunction>
double _pf(const Eigen::VectorXd& v, const AggregatorFunction& afun) const
{
double p = 1.0;
for (size_t i = 1; i < _models.size(); ++i) {
Eigen::VectorXd mu;
double sigma_sq;
std::tie(mu, sigma_sq) = _models[i].query(v);
double sigma = std::sqrt(sigma_sq);
double Z = afun(mu) / sigma;
double phi = (1.0 - std::exp(-0.5 * std::pow(Z, 2.0)) / std::sqrt(2.0 * M_PI));
p *= phi;
}
return p;
}
};
}
}
}
#endif
//| Copyright Inria May 2015
//| This project has received funding from the European Research Council (ERC) under
//| the European Union's Horizon 2020 research and innovation programme (grant
//| agreement No 637972) - see http://www.resibots.eu
//|
//| Contributor(s):
//| - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
//| - Antoine Cully (antoinecully@gmail.com)
//| - Kontantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
//| - Federico Allocati (fede.allocati@gmail.com)
//| - Vaios Papaspyros (b.papaspyros@gmail.com)
//|
//| This software is a computer library whose purpose is to optimize continuous,
//| black-box functions. It mainly implements Gaussian processes and Bayesian
//| optimization.
//| Main repository: http://github.com/resibots/limbo
//| Documentation: http://www.resibots.eu/limbo
//|
//| This software is governed by the CeCILL-C license under French law and
//| abiding by the rules of distribution of free software. You can use,
//| modify and/ or redistribute the software under the terms of the CeCILL-C
//| license as circulated by CEA, CNRS and INRIA at the following URL
//| "http://www.cecill.info".
//|
//| As a counterpart to the access to the source code and rights to copy,
//| modify and redistribute granted by the license, users are provided only
//| with a limited warranty and the software's author, the holder of the
//| economic rights, and the successive licensors have only limited
//| liability.
//|
//| In this respect, the user's attention is drawn to the risks associated
//| with loading, using, modifying and/or developing or reproducing the
//| software by the user in light of its specific status of free software,
//| that may mean that it is complicated to manipulate, and that also
//| therefore means that it is reserved for developers and experienced
//| professionals having in-depth computer knowledge. Users are therefore
//| encouraged to load and test the software's suitability as regards their
//| requirements in conditions enabling the security of their systems and/or
//| data to be ensured and, more generally, to use and operate it in the
//| same conditions as regards security.
//|
//| The fact that you are presently reading this means that you have had
//| knowledge of the CeCILL-C license and that you accept its terms.
//|
#ifndef LIMBO_BAYES_OPT_CBOPTIMIZER_HPP
#define LIMBO_BAYES_OPT_CBOPTIMIZER_HPP
#include <iostream>
#include <algorithm>
#include <iterator>
#include <boost/parameter/aux_/void.hpp>
#include <Eigen/Core>
#include <limbo/tools/macros.hpp>
#include <limbo/tools/random_generator.hpp>
#include <limbo/bayes_opt/bo_base.hpp>
#ifdef USE_NLOPT
#include <limbo/opt/nlopt_no_grad.hpp>
#elif defined USE_LIBCMAES
#include <limbo/opt/cmaes.hpp>
#else
#include <limbo/opt/grid_search.hpp>
#endif
namespace limbo {
namespace defaults {
struct cbayes_opt_boptimizer {
BO_PARAM(double, noise, 1e-6);
BO_PARAM(int, hp_period, -1);
};
}
namespace experimental {
BOOST_PARAMETER_TEMPLATE_KEYWORD(c_acquiopt)
namespace bayes_opt {
typedef boost::parameter::parameters<boost::parameter::optional<tag::c_acquiopt>,
boost::parameter::optional<limbo::tag::statsfun>,
boost::parameter::optional<limbo::tag::initfun>,
boost::parameter::optional<limbo::tag::acquifun>,
boost::parameter::optional<limbo::tag::stopcrit>,
boost::parameter::optional<limbo::tag::modelfun>>
cboptimizer_signature;
// clang-format off
/**
The classic Bayesian optimization algorithm.
\\rst
References: :cite:`brochu2010tutorial,Mockus2013`
\\endrst
This class takes the same template parameters as BoBase. It adds:
\\rst
+---------------------+------------+----------+---------------+
|type |typedef | argument | default |
+=====================+============+==========+===============+
|acqui. optimizer |acqui_opt_t | c_acquiopt | see below |
+---------------------+------------+----------+---------------+
\\endrst
The default value of acqui_opt_t is:
- ``opt::Cmaes<Params>`` if libcmaes was found in `waf configure`
- ``opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L_RAND>`` if NLOpt was found but libcmaes was not found
- ``opt::GridSearch<Params>`` otherwise (please do not use this: the algorithm will not work at all!)
*/
template <class Params,
class A1 = boost::parameter::void_,
class A2 = boost::parameter::void_,
class A3 = boost::parameter::void_,
class A4 = boost::parameter::void_,
class A5 = boost::parameter::void_,
class A6 = boost::parameter::void_>
// clang-format on
class CBOptimizer : public limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6> {
public:
// defaults
struct defaults {
#ifdef USE_NLOPT
typedef opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L_RAND> acquiopt_t;
#elif defined(USE_LIBCMAES)
typedef opt::Cmaes<Params> acquiopt_t;
#else
#warning NO NLOpt, and NO Libcmaes: the acquisition function will be optimized by a grid search algorithm (which is usually bad). Please install at least NLOpt or libcmaes to use limbo!.
typedef opt::GridSearch<Params> acquiopt_t;
#endif
};
/// link to the corresponding BoBase (useful for typedefs)
typedef limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6> base_t;
typedef typename base_t::model_t model_t;
typedef typename base_t::acquisition_function_t acquisition_function_t;
// extract the types
typedef typename cboptimizer_signature::bind<A1, A2, A3, A4, A5, A6>::type args;
typedef typename boost::parameter::binding<args, tag::c_acquiopt, typename defaults::acquiopt_t>::type acqui_optimizer_t;
/// The main function (run the Bayesian optimization algorithm)
template <typename StateFunction, typename AggregatorFunction = FirstElem>
void optimize(const StateFunction& sfun, const AggregatorFunction& afun = AggregatorFunction(), bool reset = true)
{
_models.resize(StateFunction::nb_constraints + 1); // reserving enough space for the objective function gp + constraints
_nb_constraints = StateFunction::nb_constraints;
_dim_out = StateFunction::dim_out;
this->_init(sfun, afun, reset);
if (!this->_observations.empty()) {
_split_observations();
for (size_t i = 0; i < _models.size(); ++i)
_models[i].compute(this->_samples, _obs[i], Eigen::VectorXd::Constant(_obs[i].size(), Params::cbayes_opt_boptimizer::noise()), this->_bl_samples, Eigen::VectorXd::Constant(this->_bl_samples.size(), Params::cbayes_opt_boptimizer::noise()));
}
else {
_models[0] = model_t(StateFunction::dim_in, StateFunction::dim_out);
for (size_t i = 1; i < _nb_constraints; ++i)
_models[i] = model_t(StateFunction::dim_in, 1); // explicitly to 1 for the constraints
}
acqui_optimizer_t acqui_optimizer;
while (!this->_stop(*this, afun)) {
acquisition_function_t acqui(_models, this->_current_iteration);
// we do not have gradient in our current acquisition function
auto acqui_optimization =
[&](const Eigen::VectorXd& x, bool g) { return opt::no_grad(acqui(x, afun)); };
Eigen::VectorXd starting_point = tools::random_vector(StateFunction::dim_in);
Eigen::VectorXd new_sample = acqui_optimizer(acqui_optimization, starting_point, true);
bool blacklisted = !this->eval_and_add(sfun, new_sample);
this->_update_stats(*this, afun, blacklisted);
if (blacklisted) {
for (size_t i = 0; i < _models.size(); ++i)
_models[i].add_bl_sample(this->_bl_samples.back(), Params::cbayes_opt_boptimizer::noise());
}
else {
for (size_t i = 0; i < _models.size(); ++i)
_models[i].add_sample(this->_samples.back(), _obs[i].back(), Params::cbayes_opt_boptimizer::noise());
}
if (Params::cbayes_opt_boptimizer::hp_period() > 0
&& (this->_current_iteration + 1) % Params::cbayes_opt_boptimizer::hp_period() == 0) {
for (size_t i = 0; i < _models.size(); ++i)
_models[i].optimize_hyperparams();
}
this->_current_iteration++;
this->_total_iterations++;
}
}
/// return the best observation so far (i.e. max(f(x)))
template <typename AggregatorFunction = FirstElem>
const Eigen::VectorXd& best_observation(const AggregatorFunction& afun = AggregatorFunction()) const