From a1de53f267bc813a7d161517e9bdf0452bb0fedb Mon Sep 17 00:00:00 2001 From: stevenj Date: Tue, 15 Nov 2011 17:19:29 -0500 Subject: [PATCH] added (untested) interface for specifying preconditioners Ignore-this: 246b0790dff127b54fbf72d34d2da402 darcs-hash:20111115221929-c8de0-a27654751310b4251bdf6961ca7f1f3fc13164f4.gz --- api/nlopt-internal.h | 1 + api/nlopt.h | 17 ++++++++--- api/optimize.c | 23 +++++++++----- api/options.c | 71 ++++++++++++++++++++++++++++++++++---------- mma/ccsa_quadratic.c | 48 ++++++++++++++++++++---------- mma/mma.h | 3 +- util/nlopt-util.h | 1 + 7 files changed, 120 insertions(+), 44 deletions(-) diff --git a/api/nlopt-internal.h b/api/nlopt-internal.h index 3b621d1..d41099f 100644 --- a/api/nlopt-internal.h +++ b/api/nlopt-internal.h @@ -38,6 +38,7 @@ struct nlopt_opt_s { unsigned n; /* the dimension of the problem (immutable) */ nlopt_func f; void *f_data; /* objective function to minimize */ + nlopt_precond pre; /* optional preconditioner for f (NULL if none) */ int maximize; /* nonzero if we are maximizing, not minimizing */ double *lb, *ub; /* lower and upper bounds (length n) */ diff --git a/api/nlopt.h b/api/nlopt.h index 2ab6ab6..46584e9 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -68,10 +68,10 @@ typedef void (*nlopt_mfunc)(unsigned m, double *result, double *gradient, /* NULL if not needed */ void *func_data); -/* a preconditioner, which computes the approximate Hessian H(x) at x. - In particular, it returns Hdx = H(x) * dx. [Array lengths = n.] */ -typedef void (*nlopt_precond)(unsigned n, const double *x, const double *dx, - double *Hdx, void *data); +/* A preconditioner, which preconditions v at x to return vpre. + (The meaning of "preconditioning" is algorithm-dependent.) */ +typedef void (*nlopt_precond)(unsigned n, const double *x, const double *v, + double *vpre, void *data); typedef enum { /* Naming conventions: @@ -202,6 +202,9 @@ NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, void *f_data); +NLOPT_EXTERN(nlopt_result) nlopt_set_precond_min_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data); +NLOPT_EXTERN(nlopt_result) nlopt_set_precond_max_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data); + NLOPT_EXTERN(nlopt_algorithm) nlopt_get_algorithm(const nlopt_opt opt); NLOPT_EXTERN(unsigned) nlopt_get_dimension(const nlopt_opt opt); @@ -223,6 +226,9 @@ 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_precond_inequality_constraint( + nlopt_opt opt, nlopt_func fc, nlopt_precond pre, void *fc_data, + double tol); NLOPT_EXTERN(nlopt_result) nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m, nlopt_mfunc fc, @@ -234,6 +240,9 @@ 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_precond_equality_constraint( + nlopt_opt opt, nlopt_func h, nlopt_precond pre, void *h_data, + double tol); NLOPT_EXTERN(nlopt_result) nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m, nlopt_mfunc h, diff --git a/api/optimize.c b/api/optimize.c index 236519f..75aa487 100644 --- a/api/optimize.c +++ b/api/optimize.c @@ -661,8 +661,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) lb, ub, x, minf, &stop, dual_opt); else ret = ccsa_quadratic_minimize( - n, f, f_data, opt->m, opt->fc, - NULL, NULL, NULL, NULL, + n, f, f_data, opt->m, opt->fc, opt->pre, lb, ub, x, minf, &stop, dual_opt); nlopt_destroy(dual_opt); return ret; @@ -797,6 +796,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) typedef struct { nlopt_func f; + nlopt_precond pre; void *f_data; } f_max_data; @@ -813,22 +813,31 @@ static double f_max(unsigned n, const double *x, double *grad, void *data) return -val; } +static void pre_max(unsigned n, const double *x, const double *v, + double *vpre, void *data) +{ + f_max_data *d = (f_max_data *) data; + unsigned i; + d->pre(n, x, v, vpre, data); + for (i = 0; i < n; ++i) vpre[i] = -vpre[i]; +} + nlopt_result NLOPT_STDCALL nlopt_optimize(nlopt_opt opt, double *x, double *opt_f) { - nlopt_func f; void *f_data; + nlopt_func f; void *f_data; nlopt_precond pre; f_max_data fmd; int maximize; nlopt_result ret; if (!opt || !opt_f || !opt->f) return NLOPT_INVALID_ARGS; - f = opt->f; f_data = opt->f_data; + f = opt->f; f_data = opt->f_data; pre = opt->pre; /* for maximizing, just minimize the f_max wrapper, which flips the sign of everything */ if ((maximize = opt->maximize)) { - fmd.f = f; fmd.f_data = f_data; - opt->f = f_max; opt->f_data = &fmd; + fmd.f = f; fmd.f_data = f_data; fmd.pre = pre; + opt->f = f_max; opt->f_data = &fmd; opt->pre = pre_max; opt->stopval = -opt->stopval; opt->maximize = 0; } @@ -853,7 +862,7 @@ done: if (maximize) { /* restore original signs */ opt->maximize = maximize; opt->stopval = -opt->stopval; - opt->f = f; opt->f_data = f_data; + opt->f = f; opt->f_data = f_data; opt->pre = pre; *opt_f = -*opt_f; } diff --git a/api/options.c b/api/options.c index 2e7ae03..f3a4320 100644 --- a/api/options.c +++ b/api/options.c @@ -68,7 +68,7 @@ nlopt_opt NLOPT_STDCALL nlopt_create(nlopt_algorithm algorithm, unsigned n) if (opt) { opt->algorithm = algorithm; opt->n = n; - opt->f = NULL; opt->f_data = NULL; + opt->f = NULL; opt->f_data = NULL; opt->pre = NULL; opt->maximize = 0; opt->munge_on_destroy = opt->munge_on_copy = NULL; @@ -212,12 +212,14 @@ oom: /*************************************************************************/ -nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt, - nlopt_func f, void *f_data) +nlopt_result NLOPT_STDCALL nlopt_set_precond_min_objective(nlopt_opt opt, + nlopt_func f, + nlopt_precond pre, + void *f_data) { if (opt) { if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data); - opt->f = f; opt->f_data = f_data; + opt->f = f; opt->f_data = f_data; opt->pre = pre; opt->maximize = 0; if (nlopt_isinf(opt->stopval) && opt->stopval > 0) opt->stopval = -HUGE_VAL; /* switch default from max to min */ @@ -226,12 +228,20 @@ nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt, return NLOPT_INVALID_ARGS; } -nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, +nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, void *f_data) +{ + nlopt_set_precond_min_objective(opt, f, NULL, f_data); +} + +nlopt_result NLOPT_STDCALL nlopt_set_precond_max_objective(nlopt_opt opt, + nlopt_func f, + nlopt_precond pre, + void *f_data) { if (opt) { if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data); - opt->f = f; opt->f_data = f_data; + opt->f = f; opt->f_data = f_data; opt->pre = pre; opt->maximize = 1; if (nlopt_isinf(opt->stopval) && opt->stopval < 0) opt->stopval = +HUGE_VAL; /* switch default from min to max */ @@ -240,6 +250,12 @@ nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, return NLOPT_INVALID_ARGS; } +nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, + nlopt_func f, void *f_data) +{ + nlopt_set_precond_max_objective(opt, f, NULL, f_data); +} + /*************************************************************************/ nlopt_result @@ -336,6 +352,7 @@ NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt) static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, nlopt_constraint **c, unsigned fm, nlopt_func fc, nlopt_mfunc mfc, + nlopt_precond pre, void *fc_data, const double *tol) { @@ -371,6 +388,7 @@ static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, (*c)[*m - 1].m = fm; (*c)[*m - 1].f = fc; + (*c)[*m - 1].pre = pre; (*c)[*m - 1].mf = mfc; (*c)[*m - 1].f_data = fc_data; (*c)[*m - 1].tol = tolcopy; @@ -400,26 +418,37 @@ NLOPT_STDCALL nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m, } if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS; else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc, - m, NULL, fc, fc_data, tol); + m, NULL, fc, NULL, fc_data, tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return ret; } nlopt_result -NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt, - nlopt_func fc, void *fc_data, - double tol) +NLOPT_STDCALL nlopt_add_precond_inequality_constraint(nlopt_opt opt, + nlopt_func fc, + nlopt_precond pre, + void *fc_data, + double tol) { nlopt_result ret; if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS; else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc, - 1, fc, NULL, fc_data, &tol); + 1, fc, NULL, pre, fc_data, &tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return ret; } +nlopt_result +NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt, + nlopt_func fc, void *fc_data, + double tol) +{ + return nlopt_add_precond_inequality_constraint(opt, fc, NULL, fc_data, + tol); +} + nlopt_result NLOPT_STDCALL nlopt_remove_equality_constraints(nlopt_opt opt) { @@ -460,28 +489,38 @@ NLOPT_STDCALL nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m, || nlopt_count_constraints(opt->p, opt->h) + m > opt->n) ret = NLOPT_INVALID_ARGS; else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h, - m, NULL, fc, fc_data, tol); + m, NULL, fc, NULL, fc_data, tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return ret; } nlopt_result -NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt, - nlopt_func fc, void *fc_data, - double tol) +NLOPT_STDCALL nlopt_add_precond_equality_constraint(nlopt_opt opt, + nlopt_func fc, + nlopt_precond pre, + void *fc_data, + double tol) { nlopt_result ret; if (!opt || !equality_ok(opt->algorithm) || nlopt_count_constraints(opt->p, opt->h) + 1 > opt->n) ret = NLOPT_INVALID_ARGS; else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h, - 1, fc, NULL, fc_data, &tol); + 1, fc, NULL, pre, fc_data, &tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return ret; } +nlopt_result +NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt, + nlopt_func fc, void *fc_data, + double tol) +{ + return nlopt_add_precond_equality_constraint(opt, fc, NULL, fc_data, tol); +} + /*************************************************************************/ #define SET(param, T, arg) \ diff --git a/mma/ccsa_quadratic.c b/mma/ccsa_quadratic.c index dbca8fa..eb3715b 100644 --- a/mma/ccsa_quadratic.c +++ b/mma/ccsa_quadratic.c @@ -212,8 +212,7 @@ nlopt_result ccsa_quadratic_minimize( unsigned n, nlopt_func f, void *f_data, unsigned m, nlopt_constraint *fc, - nlopt_precond pre, void *pre_data, /* for f */ - nlopt_precond *prec, void **prec_data, /* length = # constraints */ + nlopt_precond pre, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ @@ -264,14 +263,30 @@ nlopt_result ccsa_quadratic_minimize( dd.rhoc = rhoc; dd.xcur = xcur; dd.gcval = gcval; - dd.pre = pre; dd.pre_data = pre_data; - dd.prec = prec; dd.prec_data = prec_data; + dd.pre = pre; dd.pre_data = f_data; + dd.prec = NULL; dd.prec_data = NULL; dd.scratch = NULL; + if (m) { + dd.prec = (nlopt_precond *) malloc(sizeof(nlopt_precond) * m); + dd.prec_data = (void **) malloc(sizeof(void *) * m); + if (!dd.prec || !dd.prec_data) { + ret = NLOPT_OUT_OF_MEMORY; + goto done; + } + for (i = ifc = 0; ifc < mfc; ++ifc) { + unsigned inext = i + fc[ifc].m; + for (; i < inext; ++i) { + dd.prec[i] = fc[ifc].pre; + dd.prec_data[i] = fc[ifc].f_data; + } + } + } + no_precond = pre == NULL; - if (prec) + if (dd.prec) for (i = 0; i < m; ++i) - no_precond = no_precond && prec[i] == NULL; + no_precond = no_precond && dd.prec[i] == NULL; if (!no_precond) { dd.scratch = (double*) malloc(sizeof(double) * (2*n)); @@ -430,15 +445,14 @@ nlopt_result ccsa_quadratic_minimize( } 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]; - inner_done = inner_done && - (dd.gcval[i] >= fcval_cur[i]); - if (fcval_cur[i] > infeasibility_cur) - infeasibility_cur = fcval_cur[i]; - } + for (; i < inext; ++i) { + feasible_cur = feasible_cur + && fcval_cur[i] <= fc[ifc].tol[i-i0]; + inner_done = inner_done && + (dd.gcval[i] >= fcval_cur[i]); + if (fcval_cur[i] > infeasibility_cur) + infeasibility_cur = fcval_cur[i]; + } } if ((fcur < *minf && (inner_done || feasible_cur || !feasible)) @@ -525,6 +539,10 @@ nlopt_result ccsa_quadratic_minimize( done: nlopt_destroy(pre_opt); if (dd.scratch) free(dd.scratch); + if (m) { + free(dd.prec_data); + free(dd.prec); + } free(sigma); return ret; } diff --git a/mma/mma.h b/mma/mma.h index a08b043..7968dc8 100644 --- a/mma/mma.h +++ b/mma/mma.h @@ -45,8 +45,7 @@ nlopt_result ccsa_quadratic_minimize( unsigned n, nlopt_func f, void *f_data, unsigned m, nlopt_constraint *fc, - nlopt_precond pre, void *pre_data, /* for f */ - nlopt_precond *prec, void **prec_data, /* length = # constraints */ + nlopt_precond pre, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ diff --git a/util/nlopt-util.h b/util/nlopt-util.h index ddbb801..ca33df4 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -111,6 +111,7 @@ typedef struct { unsigned m; /* dimensional of constraint: mf maps R^n -> R^m */ nlopt_func f; /* one-dimensional constraint, requires m == 1 */ nlopt_mfunc mf; + nlopt_precond pre; /* preconditioner for f (NULL if none or if mf) */ void *f_data; double *tol; } nlopt_constraint; -- 2.30.2