chiark / gitweb /
added -C option to force contraints to be active
[nlopt.git] / test / testopt.cpp
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <string.h>
5
6 #include "config.h"
7
8 #ifdef HAVE_UNISTD_H
9 #  include <unistd.h>
10 #endif
11 #ifdef HAVE_GETOPT_H
12 #  include <getopt.h>
13 #endif
14
15 #include "nlopt.h"
16 #include "nlopt-util.h"
17 #include "testfuncs.h"
18
19 static nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L;
20 static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, minf_max_delta = -HUGE_VAL;
21 static int maxeval = 1000, iterations = 1, center_start = 0;
22 static double maxtime = 0.0;
23 static double xinit_tol = -1;
24 static int force_constraints = 0;
25
26 static void listalgs(FILE *f)
27 {
28   int i;
29   fprintf(f, "Available algorithms:\n");
30   for (i = 0; i < NLOPT_NUM_ALGORITHMS; ++i)
31     fprintf(f, "  %2d: %s\n", i, nlopt_algorithm_name((nlopt_algorithm) i));
32 }
33
34 static void listfuncs(FILE *f)
35 {
36   int i;
37   fprintf(f, "Available objective functions:\n");
38   for (i = 0; i < NTESTFUNCS; ++i)
39     fprintf(f, "  %2d: %s (%d dims)\n", i, testfuncs[i].name, testfuncs[i].n);
40 }
41
42 static int test_function(int ifunc)
43 {
44   testfunc func;
45   int i, iter;
46   double *x, minf, minf_max, f0, *xtabs, *lb, *ub;
47   nlopt_result ret;
48   double start = nlopt_seconds();
49   
50   if (ifunc < 0 || ifunc >= NTESTFUNCS) {
51     fprintf(stderr, "testopt: invalid function %d\n", ifunc);
52     listfuncs(stderr);
53     return 0;
54   }
55   func = testfuncs[ifunc];
56   x = (double *) malloc(sizeof(double) * func.n * 5);
57   if (!x) { fprintf(stderr, "testopt: Out of memory!\n"); return 0; }
58
59   lb = x + func.n * 3;
60   ub = lb + func.n;
61   xtabs = x + func.n * 2;
62
63   for (i = 0; i < func.n; ++i) xtabs[i] = xtol_abs;
64   minf_max = minf_max_delta > (-HUGE_VAL) ? minf_max_delta + func.minf : (-HUGE_VAL);
65   
66   printf("-----------------------------------------------------------\n");
67   printf("Optimizing %s (%d dims) using %s algorithm\n",
68          func.name, func.n, nlopt_algorithm_name(algorithm));
69   printf("lower bounds at lb = [");
70   for (i = 0; i < func.n; ++i) printf(" %g", func.lb[i]);
71   printf("]\n");
72   printf("upper bounds at ub = [");
73   for (i = 0; i < func.n; ++i) printf(" %g", func.ub[i]);
74   printf("]\n");
75   memcpy(lb, func.lb, func.n * sizeof(double));
76   memcpy(ub, func.ub, func.n * sizeof(double));
77   if (force_constraints) {
78     for (i = 0; i < func.n; ++i) {
79       if (nlopt_iurand(2) == 0)
80         ub[i] = nlopt_urand(lb[i], func.xmin[i]);
81       else
82         lb[i] = nlopt_urand(func.xmin[i], ub[i]);
83     }
84     printf("adjusted lower bounds at lb = [");
85     for (i = 0; i < func.n; ++i) printf(" %g", lb[i]);
86     printf("]\n");
87     printf("adjusted upper bounds at ub = [");
88     for (i = 0; i < func.n; ++i) printf(" %g", ub[i]);
89     printf("]\n");
90   }
91
92
93   printf("Starting guess x = [");
94   for (i = 0; i < func.n; ++i) {
95     if (center_start)
96       x[i] = (ub[i] + lb[i]) * 0.5;
97     else if (xinit_tol < 0) { /* random starting point near center of box */
98       double dx = (ub[i] - lb[i]) * 0.25;
99       double xm = 0.5 * (ub[i] + lb[i]);
100       x[i] = nlopt_urand(xm - dx, xm + dx);
101     }
102     else {
103       x[i] = nlopt_urand(-xinit_tol, xinit_tol)
104         + (1 + nlopt_urand(-xinit_tol, xinit_tol)) * func.xmin[i];
105       if (x[i] > ub[i]) x[i] = ub[i];
106       else if (x[i] < lb[i]) x[i] = lb[i];
107     }
108     printf(" %g", x[i]);
109   }
110   printf("]\n");
111   f0 = func.f(func.n, x, x + func.n, func.f_data);
112   printf("Starting function value = %g\n", f0);
113
114   if (testfuncs_verbose && func.has_gradient) {
115     printf("checking gradient:\n");
116     for (i = 0; i < func.n; ++i) {
117       double f;
118       x[i] *= 1 + 1e-6;
119       f = func.f(func.n, x, NULL, func.f_data);
120       x[i] /= 1 + 1e-6;
121       printf("  grad[%d] = %g vs. numerical derivative %g\n",
122              i, x[i + func.n], (f - f0) / (x[i] * 1e-6));
123     }
124   }
125
126   for (iter = 0; iter < iterations; ++iter) {
127     testfuncs_counter = 0;
128     ret = nlopt_minimize(algorithm,
129                          func.n, func.f, func.f_data,
130                          lb, ub,
131                          x, &minf,
132                          minf_max, ftol_rel, ftol_abs, xtol_rel, xtabs,
133                          maxeval, maxtime);
134     printf("finished after %g seconds.\n", nlopt_seconds() - start);
135     printf("return code %d from nlopt_minimize\n", ret);
136     if (ret < 0) {
137       fprintf(stderr, "testopt: error in nlopt_minimize\n");
138       return 0;
139     }
140     printf("Found minimum f = %g after %d evaluations.\n", 
141            minf, testfuncs_counter);
142     printf("Minimum at x = [");
143     for (i = 0; i < func.n; ++i) printf(" %g", x[i]);
144     printf("]\n");
145     printf("|f - minf| = %g, |f - minf| / |minf| = %e\n",
146            fabs(minf - func.minf), fabs(minf - func.minf) / fabs(func.minf));
147   }
148   printf("vs. global minimum f = %g at x = [", func.minf);
149   for (i = 0; i < func.n; ++i) printf(" %g", func.xmin[i]);
150   printf("]\n");
151   
152   free(x);
153   return 1;
154 }
155
156 static void usage(FILE *f)
157 {
158   fprintf(f, "Usage: testopt [OPTIONS]\n"
159           "Options:\n"
160           "     -h : print this help\n"
161           "     -L : list available algorithms and objective functions\n"
162           "     -v : verbose mode\n"
163           " -a <n> : use optimization algorithm <n>\n"
164           " -o <n> : use objective function <n>\n"
165           " -0 <x> : starting guess within <x> + (1+<x>) * optimum\n"
166           "     -c : starting guess at center of cell\n"
167           "     -C : put optimum outside of bound constraints\n"
168           " -e <n> : use at most <n> evals (default: %d, 0 to disable)\n"
169           " -t <t> : use at most <t> seconds (default: disabled)\n"
170           " -x <t> : relative tolerance <t> on x (default: disabled)\n"
171           " -X <t> : absolute tolerance <t> on x (default: disabled)\n"
172           " -f <t> : relative tolerance <t> on f (default: disabled)\n"
173           " -F <t> : absolute tolerance <t> on f (default: disabled)\n"
174           " -m <m> : stop when minf+<m> is reached (default: disabled)\n"
175           " -i <n> : iterate optimization <n> times (default: 1)\n"
176           " -r <s> : use random seed <s> for starting guesses\n"
177           , maxeval);
178 }
179
180 int main(int argc, char **argv)
181 {
182   int c;
183   
184   nlopt_srand_time();
185   testfuncs_verbose = 0;
186   
187   if (argc <= 1)
188     usage(stdout);
189   
190   while ((c = getopt(argc, argv, "hLvCc0:r:a:o:i:e:t:x:X:f:F:m:")) != -1)
191     switch (c) {
192     case 'h':
193       usage(stdout);
194       return EXIT_SUCCESS;
195     case 'L':
196       listalgs(stdout);
197       listfuncs(stdout);
198       return EXIT_SUCCESS;
199     case 'v':
200       testfuncs_verbose = 1;
201       break;
202     case 'C':
203       force_constraints = 1;
204       break;
205     case 'r':
206       nlopt_srand((unsigned long) atoi(optarg));
207       break;
208     case 'a':
209       c = atoi(optarg);
210       if (c < 0 || c >= NLOPT_NUM_ALGORITHMS) {
211         fprintf(stderr, "testopt: invalid algorithm %d\n", c);
212         listalgs(stderr);
213         return EXIT_FAILURE;
214       }
215       algorithm = (nlopt_algorithm) c;
216       break;
217     case 'o':
218       if (!test_function(atoi(optarg)))
219         return EXIT_FAILURE;
220       break;
221     case 'e':
222       maxeval = atoi(optarg);
223       break;
224     case 'i':
225       iterations = atoi(optarg);
226       break;
227     case 't':
228       maxtime = atof(optarg);
229       break;
230     case 'x':
231       xtol_rel = atof(optarg);
232       break;
233     case 'X':
234       xtol_abs = atof(optarg);
235       break;
236     case 'f':
237       ftol_rel = atof(optarg);
238       break;
239     case 'F':
240       ftol_abs = atof(optarg);
241       break;
242     case 'm':
243       minf_max_delta = atof(optarg);
244       break;
245     case 'c':
246       center_start = 1;
247       break;
248     case '0':
249       center_start = 0;
250       xinit_tol = atof(optarg);
251       break;
252     default:
253       fprintf(stderr, "harminv: invalid argument -%c\n", c);
254       usage(stderr);
255       return EXIT_FAILURE;
256     }
257   
258   return EXIT_SUCCESS;
259 }