From 59b1360e2f8d746fd65c8da9aabacaa150290e7c Mon Sep 17 00:00:00 2001 From: stevenj Date: Thu, 23 Aug 2007 16:17:09 -0400 Subject: [PATCH] added test program and some test objectives darcs-hash:20070823201709-c8de0-e8894aceee594f8bd61b5d13d85c7ae5d0692dfa.gz --- Makefile.am | 2 +- api/nlopt.c | 17 ++++- api/nlopt.h | 16 +++-- configure.ac | 2 + stogo/global.cc | 43 ++++++----- stogo/global.h | 7 +- stogo/stogo.cc | 2 +- stogo/stogo.h | 4 +- test/Makefile.am | 6 ++ test/testfuncs.c | 182 +++++++++++++++++++++++++++++++++++++++++++++++ test/testfuncs.h | 31 ++++++++ test/testopt.cpp | 158 ++++++++++++++++++++++++++++++++++++++++ 12 files changed, 441 insertions(+), 29 deletions(-) create mode 100644 test/Makefile.am create mode 100644 test/testfuncs.c create mode 100644 test/testfuncs.h create mode 100644 test/testopt.cpp diff --git a/Makefile.am b/Makefile.am index 817c623..0fa9377 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,7 +3,7 @@ lib_LTLIBRARIES = libnlopt.la ACLOCAL_AMFLAGS=-I ./m4 -SUBDIRS=lbfgs subplex direct stogo api +SUBDIRS=lbfgs subplex direct stogo api . test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 libnlopt_la_SOURCES = diff --git a/api/nlopt.c b/api/nlopt.c index e63e329..ef3acf9 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -4,6 +4,19 @@ #include "nlopt.h" #include "config.h" +static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { + "DIRECT (global)", + "Subplex (local)", + "StoGO (global)", + "Low-storage BFGS (LBFGS) (local)" +}; + +const char *nlopt_algorithm_name(nlopt_algorithm a) +{ + if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN"; + return nlopt_algorithm_names[a]; +} + static int my_isinf(double x) { return x == HUGE_VAL #ifdef HAVE_ISINF @@ -46,7 +59,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) #define MIN(a,b) ((a) < (b) ? (a) : (b)) nlopt_result nlopt_minimize( - nlopt_method method, + nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ @@ -59,7 +72,7 @@ nlopt_result nlopt_minimize( d.f = f; d.f_data = f_data; - switch (method) { + switch (algorithm) { case NLOPT_GLOBAL_DIRECT: switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin, maxeval, 500, ftol_rel, ftol_abs, diff --git a/api/nlopt.h b/api/nlopt.h index f8a8c30..60c21a4 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -11,14 +11,18 @@ typedef double (*nlopt_func)(int n, const double *x, void *func_data); typedef enum { - /* non-gradient methods */ - NLOPT_GLOBAL_DIRECT, + /* non-gradient algorithms */ + NLOPT_GLOBAL_DIRECT = 0, NLOPT_LOCAL_SUBPLEX, - /* gradient-based methods */ + /* gradient-based algorithms */ NLOPT_GLOBAL_STOGO, - NLOPT_LOCAL_LBFGS -} nlopt_method; + NLOPT_LOCAL_LBFGS, + + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ +} nlopt_algorithm; + +extern const char *nlopt_algorithm_name(nlopt_algorithm a); typedef enum { NLOPT_FAILURE = -1, /* generic failure code */ @@ -34,7 +38,7 @@ typedef enum { } nlopt_result; extern nlopt_result nlopt_minimize( - nlopt_method method, + nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ diff --git a/configure.ac b/configure.ac index 14580de..92fcebe 100644 --- a/configure.ac +++ b/configure.ac @@ -22,6 +22,7 @@ AC_PROG_LIBTOOL dnl Checks for typedefs, structures, and compiler characteristics. AC_HEADER_STDC AC_HEADER_TIME +AC_CHECK_HEADERS([unistd.h getopt.h]) AC_C_CONST AC_C_INLINE @@ -82,6 +83,7 @@ AC_CONFIG_FILES([ stogo/Makefile lbfgs/Makefile subplex/Makefile + test/Makefile ]) AC_OUTPUT diff --git a/stogo/global.cc b/stogo/global.cc index 9496358..80decbe 100644 --- a/stogo/global.cc +++ b/stogo/global.cc @@ -26,6 +26,8 @@ time_t StartTime; double MacEpsilon ; int FC=0, GC=0 ; +int stogo_verbose = 0; /* set to nonzero for verbose output */ + Global::Global(RTBox D, Pobj o, Pgrad g, GlobalParams P): Domain(D) { dim=Domain.GetDim(); @@ -121,8 +123,10 @@ double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) { box.AddTrial(tmpTrial) ; if (tmpTrial.objval<=fbound+mu && tmpTrial.objval<=box.fmin+mu) { - cout << "Found a candidate, x=" << tmpTrial.xvals; - cout << " F=" < +#include + +#include "testfuncs.h" +#include "config.h" + +#define UNUSED(x) (void) x + +static double sqr(double x) { return x * x; } + +int testfuncs_verbose = 0; +int testfuncs_counter = 0; + +static double testfuncs_status(int n, const double *x, double f) +{ + ++testfuncs_counter; + if (testfuncs_verbose) { + int i; + printf("f_%d (%g", testfuncs_counter, x[0]); + for (i = 1; i < n; ++i) printf(", %g", x[i]); + printf(") = %g\n", f); + } + return f; +} + +#define RETURN(f) return testfuncs_status(n, x, f); + +/****************************************************************************/ +static double rosenbrock_f(int n, const double *x, double *grad, void *data) +{ + double a = x[1] - x[0] * x[0], b = 1 - x[0]; + UNUSED(data); + if (grad) { + grad[0] = -400 * a * x[0] - 2*b; + grad[1] = -200 * a; + } + RETURN(100 * sqr(a) + sqr(b)); +} + +static const double rosenbrock_lb[2] = {-2, -2}; +static const double rosenbrock_ub[2] = {2, 2}; +static const double rosenbrock_xmin[2] = {1, 1}; + +/****************************************************************************/ +static double mccormic_f(int n, const double *x, double *grad, void *data) +{ + double a = x[0] + x[1], b = x[0] - x[1]; + UNUSED(data); + if (grad) { + grad[0] = cos(a) + 2*b - 1.5; + grad[1] = cos(a) - 2*b + 2.5; + } + RETURN(sin(a) + sqr(b) - 1.5*x[0] + 2.5*x[1] + 1); +} + +static const double mccormic_lb[2] = {-1.5, -3}; +static const double mccormic_ub[2] = {4, 4}; +static const double mccormic_xmin[2] = {-0.54719, 1.54719}; + +/****************************************************************************/ +static double boxbetts_f(int n, const double *x, double *grad, void *data) +{ + int i; + double f = 0; + UNUSED(data); + if (grad) + grad[0] = grad[1] = grad[2] = 0; + for (i = 1; i <= 10; ++i) { + double e0 = exp(-0.1*i*x[0]); + double e1 = exp(-0.1*i*x[1]); + double e2 = exp(-0.1*i) - exp((double) -i); + double g = e0 - e1 - e2 * x[2]; + f += sqr(g); + if (grad) { + grad[0] += (2 * g) * (-0.1*i*e0); + grad[1] += (2 * g) * (0.1*i*e1); + grad[2] += -(2 * g) * e2; + } + } + RETURN(f); +} + +static const double boxbetts_lb[3] = {0.9,9,0.9}; +static const double boxbetts_ub[3] = {1.2,11.2,1.2}; +static const double boxbetts_xmin[3] = {1,10,1}; + +/****************************************************************************/ +static double paviani_f(int n, const double *x, double *grad, void *data) +{ + int i; + double f = 0, prod = 1; + UNUSED(data); + if (grad) for (i = 0; i < 10; ++i) grad[i] = 0; + for (i = 0; i < 10; ++i) { + double ln1 = log(x[i] - 2); + double ln2 = log(10 - x[i]); + f += sqr(ln1) + sqr(ln2); + if (grad) + grad[i] += 2 * ln1 / (x[i] - 2) - 2 * ln2 / (10 - x[i]); + prod *= x[i]; + } + f -= (prod = pow(prod, 0.2)); + if (grad) + for (i = 0; i < 10; ++i) + grad[i] -= 0.2 * prod / x[i]; + RETURN(f); +} + +static const double paviani_lb[10] = {2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001,2.001}; +static const double paviani_ub[10] = {9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999,9.999}; +static const double paviani_xmin[10] = {9.350266,9.350266,9.350266,9.350266,9.350266,9.350266,9.350266,9.350266,9.350266,9.350266}; + +/****************************************************************************/ +static double grosenbrock_f(int n, const double *x, double *grad, void *data) +{ + int i; + double f = 0; + UNUSED(data); + if (grad) grad[0] = 0; + for (i = 0; i < 29; ++i) { + double a = x[i+1] - x[i] * x[i], b = 1 - x[i]; + if (grad) { + grad[i] += -400 * a * x[i] - 2*b; + grad[i+1] = -200 * a; + } + f += 100 * sqr(a) + sqr(b); + } + RETURN(f); +} + +static const double grosenbrock_lb[30] = {-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30}; +static const double grosenbrock_ub[30] = {30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30}; +static const double grosenbrock_xmin[30] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; + +/****************************************************************************/ +static double goldsteinprice_f(int n, const double *x, double *grad, void *data) +{ + double x0, x1, a1, a12, a2, b1, b12, b2; + UNUSED(data); + x0 = x[0]; x1 = x[1]; + a1 = x0+x1+1; a12 = sqr(a1); + a2 = 19 - 14*x0 + 3*x0*x0 - 14*x1 + 6*x0*x1 + 3*x1*x1; + b1 = 2*x0-3*x1; b12 = sqr(b1); + b2 = 18 - 32*x0 + 12*x0*x0 + 48*x1 - 36*x0*x1 + 27*x1*x1; + if (grad) { + grad[0] = (1 + a12 * a2) * (2 * b1 * 2 * b2 + + b12 * (-32 + 24*x0 - 36*x1)) + + (2 * a1 * a2 + a12 * (-14 + 6*x0 + 6*x1)) * (30 + b12 * b2); + grad[1] = (1 + a12 * a2) * (2 * b1 * (-3) * b2 + + b12 * (48 - 36*x0 + 54 * x1)) + + (2 * a1 * a2 + a12 * (-14 + 6*x0 + 6*x1)) * (30 + b12 * b2); + } + RETURN((1 + a12 * a2) * (30 + b12 * b2)); +} + +static const double goldsteinprice_lb[2] = {-2, -2}; +static const double goldsteinprice_ub[2] = {2, 2}; +static const double goldsteinprice_xmin[2] = {0, -1}; + +/****************************************************************************/ +/****************************************************************************/ + +const testfunc testfuncs[NTESTFUNCS] = { + { rosenbrock_f, NULL, 1, 2, + rosenbrock_lb, rosenbrock_ub, rosenbrock_xmin, + 0.0, "Rosenbrock function" }, + { mccormic_f, NULL, 1, 2, + mccormic_lb, mccormic_ub, mccormic_xmin, + -1.9133, "McCormic function" }, + { boxbetts_f, NULL, 1, 3, + boxbetts_lb, boxbetts_ub, boxbetts_xmin, + 0.0, "Box and Betts exponential quadratic sum" }, + { paviani_f, NULL, 1, 10, + paviani_lb, paviani_ub, paviani_xmin, + -45.778470, "Paviani function" }, + { grosenbrock_f, NULL, 1, 30, + grosenbrock_lb, grosenbrock_ub, grosenbrock_xmin, + 0.0, "Generalized Rosenbrock function" }, + { goldsteinprice_f, NULL, 1, 2, + goldsteinprice_lb, goldsteinprice_ub, goldsteinprice_xmin, + 3.0, "Goldstein and Price function" } +}; diff --git a/test/testfuncs.h b/test/testfuncs.h new file mode 100644 index 0000000..df81446 --- /dev/null +++ b/test/testfuncs.h @@ -0,0 +1,31 @@ +#ifndef TESTFUNCS_H +#define TESTFUNCS_H + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +#include "nlopt.h" + +typedef struct { + nlopt_func f; + void *f_data; + int has_gradient; + int n; + const double *lb, *ub, *xmin; + double fmin; + const char *name; +} testfunc; + +#define NTESTFUNCS 6 +extern const testfunc testfuncs[NTESTFUNCS]; + +extern int testfuncs_verbose; +extern int testfuncs_counter; + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif diff --git a/test/testopt.cpp b/test/testopt.cpp new file mode 100644 index 0000000..e2921ba --- /dev/null +++ b/test/testopt.cpp @@ -0,0 +1,158 @@ +#include +#include +#include + +#include "config.h" + +#ifdef HAVE_UNISTD_H +# include +#endif +#ifdef HAVE_GETOPT_H +# include +#endif +#if TIME_WITH_SYS_TIME +# include +# include +#else +# if HAVE_SYS_TIME_H +# include +# else +# include +# endif +#endif + +#include "nlopt.h" +#include "testfuncs.h" + +static double urand(double a, double b) +{ + return a + (rand() * (b - a) / RAND_MAX); +} + +static nlopt_algorithm algorithm = NLOPT_GLOBAL_DIRECT; +static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0; +static int maxeval = 1000; +static double maxtime = 0.0; + +static int test_function(int ifunc) +{ + testfunc func; + int i; + double *x, fmin, f0; + + if (ifunc < 0 || ifunc >= NTESTFUNCS) { + fprintf(stderr, "testopt: invalid function %d\n", ifunc); + return 0; + } + func = testfuncs[ifunc]; + if (!func.has_gradient && algorithm >= NLOPT_GLOBAL_STOGO) { + fprintf(stderr, + "testopt: A function with gradients is required for %s\n", + nlopt_algorithm_name(algorithm)); + return 0; + } + x = (double *) malloc(sizeof(double) * func.n * 2); + if (!x) { fprintf(stderr, "testopt: Out of memory!\n"); return 0; } + + + printf("-----------------------------------------------------------\n"); + printf("Optimizing %s (%d dims) using %s algorithm\n", + func.name, func.n, nlopt_algorithm_name(algorithm)); + printf("Starting guess x = ["); + for (i = 0; i < func.n; ++i) + printf(" %g", x[i] = urand(func.lb[i], func.ub[i])); + printf("]\n"); + f0 = func.f(func.n, x, x + func.n, func.f_data); + printf("Starting function value = %g\n", f0); + + if (testfuncs_verbose && func.has_gradient) { + printf("checking gradient:\n"); + for (i = 0; i < func.n; ++i) { + double f; + x[i] *= 1 + 1e-6; + f = func.f(func.n, x, NULL, func.f_data); + x[i] /= 1 + 1e-6; + printf(" grad[%d] = %g vs. numerical derivative %g\n", + i, x[i + func.n], (f - f0) / (x[i] * 1e-6)); + } + } + + testfuncs_counter = 0; + if (nlopt_minimize(algorithm, + func.n, func.f, func.f_data, + func.lb, func.ub, + x, &fmin, + HUGE_VAL, ftol_rel, ftol_abs, xtol_rel, NULL, + maxeval, maxtime) < 0) { + fprintf(stderr, "testopt: error in nlopt_minimize\n"); + return 0; + } + printf("Found minimum f = %g after %d evaluations.\n", + fmin, testfuncs_counter); + printf("Minimum at x = ["); + for (i = 0; i < func.n; ++i) printf(" %g", x[i]); + printf("]\n"); + printf("vs. global minimum f = %g at x = [", func.fmin); + for (i = 0; i < func.n; ++i) printf(" %g", func.xmin[i]); + printf("]\n"); + printf("|f - fmin| = %g, |f - fmin| / |fmin| = %e\n", + fabs(fmin - func.fmin), fabs(fmin - func.fmin) / fabs(func.fmin)); + + free(x); + return 1; +} + +static void usage(FILE *f) +{ + fprintf(f, "Usage: testopt [OPTIONS]\n" + "Options:\n" + " -h : print this help\n" + " -v : verbose mode\n" + " -r : use random seed for starting guesses\n" + " -a : use optimization algorithm \n" + " -o : use objective function \n" + " -e : use at most evals (default: %d)\n", + maxeval); +} + +int main(int argc, char **argv) +{ + int c; + + srand((unsigned) time(NULL)); + testfuncs_verbose = 0; + + while ((c = getopt(argc, argv, "hvra:o:e:")) != -1) + switch (c) { + case 'h': + usage(stdout); + return EXIT_SUCCESS; + case 'v': + testfuncs_verbose = 1; + break; + case 'r': + srand((unsigned) atoi(optarg)); + break; + case 'a': + c = atoi(optarg); + if (c < 0 || c >= NLOPT_NUM_ALGORITHMS) { + fprintf(stderr, "testopt: invalid algorithm %d\n", c); + return EXIT_FAILURE; + } + algorithm = (nlopt_algorithm) c; + break; + case 'o': + if (!test_function(atoi(optarg))) + return EXIT_FAILURE; + break; + case 'e': + maxeval = atoi(optarg); + break; + default: + fprintf(stderr, "harminv: invalid argument -%c\n", c); + usage(stderr); + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} -- 2.30.2