chiark / gitweb /
octave4.4
[nlopt.git] / test / testfuncs.c
1 #include <stdio.h>
2 #include <math.h>
3
4 #include "testfuncs.h"
5 #include "config.h"
6
7 #define UNUSED(x) (void) x
8
9 static double sqr(double x) { return x * x; }
10
11 int testfuncs_verbose = 0;
12 int testfuncs_counter = 0;
13
14 static double testfuncs_status(unsigned n, const double *x, double f)
15 {
16      ++testfuncs_counter;
17      if (testfuncs_verbose) {
18           unsigned i;
19           printf("f_%d (%g", testfuncs_counter, x[0]);
20           for (i = 1; i < n; ++i) printf(", %g", x[i]);
21           printf(") = %g\n", f);
22      }
23      return f;
24 }
25
26 #define RETURN(f) return testfuncs_status(n, x, f);
27
28 #define PI2 6.283185307179586 /* 2*pi */
29 #define PI3 9.424777960769379 /* 3*pi */
30 #define PI4 12.5663706143592 /* 4*pi */
31
32 /****************************************************************************/
33 static double rosenbrock_f(unsigned n, const double *x, double *grad, void *data)
34 {
35      double a = x[1] - x[0] * x[0], b = 1 - x[0];
36      UNUSED(data);
37      if (grad) {
38           grad[0] = -400 * a * x[0] - 2*b;
39           grad[1] = 200 * a;
40      }
41      RETURN(100 * sqr(a) + sqr(b));
42 }
43
44 static const double rosenbrock_lb[2] = {-2, -2};
45 static const double rosenbrock_ub[2] = {2, 2};
46 static const double rosenbrock_xmin[2] = {1, 1};
47
48 /****************************************************************************/
49 static double mccormic_f(unsigned n, const double *x, double *grad, void *data)
50 {
51      double a = x[0] + x[1], b = x[0] - x[1];
52      UNUSED(data);
53      if (grad) {
54           grad[0] = cos(a) + 2*b - 1.5;
55           grad[1] = cos(a) - 2*b + 2.5;
56      }
57      RETURN(sin(a) + sqr(b) - 1.5*x[0] + 2.5*x[1] + 1);
58 }
59
60 static const double mccormic_lb[2] = {-1.5, -3};
61 static const double mccormic_ub[2] = {4, 4};
62 static const double mccormic_xmin[2] = {-0.547197553, -1.54719756};
63
64 /****************************************************************************/
65 static double boxbetts_f(unsigned n, const double *x, double *grad, void *data)
66 {
67      unsigned i;
68      double f = 0;
69      UNUSED(data);
70      if (grad)
71           grad[0] = grad[1] = grad[2] = 0;
72      for (i = 1; i <= 10; ++i) {
73           double e0 = exp(-0.1*i*x[0]);
74           double e1 = exp(-0.1*i*x[1]);
75           double e2 = exp(-0.1*i) - exp((double) -i);
76           double g = e0 - e1 - e2 * x[2];
77           f += sqr(g);
78           if (grad) {
79                grad[0] += (2 * g) * (-0.1*i*e0);
80                grad[1] += (2 * g) * (0.1*i*e1);
81                grad[2] += -(2 * g) * e2;
82           }
83      }
84      RETURN(f);
85 }
86
87 static const double boxbetts_lb[3] = {0.9,9,0.9};
88 static const double boxbetts_ub[3] = {1.2,11.2,1.2};
89 static const double boxbetts_xmin[3] = {1,10,1};
90
91 /****************************************************************************/
92 static double paviani_f(unsigned n, const double *x, double *grad, void *data)
93 {
94      unsigned i;
95      double f = 0, prod = 1;
96      UNUSED(data);
97      if (grad) for (i = 0; i < 10; ++i) grad[i] = 0;
98      for (i = 0; i < 10; ++i) {
99           double ln1 = log(x[i] - 2);
100           double ln2 = log(10 - x[i]);
101           f += sqr(ln1) + sqr(ln2);
102           if (grad)
103                grad[i] += 2 * ln1 / (x[i] - 2) - 2 * ln2 / (10 - x[i]);
104           prod *= x[i];
105      }
106      f -= (prod = pow(prod, 0.2));
107      if (grad)
108           for (i = 0; i < 10; ++i)
109                grad[i] -= 0.2 * prod / x[i];
110      RETURN(f);
111 }
112
113 static const double paviani_lb[10] = {2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001};
114 static const double paviani_ub[10] = {9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999};
115 static const double paviani_xmin[10] = {9.35026583,9.35026583,9.35026583,9.35026583,9.35026583,9.35026583,9.35026583,9.35026583,9.35026583,9.35026583};
116
117 /****************************************************************************/
118 static double grosenbrock_f(unsigned n, const double *x, double *grad, void *data)
119 {
120      unsigned i;
121      double f = 0;
122      UNUSED(data);
123      if (grad) grad[0] = 0;
124      for (i = 0; i < 29; ++i) {
125           double a = x[i+1] - x[i] * x[i], b = 1 - x[i];
126           if (grad) {
127                grad[i] += -400 * a * x[i] - 2*b;
128                grad[i+1] = 200 * a;
129           }
130           f += 100 * sqr(a) + sqr(b);
131      }
132      RETURN(f);
133 }
134
135 static const double grosenbrock_lb[30] = {-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30};
136 static const double grosenbrock_ub[30] = {30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30};
137 static const double grosenbrock_xmin[30] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
138
139 /****************************************************************************/
140 static double goldsteinprice_f(unsigned n, const double *x, double *grad, void *data)
141 {
142      double x0, x1, a1, a12, a2, b1, b12, b2;
143      UNUSED(data);
144      x0 = x[0]; x1 = x[1];
145      a1 = x0+x1+1; a12 = sqr(a1);
146      a2 = 19 - 14*x0 + 3*x0*x0 - 14*x1 + 6*x0*x1 + 3*x1*x1;
147      b1 = 2*x0-3*x1; b12 = sqr(b1);
148      b2 = 18 - 32*x0 + 12*x0*x0 + 48*x1 - 36*x0*x1 + 27*x1*x1;
149      if (grad) {
150           grad[0] = (1 + a12 * a2) * (2 * b1 * 2 * b2 
151                                       + b12 * (-32 + 24*x0 - 36*x1))
152                + (2 * a1 * a2 + a12 * (-14 + 6*x0 + 6*x1)) * (30 + b12 * b2);
153           grad[1] = (1 + a12 * a2) * (2 * b1 * (-3) * b2
154                                       + b12 * (48 - 36*x0 + 54 * x1))
155                + (2 * a1 * a2 + a12 * (-14 + 6*x0 + 6*x1)) * (30 + b12 * b2);
156      }
157      RETURN((1 + a12 * a2) * (30 + b12 * b2));
158 }
159
160 static const double goldsteinprice_lb[2] = {-2, -2};
161 static const double goldsteinprice_ub[2] = {2, 2};
162 static const double goldsteinprice_xmin[2] = {0, -1};
163
164 /****************************************************************************/
165 static double shekel_f(unsigned n, const double *x, double *grad, void *data)
166 {
167      static const double A[10][4] = { {4,4,4,4},
168                                       {1,1,1,1},
169                                       {8,8,8,8},
170                                       {6,6,6,6},
171                                       {3,7,3,7},
172                                       {2,9,2,9},
173                                       {5,5,3,3},
174                                       {8,1,8,1},
175                                       {6,2,6,2},
176                                       {7,3.6,7,3.6} };
177      static const double c[10] = {.1,.2,.2,.4,.4,.6,.3,.7,.5,.5};
178      unsigned i;
179      double f = 0;
180      if (grad) for (i = 0; i < n; ++i) grad[i] = 0;
181      unsigned m = *((unsigned *) data);
182      for (i = 0; i < m; ++i) {
183           double fi = 1.0 / (c[i] 
184                              + sqr(x[0]-A[i][0])
185                              + sqr(x[1]-A[i][1])
186                              + sqr(x[2]-A[i][2])
187                              + sqr(x[3]-A[i][3]));
188           f -= fi;
189           if (grad) {
190                grad[0] += (2*fi*fi) * (x[0]-A[i][0]);
191                grad[1] += (2*fi*fi) * (x[1]-A[i][1]);
192                grad[2] += (2*fi*fi) * (x[2]-A[i][2]);
193                grad[3] += (2*fi*fi) * (x[3]-A[i][3]);
194           }
195      }
196      RETURN(f);
197 }
198
199 static unsigned shekel_m[3] = {5,7,10};
200 static const double shekel_lb[4] = {0,0,0,0};
201 static const double shekel_ub[4] = {10,10,10,10};
202 static const double shekel0_xmin[4] = {4.000037154,4.000133276,4.000037154,4.000133276};
203 static const double shekel1_xmin[4] = {4.000572917,4.000689366,3.999489709,3.999606158};
204 static const double shekel2_xmin[4] = {4.000746531,4.000592935,3.999663399,3.999509801};
205
206 /****************************************************************************/
207 static double levy_f(unsigned n, const double *x, double *grad, void *data)
208 {
209      UNUSED(data);
210      unsigned i;
211      double a = x[n-1] - 1, b = 1 + sqr(sin(PI2*x[n-1]));
212      double f = sqr(sin(PI3*x[0])) + a * b;
213      if (grad) {
214           for (i = 0; i < n; ++i) grad[i] = 0;
215           grad[0] = 2 * PI3 * sin(PI3*x[0]) * cos(PI3*x[0]);
216           grad[n-1] += b + a * 2 * PI2 * sin(PI2*x[n-1]) * cos(PI2*x[n-1]);
217      }
218      for (i = 0; i < n-1; ++i) {
219           a = x[i] - 1;
220           b = 1 + sqr(sin(PI3*x[i+1]));
221           f += sqr(a) * b;
222           if (grad) {
223                grad[i] += 2 * a * b;
224                grad[i+1] += 2*PI3 * sqr(a) * sin(PI3*x[i+1])*cos(PI3*x[i+1]);
225           }
226      }
227      RETURN(f);
228 }
229
230 static const double levy_lb[7] = {-5,-5,-5,-5,-5,-5,-5};
231 static const double levy_ub[7] = {5,5,5,5,5,5,5};
232 static const double levy_xmin[7] = {1,1,1,1,1,1,-4.75440246};
233 static const double levy4_lb[4] = {-10,-10,-10,-10};
234 static const double levy4_ub[4] = {10,10,10,10};
235 static const double levy4_xmin[4] = {1,1,1,-9.75235596};
236
237 /****************************************************************************/
238 static double griewank_f(unsigned n, const double *x, double *grad, void *data)
239 {
240      unsigned i;
241      double f = 1, p = 1;
242      UNUSED(data);
243      for (i = 0; i < n; ++i) {
244           f += sqr(x[i]) * 0.00025;
245           p *= cos(x[i] / sqrt(i + 1.));
246           if (grad) grad[i] = x[i] * 0.0005;
247      }
248      f -= p;
249      if (grad)
250           for (i = 0; i < n; ++i)
251                grad[i] += p * tan(x[i] / sqrt(i + 1.)) / sqrt(i + 1.);
252      RETURN(f);
253 }
254
255 static const double griewank_lb[10] = {-500,-500,-500,-500,-500,-500,-500,-500,-500,-500};
256 static const double griewank_ub[10] = {600,600,600,600,600,600,600,600,600,600};
257 static const double griewank_xmin[10] = {0,0,0,0,0,0,0,0,0,0};
258
259 /****************************************************************************/
260 static double sixhumpcamel_f(unsigned n, const double *x, double *grad, void *data)
261 {
262      UNUSED(data);
263      if (grad) {
264           grad[0] = 8*x[0] - 2.1*4*pow(x[0],3.) + 2*pow(x[0],5.) + x[1];
265           grad[1] = x[0] - 8*x[1] + 16*pow(x[1],3.);
266      }
267      RETURN(4*sqr(x[0]) - 2.1 * pow(x[0],4.) + pow(x[0],6.)/3. 
268             + x[0]*x[1] - 4*sqr(x[1]) + 4*pow(x[1],4.));
269 }
270
271 static const double sixhumpcamel_lb[2] = {-5,-5};
272 static const double sixhumpcamel_ub[2] = {5,5};
273 static const double sixhumpcamel_xmin[2] = {0.08984201317, -0.7126564032};
274
275 /****************************************************************************/
276 static double convexcosh_f(unsigned n, const double *x, double *grad, void *data)
277 {
278      unsigned i;
279      double f = 1;
280      UNUSED(data);
281      for (i = 0; i < n; ++i)
282           f *= cosh((x[i] - i) * (i+1));
283      if (grad)
284           for (i = 0; i < n; ++i)
285                grad[i] = f * tanh((x[i] - i) * (i+1)) * (i+1);
286      RETURN(f);
287 }
288
289 static const double convexcosh_lb[10] = {-1,0,0,0,0,0,0,0,0,0};
290 static const double convexcosh_ub[10] = {2,3,6,7,8,10,11,13,14,16};
291 static const double convexcosh_xmin[10] = {0,1,2,3,4,5,6,7,8,9};
292
293 /****************************************************************************/
294 static double branin_f(unsigned n, const double *x, double *grad, void *data)
295 {
296      double a = 1 - 2*x[1] + 0.05 * sin(PI4 * x[1]) - x[0];
297      double b = x[1] - 0.5 * sin(PI2 * x[0]);
298      UNUSED(data);
299      if (grad) {
300           grad[0] = -2*a - cos(PI2 * x[0]) * PI2 * b;
301           grad[1] = 2*a*(0.05 * PI4 * cos(PI4*x[1]) - 2) + 2*b;
302      }
303      RETURN(sqr(a) + sqr(b));
304 }
305
306 static const double branin_lb[2] = {-10,-10};
307 static const double branin_ub[2] = {10,10};
308 static const double branin_xmin[2] = {1,0};
309
310 /****************************************************************************/
311 static double shubert_f(unsigned n, const double *x, double *grad, void *data)
312 {
313      UNUSED(data);
314      unsigned i, j;
315      double f = 0;
316      for (j = 1; j <= 5; ++j)
317           for (i = 0; i < n; ++i)
318                f -= j * sin((j+1) * x[i] + j);
319      if (grad) {
320           for (i = 0; i < n; ++i) {
321                grad[i] = 0;
322                for (j = 1; j <= 5; ++j)
323                     grad[i] -= j * (j+1) * cos((j+1) * x[i] + j);
324           }
325      }
326      RETURN(f);
327 }
328
329 static const double shubert_lb[2] = {-10,-10};
330 static const double shubert_ub[2] = {10,10};
331 static const double shubert_xmin[2] = {-6.774576143, -6.774576143};
332
333 /****************************************************************************/
334 static double hansen_f(unsigned n, const double *x, double *grad, void *data)
335 {
336      unsigned i;
337      double a = 0, b = 0;
338      UNUSED(data);
339      for (i = 1; i <= 5; ++i)
340           a += i * cos((i-1) * x[0] + i);
341      for (i = 1; i <= 5; ++i)
342           b += i * cos((i+1) * x[1] + i);
343      if (grad) {
344           grad[0] = 0;
345           for (i = 1; i <= 5; ++i)
346                grad[0] -= i * (i-1) * sin((i-1) * x[0] + i);
347           grad[0] *= b;
348           grad[1] = 0;
349           for (i = 1; i <= 5; ++i)
350                grad[1] -= i * (i+1) * sin((i+1) * x[1] + i);
351           grad[1] *= a;
352      }
353      RETURN(a*b);
354 }
355
356 static const double hansen_lb[2] = {-10,-10};
357 static const double hansen_ub[2] = {10,10};
358 static const double hansen_xmin[2] = {-1.306707704,-1.425128429};
359
360 /****************************************************************************/
361 static double osc1d_f(unsigned n, const double *x, double *grad, void *data)
362 {
363      double y = *x - 1.23456789;
364      UNUSED(data);
365      if (grad) grad[0] = y*0.02 + sin(y - 2*sin(3*y)) * (1 - 6*cos(3*y));
366      RETURN(sqr(y*0.1) - cos(y - 2*sin(3*y)));
367 }
368
369 static const double osc1d_lb[1] = {-5};
370 static const double osc1d_ub[1] = {5};
371 static const double osc1d_xmin[1] = {1.23456789};
372
373 /****************************************************************************/
374 static double corner4d_f(unsigned n, const double *x, double *grad, void *data)
375 {
376      UNUSED(data);
377      UNUSED(n);
378      double u = x[0] + x[1] * x[2] * sin(2 * x[3]);
379      double v = x[0] + 2*sin(u);
380      if (grad) {
381           grad[0] = 2*v * (1 + 2*cos(u));
382           grad[1] = 2*v * 2*cos(u) * x[2] * sin(2*x[3]) + 0.1;
383           grad[2] = 2*v * 2*cos(u) * x[1] * sin(2*x[3]) + 0.1;
384           grad[3] = 2*v * 2*cos(u) * x[1]*x[2] * cos(2*x[3]) * 2 + 0.1;
385      }
386      RETURN(1 + v*v + 0.1*(x[1]+x[2]+x[3]));
387 }
388
389 static const double corner4d_lb[4] = {0,0,0,0};
390 static const double corner4d_ub[4] = {1,1,1,1};
391 static const double corner4d_xmin[4] = {0,0,0,0};
392
393 /****************************************************************************/
394 static double side4d_f(unsigned n, const double *x, double *grad, void *data)
395 {
396      UNUSED(data);
397      UNUSED(n);
398      double x0, x1, x2, x3, d0,d1,d2,d3;
399      const double w0 = 0.1, w1 = 0.2, w2 = 0.3, w3 = 0.4;
400      x0 = +0.4977 * x[0] - 0.3153 * x[1] - 0.5066 * x[2] - 0.4391 * x[3];
401      x1 = -0.3153 * x[0] + 0.3248 * x[1] - 0.4382 * x[2] - 0.4096 * x[3];
402      x2 = -0.5066 * x[0] - 0.4382 * x[1] + 0.3807 * x[2] - 0.4543 * x[3];
403      x3 = -0.4391 * x[0] - 0.4096 * x[1] - 0.4543 * x[2] + 0.5667 * x[3];
404
405      d0 = -1. / (x0*x0 + w0*w0);
406      d1 = -1. / (x1*x1 + w1*w1);
407      d2 = -1. / (x2*x2 + w2*w2);
408      d3 = -1. / (x3*x3 + w3*w3);
409
410      if (grad) {
411           grad[0] = 2 * (x0*d0*d0 * +0.4977 +
412                          x1*d1*d1 * -0.3153 +
413                          x2*d2*d2 * -0.5066 +
414                          x3*d3*d3 * -0.4391);
415           grad[1] = 2 * (x0*d0*d0 * -0.3153 +
416                          x1*d1*d1 * +0.3248 +
417                          x2*d2*d2 * -0.4382 +
418                          x3*d3*d3 * -0.4096);
419           grad[2] = 2 * (x0*d0*d0 * -0.5066 +
420                          x1*d1*d1 * -0.4382 +
421                          x2*d2*d2 * +0.3807 +
422                          x3*d3*d3 * -0.4543);
423           grad[3] = 2 * (x0*d0*d0 * -0.4391 +
424                          x1*d1*d1 * -0.4096 +
425                          x2*d2*d2 * -0.4543 +
426                          x3*d3*d3 * +0.5667);
427      }
428      RETURN(d0 + d1 + d2 + d3);
429 }
430
431 static const double side4d_lb[4] = {0.1,-1,-1,-1};
432 static const double side4d_ub[4] = {1,1,1,1};
433 static const double side4d_xmin[4] = {0.1,0.102971169,0.0760520641,-0.0497098571};
434
435 /****************************************************************************/
436 /****************************************************************************/
437
438 const testfunc testfuncs[NTESTFUNCS] = {
439      { rosenbrock_f, NULL, 1, 2,
440        rosenbrock_lb, rosenbrock_ub, rosenbrock_xmin,
441        0.0, "Rosenbrock function" },
442      { mccormic_f, NULL, 1, 2,
443        mccormic_lb, mccormic_ub, mccormic_xmin,
444        -1.91322295, "McCormic function" },
445      { boxbetts_f, NULL, 1, 3,
446        boxbetts_lb, boxbetts_ub, boxbetts_xmin,
447        0.0, "Box and Betts exponential quadratic sum" },
448      { paviani_f, NULL, 1, 10,
449        paviani_lb, paviani_ub, paviani_xmin,
450        -45.7784697, "Paviani function" },
451      { grosenbrock_f, NULL, 1, 30,
452        grosenbrock_lb, grosenbrock_ub, grosenbrock_xmin,
453        0.0, "Generalized Rosenbrock function" },
454      { goldsteinprice_f, NULL, 1, 2,
455        goldsteinprice_lb, goldsteinprice_ub, goldsteinprice_xmin,
456        3.0, "Goldstein and Price function" },
457      { shekel_f, shekel_m + 0, 1, 4,
458        shekel_lb, shekel_ub, shekel0_xmin,
459        -10.15319968, "Shekel m=5 function" },
460      { shekel_f, shekel_m + 1, 1, 4,
461        shekel_lb, shekel_ub, shekel1_xmin,
462        -10.40294057, "Shekel m=7 function" },
463      { shekel_f, shekel_m + 2, 1, 4,
464        shekel_lb, shekel_ub, shekel2_xmin,
465        -10.53640982, "Shekel m=10 function" },
466      { levy_f, NULL, 1, 4,
467        levy4_lb, levy4_ub, levy4_xmin,
468        -21.50235596, "Levy n=4 function" },
469      { levy_f, NULL, 1, 5,
470        levy_lb, levy_ub, levy_xmin+2,
471        -11.50440302, "Levy n=5 function" },
472      { levy_f, NULL, 1, 6,
473        levy_lb, levy_ub, levy_xmin+1,
474        -11.50440302, "Levy n=6 function" },
475      { levy_f, NULL, 1, 7,
476        levy_lb, levy_ub, levy_xmin,
477        -11.50440302, "Levy n=7 function" },
478      { griewank_f, NULL, 1, 10,
479        griewank_lb, griewank_ub, griewank_xmin,
480        0.0, "Griewank function" },
481      { sixhumpcamel_f, NULL, 1, 2,
482        sixhumpcamel_lb, sixhumpcamel_ub, sixhumpcamel_xmin,
483        -1.031628453, "Six-hump camel back function" },
484      { convexcosh_f, NULL, 1, 10,
485        convexcosh_lb, convexcosh_ub, convexcosh_xmin,
486        1.0, "Convex product of cosh functions" },
487      { branin_f, NULL, 1, 2,
488        branin_lb, branin_ub, branin_xmin,
489        -.0, "Branin function" },
490      { shubert_f, NULL, 1, 2,
491        shubert_lb, shubert_ub, shubert_xmin,
492        -24.06249888, "Shubert function" },
493      { hansen_f, NULL, 1, 2,
494        hansen_lb, hansen_ub, hansen_xmin,
495        -176.5417931367, "Hansen function" },
496      { osc1d_f, NULL, 1, 1,
497        osc1d_lb, osc1d_ub, osc1d_xmin,
498        -1.0, "1d oscillating function with a single minimum" },
499      { corner4d_f, NULL, 1, 4,
500        corner4d_lb, corner4d_ub, corner4d_xmin,
501        1.0, "4d function with minimum at corner" },
502      { side4d_f, NULL, 1, 4,
503        side4d_lb, side4d_ub, side4d_xmin,
504        -141.285020472, "4d function with minimum at side" }
505 };