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 "config.h"
29
30 #ifdef HAVE_UNISTD_H
31 # include <unistd.h>
32 #endif
33 #ifdef HAVE_GETOPT_H
34 # include <getopt.h>
35 #endif
36
37 #define USE_FEENABLEEXCEPT 0
38 #if USE_FEENABLEEXCEPT
39 # include <fenv.h>
40 extern "C" int feenableexcept (int EXCEPTS);
41 #endif
42
43
44 #include "nlopt.h"
45 #include "nlopt-util.h"
46 #include "testfuncs.h"
47
48 static nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L;
49 static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, minf_max_delta = -HUGE_VAL;
50 static int maxeval = 1000, iterations = 1, center_start = 0;
51 static double maxtime = 0.0;
52 static double xinit_tol = -1;
53 static int force_constraints = 0;
54 static int fix_bounds[100] = {0,0,0,0,0,0,0,0,0,0,
55 0,0,0,0,0,0,0,0,0,0,
56 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
listalgs(FILE * f)65 static void listalgs(FILE *f)
66 {
67 int i;
68 fprintf(f, "Available algorithms:\n");
69 for (i = 0; i < NLOPT_NUM_ALGORITHMS; ++i)
70 fprintf(f, " %2d: %s\n", i, nlopt_algorithm_name((nlopt_algorithm) i));
71 }
72
listfuncs(FILE * f)73 static void listfuncs(FILE *f)
74 {
75 int i;
76 fprintf(f, "Available objective functions:\n");
77 for (i = 0; i < NTESTFUNCS; ++i)
78 fprintf(f, " %2d: %s (%d dims)\n", i, testfuncs[i].name, testfuncs[i].n);
79 }
80
81 typedef struct {
82 const double *lb, *ub;
83 nlopt_func f;
84 void *f_data;
85 } bounds_wrap_data;
86
bounds_wrap_func(int n,const double * x,double * grad,void * d_)87 static double bounds_wrap_func(int n, const double *x, double *grad, void *d_)
88 {
89 bounds_wrap_data *d = (bounds_wrap_data *) d_;
90 int i;
91 double b = 0;
92 for (i = 0; i < n; ++i) {
93 if (x[i] < d->lb[i]) {
94 b = d->lb[i];
95 break;
96 }
97 else if (x[i] > d->ub[i]) {
98 b = d->ub[i];
99 break;
100 }
101 }
102 if (i < n)
103 fprintf(stderr, "WARNING: bounds violated by x[%d] = %g = %g + %g\n",
104 i, x[i], b, x[i] - b);
105 return d->f(n, x, grad, d->f_data);
106 }
107
test_function(int ifunc)108 static int test_function(int ifunc)
109 {
110 testfunc func;
111 int i, iter;
112 double *x, minf, minf_max, f0, *xtabs, *lb, *ub;
113 nlopt_result ret;
114 double start = nlopt_seconds();
115 int total_count = 0, max_count = 0, min_count = 1<<30;
116 double total_err = 0, max_err = 0;
117 bounds_wrap_data bw;
118
119 if (ifunc < 0 || ifunc >= NTESTFUNCS) {
120 fprintf(stderr, "testopt: invalid function %d\n", ifunc);
121 listfuncs(stderr);
122 return 0;
123 }
124 func = testfuncs[ifunc];
125 x = (double *) malloc(sizeof(double) * func.n * 5);
126 if (!x) { fprintf(stderr, "testopt: Out of memory!\n"); return 0; }
127
128 lb = x + func.n * 3;
129 ub = lb + func.n;
130 xtabs = x + func.n * 2;
131 bw.lb = lb;
132 bw.ub = ub;
133 bw.f = func.f;
134 bw.f_data = func.f_data;
135
136 for (i = 0; i < func.n; ++i) xtabs[i] = xtol_abs;
137 minf_max = minf_max_delta > (-HUGE_VAL) ? minf_max_delta + func.minf : (-HUGE_VAL);
138
139 printf("-----------------------------------------------------------\n");
140 printf("Optimizing %s (%d dims) using %s algorithm\n",
141 func.name, func.n, nlopt_algorithm_name(algorithm));
142 printf("lower bounds at lb = [");
143 for (i = 0; i < func.n; ++i) printf(" %g", func.lb[i]);
144 printf("]\n");
145 printf("upper bounds at ub = [");
146 for (i = 0; i < func.n; ++i) printf(" %g", func.ub[i]);
147 printf("]\n");
148 memcpy(lb, func.lb, func.n * sizeof(double));
149 memcpy(ub, func.ub, func.n * sizeof(double));
150 for (i = 0; i < func.n; ++i) if (fix_bounds[i]) {
151 printf("fixing bounds for dim[%d] to xmin[%d]=%g\n",
152 i, i, func.xmin[i]);
153 lb[i] = ub[i] = func.xmin[i];
154 }
155 if (force_constraints) {
156 for (i = 0; i < func.n; ++i) {
157 if (nlopt_iurand(2) == 0)
158 ub[i] = nlopt_urand(lb[i], func.xmin[i]);
159 else
160 lb[i] = nlopt_urand(func.xmin[i], ub[i]);
161 }
162 printf("adjusted lower bounds at lb = [");
163 for (i = 0; i < func.n; ++i) printf(" %g", lb[i]);
164 printf("]\n");
165 printf("adjusted upper bounds at ub = [");
166 for (i = 0; i < func.n; ++i) printf(" %g", ub[i]);
167 printf("]\n");
168 }
169
170 if (fabs(func.f(func.n, func.xmin, 0, func.f_data) - func.minf) > 1e-8) {
171 fprintf(stderr, "BUG: function does not achieve given lower bound!\n");
172 fprintf(stderr, "f(%g", func.xmin[0]);
173 for (i = 1; i < func.n; ++i) fprintf(stderr, ", %g", func.xmin[i]);
174 fprintf(stderr, ") = %0.16g instead of %0.16g, |diff| = %g\n",
175 func.f(func.n, func.xmin, 0, func.f_data), func.minf,
176 fabs(func.f(func.n, func.xmin, 0, func.f_data) - func.minf));
177 return 0;
178 }
179
180 for (iter = 0; iter < iterations; ++iter) {
181 double val;
182 testfuncs_counter = 0;
183
184 printf("Starting guess x = [");
185 for (i = 0; i < func.n; ++i) {
186 if (center_start)
187 x[i] = (ub[i] + lb[i]) * 0.5;
188 else if (xinit_tol < 0) { /* random starting point near center of box */
189 double dx = (ub[i] - lb[i]) * 0.25;
190 double xm = 0.5 * (ub[i] + lb[i]);
191 x[i] = nlopt_urand(xm - dx, xm + dx);
192 }
193 else {
194 x[i] = nlopt_urand(-xinit_tol, xinit_tol)
195 + (1 + nlopt_urand(-xinit_tol, xinit_tol)) * func.xmin[i];
196 if (x[i] > ub[i]) x[i] = ub[i];
197 else if (x[i] < lb[i]) x[i] = lb[i];
198 }
199 printf(" %g", x[i]);
200 }
201 printf("]\n");
202 f0 = func.f(func.n, x, x + func.n, func.f_data);
203 printf("Starting function value = %g\n", f0);
204
205 if (iter == 0 && testfuncs_verbose && func.has_gradient) {
206 printf("checking gradient:\n");
207 for (i = 0; i < func.n; ++i) {
208 double f;
209 x[i] *= 1 + 1e-6;
210 f = func.f(func.n, x, NULL, func.f_data);
211 x[i] /= 1 + 1e-6;
212 printf(" grad[%d] = %g vs. numerical derivative %g\n",
213 i, x[i + func.n], (f - f0) / (x[i] * 1e-6));
214 }
215 }
216
217 testfuncs_counter = 0;
218 ret = nlopt_minimize(algorithm,
219 func.n, bounds_wrap_func, &bw,
220 lb, ub,
221 x, &minf,
222 minf_max, ftol_rel, ftol_abs, xtol_rel, xtabs,
223 maxeval, maxtime);
224 printf("finished after %g seconds.\n", nlopt_seconds() - start);
225 printf("return code %d from nlopt_minimize\n", ret);
226 if (ret < 0 && ret != NLOPT_ROUNDOFF_LIMITED
227 && ret != NLOPT_FORCED_STOP) {
228 fprintf(stderr, "testopt: error in nlopt_minimize\n");
229 free(x);
230 return 0;
231 }
232 printf("Found minimum f = %g after %d evaluations.\n",
233 minf, testfuncs_counter);
234 total_count += testfuncs_counter;
235 if (testfuncs_counter > max_count) max_count = testfuncs_counter;
236 if (testfuncs_counter < min_count) min_count = testfuncs_counter;
237 printf("Minimum at x = [");
238 for (i = 0; i < func.n; ++i) printf(" %g", x[i]);
239 printf("]\n");
240 if (func.minf == 0)
241 printf("|f - minf| = %g\n", fabs(minf - func.minf));
242 else
243 printf("|f - minf| = %g, |f - minf| / |minf| = %e\n",
244 fabs(minf - func.minf), fabs(minf - func.minf) / fabs(func.minf));
245 total_err += fabs(minf - func.minf);
246 if (fabs(minf - func.minf) > max_err)
247 max_err = fabs(minf - func.minf);
248 printf("vs. global minimum f = %g at x = [", func.minf);
249 for (i = 0; i < func.n; ++i) printf(" %g", func.xmin[i]);
250 printf("]\n");
251
252 val = func.f(func.n, x, NULL, func.f_data);
253 if (val != minf) {
254 fprintf(stderr, "Mismatch %g between returned minf=%g and f(x) = %g\n",
255 minf - val, minf, val);
256 free(x);
257 return 0;
258 }
259 }
260 if (iterations > 1)
261 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);
262
263 free(x);
264 return 1;
265 }
266
usage(FILE * f)267 static void usage(FILE *f)
268 {
269 fprintf(f, "Usage: testopt [OPTIONS]\n"
270 "Options:\n"
271 " -h : print this help\n"
272 " -L : list available algorithms and objective functions\n"
273 " -v : verbose mode\n"
274 " -a <n> : use optimization algorithm <n>\n"
275 " -o <n> : use objective function <n>\n"
276 " -0 <x> : starting guess within <x> + (1+<x>) * optimum\n"
277 " -b <dim0,dim1,...>: eliminate given dims by equating bounds\n"
278 " -c : starting guess at center of cell\n"
279 " -C : put optimum outside of bound constraints\n"
280 " -e <n> : use at most <n> evals (default: %d, 0 to disable)\n"
281 " -t <t> : use at most <t> seconds (default: disabled)\n"
282 " -x <t> : relative tolerance <t> on x (default: disabled)\n"
283 " -X <t> : absolute tolerance <t> on x (default: disabled)\n"
284 " -f <t> : relative tolerance <t> on f (default: disabled)\n"
285 " -F <t> : absolute tolerance <t> on f (default: disabled)\n"
286 " -m <m> : stop when minf+<m> is reached (default: disabled)\n"
287 " -i <n> : iterate optimization <n> times (default: 1)\n"
288 " -r <s> : use random seed <s> for starting guesses\n"
289 , maxeval);
290 }
291
main(int argc,char ** argv)292 int main(int argc, char **argv)
293 {
294 int c;
295
296 nlopt_srand_time();
297 testfuncs_verbose = 0;
298
299 if (argc <= 1)
300 usage(stdout);
301
302 #if USE_FEENABLEEXCEPT
303 feenableexcept(FE_INVALID);
304 #endif
305
306 while ((c = getopt(argc, argv, "hLvCc0:r:a:o:i:e:t:x:X:f:F:m:b:")) != -1)
307 switch (c) {
308 case 'h':
309 usage(stdout);
310 return EXIT_SUCCESS;
311 case 'L':
312 listalgs(stdout);
313 listfuncs(stdout);
314 return EXIT_SUCCESS;
315 case 'v':
316 testfuncs_verbose = 1;
317 break;
318 case 'C':
319 force_constraints = 1;
320 break;
321 case 'r':
322 nlopt_srand((unsigned long) atoi(optarg));
323 break;
324 case 'a':
325 c = atoi(optarg);
326 if (c < 0 || c >= NLOPT_NUM_ALGORITHMS) {
327 fprintf(stderr, "testopt: invalid algorithm %d\n", c);
328 listalgs(stderr);
329 return EXIT_FAILURE;
330 }
331 algorithm = (nlopt_algorithm) c;
332 break;
333 case 'o':
334 if (!test_function(atoi(optarg)))
335 return EXIT_FAILURE;
336 break;
337 case 'e':
338 maxeval = atoi(optarg);
339 break;
340 case 'i':
341 iterations = atoi(optarg);
342 break;
343 case 't':
344 maxtime = atof(optarg);
345 break;
346 case 'x':
347 xtol_rel = atof(optarg);
348 break;
349 case 'X':
350 xtol_abs = atof(optarg);
351 break;
352 case 'f':
353 ftol_rel = atof(optarg);
354 break;
355 case 'F':
356 ftol_abs = atof(optarg);
357 break;
358 case 'm':
359 minf_max_delta = atof(optarg);
360 break;
361 case 'c':
362 center_start = 1;
363 break;
364 case '0':
365 center_start = 0;
366 xinit_tol = atof(optarg);
367 break;
368 case 'b': {
369 const char *s = optarg;
370 while (s && *s) {
371 int b = atoi(s);
372 if (b < 0 || b >= 100) {
373 fprintf(stderr, "invalid -b argument");
374 return EXIT_FAILURE;
375 }
376 fix_bounds[b] = 1;
377 s = strchr(s, ','); if (s) ++s;
378 }
379 break;
380 }
381 default:
382 fprintf(stderr, "harminv: invalid argument -%c\n", c);
383 usage(stderr);
384 return EXIT_FAILURE;
385 }
386
387 return EXIT_SUCCESS;
388 }
389