test/testopt.cpp

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