bo_functions.cpp 10.5 KB
Newer Older
1
2
#include <cmath>
#include <algorithm>
3
4
5
6
#include <string>
#include <vector>
#include <utility>
#include <iostream>
7
8
9
10
#ifdef USE_TBB
#include <tbb/task_scheduler_init.h>
#include <tbb/parallel_for.h>
#include <tbb/concurrent_hash_map.h>
11
12
#include <map>
#else
13
14
#endif

15
#include <Eigen/Core>
16

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

using namespace limbo;

static constexpr int nb_replicates = 4;

25
namespace colors {
26
27
28
29
30
    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";
31
32
33
}

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

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

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

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

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

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

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

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

90
struct Ellipsoid {
91
92
93
94
95
96
97
98
99
100
101
102
    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);
    }
103
104
};

105
struct Rastrigin {
106
107
108
109
110
111
112
113
114
115
116
    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);
    }
};
117
118

// see : http://www.sfu.ca/~ssurjano/hart3.html
119
struct Hartman3 {
120
121
122
123
124
125
126
127
128
129
130
131
132
    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;
133
        for (int i = 0; i < 4; i++) {
134
            double s = 0.0f;
135
            for (size_t j = 0; j < 3; j++) {
136
137
138
139
140
                s += a(i, j) * (x(j) - p(i, j)) * (x(j) - p(i, j));
            }
            res += alpha(i) * exp(-s);
        }
        return make_v1(res);
141
142
143
144
    }
};

// see : http://www.sfu.ca/~ssurjano/hart6.html
145
struct Hartman6 {
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    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;
162
        for (int i = 0; i < 4; i++) {
163
            double s = 0.0f;
164
            for (size_t j = 0; j < 6; j++) {
165
166
167
168
169
                s += a(i, j) * (x(j) - p(i, j)) * (x(j) - p(i, j));
            }
            res += alpha(i) * exp(-s);
        }
        return make_v1(res);
170
171
172
173
174
    }
};

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

179
180
181
182
    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)));
183

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

188
189
struct Params {
    struct boptimizer {
190
191
192
193
        BO_PARAM(double, noise, 0.0);
        BO_PARAM(int, dump_period, -1);
    };

194
    struct init {
195
196
197
        BO_PARAM(int, nb_samples, 10);
    };

198
    struct maxiterations {
199
200
201
        BO_PARAM(int, n_iterations, 40);
    };

202
    struct gp_ucb : public defaults::gp_ucb {
203
204
    };

205
    struct gp_auto : public defaults::gp_auto {
206
207
    };

208
    struct cmaes : public defaults::cmaes {
209
210
    };
};
211

212
213
214
215
template <typename T>
void print_res(const T& r)
{
    std::cout << "====== RESULTS ======" << std::endl;
216
217
    for (auto x : r) {
        for (auto y : x.second) {
218
            std::cout << x.first << "\t =>"
219
                      << " found :" << y.second << " expected " << y.first
220
221
                      << std::endl;
        }
222
        std::vector<std::pair<double, double>>& v = x.second;
223
        std::sort(v.begin(), v.end(),
224
225
            [](const std::pair<double, double>& x1,
                      const std::pair<double, double>& x2) {
226
                // clang-format off
227
                return x1.second < x2.second;
228
229
                // clang-format on
            });
230
231
        double med = v[v.size() / 2].second;
        if (fabs(v[0].first - med) < 0.05)
232
233
234
235
236
            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 << " ";
237
        std::cout << "Median: " << med << " error :" << fabs(v[0].first - med)
238
                  << std::endl;
239
240
241
    }
}

242
243
244
bool is_in_argv(int argc, char** argv, const char* needle)
{
    auto it = std::find_if(argv, argv + argc,
245
        [=](const char* s) {
246
247
248
249
250
            // clang-format off
            return strcmp(needle, s) == 0;
            // clang-format on
        });
    return !(it == argv + argc);
251
252
}

253
254
255
256
257
258
259
template <typename T1, typename T2>
void add_to_results(const char* key, T1& map, const T2& p)
{
    typename T1::accessor a;
    if (!map.find(a, key))
        map.insert(a, key);
    a->second.push_back(p);
260
261
}

262
int main(int argc, char** argv)
263
{    
264
    tools::par::init();
265
266

#ifdef USE_TBB
267
    typedef tbb::concurrent_hash_map<std::string, std::vector<std::pair<double, double>>>
268
       res_t;
269
#else
270
    typedef std::map<std::string,
271
        std::vector<std::pair<double, double>>>
272
        res_t;
273
#endif
274
275
    res_t results;

276
277
    typedef bayes_opt::BOptimizer<Params> Opt_t;

278
    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "sphere"))
279
        tools::par::replicate(nb_replicates, [&]() {
280
281
282
283
                // clang-format off
                Opt_t opt;
                opt.optimize(Sphere());
                Eigen::Vector2d s_val(0.5, 0.5);
284
                double x_opt = FirstElem()(Sphere()(s_val));
285
                add_to_results("Sphere", results, std::make_pair(x_opt, opt.best_observation()));
286
287
            // clang-format on
        });
288
289

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "ellipsoid"))
290
        tools::par::replicate(nb_replicates, [&]() {
291
292
293
294
                // clang-format off
                Opt_t opt;
                opt.optimize(Ellipsoid());
                Eigen::Vector2d s_val(0.5, 0.5);
295
                double x_opt = FirstElem()(Ellipsoid()(s_val));
296
                add_to_results("Ellipsoid", results, std::make_pair(x_opt, opt.best_observation()));
297
298
            // clang-format on
        });
299
300

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

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "hartman3"))
312
        tools::par::replicate(nb_replicates, [&]() {
313
314
315
316
317
                // clang-format off
                Opt_t opt;
                opt.optimize(Hartman3());
                // double s_max = 3.86278;
                Eigen::Vector3d s_val(0.114614, 0.555549, 0.852547);
318
                double x_opt = FirstElem()(Hartman3()(s_val));
319
                add_to_results("Hartman 3", results, std::make_pair(x_opt, opt.best_observation()));
320
321
            // clang-format on
        });
322
323

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "hartman6"))
324
        tools::par::replicate(nb_replicates, [&]() {
325
326
327
328
329
330
                // 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;
331
                double x_opt = FirstElem()(Hartman6()(s_val));
332
                add_to_results("Hartman 6", results, std::make_pair(x_opt, opt.best_observation()));
333
334
            // clang-format on
        });
335
336

    if (!is_in_argv(argc, argv, "--only") || is_in_argv(argc, argv, "golden_price"))
337
        tools::par::replicate(nb_replicates, [&]() {
338
339
340
341
342
                // clang-format off
                Opt_t opt;
                opt.optimize(GoldenPrice());
                //    double s_max = -log(3);
                Eigen::Vector2d s_val(0.5, 0.25);
343
                double x_opt = FirstElem()(GoldenPrice()(s_val));
344
                add_to_results("Golden Price", results, std::make_pair(x_opt, opt.best_observation()));
345
346
            // clang-format on
        });
347
348
349
350
351

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

    print_res(results);
    return 0;
352
}