From: stevenj Date: Fri, 15 Aug 2008 00:31:21 +0000 (-0400) Subject: initial (untested) stab at ~full MMA X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=b5092ed269150c48df1151159ad65f9de4cb3c8e;p=nlopt.git initial (untested) stab at ~full MMA darcs-hash:20080815003121-c8de0-638a9763e951bdab04146be55d5fc80125b8e5df.gz --- diff --git a/api/nlopt.c b/api/nlopt.c index 71570bd..fb58dd1 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -17,7 +17,7 @@ # define MY_INF HUGE_VAL #endif -static int my_isinf(double x) { +int nlopt_isinf(double x) { return fabs(x) >= HUGE_VAL * 0.99 #ifdef HAVE_ISINF || isinf(x) @@ -131,7 +131,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) nlopt_data *data = (nlopt_data *) data_; double f; f = data->f(n, x, NULL, data->f_data); - *undefined = isnan(f) || my_isinf(f); + *undefined = isnan(f) || nlopt_isinf(f); return f; } @@ -202,7 +202,12 @@ static nlopt_result nlopt_minimize_( nlopt_stopping stop; /* some basic argument checks */ - if (n <= 0 || !f || !lb || !ub || !x || !minf) + if (!minf || !f) return NLOPT_INVALID_ARGS; + if (n == 0) { /* trivial case: no degrees of freedom */ + *minf = f(n, x, NULL, f_data); + return NLOPT_SUCCESS; + } + else if (n < 0 || !lb || !ub || !x) return NLOPT_INVALID_ARGS; d.f = f; @@ -220,7 +225,8 @@ static nlopt_result nlopt_minimize_( return NLOPT_INVALID_ARGS; stop.n = n; - stop.minf_max = (isnan(minf_max) || (my_isinf(minf_max) && minf_max < 0)) + stop.minf_max = (isnan(minf_max) + || (nlopt_isinf(minf_max) && minf_max < 0)) ? -MY_INF : minf_max; stop.ftol_rel = ftol_rel; stop.ftol_abs = ftol_abs; @@ -297,11 +303,11 @@ static nlopt_result nlopt_minimize_( double *scale = (double *) malloc(sizeof(double) * n); if (!scale) return NLOPT_OUT_OF_MEMORY; for (i = 0; i < n; ++i) { - if (!my_isinf(ub[i]) && !my_isinf(lb[i])) + if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i])) scale[i] = (ub[i] - lb[i]) * 0.01; - else if (!my_isinf(lb[i]) && x[i] > lb[i]) + else if (!nlopt_isinf(lb[i]) && x[i] > lb[i]) scale[i] = (x[i] - lb[i]) * 0.01; - else if (!my_isinf(ub[i]) && x[i] < ub[i]) + else if (!nlopt_isinf(ub[i]) && x[i] < ub[i]) scale[i] = (ub[i] - x[i]) * 0.01; else scale[i] = 0.01 * x[i] + 0.0001; @@ -325,11 +331,11 @@ static nlopt_result nlopt_minimize_( case NLOPT_LN_PRAXIS: { double h0 = HUGE_VAL; for (i = 0; i < n; ++i) { - if (!my_isinf(ub[i]) && !my_isinf(lb[i])) + if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i])) h0 = MIN(h0, (ub[i] - lb[i]) * 0.01); - else if (!my_isinf(lb[i]) && x[i] > lb[i]) + else if (!nlopt_isinf(lb[i]) && x[i] > lb[i]) h0 = MIN(h0, (x[i] - lb[i]) * 0.01); - else if (!my_isinf(ub[i]) && x[i] < ub[i]) + else if (!nlopt_isinf(ub[i]) && x[i] < ub[i]) h0 = MIN(h0, (ub[i] - x[i]) * 0.01); else h0 = MIN(h0, 0.01 * x[i] + 0.0001); @@ -343,8 +349,8 @@ static nlopt_result nlopt_minimize_( int iret, *nbd = (int *) malloc(sizeof(int) * n); if (!nbd) return NLOPT_OUT_OF_MEMORY; for (i = 0; i < n; ++i) { - int linf = my_isinf(lb[i]) && lb[i] < 0; - int uinf = my_isinf(ub[i]) && ub[i] > 0; + int linf = nlopt_isinf(lb[i]) && lb[i] < 0; + int uinf = nlopt_isinf(ub[i]) && ub[i] > 0; nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2)); } iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub, diff --git a/mma/mma.c b/mma/mma.c index 891d83d..280ee3f 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -13,33 +13,188 @@ int mma_verbose = 0; /* > 0 for verbose output */ /* magic minimum value for rho in MMA ... the 2002 paper says it should be a "fixed, strictly positive `small' number, e.g. 1e-5" ... grrr, I hate these magic numbers, which seem like they - should depend on the objective function in some way */ + should depend on the objective function in some way ... in particular, + note that rho is dimensionful (= dimensions of objective function) */ #define MMA_RHOMIN 1e-5 +/***********************************************************************/ +/* function for MMA's dual optimization of the approximant function */ + +typedef struct { + int n; /* must be set on input to dimension of x */ + const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */ + const double *dfcdx; /* m-by-n array of fc gradients */ + double fval, rho; /* must be set on input */ + const double *fcval, *rhoc; /* arrays of length m */ + double *xcur; /* array of length n, output each time */ + double gval, wval, *gcval; /* output each time (array length m) */ +} dual_data; + +static double dual_func(int m, const double *y, double *grad, void *d_) +{ + dual_data *d = (dual_data *) d_; + int n = d->n; + const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma, + *dfdx = d->dfdx; + const double *dfcdx = d->dfcdx; + double rho = d->rho, fval = d->fval; + const double *rhoc = d->rhoc, *fcval = d->fcval; + double *xcur = d->xcur; + double *gcval = d->gcval; + int i, j; + double val; + + val = d->gval = fval; + d->wval = 0; + for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]); + if (grad) { for (i = 0; i < m; ++i) grad[i] = fcval[i]; } + + for (j = 0; j < n; ++j) { + int x_pinned; + double a, b, dx, denom, denominv, ca, cb, sigma2; + + /* first, compute xcur[j] for y */ + a = dfdx[j]; + b = fabs(dfdx[j]) * sigma[j] + 0.5 * rho; + for (i = 0; i < m; ++i) { + a += dfcdx[i*n + j] * y[i]; + b += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i]; + } + sigma2 = sigma[j] * sigma[j]; + dx = -a * sigma2 / b; + xcur[j] = x[j] + dx; + if (xcur[j] > ub[j]) xcur[j] = ub[j]; + if (xcur[j] < lb[j]) xcur[j] = lb[j]; + if (xcur[j] > x[j] + 0.9*sigma[j]) xcur[j] = x[j] + 0.9*sigma[j]; + if (xcur[j] < x[j] - 0.9*sigma[j]) xcur[j] = x[j] - 0.9*sigma[j]; + x_pinned = xcur[j] != x[j] + dx; + dx = xcur[j] - x[j]; + + /* function value: */ + denom = sigma2 - dx*dx; + denominv = 1.0 / denom; + ca = sigma2 * dx; + cb = dx*dx; + val += (a*ca + b*cb) * denominv; + + /* update gval, wval, gcval */ + d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb) + * denominv; + d->wval += 0.5 * cb * denominv; + for (i = 0; i < m; ++i) + gcval[i] += (dfcdx[i*n+j] * ca + + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[j]) * cb) + * denominv; + + /* gradient */ + if (grad) + for (i = 0; i < m; ++i) { + double dady, dbdy; + dady = dfcdx[i*n + j]; + dbdy = fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]; + grad[i] += (dady*ca + dbdy*cb) * denominv; + if (!x_pinned) { /* dx/dy nonzero */ + int k; + double dxdy; + dxdy = (a*dbdy - dady*b) * (sigma2 / (b*b)); + grad[i] += + ((sigma2 * dfdx[j] + 2 * (fabs(dfdx[j])*sigma[j] + + 0.5*rho) * dx) + * denom + + 2*dx * (dfdx[j]*ca + (fabs(dfdx[j])*sigma[j] + + 0.5*rho) * cb)) + * (denominv*denominv * dxdy); + for (k = 0; k < m; ++k) + grad[i] += + ((sigma2 * dfcdx[k*n + j] + + 2 * (fabs(dfcdx[k*n + j])*sigma[j] + + 0.5*rhoc[k]) * dx) + * denom + + 2*dx * (dfcdx[k*n + j]*ca + + (fabs(dfcdx[k*n + j])*sigma[j] + + 0.5*rhoc[k]) * cb)) + * (denominv*denominv * dxdy) * y[k]; + } + } + } + + /* negate because we are maximizing */ + if (grad) { for (i = 0; i < m; ++i) grad[i] = -grad[i]; } + return -val; +} + +/***********************************************************************/ + nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, + void *fc_data_, ptrdiff_t fc_data_size, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, - nlopt_stopping *stop) + nlopt_stopping *stop, + nlopt_algorithm dual_alg, + double dual_tolrel, int dual_maxeval) { nlopt_result ret = NLOPT_SUCCESS; double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur; - int j, k = 0; + double *dfcdx, *dfcdx_cur, *fcval, *rhoc, *gcval, *y, *dual_lb, *dual_ub; + int i, j, k = 0; + char *fc_data = (char *) fc_data_; + dual_data dd; + int feasible; - sigma = (double *) malloc(sizeof(double) * 6*n); + sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*6)); if (!sigma) return NLOPT_OUT_OF_MEMORY; dfdx = sigma + n; dfdx_cur = dfdx + n; xcur = dfdx_cur + n; xprev = xcur + n; xprevprev = xprev + n; - for (j = 0; j < n; ++j) - sigma[j] = 0.5 * (ub[j] - lb[j]); + fcval = xprevprev + n; + rhoc = fcval + m; + gcval = rhoc + m; + dual_lb = gcval + m; + dual_ub = dual_lb + m; + y = dual_ub + m; + dfcdx = y + m; + dfcdx_cur = dfcdx + m*n; + + dd.n = n; + dd.x = x; + dd.lb = lb; + dd.ub = ub; + dd.sigma = sigma; + dd.dfdx = dfdx; + dd.dfcdx = dfcdx; + dd.fcval = fcval; + dd.rhoc = rhoc; + dd.xcur = xcur; + dd.gcval = gcval; + + for (j = 0; j < n; ++j) { + if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j])) + sigma[j] = 1.0; /* arbitrary default */ + else + sigma[j] = 0.5 * (ub[j] - lb[j]); + } rho = 1.0; + for (i = 0; i < m; ++i) { + rhoc[i] = 1.0; + dual_lb[i] = y[i] = 0.0; + dual_ub[i] = HUGE_VAL; + } - fcur = *minf = f(n, x, dfdx, f_data); - memcpy(xcur, x, sizeof(double) * n); + dd.fval = fcur = *minf = f(n, x, dfdx, f_data); stop->nevals++; + feasible = 1; + for (i = 0; i < m; ++i) { + fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_data_size * i); + feasible = feasible && (fcval[i] <= 0); + } + memcpy(xcur, x, sizeof(double) * n); + + if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */ + while (1) { /* outer iterations */ double fprev = fcur; if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; @@ -50,40 +205,52 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, memcpy(xprev, xcur, sizeof(double) * n); while (1) { /* inner iterations */ - double gcur = *minf, w = 0; - for (j = 0; j < n; ++j) { - /* because we have no constraint functions, minimizing - the MMA approximate function can be done analytically */ - double dx = -sigma[j]*sigma[j]*dfdx[j] - / (2*sigma[j]*fabs(dfdx[j]) + rho); - xcur[j] = x[j] + dx; - if (xcur[j] > x[j] + 0.9*sigma[j]) - xcur[j] = x[j] + 0.9*sigma[j]; - else if (xcur[j] < x[j] - 0.9*sigma[j]) - xcur[j] = x[j] - 0.9*sigma[j]; - if (xcur[j] > ub[j]) xcur[j] = ub[j]; - else if (xcur[j] < lb[j]) xcur[j] = lb[j]; - dx = xcur[j] - x[j]; - gcur += (sigma[j]*sigma[j]*dfdx[j]*dx - + sigma[j]*fabs(dfdx[j])*dx*dx - + 0.5*rho*dx*dx) / (sigma[j]*sigma[j]-dx*dx); - w += 0.5*rho*dx*dx / (sigma[j]*sigma[j]-dx*dx); - } + double min_dual; + int feasible_cur, inner_done; + + /* solve dual problem */ + dd.rho = rho; + nlopt_minimize(dual_alg, m, dual_func, &dd, + dual_lb, dual_ub, y, &min_dual, + -HUGE_VAL, dual_tolrel, 0, 0, NULL, dual_maxeval, + stop->maxtime - (nlopt_seconds() - stop->start)); + dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */ fcur = f(n, xcur, dfdx_cur, f_data); stop->nevals++; - if (fcur < *minf) { - *minf = fcur; + inner_done = dd.gval >= fcur; + feasible_cur = 1; + for (i = 0; i < m; ++i) { + fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n, + fc_data + fc_data_size * i); + feasible_cur = feasible_cur && (fcval_cur[i] <= 0); + inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]); + } + + if (fcur < *minf && (feasible_cur || !feasible)) { + feasible = feasible_cur; + dd.fval = *minf = fcur; + memcpy(fcval, fcval_cur, sizeof(double)*m); memcpy(x, xcur, sizeof(double)*n); memcpy(dfdx, dfdx_cur, sizeof(double)*n); + memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m); } if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED; else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; if (ret != NLOPT_SUCCESS) goto done; - if (gcur >= fcur) break; - rho = MIN(10*rho, 1.1 * (rho + (fcur - gcur) / w)); + if (inner_done) break; + + if (fcur > dd.gval) + rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval)); + for (i = 0; i < m; ++i) + if (fcval_cur[i] > dd.gcval[i]) + rhoc[i] = + MIN(10*rhoc[i], + 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i]) + / dd.wval)); + if (mma_verbose) printf("MMA inner iteration: rho -> %g\n", rho); } @@ -98,6 +265,8 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, rho = MAX(0.1 * rho, MMA_RHOMIN); if (mma_verbose) printf("MMA outer iteration: rho -> %g\n", rho); + for (i = 0; i < m; ++i) + rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN); if (k > 1) for (j = 0; j < n; ++j) { double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]); diff --git a/util/nlopt-util.h b/util/nlopt-util.h index e0eb134..fffb750 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -6,6 +6,8 @@ extern "C" { #endif /* __cplusplus */ +int nlopt_isinf(double x); + /* seconds timer */ extern double nlopt_seconds(void); extern unsigned long nlopt_time_seed(void);