1 /* Copyright (c) 2007-2011 Massachusetts Institute of Technology
2 *
3 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including
6 * without limitation the rights to use, copy, modify, merge, publish,
7 * distribute, sublicense, and/or sell copies of the Software, and to
8 * permit persons to whom the Software is furnished to do so, subject to
9 * the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21 */
22
23 #include <stdlib.h>
24 #include <stdio.h>
25 #include <math.h>
26 #include <string.h>
27
28 #include "nlopt_config.h"
29
30 #ifdef HAVE_UNISTD_H
31 # include <unistd.h>
32 #endif
33 #if defined(HAVE_GETOPT_H) && defined(HAVE_GETOPT)
34 # include <getopt.h>
35 #else
36 # include "nlopt-getopt.h"
37 #endif
38
39 #define USE_FEENABLEEXCEPT 0
40 #if USE_FEENABLEEXCEPT
41 # include <fenv.h>
42 extern "C" int feenableexcept(int EXCEPTS);
43 #endif
44
45
46 #include "nlopt.h"
47 #include "nlopt-util.h"
48 #include "testfuncs.h"
49
50 static nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L;
51 static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, minf_max_delta;
52 static int maxeval = 1000, iterations = 1, center_start = 0;
53 static double maxtime = 0.0;
54 static double xinit_tol = -1;
55 static int force_constraints = 0;
56 static int fix_bounds[100] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
57 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
58 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
59 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
60 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
61 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
62 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
63 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
64 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
65 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
66 };
67
listalgs(FILE * f)68 static void listalgs(FILE * f)
69 {
70 int i;
71 fprintf(f, "Available algorithms:\n");
72 for (i = 0; i < NLOPT_NUM_ALGORITHMS; ++i)
73 fprintf(f, " %2d: %s\n", i, nlopt_algorithm_name((nlopt_algorithm) i));
74 }
75
listfuncs(FILE * f)76 static void listfuncs(FILE * f)
77 {
78 int i;
79 fprintf(f, "Available objective functions:\n");
80 for (i = 0; i < NTESTFUNCS; ++i)
81 fprintf(f, " %2d: %s (%d dims)\n", i, testfuncs[i].name, testfuncs[i].n);
82 }
83
84 typedef struct {
85 const double *lb, *ub;
86 nlopt_func f;
87 void *f_data;
88 } bounds_wrap_data;
89
bounds_wrap_func(unsigned n,const double * x,double * grad,void * d_)90 static double bounds_wrap_func(unsigned n, const double *x, double *grad, void *d_)
91 {
92 bounds_wrap_data *d = (bounds_wrap_data *) d_;
93 unsigned i;
94 double b = 0;
95 for (i = 0; i < n; ++i) {
96 if (x[i] < d->lb[i]) {
97 b = d->lb[i];
98 break;
99 } else if (x[i] > d->ub[i]) {
100 b = d->ub[i];
101 break;
102 }
103 }
104 if (i < n)
105 fprintf(stderr, "WARNING: bounds violated by x[%u] = %g = %g + %g\n", i, x[i], b, x[i] - b);
106 return d->f(n, x, grad, d->f_data);
107 }
108
test_function(int ifunc)109 static int test_function(int ifunc)
110 {
111 nlopt_opt opt;
112 testfunc func;
113 int i, iter;
114 double *x, minf, minf_max, f0, *xtabs, *lb, *ub;
115 nlopt_result ret;
116 double start = nlopt_seconds();
117 int total_count = 0, max_count = 0, min_count = 1 << 30;
118 double total_err = 0, max_err = 0;
119 bounds_wrap_data bw;
120
121 if (ifunc < 0 || ifunc >= NTESTFUNCS) {
122 fprintf(stderr, "testopt: invalid function %d\n", ifunc);
123 listfuncs(stderr);
124 return 0;
125 }
126 func = testfuncs[ifunc];
127 x = (double *) malloc(sizeof(double) * func.n * 5);
128 if (!x) {
129 fprintf(stderr, "testopt: Out of memory!\n");
130 return 0;
131 }
132
133 lb = x + func.n * 3;
134 ub = lb + func.n;
135 xtabs = x + func.n * 2;
136 bw.lb = lb;
137 bw.ub = ub;
138 bw.f = func.f;
139 bw.f_data = func.f_data;
140
141 for (i = 0; i < func.n; ++i)
142 xtabs[i] = xtol_abs;
143 minf_max = minf_max_delta > (-HUGE_VAL) ? minf_max_delta + func.minf : (-HUGE_VAL);
144
145 printf("-----------------------------------------------------------\n");
146 printf("Optimizing %s (%d dims) using %s algorithm\n", func.name, func.n, nlopt_algorithm_name(algorithm));
147 printf("lower bounds at lb = [");
148 for (i = 0; i < func.n; ++i)
149 printf(" %g", func.lb[i]);
150 printf("]\n");
151 printf("upper bounds at ub = [");
152 for (i = 0; i < func.n; ++i)
153 printf(" %g", func.ub[i]);
154 printf("]\n");
155 memcpy(lb, func.lb, func.n * sizeof(double));
156 memcpy(ub, func.ub, func.n * sizeof(double));
157 for (i = 0; i < func.n; ++i)
158 if (fix_bounds[i]) {
159 printf("fixing bounds for dim[%d] to xmin[%d]=%g\n", i, i, func.xmin[i]);
160 lb[i] = ub[i] = func.xmin[i];
161 }
162 if (force_constraints) {
163 for (i = 0; i < func.n; ++i) {
164 if (nlopt_iurand(2) == 0)
165 ub[i] = nlopt_urand(lb[i], func.xmin[i]);
166 else
167 lb[i] = nlopt_urand(func.xmin[i], ub[i]);
168 }
169 printf("adjusted lower bounds at lb = [");
170 for (i = 0; i < func.n; ++i)
171 printf(" %g", lb[i]);
172 printf("]\n");
173 printf("adjusted upper bounds at ub = [");
174 for (i = 0; i < func.n; ++i)
175 printf(" %g", ub[i]);
176 printf("]\n");
177 }
178
179 if (fabs(func.f(func.n, func.xmin, 0, func.f_data) - func.minf) > 1e-8) {
180 fprintf(stderr, "BUG: function does not achieve given lower bound!\n");
181 fprintf(stderr, "f(%g", func.xmin[0]);
182 for (i = 1; i < func.n; ++i)
183 fprintf(stderr, ", %g", func.xmin[i]);
184 fprintf(stderr, ") = %0.16g instead of %0.16g, |diff| = %g\n", func.f(func.n, func.xmin, 0, func.f_data), func.minf, fabs(func.f(func.n, func.xmin, 0, func.f_data) - func.minf));
185 free(x);
186 return 0;
187 }
188
189 for (iter = 0; iter < iterations; ++iter) {
190 double val;
191 testfuncs_counter = 0;
192
193 printf("Starting guess x = [");
194 for (i = 0; i < func.n; ++i) {
195 if (center_start)
196 x[i] = (ub[i] + lb[i]) * 0.5;
197 else if (xinit_tol < 0) { /* random starting point near center of box */
198 double dx = (ub[i] - lb[i]) * 0.25;
199 double xm = 0.5 * (ub[i] + lb[i]);
200 x[i] = nlopt_urand(xm - dx, xm + dx);
201 } else {
202 x[i] = nlopt_urand(-xinit_tol, xinit_tol)
203 + (1 + nlopt_urand(-xinit_tol, xinit_tol)) * func.xmin[i];
204 if (x[i] > ub[i])
205 x[i] = ub[i];
206 else if (x[i] < lb[i])
207 x[i] = lb[i];
208 }
209 printf(" %g", x[i]);
210 }
211 printf("]\n");
212 f0 = func.f(func.n, x, x + func.n, func.f_data);
213 printf("Starting function value = %g\n", f0);
214
215 if (iter == 0 && testfuncs_verbose && func.has_gradient) {
216 printf("checking gradient:\n");
217 for (i = 0; i < func.n; ++i) {
218 double f;
219 x[i] *= 1 + 1e-6;
220 f = func.f(func.n, x, NULL, func.f_data);
221 x[i] /= 1 + 1e-6;
222 printf(" grad[%d] = %g vs. numerical derivative %g\n", i, x[i + func.n], (f - f0) / (x[i] * 1e-6));
223 }
224 }
225
226 testfuncs_counter = 0;
227 opt = nlopt_create(algorithm, func.n);
228 nlopt_set_min_objective(opt, bounds_wrap_func, &bw);
229 nlopt_set_lower_bounds(opt, lb);
230 nlopt_set_upper_bounds(opt, ub);
231 nlopt_set_stopval(opt, minf_max);
232 nlopt_set_ftol_rel(opt, ftol_rel);
233 nlopt_set_ftol_abs(opt, ftol_abs);
234 nlopt_set_xtol_rel(opt, xtol_rel);
235 nlopt_set_xtol_abs(opt, xtabs);
236 nlopt_set_maxeval(opt, maxeval);
237 nlopt_set_maxtime(opt, maxtime);
238 ret = nlopt_optimize(opt, x, &minf);
239 printf("finished after %g seconds.\n", nlopt_seconds() - start);
240 printf("return code %d from nlopt_minimize\n", ret);
241 if (ret < 0 && ret != NLOPT_ROUNDOFF_LIMITED && ret != NLOPT_FORCED_STOP) {
242 fprintf(stderr, "testopt: error in nlopt_minimize\n");
243 free(x);
244 return 0;
245 }
246 printf("Found minimum f = %g after %d evaluations (numevals = %d).\n", minf, testfuncs_counter, nlopt_get_numevals(opt));
247 nlopt_destroy(opt);
248 total_count += testfuncs_counter;
249 if (testfuncs_counter > max_count)
250 max_count = testfuncs_counter;
251 if (testfuncs_counter < min_count)
252 min_count = testfuncs_counter;
253 printf("Minimum at x = [");
254 for (i = 0; i < func.n; ++i)
255 printf(" %g", x[i]);
256 printf("]\n");
257 if (func.minf == 0)
258 printf("|f - minf| = %g\n", fabs(minf - func.minf));
259 else
260 printf("|f - minf| = %g, |f - minf| / |minf| = %e\n", fabs(minf - func.minf), fabs(minf - func.minf) / fabs(func.minf));
261 total_err += fabs(minf - func.minf);
262 if (fabs(minf - func.minf) > max_err)
263 max_err = fabs(minf - func.minf);
264 printf("vs. global minimum f = %g at x = [", func.minf);
265 for (i = 0; i < func.n; ++i)
266 printf(" %g", func.xmin[i]);
267 printf("]\n");
268
269 val = func.f(func.n, x, NULL, func.f_data);
270 if (fabs(val - minf) > 1e-12) {
271 fprintf(stderr, "Mismatch %g between returned minf=%g and f(x) = %g\n", minf - val, minf, val);
272 free(x);
273 return 0;
274 }
275 }
276 if (iterations > 1)
277 printf("average #evaluations = %g (%d-%d)\naverage |f-minf| = %g, max |f-minf| = %g\n", total_count * 1.0 / iterations, min_count, max_count, total_err / iterations, max_err);
278
279 free(x);
280 return 1;
281 }
282
usage(FILE * f)283 static void usage(FILE * f)
284 {
285 fprintf(f, "Usage: testopt [OPTIONS]\n"
286 "Options:\n"
287 " -h : print this help\n"
288 " -L : list available algorithms and objective functions\n"
289 " -v : verbose mode\n"
290 " -a <n> : use optimization algorithm <n>\n"
291 " -o <n> : use objective function <n>\n"
292 " -0 <x> : starting guess within <x> + (1+<x>) * optimum\n"
293 " -b <dim0,dim1,...>: eliminate given dims by equating bounds\n");
294 fprintf(f,
295 " -c : starting guess at center of cell\n"
296 " -C : put optimum outside of bound constraints\n"
297 " -e <n> : use at most <n> evals (default: %d, 0 to disable)\n"
298 " -t <t> : use at most <t> seconds (default: disabled)\n"
299 " -x <t> : relative tolerance <t> on x (default: disabled)\n"
300 " -X <t> : absolute tolerance <t> on x (default: disabled)\n", maxeval);
301 fprintf(f,
302 " -f <t> : relative tolerance <t> on f (default: disabled)\n"
303 " -F <t> : absolute tolerance <t> on f (default: disabled)\n"
304 " -m <m> : stop when minf+<m> is reached (default: disabled)\n"
305 " -i <n> : iterate optimization <n> times (default: 1)\n"
306 " -r <s> : use random seed <s> for starting guesses\n");
307 }
308
main(int argc,char ** argv)309 int main(int argc, char **argv)
310 {
311 int c;
312
313 nlopt_srand_time();
314 testfuncs_verbose = 0;
315 minf_max_delta = -HUGE_VAL;
316
317 if (argc <= 1)
318 usage(stdout);
319
320 #if USE_FEENABLEEXCEPT
321 feenableexcept(FE_INVALID);
322 #endif
323
324 while ((c = getopt(argc, argv, "hLvVCc0:r:a:o:i:e:t:x:X:f:F:m:b:")) != -1)
325 switch (c) {
326 case 'h':
327 usage(stdout);
328 return EXIT_SUCCESS;
329 case 'L':
330 listalgs(stdout);
331 listfuncs(stdout);
332 return EXIT_SUCCESS;
333 case 'v':
334 testfuncs_verbose = 1;
335 break;
336 case 'V': {
337 int major, minor, patch;
338 nlopt_version(&major, &minor, &patch);
339 printf("NLopt version %d.%d.%d\n", major, minor, patch);
340 return EXIT_SUCCESS;
341 }
342 case 'C':
343 force_constraints = 1;
344 break;
345 case 'r':
346 nlopt_srand((unsigned long) atoi(optarg));
347 break;
348 case 'a':
349 c = atoi(optarg);
350 if (c < 0 || c >= NLOPT_NUM_ALGORITHMS) {
351 fprintf(stderr, "testopt: invalid algorithm %d\n", c);
352 listalgs(stderr);
353 return EXIT_FAILURE;
354 }
355 algorithm = (nlopt_algorithm) c;
356 break;
357 case 'o':
358 if (!test_function(atoi(optarg)))
359 return EXIT_FAILURE;
360 break;
361 case 'e':
362 maxeval = atoi(optarg);
363 break;
364 case 'i':
365 iterations = atoi(optarg);
366 break;
367 case 't':
368 maxtime = atof(optarg);
369 break;
370 case 'x':
371 xtol_rel = atof(optarg);
372 break;
373 case 'X':
374 xtol_abs = atof(optarg);
375 break;
376 case 'f':
377 ftol_rel = atof(optarg);
378 break;
379 case 'F':
380 ftol_abs = atof(optarg);
381 break;
382 case 'm':
383 minf_max_delta = atof(optarg);
384 break;
385 case 'c':
386 center_start = 1;
387 break;
388 case '0':
389 center_start = 0;
390 xinit_tol = atof(optarg);
391 break;
392 case 'b':{
393 const char *s = optarg;
394 while (s && *s) {
395 int b = atoi(s);
396 if (b < 0 || b >= 100) {
397 fprintf(stderr, "invalid -b argument");
398 return EXIT_FAILURE;
399 }
400 fix_bounds[b] = 1;
401 s = strchr(s, ',');
402 if (s)
403 ++s;
404 }
405 break;
406 }
407 default:
408 fprintf(stderr, "harminv: invalid argument -%c\n", c);
409 usage(stderr);
410 return EXIT_FAILURE;
411 }
412
413 return EXIT_SUCCESS;
414 }
415