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