cboptimizer.hpp 12.7 KB
Newer Older
1
2
3
4
5
6
7
8
//| 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)
9
//|   - Konstantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
10
11
//|   - Federico Allocati (fede.allocati@gmail.com)
//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
12
//|   - Roberto Rama (bertoski@gmail.com)
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
//|
//| 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 <algorithm>
50
#include <iostream>
51
52
53
54
55
56
#include <iterator>

#include <boost/parameter/aux_/void.hpp>

#include <Eigen/Core>

57
#include <limbo/bayes_opt/bo_base.hpp>
58
59
60
61
62
63
64
65
66
67
68
69
#include <limbo/tools/macros.hpp>
#include <limbo/tools/random_generator.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 {
70
        struct bayes_opt_cboptimizer {
71
            BO_PARAM(int, hp_period, -1);
72
            BO_PARAM(bool, bounded, true);
73
74
75
76
        };
    }

    namespace experimental {
77
        BOOST_PARAMETER_TEMPLATE_KEYWORD(constraint_modelfun)
78
79
80

        namespace bayes_opt {

81
            using cboptimizer_signature = boost::parameter::parameters<boost::parameter::optional<limbo::tag::acquiopt>,
82
83
84
85
                boost::parameter::optional<limbo::tag::statsfun>,
                boost::parameter::optional<limbo::tag::initfun>,
                boost::parameter::optional<limbo::tag::acquifun>,
                boost::parameter::optional<limbo::tag::stopcrit>,
86
                boost::parameter::optional<limbo::tag::modelfun>,
87
                boost::parameter::optional<limbo::experimental::tag::constraint_modelfun>>;
88
89
90
91
92
93
94
95
96
97
98
99
100
101

            // 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       |
        +=====================+============+==========+===============+
Vaios Papaspyros's avatar
Vaios Papaspyros committed
102
        |acqui. optimizer     |acqui_opt_t | acquiopt | see below     |
103
104
105
106
107
108
109
110
111
112
113
114
115
116
        +---------------------+------------+----------+---------------+
        \\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_,
117
118
          class A6 = boost::parameter::void_,
          class A7 = boost::parameter::void_>
119
120
121
122
123
124
            // clang-format on
            class CBOptimizer : public limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6> {
            public:
                // defaults
                struct defaults {
#ifdef USE_NLOPT
125
                    using acquiopt_t = limbo::opt::NLOptNoGrad<Params, nlopt::GN_DIRECT_L_RAND>;
126
#elif defined(USE_LIBCMAES)
127
                    using acquiopt_t = limbo::opt::Cmaes<Params>;
128
129
#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!.
130
                    using acquiopt_t = limbo::opt::GridSearch<Params>;
131
#endif
132
133
134
                    using kf_t = limbo::kernel::Exp<Params>;
                    using mean_t = limbo::mean::Constant<Params>;
                    using constraint_model_t = limbo::model::GP<Params, kf_t, mean_t>;
135
136
                };
                /// link to the corresponding BoBase (useful for typedefs)
137
138
                using base_t = limbo::bayes_opt::BoBase<Params, A1, A2, A3, A4, A5, A6>;
                using model_t = typename base_t::model_t;
139

140
                using acquisition_function_t = typename base_t::acquisition_function_t;
141
                // extract the types
142
143
                using args = typename cboptimizer_signature::bind<A1, A2, A3, A4, A5, A6, A7>::type;
                using acqui_optimizer_t = typename boost::parameter::binding<args, limbo::tag::acquiopt, typename defaults::acquiopt_t>::type;
144

145
                using constraint_model_t = typename boost::parameter::binding<args, limbo::experimental::tag::constraint_modelfun, typename defaults::constraint_model_t>::type;
146

147
148
149
150
                /// 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)
                {
151
152
                    _nb_constraints = StateFunction::nb_constraints();
                    _dim_out = StateFunction::dim_out();
153
154
155
156
157

                    this->_init(sfun, afun, reset);

                    if (!this->_observations.empty()) {
                        _split_observations();
158
                        _model.compute(this->_samples, _obs[0]);
159
                        if (_nb_constraints > 0)
160
                            _constraint_model.compute(this->_samples, _obs[1]);
161
162
                    }
                    else {
163
                        _model = model_t(StateFunction::dim_in(), StateFunction::dim_out());
164
                        if (_nb_constraints > 0)
165
                            _constraint_model = constraint_model_t(StateFunction::dim_in(), _nb_constraints);
166
167
168
169
170
                    }

                    acqui_optimizer_t acqui_optimizer;

                    while (!this->_stop(*this, afun)) {
171
                        acquisition_function_t acqui(_model, _constraint_model, this->_current_iteration);
172
173

                        auto acqui_optimization =
174
                            [&](const Eigen::VectorXd& x, bool g) { return acqui(x, afun, g); };
175
                        Eigen::VectorXd starting_point = tools::random_vector(StateFunction::dim_in(), Params::bayes_opt_cboptimizer::bounded());
176
                        Eigen::VectorXd new_sample = acqui_optimizer(acqui_optimization, starting_point, Params::bayes_opt_cboptimizer::bounded());
177
178
179
180
                        this->eval_and_add(sfun, new_sample);

                        this->_update_stats(*this, afun);

181
                        _model.add_sample(this->_samples.back(), _obs[0].back());
182
                        if (_nb_constraints > 0)
183
                            _constraint_model.add_sample(this->_samples.back(), _obs[1].back());
184

185
186
                        if (Params::bayes_opt_cboptimizer::hp_period() > 0
                            && (this->_current_iteration + 1) % Params::bayes_opt_cboptimizer::hp_period() == 0) {
187
188
189
                            _model.optimize_hyperparams();
                            if (_nb_constraints > 0)
                                _constraint_model.optimize_hyperparams();
190
191
192
193
194
195
196
197
198
199
200
                        }

                        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
                {
201
                    _dim_out = _model.dim_out();
202
203
                    _split_observations();

204
205
206
207
                    std::vector<Eigen::VectorXd> obs = _feasible_observations();
                    if (obs.size() > 0)
                        _obs[0] = obs;

208
209
210
211
212
213
214
215
216
217
                    auto rewards = std::vector<double>(_obs[0].size());
                    std::transform(_obs[0].begin(), _obs[0].end(), rewards.begin(), afun);
                    auto max_e = std::max_element(rewards.begin(), rewards.end());
                    return _obs[0][std::distance(rewards.begin(), max_e)];
                }

                /// return the best sample so far (i.e. the argmax(f(x)))
                template <typename AggregatorFunction = FirstElem>
                const Eigen::VectorXd& best_sample(const AggregatorFunction& afun = AggregatorFunction()) const
                {
218
                    _dim_out = _model.dim_out();
219
220
                    _split_observations();

221
222
223
224
                    std::vector<Eigen::VectorXd> obs = _feasible_observations();
                    if (obs.size() > 0)
                        _obs[0] = obs;

225
226
227
228
229
230
                    auto rewards = std::vector<double>(_obs[0].size());
                    std::transform(_obs[0].begin(), _obs[0].end(), rewards.begin(), afun);
                    auto max_e = std::max_element(rewards.begin(), rewards.end());
                    return this->_samples[std::distance(rewards.begin(), max_e)];
                }

231
                const model_t& model() const { return _model; } // returns model for the objective function
232
            protected:
233
234
235
                model_t _model;
                constraint_model_t _constraint_model;

236
237
238
239
                size_t _nb_constraints;
                mutable size_t _dim_out;
                mutable std::vector<std::vector<Eigen::VectorXd>> _obs;

240
241
242
243
244
245
246
247
248
249
250
                std::vector<Eigen::VectorXd> _feasible_observations() const
                {
                    std::vector<Eigen::VectorXd> feasible_obs;
                    for (size_t i = 0; i < _obs[0].size(); ++i) {
                        if (_obs[1][i].prod() > 0)
                            feasible_obs.push_back(_obs[0][i]);
                    }

                    return feasible_obs;
                }

251
252
253
                void _split_observations() const
                {
                    _obs.clear();
254
                    _obs.resize(2);
255
256

                    for (size_t i = 0; i < this->_observations.size(); ++i) {
257
                        assert(size_t(this->_observations[i].size()) == _dim_out + _nb_constraints);
258

259
                        Eigen::VectorXd vec_obj(_dim_out);
260
                        for (size_t j = 0; j < _dim_out; ++j)
261
262
263
264
265
                            vec_obj(j) = this->_observations[i](j);
                        _obs[0].push_back(vec_obj);

                        if (_nb_constraints > 0) {
                            Eigen::VectorXd vec_con(_nb_constraints);
266
                            for (int j = _dim_out, ind = 0; j < this->_observations[i].size(); ++j, ++ind)
267
268
269
                                vec_con(ind) = this->_observations[i](j);
                            _obs[1].push_back(vec_con);
                        }
270
271
272
273
274
275
276
277
                    }
                }
            };
        }
    }
}

#endif