From: stevenj Date: Thu, 1 Jul 2010 22:24:52 +0000 (-0400) Subject: add mconstraint to C API (not yet tested), thanks to Dmitrey of OpenOpt for the sugge... X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=8a840814c441600a523ff3751f31e5f947b07627;p=nlopt.git add mconstraint to C API (not yet tested), thanks to Dmitrey of OpenOpt for the suggestion darcs-hash:20100701222452-c8de0-c653c88d49348e7d42440283baf3e0dd90d1588d.gz --- diff --git a/api/nlopt-internal.h b/api/nlopt-internal.h index e8a80ff..7de7c55 100644 --- a/api/nlopt-internal.h +++ b/api/nlopt-internal.h @@ -69,6 +69,8 @@ struct nlopt_opt_s { nlopt_opt local_opt; /* local optimizer */ unsigned stochastic_population; /* population size for stochastic algs */ double *dx; /* initial step sizes (length n) for nonderivative algs */ + + double *work; /* algorithm-specific workspace during optimization */ }; /*********************************************************************/ diff --git a/api/nlopt.h b/api/nlopt.h index 29b7d18..c1ebbc3 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -59,6 +59,11 @@ typedef double (*nlopt_func)(unsigned n, const double *x, double *gradient, /* NULL if not needed */ void *func_data); +typedef void (*nlopt_mfunc)(unsigned m, double *result, + unsigned n, const double *x, + double *gradient, /* NULL if not needed */ + void *func_data); + typedef enum { /* Naming conventions: @@ -205,12 +210,22 @@ NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_constraint(nlopt_opt opt, nlopt_func fc, void *fc_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_mconstraint(nlopt_opt opt, + unsigned m, + nlopt_mfunc fc, + void *fc_data, + const double *tol); NLOPT_EXTERN(nlopt_result) nlopt_remove_equality_constraints(nlopt_opt opt); NLOPT_EXTERN(nlopt_result) nlopt_add_equality_constraint(nlopt_opt opt, nlopt_func h, void *h_data, double tol); +NLOPT_EXTERN(nlopt_result) nlopt_add_equality_mconstraint(nlopt_opt opt, + unsigned m, + nlopt_mfunc h, + void *h_data, + const double *tol); /* stopping criteria: */ diff --git a/api/optimize.c b/api/optimize.c index fe2164f..92e0797 100644 --- a/api/optimize.c +++ b/api/optimize.c @@ -89,12 +89,15 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) { nlopt_opt data = (nlopt_opt) data_; double f; - unsigned i; + unsigned i, j; f = data->f((unsigned) n, x, NULL, data->f_data); *undefined = isnan(f) || nlopt_isinf(f); - for (i = 0; i < data->m && !*undefined; ++i) - if (data->fc[i].f((unsigned) n, x, NULL, data->fc[i].f_data) > 0) - *undefined = 1; + for (i = 0; i < data->m && !*undefined; ++i) { + nlopt_eval_constraint(data->work, NULL, data->fc+i, (unsigned) n, x); + for (j = 0; j < data->fc[i].m; ++j) + if (data->work[j] > 0) + *undefined = 1; + } return f; } @@ -213,16 +216,23 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); case NLOPT_GN_ORIG_DIRECT: - case NLOPT_GN_ORIG_DIRECT_L: + case NLOPT_GN_ORIG_DIRECT_L: { + direct_return_code dret; if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS; - switch (direct_optimize(f_direct, opt, ni, lb, ub, x, minf, - stop.maxeval, -1, 0.0, 0.0, - pow(stop.xtol_rel, (double) n), -1.0, - stop.minf_max, 0.0, - NULL, - algorithm == NLOPT_GN_ORIG_DIRECT - ? DIRECT_ORIGINAL - : DIRECT_GABLONSKY)) { + opt->work = (double*) malloc(sizeof(double) * + nlopt_max_constraint_dim(opt->m, + opt->fc)); + if (!opt->work) return NLOPT_OUT_OF_MEMORY; + dret = direct_optimize(f_direct, opt, ni, lb, ub, x, minf, + stop.maxeval, -1, 0.0, 0.0, + pow(stop.xtol_rel, (double) n), -1.0, + stop.minf_max, 0.0, + NULL, + algorithm == NLOPT_GN_ORIG_DIRECT + ? DIRECT_ORIGINAL + : DIRECT_GABLONSKY); + free(opt->work); opt->work = NULL; + switch (dret) { case DIRECT_INVALID_BOUNDS: case DIRECT_MAXFEVAL_TOOBIG: case DIRECT_INVALID_ARGS: @@ -241,6 +251,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) return NLOPT_XTOL_REACHED; case DIRECT_OUT_OF_MEMORY: return NLOPT_OUT_OF_MEMORY; + } break; } diff --git a/api/options.c b/api/options.c index 50a376e..db955cc 100644 --- a/api/options.c +++ b/api/options.c @@ -33,21 +33,26 @@ void NLOPT_STDCALL nlopt_destroy(nlopt_opt opt) { if (opt) { + unsigned i; if (opt->munge_on_destroy) { nlopt_munge munge = opt->munge_on_destroy; - unsigned i; munge(opt->f_data); for (i = 0; i < opt->m; ++i) munge(opt->fc[i].f_data); for (i = 0; i < opt->p; ++i) munge(opt->h[i].f_data); } + for (i = 0; i < opt->m; ++i) + free(opt->fc[i].tol); + for (i = 0; i < opt->p; ++i) + free(opt->h[i].tol); free(opt->lb); free(opt->ub); free(opt->xtol_abs); free(opt->fc); free(opt->h); nlopt_destroy(opt->local_opt); free(opt->dx); + free(opt->work); free(opt); } } @@ -84,6 +89,7 @@ nlopt_opt NLOPT_STDCALL nlopt_create(nlopt_algorithm algorithm, unsigned n) opt->local_opt = NULL; opt->stochastic_population = 0; opt->dx = NULL; + opt->work = NULL; if (n > 0) { opt->lb = (double *) malloc(sizeof(double) * (n)); @@ -108,6 +114,7 @@ oom: nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt) { nlopt_opt nopt = NULL; + unsigned i; if (opt) { nopt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s)); *nopt = *opt; @@ -116,6 +123,7 @@ nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt) nopt->m_alloc = nopt->p_alloc = 0; nopt->local_opt = NULL; nopt->dx = NULL; + nopt->work = NULL; opt->force_stop_child = NULL; nlopt_munge munge = nopt->munge_on_copy; @@ -141,12 +149,21 @@ nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt) * (opt->m)); if (!nopt->fc) goto oom; memcpy(nopt->fc, opt->fc, sizeof(nlopt_constraint) * (opt->m)); + for (i = 0; i < opt->m; ++i) nopt->fc[i].tol = NULL; if (munge) - for (unsigned i = 0; i < opt->m; ++i) + for (i = 0; i < opt->m; ++i) if (nopt->fc[i].f_data && !(nopt->fc[i].f_data = munge(nopt->fc[i].f_data))) goto oom; + for (i = 0; i < opt->m; ++i) + if (opt->fc[i].tol) { + nopt->fc[i].tol = (double *) malloc(sizeof(double) + * nopt->fc[i].m); + if (!nopt->fc[i].tol) goto oom; + memcpy(nopt->fc[i].tol, opt->fc[i].tol, + sizeof(double) * nopt->fc[i].m); + } } if (opt->p) { @@ -155,12 +172,21 @@ nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt) * (opt->p)); if (!nopt->h) goto oom; memcpy(nopt->h, opt->h, sizeof(nlopt_constraint) * (opt->p)); + for (i = 0; i < opt->p; ++i) nopt->h[i].tol = NULL; if (munge) - for (unsigned i = 0; i < opt->p; ++i) + for (i = 0; i < opt->p; ++i) if (nopt->h[i].f_data && !(nopt->h[i].f_data = munge(nopt->h[i].f_data))) goto oom; + for (i = 0; i < opt->p; ++i) + if (opt->h[i].tol) { + nopt->h[i].tol = (double *) malloc(sizeof(double) + * nopt->h[i].m); + if (!nopt->h[i].tol) goto oom; + memcpy(nopt->h[i].tol, opt->h[i].tol, + sizeof(double) * nopt->h[i].m); + } } if (opt->local_opt) { @@ -288,12 +314,15 @@ NLOPT_STDCALL nlopt_get_upper_bounds(nlopt_opt opt, double *ub) nlopt_result NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt) { + unsigned i; if (!opt) return NLOPT_INVALID_ARGS; if (opt->munge_on_destroy) { nlopt_munge munge = opt->munge_on_destroy; - for (unsigned i = 0; i < opt->m; ++i) + for (i = 0; i < opt->m; ++i) munge(opt->fc[i].f_data); } + for (i = 0; i < opt->m; ++i) + free(opt->fc[i].tol); free(opt->fc); opt->fc = NULL; opt->m = opt->m_alloc = 0; @@ -302,9 +331,19 @@ NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt) static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, nlopt_constraint **c, - nlopt_func fc, void *fc_data, - double tol) + unsigned fm, nlopt_func fc, nlopt_mfunc mfc, + void *fc_data, + const double *tol) { + double *tolcopy; + + if ((fc && mfc) || (fc && fm != 1) || (!fc && !mfc) || tol < 0) + return NLOPT_INVALID_ARGS; + + tolcopy = (double *) malloc(sizeof(double) * fm); + if (fm && !tolcopy) return NLOPT_OUT_OF_MEMORY; + memcpy(tolcopy, tol, sizeof(double) * fm); + *m += 1; if (*m > *m_alloc) { /* allocate by repeated doubling so that @@ -315,65 +354,94 @@ static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, * (*m_alloc)); if (!*c) { *m_alloc = *m = 0; + free(tolcopy); return NLOPT_OUT_OF_MEMORY; } } + (*c)[*m - 1].m = fm; (*c)[*m - 1].f = fc; + (*c)[*m - 1].mf = mfc; (*c)[*m - 1].f_data = fc_data; - (*c)[*m - 1].tol = tol; + (*c)[*m - 1].tol = tolcopy; return NLOPT_SUCCESS; } +static int inequality_ok(nlopt_algorithm algorithm) { + /* nonlinear constraints are only supported with some algorithms */ + return (algorithm == NLOPT_LD_MMA + || algorithm == NLOPT_LN_COBYLA + || AUGLAG_ALG(algorithm) + || algorithm == NLOPT_GN_ISRES + || algorithm == NLOPT_GN_ORIG_DIRECT + || algorithm == NLOPT_GN_ORIG_DIRECT_L); +} + +nlopt_result +NLOPT_STDCALL nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m, + nlopt_mfunc fc, void *fc_data, + const double *tol) +{ + if (!m) return NLOPT_SUCCESS; /* empty constraints are always ok */ + if (!opt || !inequality_ok(opt->algorithm)) return NLOPT_INVALID_ARGS; + return add_constraint(&opt->m, &opt->m_alloc, &opt->fc, + m, NULL, fc, fc_data, tol); +} + nlopt_result NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt, nlopt_func fc, void *fc_data, double tol) { - if (opt && fc && tol >= 0) { - /* nonlinear constraints are only supported with some algorithms */ - if (opt->algorithm != NLOPT_LD_MMA - && opt->algorithm != NLOPT_LN_COBYLA - && !AUGLAG_ALG(opt->algorithm) - && opt->algorithm != NLOPT_GN_ISRES - && opt->algorithm != NLOPT_GN_ORIG_DIRECT - && opt->algorithm != NLOPT_GN_ORIG_DIRECT_L) - return NLOPT_INVALID_ARGS; - return add_constraint(&opt->m, &opt->m_alloc, &opt->fc, - fc, fc_data, tol); - } - return NLOPT_INVALID_ARGS; + if (!opt || !inequality_ok(opt->algorithm)) return NLOPT_INVALID_ARGS; + return add_constraint(&opt->m, &opt->m_alloc, &opt->fc, + 1, fc, NULL, fc_data, &tol); } nlopt_result NLOPT_STDCALL nlopt_remove_equality_constraints(nlopt_opt opt) { + unsigned i; if (!opt) return NLOPT_INVALID_ARGS; if (opt->munge_on_destroy) { nlopt_munge munge = opt->munge_on_destroy; - for (unsigned i = 0; i < opt->p; ++i) + for (i = 0; i < opt->p; ++i) munge(opt->h[i].f_data); } + for (i = 0; i < opt->p; ++i) + free(opt->h[i].tol); free(opt->h); opt->h = NULL; opt->p = opt->p_alloc = 0; return NLOPT_SUCCESS; } +static int equality_ok(nlopt_algorithm algorithm) { + /* equality constraints (h(x) = 0) only via some algorithms */ + return (AUGLAG_ALG(algorithm) + || algorithm == NLOPT_GN_ISRES + || algorithm == NLOPT_LN_COBYLA); +} + +nlopt_result +NLOPT_STDCALL nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m, + nlopt_mfunc fc, void *fc_data, + const double *tol) +{ + if (!m) return NLOPT_SUCCESS; /* empty constraints are always ok */ + if (!opt || !equality_ok(opt->algorithm)) return NLOPT_INVALID_ARGS; + return add_constraint(&opt->p, &opt->p_alloc, &opt->h, + m, NULL, fc, fc_data, tol); +} + nlopt_result NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt, - nlopt_func h, void *h_data, + nlopt_func fc, void *fc_data, double tol) { - if (opt && h && tol >= 0) { - /* equality constraints (h(x) = 0) only via some algorithms */ - if (!AUGLAG_ALG(opt->algorithm) && opt->algorithm != NLOPT_GN_ISRES - && opt->algorithm != NLOPT_LN_COBYLA) - return NLOPT_INVALID_ARGS; - return add_constraint(&opt->p, &opt->p_alloc, &opt->h, - h, h_data, tol); - } - return NLOPT_INVALID_ARGS; + if (!opt || !equality_ok(opt->algorithm)) return NLOPT_INVALID_ARGS; + return add_constraint(&opt->p, &opt->p_alloc, &opt->h, + 1, fc, NULL, fc_data, &tol); } /*************************************************************************/ diff --git a/auglag/auglag.c b/auglag/auglag.c index a9928cc..b86b9ba 100644 --- a/auglag/auglag.c +++ b/auglag/auglag.c @@ -14,10 +14,10 @@ int auglag_verbose = 0; typedef struct { nlopt_func f; void *f_data; - int m; nlopt_constraint *fc; - int p; nlopt_constraint *h; + int m, mm; nlopt_constraint *fc; + int p, pp; nlopt_constraint *h; double rho, *lambda, *mu; - double *gradtmp; + double *restmp, *gradtmp; nlopt_stopping *stop; } auglag_data; @@ -26,28 +26,34 @@ static double auglag(unsigned n, const double *x, double *grad, void *data) { auglag_data *d = (auglag_data *) data; double *gradtmp = grad ? d->gradtmp : NULL; + double *restmp = d->restmp; double rho = d->rho; const double *lambda = d->lambda, *mu = d->mu; double L; - int i; - unsigned j; + int i, ii; + unsigned j, k; L = d->f(n, x, grad, d->f_data); - for (i = 0; i < d->p; ++i) { - double h; - h = d->h[i].f(n, x, gradtmp, d->h[i].f_data) + lambda[i] / rho; - L += 0.5 * rho * h*h; - if (grad) for (j = 0; j < n; ++j) grad[j] += (rho * h) * gradtmp[j]; + for (ii = i = 0; i < d->p; ++i) { + nlopt_eval_constraint(restmp, gradtmp, d->h + i, n, x); + for (k = 0; k < d->h[i].m; ++k) { + double h = restmp[k] + lambda[ii++] / rho; + L += 0.5 * rho * h*h; + if (grad) for (j = 0; j < n; ++j) + grad[j] += (rho * h) * gradtmp[k*n + j]; + } } - for (i = 0; i < d->m; ++i) { - double fc; - fc = d->fc[i].f(n, x, gradtmp, d->fc[i].f_data) + mu[i] / rho; - if (fc > 0) { - L += 0.5 * rho * fc*fc; - if (grad) for (j = 0; j < n; ++j) - grad[j] += (rho * fc) * gradtmp[j]; + for (ii = i = 0; i < d->m; ++i) { + nlopt_eval_constraint(restmp, gradtmp, d->fc + i, n, x); + for (k = 0; k < d->h[i].m; ++k) { + double fc = restmp[k] + mu[ii++] / rho; + if (fc > 0) { + L += 0.5 * rho * fc*fc; + if (grad) for (j = 0; j < n; ++j) + grad[j] += (rho * fc) * gradtmp[k*n + j]; + } } } @@ -71,8 +77,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, nlopt_result ret = NLOPT_SUCCESS; double ICM = HUGE_VAL, minf_penalty = HUGE_VAL, penalty; double *xcur = NULL, fcur; - int i, feasible, minf_feasible = 0; + int i, ii, k, feasible, minf_feasible = 0; int auglag_iters = 0; + int max_constraint_dim; /* magic parameters from Birgin & Martinez */ const double tau = 0.5, gam = 10; @@ -83,6 +90,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, d.p = p; d.h = h; d.stop = stop; + max_constraint_dim = MAX(nlopt_max_constraint_dim(m, fc), + nlopt_max_constraint_dim(p, h)); + /* whether we handle inequality constraints via the augmented Lagrangian penalty function, or directly in the sub-algorithm */ if (sub_has_fc) @@ -90,6 +100,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, else m = 0; + d.mm = nlopt_count_constraints(d.m, fc); + d.pp = nlopt_count_constraints(d.p, fc); + ret = nlopt_set_min_objective(sub_opt, auglag, &d); if (ret<0) return ret; ret = nlopt_set_lower_bounds(sub_opt, lb); if (ret<0) return ret; ret = nlopt_set_upper_bounds(sub_opt, ub); if (ret<0) return ret; @@ -97,19 +110,28 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, ret = nlopt_remove_inequality_constraints(sub_opt); if (ret<0) return ret; ret = nlopt_remove_equality_constraints(sub_opt); if (ret<0) return ret; for (i = 0; i < m; ++i) { - ret = nlopt_add_inequality_constraint(sub_opt, fc[i].f, fc[i].f_data, - fc[i].tol); + if (fc[i].f) + ret = nlopt_add_inequality_constraint(sub_opt, + fc[i].f, fc[i].f_data, + fc[i].tol[0]); + else + ret = nlopt_add_inequality_mconstraint(sub_opt, fc[i].m, + fc[i].mf, fc[i].f_data, + fc[i].tol); if (ret < 0) return ret; } - xcur = (double *) malloc(sizeof(double) * (2*n + d.p + d.m)); + xcur = (double *) malloc(sizeof(double) * (n + + max_constraint_dim * (1 + n) + + d.pp + d.mm)); if (!xcur) return NLOPT_OUT_OF_MEMORY; memcpy(xcur, x, sizeof(double) * n); - d.gradtmp = xcur + n; - memset(d.gradtmp, 0, sizeof(double) * (n + d.p + d.m)); - d.lambda = d.gradtmp + n; - d.mu = d.lambda + d.p; + d.restmp = xcur + n; + d.gradtmp = d.restmp + max_constraint_dim; + memset(d.gradtmp, 0, sizeof(double) * (n*max_constraint_dim + d.pp+d.mm)); + d.lambda = d.gradtmp + n * max_constraint_dim; + d.mu = d.lambda + d.pp; *minf = HUGE_VAL; @@ -121,16 +143,22 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, penalty = 0; feasible = 1; for (i = 0; i < d.p; ++i) { - double hi = h[i].f(n, xcur, NULL, d.h[i].f_data); - penalty += fabs(hi); - feasible = feasible && fabs(hi) <= h[i].tol; - con2 += hi * hi; + nlopt_eval_constraint(d.restmp, NULL, d.h + i, n, xcur); + for (k = 0; k < d.h[i].m; ++k) { + double hi = d.restmp[k]; + penalty += fabs(hi); + feasible = feasible && fabs(hi) <= h[i].tol[k]; + con2 += hi * hi; + } } for (i = 0; i < d.m; ++i) { - double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data); - penalty += fci > 0 ? fci : 0; - feasible = feasible && fci <= fc[i].tol; - if (fci > 0) con2 += fci * fci; + nlopt_eval_constraint(d.restmp, NULL, d.fc + i, n, xcur); + for (k = 0; k < d.fc[i].m; ++k) { + double fci = d.restmp[k]; + penalty += fci > 0 ? fci : 0; + feasible = feasible && fci <= fc[i].tol[k]; + if (fci > 0) con2 += fci * fci; + } } *minf = fcur; minf_penalty = penalty; @@ -142,9 +170,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, if (auglag_verbose) { printf("auglag: initial rho=%g\nauglag initial lambda=", d.rho); - for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]); + for (i = 0; i < d.pp; ++i) printf(" %g", d.lambda[i]); printf("\nauglag initial mu = "); - for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]); + for (i = 0; i < d.mm; ++i) printf(" %g", d.mu[i]); printf("\n"); } @@ -163,21 +191,27 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, ICM = 0; penalty = 0; feasible = 1; - for (i = 0; i < d.p; ++i) { - double hi = h[i].f(n, xcur, NULL, d.h[i].f_data); - double newlam = d.lambda[i] + d.rho * hi; - penalty += fabs(hi); - feasible = feasible && fabs(hi) <= h[i].tol; - ICM = MAX(ICM, fabs(hi)); - d.lambda[i] = MIN(MAX(lam_min, newlam), lam_max); + for (i = ii = 0; i < d.p; ++i) { + nlopt_eval_constraint(d.restmp, NULL, d.h + i, n, xcur); + for (k = 0; k < d.h[i].m; ++k) { + double hi = d.restmp[k]; + double newlam = d.lambda[ii] + d.rho * hi; + penalty += fabs(hi); + feasible = feasible && fabs(hi) <= h[i].tol[k]; + ICM = MAX(ICM, fabs(hi)); + d.lambda[ii++] = MIN(MAX(lam_min, newlam), lam_max); + } } - for (i = 0; i < d.m; ++i) { - double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data); - double newmu = d.mu[i] + d.rho * fci; - penalty += fci > 0 ? fci : 0; - feasible = feasible && fci <= fc[i].tol; - ICM = MAX(ICM, fabs(MAX(fci, -d.mu[i] / d.rho))); - d.mu[i] = MIN(MAX(0.0, newmu), mu_max); + for (i = ii = 0; i < d.m; ++i) { + nlopt_eval_constraint(d.restmp, NULL, d.fc + i, n, xcur); + for (k = 0; k < d.fc[i].m; ++k) { + double fci = d.restmp[k]; + double newmu = d.mu[ii] + d.rho * fci; + penalty += fci > 0 ? fci : 0; + feasible = feasible && fci <= fc[i].tol[k]; + ICM = MAX(ICM, fabs(MAX(fci, -d.mu[ii] / d.rho))); + d.mu[ii++] = MIN(MAX(0.0, newmu), mu_max); + } } if (ICM > tau * prev_ICM) { d.rho *= gam; @@ -188,9 +222,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, if (auglag_verbose) { printf("auglag %d: ICM=%g (%sfeasible), rho=%g\nauglag lambda=", auglag_iters, ICM, feasible ? "" : "not ", d.rho); - for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]); + for (i = 0; i < d.pp; ++i) printf(" %g", d.lambda[i]); printf("\nauglag %d: mu = ", auglag_iters); - for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]); + for (i = 0; i < d.mm; ++i) printf(" %g", d.mu[i]); printf("\n"); } @@ -223,7 +257,9 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, Besides the fact that these kinds of absolute tolerances (non-scale-invariant) are unsatisfying and it is not clear how the user should specify it, the ICM <= epsilon - condition seems not too different from requiring feasibility. */ + condition seems not too different from requiring feasibility, + especially now that the user can provide constraint-specific + tolerances analogous to epsilon. */ if (ICM == 0) return NLOPT_FTOL_REACHED; } while (1); diff --git a/cobyla/cobyla.c b/cobyla/cobyla.c index fe705be..ea981cd 100644 --- a/cobyla/cobyla.c +++ b/cobyla/cobyla.c @@ -70,12 +70,13 @@ typedef struct { nlopt_constraint *h; double *xtmp; const double *lb, *ub; + double *con_tol; } func_wrap_state; static int func_wrap(int n, int m, double *x, double *f, double *con, func_wrap_state *s) { - int i, j; + int i, j, k; double *xtmp = s->xtmp; const double *lb = s->lb, *ub = s->ub; @@ -91,12 +92,18 @@ static int func_wrap(int n, int m, double *x, double *f, double *con, } *f = s->f(n, xtmp, NULL, s->f_data); - for (i = 0; i < s->m_orig; ++i) - con[i] = -s->fc[i].f(n, xtmp, NULL, s->fc[i].f_data); + i = 0; + for (j = 0; j < s->m_orig; ++j) { + nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp); + for (k = 0; k < s->fc[j].m; ++k) + con[i + k] = -con[i + k]; + i += s->fc[j].m; + } for (j = 0; j < s->p; ++j) { - double h = s->h[j].f(n, xtmp, NULL, s->h[j].f_data); - con[i++] = h; - con[i++] = -h; + nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp); + for (k = 0; k < s->h[j].m; ++k) + con[(i + s->h[j].m) + k] = -con[i + k]; + i += 2 * s->h[j].m; } for (j = 0; j < n; ++j) { if (!nlopt_isinf(lb[j])) @@ -168,7 +175,7 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, nlopt_stopping *stop, double rhobegin) { - int j; + int i, j; func_wrap_state s; nlopt_result ret; @@ -182,7 +189,7 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, if (!s.xtmp) return NLOPT_OUT_OF_MEMORY; /* each equality constraint gives two inequality constraints */ - m += 2*p; + m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h); /* add constraints for lower/upper bounds (if any) */ for (j = 0; j < n; ++j) { @@ -192,6 +199,23 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, ++m; } + s.con_tol = (double *) malloc(sizeof(double) * m); + if (m && !s.con_tol) { + free(s.xtmp); + return NLOPT_OUT_OF_MEMORY; + } + for (j = 0; j < m; ++j) s.con_tol[j] = 0; + for (j = i = 0; i < s.m_orig; ++i) { + int j0 = j, jnext = j + fc[i].m; + for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - j0]; + } + for (i = 0; i < s.p; ++i) { + int j0 = j, jnext = j + h[i].m; + for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0]; + j0 = j; jnext = j + h[i].m; + for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0]; + } + ret = cobyla(n, m, x, minf, rhobegin, stop, lb, ub, COBYLA_MSG_NONE, func_wrap, &s); @@ -201,6 +225,7 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, if (x[j] > ub[j]) x[j] = ub[j]; } + free(s.con_tol); free(s.xtmp); return ret; } @@ -534,7 +559,7 @@ L40: for (k = 1; k <= i__1; ++k) { d__1 = resmax, d__2 = -con[k]; resmax = max(d__1,d__2); - if (d__2 > (k <= state->m_orig ? state->fc[k-1].tol : 0)) + if (d__2 > state->con_tol[k-1]) feasible = 0; /* SGJ, 2010 */ } } diff --git a/isres/isres.c b/isres/isres.c index 47cc65c..5b4fec3 100644 --- a/isres/isres.c +++ b/isres/isres.c @@ -53,6 +53,8 @@ static int key_compare(void *keys_, const void *a_, const void *b_) return keys[*a] < keys[*b] ? -1 : (keys[*a] > keys[*b] ? +1 : 0); } +static unsigned imax2(unsigned a, unsigned b) { return (a > b ? a : b); } + nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, int m, nlopt_constraint *fc, /* fc <= 0 */ int p, nlopt_constraint *h, /* h == 0 */ @@ -78,6 +80,8 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, int mp = m + p; double minf_penalty = HUGE_VAL, minf_gpenalty = HUGE_VAL; double taup, tau; + double *results = 0; /* scratch space for mconstraint results */ + unsigned ires; *minf = HUGE_VAL; @@ -92,11 +96,16 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, for (j = 0; j < n; ++j) if (nlopt_isinf(lb[j]) || nlopt_isinf(ub[j])) return NLOPT_INVALID_ARGS; + ires = imax2(nlopt_max_constraint_dim(m, fc), + nlopt_max_constraint_dim(p, h)); + results = (double *) malloc(ires * sizeof(double)); + if (ires > 0 && !results) return NLOPT_OUT_OF_MEMORY; + sigmas = (double*) malloc(sizeof(double) * (population*n*2 + population + population + n)); - if (!sigmas) return NLOPT_OUT_OF_MEMORY; + if (!sigmas) { free(results); return NLOPT_OUT_OF_MEMORY; } xs = sigmas + population*n; fval = xs + population*n; penalty = fval + population; @@ -124,16 +133,24 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, fval[k] = f(n, xs + k*n, NULL, f_data); penalty[k] = 0; for (c = 0; c < m; ++c) { /* inequality constraints */ - double gval = fc[c].f(n, xs + k*n, NULL, fc[c].f_data); - if (gval > fc[c].tol) feasible = 0; - if (gval < 0) gval = 0; - penalty[k] += gval*gval; + nlopt_eval_constraint(results, NULL, + fc + c, n, xs + k*n); + for (ires = 0; ires < fc[c].m; ++ires) { + double gval = results[ires]; + if (gval > fc[c].tol[ires]) feasible = 0; + if (gval < 0) gval = 0; + penalty[k] += gval*gval; + } } gpenalty = penalty[k]; for (c = m; c < mp; ++c) { /* equality constraints */ - double hval = h[c-m].f(n, xs + k*n, NULL, h[c-m].f_data); - if (fabs(hval) > h[c-m].tol) feasible = 0; - penalty[k] += hval*hval; + nlopt_eval_constraint(results, NULL, + h + (c-m), n, xs + k*n); + for (ires = 0; ires < h[c-m].m; ++ires) { + double hval = results[ires]; + if (fabs(hval) > h[c-m].tol[ires]) feasible = 0; + penalty[k] += hval*hval; + } } if (penalty[k] > 0) all_feasible = 0; @@ -255,5 +272,6 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, done: if (irank) free(irank); if (sigmas) free(sigmas); + if (results) free(results); return ret; } diff --git a/mma/mma.c b/mma/mma.c index 9de104f..7d55335 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -157,11 +157,13 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur; double *dfcdx, *dfcdx_cur; double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub; - unsigned i, j, k = 0; + unsigned i, ifc, j, k = 0; dual_data dd; int feasible; double infeasibility; - + unsigned mfc; + + m = nlopt_count_constraints(mfc = m, fc); sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7)); if (!sigma) return NLOPT_OUT_OF_MEMORY; dfdx = sigma + n; @@ -209,8 +211,12 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, memcpy(xcur, x, sizeof(double) * n); feasible = 1; infeasibility = 0; + for (i = ifc = 0; ifc < mfc; ++ifc) { + nlopt_eval_constraint(fcval + i, dfcdx + i*n, + fc + ifc, n, x); + i += fc[ifc].m; + } for (i = 0; i < m; ++i) { - fcval[i] = fc[i].f(n, x, dfcdx + i*n, fc[i].f_data); feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i])); if (fcval[i] > infeasibility) infeasibility = fcval[i]; } @@ -279,20 +285,25 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, feasible_cur = 1; infeasibility_cur = 0; new_infeasible_constraint = 0; inner_done = dd.gval >= fcur; - for (i = 0; i < m; ++i) { - fcval_cur[i] = fc[i].f(n, xcur, dfcdx_cur + i*n, - fc[i].f_data); - if (!isnan(fcval_cur[i])) { - feasible_cur = feasible_cur - && (fcval_cur[i] <= fc[i].tol); - if (!isnan(fcval[i])) - inner_done = inner_done && - (dd.gcval[i] >= fcval_cur[i]); - else if (fcval_cur[i] > 0) - new_infeasible_constraint = 1; - if (fcval_cur[i] > infeasibility_cur) - infeasibility_cur = fcval_cur[i]; - } + for (i = ifc = 0; ifc < mfc; ++ifc) { + nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n, + fc + ifc, n, xcur); + i += fc[ifc].m; + } + for (i = ifc = 0; ifc < mfc; ++ifc) { + unsigned i0 = i, inext = i + fc[ifc].m; + for (; i < inext; ++i) + if (!isnan(fcval_cur[i])) { + feasible_cur = feasible_cur + && (fcval_cur[i] <= fc[ifc].tol[i-i0]); + if (!isnan(fcval[i])) + inner_done = inner_done && + (dd.gcval[i] >= fcval_cur[i]); + else if (fcval_cur[i] > 0) + new_infeasible_constraint = 1; + if (fcval_cur[i] > infeasibility_cur) + infeasibility_cur = fcval_cur[i]; + } } if ((fcur < *minf && (inner_done || feasible_cur || !feasible)) diff --git a/util/nlopt-util.h b/util/nlopt-util.h index 715e669..484c3b7 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -102,11 +102,19 @@ extern nlopt_result nlopt_optimize_limited(nlopt_opt opt, "feasible" for purposes of stopping of the constraint is violated by at most tol. */ typedef struct { - nlopt_func f; + unsigned m; /* dimensional of constraint: mf maps R^n -> R^m */ + nlopt_func f; /* one-dimensional constraint, requires m == 1 */ + nlopt_mfunc mf; void *f_data; - double tol; + double *tol; } nlopt_constraint; +extern unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint *c); +extern unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint *c); +extern void nlopt_eval_constraint(double *result, double *grad, + const nlopt_constraint *c, + unsigned n, const double *x); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/util/stop.c b/util/stop.c index eac93dc..b9cc05d 100644 --- a/util/stop.c +++ b/util/stop.c @@ -100,3 +100,30 @@ int nlopt_stop_forced(const nlopt_stopping *stop) { return stop->force_stop && *(stop->force_stop); } + +unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint *c) +{ + unsigned i, count = 0; + for (i = 0; i < p; ++i) + count += c[i].m; + return count; +} + +unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint *c) +{ + unsigned i, max_dim = 0; + for (i = 0; i < p; ++i) + if (c[i].m > max_dim) + max_dim = c[i].m; + return max_dim; +} + +void nlopt_eval_constraint(double *result, double *grad, + const nlopt_constraint *c, + unsigned n, const double *x) +{ + if (c->f) + result[0] = c->f(n, x, grad, c->f_data); + else + c->mf(c->m, result, n, x, grad, c->f_data); +}