bo_functions.cpp 11 KB
Newer Older
1
2
#include <cmath>
#include <algorithm>
3
4
5
6
#include <string>
#include <vector>
#include <utility>
#include <iostream>
7
8
9

#include <Eigen/Core>

10
11
12
13
#ifdef USE_TBB
#include <tbb/task_scheduler_init.h>
#include <tbb/parallel_for.h>
#include <tbb/concurrent_hash_map.h>
14
#else
15
#include <map>
16
17
#endif

18
19
20
#include <limbo/tools/macros.hpp>
#include <limbo/bayes_opt/boptimizer.hpp>
#include <limbo/tools/parallel.hpp>
21
22
23
24
25

using namespace limbo;

static constexpr int nb_replicates = 4;

26
namespace colors {
27
28
29
30
31
    static const char* red = "\33[31m";
    static const char* green = "\33[32m";
    static const char* yellow = "\33[33m";
    static const char* reset = "\33[0m";
    static const char* bold = "\33[1m";
32
33
34
}

// support functions
35
36
37
38
39
40
41
inline double sign(double x)
{
    if (x < 0)
        return -1;
    if (x > 0)
        return 1;
    return 0;
42
43
}

44
45
46
47
48
inline double hat(double x)
{
    if (x != 0)
        return log(fabs(x));
    return 0;
49
50
}

51
52
53
54
55
inline double c1(double x)
{
    if (x > 0)
        return 10;
    return 5.5;
56
57
}

58
59
60
61
62
inline double c2(double x)
{
    if (x > 0)
        return 7.9;
    return 3.1;
63
64
}

65
66
67
68
69
70
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;
71
72
}

73
74
75
76
77
Eigen::VectorXd make_v1(double x)
{
    Eigen::VectorXd v1(1);
    v1 << x;
    return v1;
78
79
}

80
struct Sphere {
81
82
    static constexpr size_t dim_in = 2;
    static constexpr size_t dim_out = 1;
83

84
85
86
87
88
    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
    {
        Eigen::Vector2d opt(0.5, 0.5);
        return make_v1(-(x - opt).squaredNorm());
    }
89
90
};

91
struct Ellipsoid {
92
93
94
95
96
97
98
99
100
101
102
103
    static constexpr size_t dim_in = 2;
    static constexpr size_t dim_out = 1;

    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
    {
        Eigen::Vector2d opt(0.5, 0.5);
        Eigen::Vector2d 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 make_v1(-r);
    }
104
105
};

106
struct Rastrigin {
107
108
109
110
111
112
113
114
115
116
117
    static constexpr size_t dim_in = 4;
    static constexpr size_t dim_out = 1;

    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
    {
        double f = 10 * x.size();
        for (int i = 0; i < x.size(); ++i)
            f += x(i) * x(i) - 10 * cos(2 * M_PI * x(i));
        return make_v1(-f);
    }
};
118
119

// see : http://www.sfu.ca/~ssurjano/hart3.html
120
struct Hartman3 {
121
122
123
124
125
126
127
128
129
130
131
132
133
    static constexpr size_t dim_in = 3;
    static constexpr size_t dim_out = 1;

    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
    {
        Eigen::Matrix<double, 4, 3> a, p;
        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;
        Eigen::Vector4d alpha;
        alpha << 1.0, 1.2, 3.0, 3.2;

        double res = 0;
134
        for (int i = 0; i < 4; i++) {
135
            double s = 0.0f;
136
            for (size_t j = 0; j < 3; j++) {
137
138
139
140
141
                s += a(i, j) * (x(j) - p(i, j)) * (x(j) - p(i, j));
            }
            res += alpha(i) * exp(-s);
        }
        return make_v1(res);
142
143
144
145
    }
};

// see : http://www.sfu.ca/~ssurjano/hart6.html
146
struct Hartman6 {
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    static constexpr size_t dim_in = 6;
    static constexpr size_t dim_out = 1;

    Eigen::VectorXd operator()(const Eigen::VectorXd& x) const
    {
        Eigen::Matrix<double, 4, 6> a, p;
        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;

        Eigen::Vector4d alpha;
        alpha << 1.0, 1.2, 3.0, 3.2;

        double res = 0;
163
        for (int i = 0; i < 4; i++) {
164
            double s = 0.0f;
165
            for (size_t j = 0; j < 6; j++) {
166
167
168
169
170
                s += a(i, j) * (x(j) - p(i, j)) * (x(j) - p(i, j));
            }
            res += alpha(i) * exp(-s);
        }
        return make_v1(res);
171
172
173
174
175
    }
};

// see : http://www.sfu.ca/~ssurjano/goldpr.html
// (with ln, as suggested in Jones et al.)
176
struct GoldenPrice {
177
178
    static constexpr size_t dim_in = 2;
    static constexpr size_t dim_out = 1;
179

180
181
182
183
    Eigen::VectorXd operator()(const Eigen::VectorXd& xx) const
    {
        Eigen::VectorXd x = (4.0 * xx).array() - 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)));
184

185
186
        return make_v1(-log(r) + 5);
    }
187
188
};

189
struct Params {
190
    struct bayes_opt_bobase {
191
        BO_PARAM(bool, stats_enabled, false);
192
193
    };

194
195
196
197
198
199
    struct bayes_opt_boptimizer {
        BO_PARAM(double, noise, 0.0);
    };

    struct init_randomsampling {
        BO_PARAM(int, samples, 10);
200
201
    };

202
203
    struct stop_maxiterations {
        BO_PARAM(int, iterations, 40);
204
205
    };

206
    struct acqui_gpucb : public defaults::acqui_gpucb {
207
208
    };

209
#ifdef USE_LIBCMAES
210
    struct opt_cmaes : public defaults::opt_cmaes {
211
    };
212
213
214
215
216
217
218
#elif defined(USE_NLOPT)
    struct opt_nloptnograd : public defaults::opt_nloptnograd {
    };
#else
    struct opt_gridsearch : public defaults::opt_gridsearch {
    };
#endif
219
    struct opt_rprop : public defaults::opt_rprop {
220
221
    };

222
    struct opt_parallelrepeater : public defaults::opt_parallelrepeater {
223
    };
224
};
225

226
227
228
229
template <typename T>
void print_res(const T& r)
{
    std::cout << "====== RESULTS ======" << std::endl;
230
231
    for (auto x : r) {
        for (auto y : x.second) {
232
            std::cout << x.first << "\t =>"
233
                      << " found :" << y.second << " expected " << y.first
234
235
                      << std::endl;
        }
236
        std::vector<std::pair<double, double>>& v = x.second;
237
        std::sort(v.begin(), v.end(),
238
239
            [](const std::pair<double, double>& x1,
                      const std::pair<double, double>& x2) {
240
                // clang-format off
241
                return x1.second < x2.second;
242
243
                // clang-format on
            });
244
245
        double med = v[v.size() / 2].second;
        if (fabs(v[0].first - med) < 0.05)
246
247
248
249
250
            std::cout << "[" << colors::green << "OK" << colors::reset << "] ";
        else
            std::cout << "[" << colors::red << "ERROR" << colors::reset << "] ";
        std::cout << colors::yellow << colors::bold << " -- " << x.first
                  << colors::reset << " ";
251
        std::cout << "Median: " << med << " error :" << fabs(v[0].first - med)
252
                  << std::endl;
253
254
255
    }
}

256
257
258
bool is_in_argv(int argc, char** argv, const char* needle)
{
    auto it = std::find_if(argv, argv + argc,
259
        [=](const char* s) {
260
261
262
263
264
            // clang-format off
            return strcmp(needle, s) == 0;
            // clang-format on
        });
    return !(it == argv + argc);
265
266
}

267
268
269
template <typename T1, typename T2>
void add_to_results(const char* key, T1& map, const T2& p)
{
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
270
#ifdef USE_TBB
271
272
273
    typename T1::accessor a;
    if (!map.find(a, key))
        map.insert(a, key);
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
274
275
276
277
278
279
#else
    typename T1::iterator a;
    a = map.find(key);
    if (a == map.end())
        map[key] = std::vector<std::pair<double, double>>();
#endif
280
    a->second.push_back(p);
281
282
}

283
int main(int argc, char** argv)
Federico Allocati's avatar
Federico Allocati committed
284
{
285
    tools::par::init();
286
287

#ifdef USE_TBB
288
    typedef tbb::concurrent_hash_map<std::string, std::vector<std::pair<double, double>>>
Federico Allocati's avatar
Federico Allocati committed
289
        res_t;
290
#else
291
    typedef std::map<std::string,
292
        std::vector<std::pair<double, double>>>
293
        res_t;
294
#endif
295
296
    res_t results;

297
298
    typedef bayes_opt::BOptimizer<Params> Opt_t;

299
    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "sphere"))
300
        tools::par::replicate(nb_replicates, [&]() {
301
302
303
304
                // clang-format off
                Opt_t opt;
                opt.optimize(Sphere());
                Eigen::Vector2d s_val(0.5, 0.5);
305
                double x_opt = FirstElem()(Sphere()(s_val));
306
                add_to_results("Sphere", results, std::make_pair(x_opt, opt.best_observation()(0)));
307
308
            // clang-format on
        });
309
310

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "ellipsoid"))
311
        tools::par::replicate(nb_replicates, [&]() {
312
313
314
315
                // clang-format off
                Opt_t opt;
                opt.optimize(Ellipsoid());
                Eigen::Vector2d s_val(0.5, 0.5);
316
                double x_opt = FirstElem()(Ellipsoid()(s_val));
317
                add_to_results("Ellipsoid", results, std::make_pair(x_opt, opt.best_observation()(0)));
318
319
            // clang-format on
        });
320
321

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "rastrigin"))
322
        tools::par::replicate(nb_replicates, [&]() {
323
324
325
326
                // clang-format off
                Opt_t opt;
                opt.optimize(Rastrigin());
                Eigen::Vector4d s_val(0, 0, 0, 0);
327
                double x_opt = FirstElem()(Rastrigin()(s_val));
328
                add_to_results("Rastrigin", results, std::make_pair(x_opt, opt.best_observation()(0)));
329
330
            // clang-format on
        });
331
332

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "hartman3"))
333
        tools::par::replicate(nb_replicates, [&]() {
334
335
336
337
338
                // clang-format off
                Opt_t opt;
                opt.optimize(Hartman3());
                // double s_max = 3.86278;
                Eigen::Vector3d s_val(0.114614, 0.555549, 0.852547);
339
                double x_opt = FirstElem()(Hartman3()(s_val));
340
                add_to_results("Hartman 3", results, std::make_pair(x_opt, opt.best_observation()(0)));
341
342
            // clang-format on
        });
343
344

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "hartman6"))
345
        tools::par::replicate(nb_replicates, [&]() {
346
347
348
349
350
351
                // clang-format off
                Opt_t opt;
                opt.optimize(Hartman6());
                Eigen::Matrix<double, 6, 1> s_val;
                s_val << 0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573;
                //double s_max = 3.32237;
352
                double x_opt = FirstElem()(Hartman6()(s_val));
353
                add_to_results("Hartman 6", results, std::make_pair(x_opt, opt.best_observation()(0)));
354
355
            // clang-format on
        });
356
357

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "golden_price"))
358
        tools::par::replicate(nb_replicates, [&]() {
359
360
361
362
363
                // clang-format off
                Opt_t opt;
                opt.optimize(GoldenPrice());
                //    double s_max = -log(3);
                Eigen::Vector2d s_val(0.5, 0.25);
364
                double x_opt = FirstElem()(GoldenPrice()(s_val));
365
                add_to_results("Golden Price", results, std::make_pair(x_opt, opt.best_observation()(0)));
366
367
            // clang-format on
        });
368
369
370
371
372

    std::cout << "Benchmark finished." << std::endl;

    print_res(results);
    return 0;
373
}