testfunctions.hpp 11 KB
Newer Older
1
2
3
4
//| 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
5
//|
6
7
8
//| 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
//| 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
19
//|
20
21
22
23
24
//| 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".
25
//|
26
27
28
29
30
//| 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.
31
//|
32
33
34
35
36
37
38
39
40
41
//| 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.
42
//|
43
44
//| 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.
45
//|
46
47
48
#define _USE_MATH_DEFINES
#include <Eigen/Core>

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#ifdef BAYES_OPT
#include <bayesopt/bayesopt.hpp>
#include <boost/numeric/ublas/assignment.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <limbo/tools/macros.hpp>

using namespace bayesopt;

using vec_t = vectord;
using mat_t = matrixd;
#define ASSIGNMENT_OP <<=
#else
using vec_t = Eigen::VectorXd;
using mat_t = Eigen::MatrixXd;
#define ASSIGNMENT_OP <<
#endif

66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
// 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)
84
        return std::log(std::abs(x));
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
    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;
}

102
inline vec_t t_osz(const vec_t& x)
103
{
104
    vec_t r = x;
105
    for (int i = 0; i < x.size(); i++)
106
        r(i) = sign(x(i)) * std::exp(hat(x(i)) + 0.049 * std::sin(c1(x(i)) * hat(x(i))) + std::sin(c2(x(i)) * hat(x(i))));
107
108
109
110
    return r;
}

struct Sphere {
111
112
    BO_PARAM(size_t, dim_in, 2);
    BO_PARAM(size_t, dim_out, 1);
113

114
    double operator()(const vec_t& x) const
115
    {
116
117
        vec_t opt(2);
        opt ASSIGNMENT_OP 0.5, 0.5;
118

119
#ifndef BAYES_OPT
120
        return (x - opt).squaredNorm();
121
122
123
#else
        return sqr(norm_2(x - opt));
#endif
124
125
    }

126
    mat_t solutions() const
127
    {
128
129
        mat_t sols(1, 2);
        sols ASSIGNMENT_OP 0.5, 0.5;
130
131
132
133
134
        return sols;
    }
};

struct Ellipsoid {
135
136
    BO_PARAM(size_t, dim_in, 2);
    BO_PARAM(size_t, dim_out, 1);
137

138
    double operator()(const vec_t& x) const
139
    {
140
141
142
        vec_t opt(2);
        opt ASSIGNMENT_OP 0.5, 0.5;
        vec_t z = t_osz(x - opt);
143
        double r = 0;
144
        for (size_t i = 0; i < dim_in(); ++i)
145
            r += std::pow(10., ((double)i) / (dim_in() - 1.0)) * z(i) * z(i) + 1;
146
147
148
        return r;
    }

149
    mat_t solutions() const
150
    {
151
152
        mat_t sols(1, 2);
        sols ASSIGNMENT_OP 0.5, 0.5;
153
154
155
156
157
        return sols;
    }
};

struct Rastrigin {
158
159
    BO_PARAM(size_t, dim_in, 4);
    BO_PARAM(size_t, dim_out, 1);
160

161
    double operator()(const vec_t& xx) const
162
    {
163
164
165
166
        vec_t x = xx;
        for (int i = 0; i < x.size(); i++)
            x(i) = 2. * xx(i) - 1.;
        double f = 10. * dim_in();
167
        for (size_t i = 0; i < dim_in(); ++i)
168
            f += x(i) * x(i) - 10. * std::cos(2 * M_PI * x(i));
169
        return f;
170
171
    }

172
    mat_t solutions() const
173
    {
174
175
176
177
178
179
180
#ifndef BAYES_OPT
        mat_t sols = Eigen::MatrixXd::Zero(1, 4);
#else
        mat_t sols = boost::numeric::ublas::zero_matrix<double>(1, 4);
#endif
        for (int i = 0; i < 4; i++)
            sols(0, i) = (sols(0, i) + 1.) / 2.;
181
        return sols;
182
183
184
185
    }
};

// see : http://www.sfu.ca/~ssurjano/hart3.html
186
struct Hartmann3 {
187
188
    BO_PARAM(size_t, dim_in, 3);
    BO_PARAM(size_t, dim_out, 1);
189

190
    double operator()(const vec_t& x) const
191
    {
192
193
194
195
196
197
198
199
200
201
202
203
204
205
        mat_t a(4, 3);
        mat_t p(4, 3);
        a ASSIGNMENT_OP 3.0, 10., 30.,
            0.1, 10., 35.,
            3.0, 10., 30.,
            0.1, 10., 35.;
        p ASSIGNMENT_OP 0.3689, 0.1170, 0.2673,
            0.4699, 0.4387, 0.7470,
            0.1091, 0.8732, 0.5547,
            0.0381, 0.5743, 0.8828;
        vec_t alpha(4);
        alpha ASSIGNMENT_OP 1.0, 1.2, 3.0, 3.2;

        double res = 0.;
206
        for (int i = 0; i < 4; i++) {
207
            double s = 0.;
208
            for (size_t j = 0; j < 3; j++) {
209
                s += a(i, j) * sqr(x(j) - p(i, j));
210
            }
211
            res += alpha(i) * std::exp(-s);
212
213
214
215
        }
        return -res;
    }

216
    mat_t solutions() const
217
    {
218
219
        mat_t sols(1, 3);
        sols ASSIGNMENT_OP 0.114614, 0.555649, 0.852547;
220
221
222
223
224
        return sols;
    }
};

// see : http://www.sfu.ca/~ssurjano/hart6.html
225
struct Hartmann6 {
226
227
    BO_PARAM(size_t, dim_in, 6);
    BO_PARAM(size_t, dim_out, 1);
228

229
    double operator()(const vec_t& x) const
230
    {
231
232
233
234
235
236
237
238
239
        mat_t a(4, 6);
        mat_t p(4, 6);
        a ASSIGNMENT_OP 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 ASSIGNMENT_OP 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.6650,
240
241
            0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381;

242
243
        vec_t alpha(4);
        alpha ASSIGNMENT_OP 1.0, 1.2, 3.0, 3.2;
244

245
        double res = 0.;
246
        for (int i = 0; i < 4; i++) {
247
            double s = 0.;
248
249
250
            for (size_t j = 0; j < 6; j++) {
                s += a(i, j) * sqr(x(j) - p(i, j));
            }
251
            res += alpha(i) * std::exp(-s);
252
253
254
255
        }
        return -res;
    }

256
    mat_t solutions() const
257
    {
258
259
        mat_t sols(1, 6);
        sols ASSIGNMENT_OP 0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573;
260
261
262
263
264
265
        return sols;
    }
};

// see : http://www.sfu.ca/~ssurjano/goldpr.html
// (with ln, as suggested in Jones et al.)
266
struct GoldsteinPrice {
267
268
    BO_PARAM(size_t, dim_in, 2);
    BO_PARAM(size_t, dim_out, 1);
269

270
    double operator()(const vec_t& xx) const
271
    {
272
273
274
275
276
277
278
279
280
281
282
        vec_t x = xx;
        for (int i = 0; i < x.size(); i++)
            x(i) = 4. * xx(i) - 2.;

        double fact1a = sqr(x(0) + x(1) + 1.);
        double fact1b = 19. - 14. * x(0) + 3. * sqr(x(0)) - 14. * x(1) + 6. * x(0) * x(1) + 3. * sqr(x(1));
        double fact1 = 1. + fact1a * fact1b;

        double fact2a = sqr(2. * x(0) - 3. * x(1));
        double fact2b = 18. - 32. * x(0) + 12. * sqr(x(0)) + 48. * x(1) - 36. * x(0) * x(1) + 27. * sqr(x(1));
        double fact2 = 30. + fact2a * fact2b;
283

284
285
286
287
        double r = fact1 * fact2;

        return (std::log(r) - 8.693) / 2.427;
        // return std::log(r) - 5.;
288
289
    }

290
    mat_t solutions() const
291
    {
292
293
        mat_t sols(1, 2);
        sols ASSIGNMENT_OP 0.5, 0.25;
294
295
296
297
298
        return sols;
    }
};

struct BraninNormalized {
299
300
    BO_PARAM(size_t, dim_in, 2);
    BO_PARAM(size_t, dim_out, 1);
301

302
    double operator()(const vec_t& x) const
303
    {
304
305
306
307
308
309
310
        double x1 = x(0) * 15 - 5;
        double x2 = x(1) * 15;

        double term1 = sqr(x2 - (5.1 * sqr(x1) / (4. * sqr(M_PI))) + 5. * x1 / M_PI - 6);
        double term2 = (10. - 10. / (8. * M_PI)) * std::cos(x1);

        return (term1 + term2 - 44.81) / 51.95;
311
312
    }

313
    mat_t solutions() const
314
    {
315
316
317
318
319
320
321
322
323
324
325
        mat_t sols(3, 2);
        sols ASSIGNMENT_OP - M_PI, 12.275,
            M_PI, 2.275,
            9.42478, 2.475;

        sols(0, 0) = (sols(0, 0) + 5.) / 15.;
        sols(1, 0) = (sols(1, 0) + 5.) / 15.;
        sols(2, 0) = (sols(2, 0) + 5.) / 15.;
        sols(0, 1) = sols(0, 1) / 15.;
        sols(1, 1) = sols(1, 1) / 15.;
        sols(2, 1) = sols(2, 1) / 15.;
326
327
328
329
330
        return sols;
    }
};

struct SixHumpCamel {
331
332
333
    BO_PARAM(size_t, dim_in, 2);
    BO_PARAM(size_t, dim_out, 1);

334
    double operator()(const vec_t& x) const
335
    {
336
337
        double x1 = -3 + 6 * x(0);
        double x2 = -2 + 4 * x(1);
338
339
        double x1_2 = sqr(x1);
        double x2_2 = sqr(x2);
340

341
        double tmp1 = (4 - 2.1 * x1_2 + sqr(x1_2) / 3.) * x1_2;
342
        double tmp2 = x1 * x2;
343
344
345
346
        double tmp3 = (-4 + 4 * x2_2) * x2_2;
        return tmp1 + tmp2 + tmp3;
    }

347
    mat_t solutions() const
348
    {
349
350
        mat_t sols(2, 2);
        sols ASSIGNMENT_OP 0.0898, -0.7126,
351
            -0.0898, 0.7126;
352
353
354
355
        sols(0, 0) = (sols(0, 0) + 3.) / 6.;
        sols(1, 0) = (sols(1, 0) + 3.) / 6.;
        sols(0, 1) = (sols(0, 1) + 2.) / 4.;
        sols(1, 1) = (sols(1, 1) + 2.) / 4.;
356
357
358
359
        return sols;
    }
};

360
#ifndef BAYES_OPT
361
362
363
template <typename Function>
class Benchmark {
public:
364
365
    BO_PARAM(size_t, dim_in, Function::dim_in());
    BO_PARAM(size_t, dim_out, Function::dim_out());
366

367
    vec_t operator()(const vec_t& x) const
368
    {
369
        vec_t res(1);
370
371
372
        res(0) = -f(x);
        return res;
    }
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
#else
template <typename Function>
class Benchmark : public bayesopt::ContinuousModel {
public:
    Benchmark(bopt_params par) : ContinuousModel(Function::dim_in(), par) {}

    double evaluateSample(const vec_t& xin)
    {
        return f(xin);
    }

    bool checkReachability(const vec_t& query)
    {
        return true;
    };
#endif
389

390
    double accuracy(double x)
391
    {
392
393
        mat_t sols = f.solutions();
#ifndef BAYES_OPT
394
        double diff = std::abs(x + f(sols.row(0)));
395
396
397
#else
        double diff = std::abs(x - f(row(sols, 0)));
#endif
398
399
        double min_diff = diff;

400
#ifndef BAYES_OPT
401
        for (int i = 1; i < sols.rows(); i++) {
402
            diff = std::abs(x + f(sols.row(i)));
403
404
405
406
#else
        for (size_t i = 1; i < sols.size1(); i++) {
            diff = std::abs(x - f(row(sols, i)));
#endif
407
408
409
410
411
412
413
414
415
            if (diff < min_diff)
                min_diff = diff;
        }

        return min_diff;
    }

    Function f;
};