From 968c07739b9ae38afb87fdf5d14ef300b9c85bca Mon Sep 17 00:00:00 2001 From: stevenj Date: Sun, 31 Aug 2008 23:27:52 -0400 Subject: [PATCH] first start at replacing original newuoa routines with nlopt versions that will allow bound constraints darcs-hash:20080901032752-c8de0-2f993fc0527b8c9f04c2e7dccc00ccd281ef0bc5.gz --- newuoa/newuoa.c | 159 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) diff --git a/newuoa/newuoa.c b/newuoa/newuoa.c index 0547c77..58f8fd1 100644 --- a/newuoa/newuoa.c +++ b/newuoa/newuoa.c @@ -28,6 +28,7 @@ #include #include #include +#include #include "newuoa.h" @@ -37,6 +38,64 @@ /*************************************************************************/ /* trsapp.f */ +typedef struct { + int npt; + double *xpt, *pq, *hq, *gq, *xopt; + double *hd; + int iter; +} quad_model_data; + +static double quad_model(int n, const double *x, double *grad, void *data) +{ + quad_model_data *d = (quad_model_data *) data; + const double *xpt = d->xpt, *pq = d->pq, *hq = d->hq, *gq = d->gq, *xopt = d->xopt; + double *hd = d->hd; + int npt = d->npt; + int i, j, k; + double val = 0; + + /* first, set hd to be Hessian matrix times x */ + memset(hd, 0, sizeof(double) * n); + /* implicit Hessian terms (a sum of outer products of xpt vectors) */ + for (k = 0; k < npt; ++k) { + double temp = 0; + for (j = 0; j < n; ++j) + temp += xpt[k + j * npt] * (xopt[j] + x[j]); + temp *= pq[k]; + for (i = 0; i < n; ++i) + hd[i] += temp * xpt[k + i * npt]; + } + /* explicit Hessian terms (stored as compressed lower triangle hq) */ + k = 0; + for (j = 0; j < n; ++j) { + for (i = 0; i < j; ++i) { + hd[j] += hq[k] * (xopt[i] + x[i]); + hd[i] += hq[k] * (xopt[j] + x[j]); + ++k; + } + hd[j] += hq[k++] * (xopt[j] + x[j]); + } + + for (i = 0; i < n; ++i) { + val += (gq[i] + 0.5 * hd[i]) * (xopt[i] + x[i]); + if (grad) grad[i] = gq[i] + hd[i]; + } + d->iter++; + return val; +} + +/* constraint function to enforce |x|^2 <= rho*rho */ +static double rho_constraint(int n, const double *x, double *grad, void *data) +{ + double rho = *((double *) data); + double val = - rho*rho; + int i; + for (i = 0; i < n; ++i) + val += x[i] * x[i]; + if (grad) for (i = 0; i < n; ++i) grad[i] = 2*x[i]; + return val; +} + static void trsapp_(int *n, int *npt, double *xopt, double *xpt, double *gq, double *hq, double *pq, double *delta, double *step, double *d__, double *g, @@ -323,6 +382,46 @@ L120: goto L90; } L160: + if (0) { + double *lb, *ub, minf, crv; + quad_model_data qmd; + qmd.npt = *npt; + qmd.xpt = &xpt[1 + 1 * xpt_dim1]; + qmd.pq = &pq[1]; + qmd.hq = &hq[1]; + qmd.gq = &gq[1]; + qmd.xopt = &xopt[1]; + qmd.hd = &hd[1]; + qmd.iter = 0; + lb = &g[1]; + ub = &hs[1]; + for (j = 0; j < *n; ++j) lb[j] = -(ub[j] = *delta); + memset(&d__[1], 0, sizeof(double) * *n); + nlopt_minimize_constrained(NLOPT_LD_MMA, *n, quad_model, &qmd, + 1, rho_constraint, delta, sizeof(double), + lb, ub, &d__[1], &minf, -HUGE_VAL, + 0., 0., 1e-8, 0, 1000, 0.); + printf("trsapp = %g vs. nlopt (%d iters) = %g (delta = %g):\n", + quad_model(*n, &step[1], 0, &qmd), qmd.iter,minf,*delta); + for (j = 1; j <= *n; ++j) + printf(" xopt[%d] = %g, step[%d] = %g vs. %g\n", j, xopt[j], j, step[j], d__[j]); + printf("|d|^2 - delta^2 = %g, |step|^2 - delta^2 = %g\n", + rho_constraint(*n, &d__[1], 0, delta), + rho_constraint(*n, &step[1], 0, delta)); + if (rho_constraint(*n, &d__[1], 0, delta) > -1e-6 * (*delta)*(*delta)) + crv = 0; + else { + for (j = 1; j <= *n; ++j) d__[j] = step[j] - xopt[j]; + quad_model(*n, &d__[1], &g[1], &qmd); + crv = gg = 0; + for (j = 1; j <= *n; ++j) { + crv += step[j] * (g[j] - gq[j]); + gg += step[j] * step[j]; + } + crv = crv / gg; + } + printf("crvmin = %g, crv = %g\n", *crvmin, crv); + } return; /* The following instructions act as a subroutine for setting the vector */ @@ -890,6 +989,40 @@ L340: /*************************************************************************/ /* biglag.f */ +typedef struct { + int npt, ndim, iter; + double *hcol, *xpt, *bmat, *xopt; + int flipsign; +} lag_data; + +/* the Lagrange function, whose absolute value biglag maximizes */ +static double lag(int n, const double *dx, double *grad, void *data) +{ + lag_data *d = (lag_data *) data; + int i, j, npt = d->npt, ndim = d->ndim; + const double *hcol = d->hcol, *xpt = d->xpt, *bmat = d->bmat, *xopt = d->xopt; + double val = 0; + for (j = 0; j < n; ++j) { + val += bmat[j * ndim] * (xopt[j] + dx[j]); + if (grad) grad[j] = bmat[j * ndim]; + } + for (i = 0; i < npt; ++i) { + double dot = 0; + for (j = 0; j < n; ++j) + dot += xpt[i + j * npt] * (xopt[j] + dx[j]); + val += 0.5 * hcol[i] * (dot * dot); + dot *= hcol[i]; + if (grad) for (j = 0; j < n; ++j) + grad[j] += dot * xpt[i + j * npt]; + } + if (d->flipsign) { + val = -val; + if (grad) for (j = 0; j < n; ++j) grad[j] = -grad[j]; + } + d->iter++; + return val; +} + static void biglag_(int *n, int *npt, double *xopt, double *xpt, double *bmat, double *zmat, int *idz, int *ndim, int *knew, double *delta, double *d__, @@ -1050,6 +1183,32 @@ static void biglag_(int *n, int *npt, double *xopt, s[i__] = gc[i__] + temp * gd[i__]; } + if (0) { + double minf, *lb, *ub; + lag_data ld; + ld.npt = *npt; + ld.ndim = *ndim; + ld.iter = 0; + ld.hcol = &hcol[1]; + ld.xpt = &xpt[1 + 1 * xpt_dim1]; + ld.bmat = &bmat[*knew + 1 * bmat_dim1]; + ld.xopt = &xopt[1]; + ld.flipsign = 0; + lb = &gc[1]; ub = &gd[1]; + for (j = 0; j < *n; ++j) + lb[j] = -(ub[j] = *delta); + ld.flipsign = lag(*n, &d__[1], 0, &ld) > 0; /* maximize if > 0 */ + nlopt_minimize_constrained(NLOPT_LD_MMA, *n, lag, &ld, + 1, rho_constraint, delta, sizeof(double), + lb, ub, &d__[1], &minf, -HUGE_VAL, + 0., 0., 1e-5, 0, 1000, 0.); + minf = -minf; + printf("biglag nlopt converged to %g in %d iters\n", minf, ld.iter); + for (j = 1; j <= *n; ++j) + printf(" d[%d] = %g\n", j, d__[j]); + return; + } + /* Begin the iteration by overwriting S with a vector that has the */ /* required length and direction, except that termination occurs if */ /* the given D and S are nearly parallel. */ -- 2.30.2