test/testopt.c

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