test/testopt.cpp

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