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 */
};
/*********************************************************************/
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:
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: */
{
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;
}
+ 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:
return NLOPT_XTOL_REACHED;
case DIRECT_OUT_OF_MEMORY:
return NLOPT_OUT_OF_MEMORY;
+ }
break;
}
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);
}
}
opt->local_opt = NULL;
opt->stochastic_population = 0;
opt->dx = NULL;
+ opt->work = NULL;
if (n > 0) {
opt->lb = (double *) malloc(sizeof(double) * (n));
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;
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;
* (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) {
* (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) {
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;
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
* (*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);
}
/*************************************************************************/
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;
{
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];
+ }
}
}
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;
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)
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;
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;
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;
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");
}
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;
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");
}
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);
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;
}
*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]))
nlopt_stopping *stop,
double rhobegin)
{
- int j;
+ int i, j;
func_wrap_state s;
nlopt_result ret;
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) {
++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);
if (x[j] > ub[j]) x[j] = ub[j];
}
+ free(s.con_tol);
free(s.xtmp);
return ret;
}
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 */
}
}
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 */
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;
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;
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;
done:
if (irank) free(irank);
if (sigmas) free(sigmas);
+ if (results) free(results);
return ret;
}
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;
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];
}
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))
"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 */
{
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);
+}
~ianmdlvl / nlopt.git / commitdiff