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