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