From b265b05006a3e6c90238634e0af74bf36f345f6b Mon Sep 17 00:00:00 2001 From: Nicola Piccinelli Date: Tue, 21 Jan 2020 13:17:56 +0100 Subject: [PATCH 1/3] Added matlab binding for vector valued constraints --- src/octave/nlopt_optimize-mex.c | 754 ++++++++++++++++++-------------- 1 file changed, 435 insertions(+), 319 deletions(-) diff --git a/src/octave/nlopt_optimize-mex.c b/src/octave/nlopt_optimize-mex.c index 499ac76e..2b094e37 100644 --- a/src/octave/nlopt_optimize-mex.c +++ b/src/octave/nlopt_optimize-mex.c @@ -22,358 +22,474 @@ /* Matlab MEX interface to NLopt, and in particular to nlopt_optimize */ +#include +#include #include #include #include -#include -#include #include "nlopt.h" -#define CHECK0(cond, msg) if (!(cond)) mexErrMsgTxt(msg); +#define CHECK0(cond, msg) \ + if (!(cond)) \ + mexErrMsgTxt(msg); -static double struct_val_default(const mxArray *s, const char *name, double dflt) +static double struct_val_default(const mxArray* s, const char* name, double dflt) { - mxArray *val = mxGetField(s, 0, name); - if (val) { - CHECK0(mxIsNumeric(val) && !mxIsComplex(val) - && mxGetM(val) * mxGetN(val) == 1, - "opt fields, other than xtol_abs, must be real scalars"); - return mxGetScalar(val); - } - return dflt; + mxArray* val = mxGetField(s, 0, name); + if (val) { + CHECK0(mxIsNumeric(val) && !mxIsComplex(val) + && mxGetM(val) * mxGetN(val) == 1, + "opt fields, other than xtol_abs, must be real scalars"); + return mxGetScalar(val); + } + return dflt; } -static double *struct_arrval(const mxArray *s, const char *name, unsigned n, - double *dflt) +static double* struct_arrval(const mxArray* s, const char* name, unsigned n, + double* dflt) { - mxArray *val = mxGetField(s, 0, name); - if (val) { - CHECK0(mxIsNumeric(val) && !mxIsComplex(val) - && mxGetM(val) * mxGetN(val) == n, - "opt vector field is not of length n"); - return mxGetPr(val); - } - return dflt; + mxArray* val = mxGetField(s, 0, name); + if (val) { + CHECK0(mxIsNumeric(val) && !mxIsComplex(val) + && mxGetM(val) * mxGetN(val) == n, + "opt vector field is not of length n"); + return mxGetPr(val); + } + return dflt; } -static mxArray *struct_funcval(const mxArray *s, const char *name) +static mxArray* struct_funcval(const mxArray* s, const char* name) { - mxArray *val = mxGetField(s, 0, name); - if (val) { - CHECK0(mxIsChar(val) || mxIsFunctionHandle(val), - "opt function field is not a function handle/name"); - return val; - } - return NULL; + mxArray* val = mxGetField(s, 0, name); + if (val) { + CHECK0(mxIsChar(val) || mxIsFunctionHandle(val), + "opt function field is not a function handle/name"); + return val; + } + return NULL; } -static double *fill(double *arr, unsigned n, double val) +static double* fill(double* arr, unsigned n, double val) { - unsigned i; - for (i = 0; i < n; ++i) arr[i] = val; - return arr; + unsigned i; + for (i = 0; i < n; ++i) + arr[i] = val; + return arr; } #define FLEN 128 /* max length of user function name */ #define MAXRHS 3 /* max nrhs for user function */ typedef struct user_function_data_s { - char f[FLEN]; - mxArray *plhs[2]; - mxArray *prhs[MAXRHS]; - int xrhs, nrhs; - int verbose, neval; - struct user_function_data_s *dpre; - nlopt_opt opt; + char f[FLEN]; + mxArray* plhs[2]; + mxArray* prhs[MAXRHS]; + int xrhs, nrhs; + int verbose, neval; + struct user_function_data_s* dpre; + nlopt_opt opt; } user_function_data; -static double user_function(unsigned n, const double *x, - double *gradient, /* NULL if not needed */ - void *d_) +static double user_function(unsigned n, const double* x, + double* gradient, /* NULL if not needed */ + void* d_) { - user_function_data *d = (user_function_data *) d_; - double f; - - d->plhs[0] = d->plhs[1] = NULL; - memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double)); - - CHECK0(0 == mexCallMATLAB(gradient ? 2 : 1, d->plhs, - d->nrhs, d->prhs, d->f), - "error calling user function"); - - CHECK0(mxIsNumeric(d->plhs[0]) && !mxIsComplex(d->plhs[0]) - && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == 1, - "user function must return real scalar"); - f = mxGetScalar(d->plhs[0]); - mxDestroyArray(d->plhs[0]); - if (gradient) { - CHECK0(mxIsDouble(d->plhs[1]) && !mxIsComplex(d->plhs[1]) - && (mxGetM(d->plhs[1]) == 1 || mxGetN(d->plhs[1]) == 1) - && mxGetM(d->plhs[1]) * mxGetN(d->plhs[1]) == n, - "gradient vector from user function is the wrong size"); - memcpy(gradient, mxGetPr(d->plhs[1]), n * sizeof(double)); - mxDestroyArray(d->plhs[1]); - } - d->neval++; - if (d->verbose) mexPrintf("nlopt_optimize eval #%d: %g\n", d->neval, f); - if (mxIsNaN(f)) nlopt_force_stop(d->opt); - return f; + user_function_data* d = (user_function_data*)d_; + double f; + + d->plhs[0] = d->plhs[1] = NULL; + memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double)); + + CHECK0(0 == mexCallMATLAB(gradient ? 2 : 1, d->plhs, d->nrhs, d->prhs, d->f), + "error calling user function"); + + CHECK0(mxIsNumeric(d->plhs[0]) && !mxIsComplex(d->plhs[0]) + && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == 1, + "user function must return real scalar"); + f = mxGetScalar(d->plhs[0]); + mxDestroyArray(d->plhs[0]); + if (gradient) { + CHECK0(mxIsDouble(d->plhs[1]) && !mxIsComplex(d->plhs[1]) + && (mxGetM(d->plhs[1]) == 1 || mxGetN(d->plhs[1]) == 1) + && mxGetM(d->plhs[1]) * mxGetN(d->plhs[1]) == n, + "gradient vector from user function is the wrong size"); + memcpy(gradient, mxGetPr(d->plhs[1]), n * sizeof(double)); + mxDestroyArray(d->plhs[1]); + } + d->neval++; + if (d->verbose) + mexPrintf("nlopt_optimize eval #%d: %g\n", d->neval, f); + if (mxIsNaN(f)) + nlopt_force_stop(d->opt); + return f; } -static void user_pre(unsigned n, const double *x, const double *v, - double *vpre, void *d_) +static void user_mfunction(unsigned m, double* result, unsigned n, const double* x, + double* gradient, /* NULL if not needed */ + void* d_) { - user_function_data *d = ((user_function_data *) d_)->dpre; - d->plhs[0] = d->plhs[1] = NULL; - memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double)); - memcpy(mxGetPr(d->prhs[d->xrhs + 1]), v, n * sizeof(double)); - - CHECK0(0 == mexCallMATLAB(1, d->plhs, - d->nrhs, d->prhs, d->f), - "error calling user function"); - - CHECK0(mxIsDouble(d->plhs[0]) && !mxIsComplex(d->plhs[0]) - && (mxGetM(d->plhs[0]) == 1 || mxGetN(d->plhs[0]) == 1) - && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == n, - "vpre vector from user function is the wrong size"); - memcpy(vpre, mxGetPr(d->plhs[0]), n * sizeof(double)); - mxDestroyArray(d->plhs[0]); - d->neval++; - if (d->verbose) mexPrintf("nlopt_optimize precond eval #%d\n", d->neval); + user_function_data* d = (user_function_data*)d_; + + d->plhs[0] = d->plhs[1] = NULL; + memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double)); + + CHECK0(0 == mexCallMATLAB(gradient ? 2 : 1, d->plhs, d->nrhs, d->prhs, d->f), + "error calling user mfunction"); + + CHECK0(mxIsNumeric(d->plhs[0]) && !mxIsComplex(d->plhs[0]) + && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == m, + "user mfunction must return an array of size equal to mfc_count parameter"); + memcpy(result, mxGetPr(d->plhs[0]), n * sizeof(double)); + mxDestroyArray(d->plhs[0]); + + if (gradient) { + CHECK0(mxIsDouble(d->plhs[1]) && !mxIsComplex(d->plhs[1]) + && (mxGetM(d->plhs[1]) == m || mxGetN(d->plhs[1]) == m) + && mxGetM(d->plhs[1]) * mxGetN(d->plhs[1]) == m * n, + "gradient vector from user mfunction is the wrong size"); + memcpy(gradient, mxGetPr(d->plhs[1]), m * n * sizeof(double)); + mxDestroyArray(d->plhs[1]); + } + + d->neval++; + + if (d->verbose) { + mexPrintf("nlopt_optimize eval #%d: ", d->neval); + for (size_t i = 0; i < n - 1; i++) { + mexPrintf("%g ", result[i]); + } + mexPrintf("%g\n", result[n - 1]); + } + + for (size_t i = 0; i < n - 1; i++) { + if (mxIsNaN(result[i])) { + nlopt_force_stop(d->opt); + break; + } + } } -#define CHECK1(cond, msg) if (!(cond)) { mxFree(tmp); nlopt_destroy(opt); nlopt_destroy(local_opt); mexWarnMsgTxt(msg); return NULL; }; - -nlopt_opt make_opt(const mxArray *opts, unsigned n) +static void user_pre(unsigned n, const double* x, const double* v, + double* vpre, void* d_) { - nlopt_opt opt = NULL, local_opt = NULL; - nlopt_algorithm algorithm; - double *tmp = NULL; - unsigned i; - - algorithm = (nlopt_algorithm) - struct_val_default(opts, "algorithm", NLOPT_NUM_ALGORITHMS); - CHECK1(((int)algorithm) >= 0 && algorithm < NLOPT_NUM_ALGORITHMS, - "invalid opt.algorithm"); - - tmp = (double *) mxCalloc(n, sizeof(double)); - opt = nlopt_create(algorithm, n); - CHECK1(opt, "nlopt: out of memory"); - - nlopt_set_lower_bounds(opt, struct_arrval(opts, "lower_bounds", n, - fill(tmp, n, -HUGE_VAL))); - nlopt_set_upper_bounds(opt, struct_arrval(opts, "upper_bounds", n, - fill(tmp, n, +HUGE_VAL))); - - nlopt_set_stopval(opt, struct_val_default(opts, "stopval", -HUGE_VAL)); - nlopt_set_ftol_rel(opt, struct_val_default(opts, "ftol_rel", 0.0)); - nlopt_set_ftol_abs(opt, struct_val_default(opts, "ftol_abs", 0.0)); - nlopt_set_xtol_rel(opt, struct_val_default(opts, "xtol_rel", 0.0)); - nlopt_set_xtol_abs(opt, struct_arrval(opts, "xtol_abs", n, - fill(tmp, n, 0.0))); - nlopt_set_x_weights(opt, struct_arrval(opts, "x_weights", n, - fill(tmp, n, 1.0))); - nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0.0) < 0 ? - 0 : struct_val_default(opts, "maxeval", 0.0)); - nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0)); - - nlopt_set_population(opt, struct_val_default(opts, "population", 0)); - nlopt_set_vector_storage(opt, struct_val_default(opts, "vector_storage", 0)); - - if (struct_arrval(opts, "initial_step", n, NULL)) - nlopt_set_initial_step(opt, - struct_arrval(opts, "initial_step", n, NULL)); - - if (mxGetField(opts, 0, "local_optimizer")) { - const mxArray *local_opts = mxGetField(opts, 0, "local_optimizer"); - CHECK1(mxIsStruct(local_opts), - "opt.local_optimizer must be a structure"); - CHECK1(local_opt = make_opt(local_opts, n), - "error initializing local optimizer"); - nlopt_set_local_optimizer(opt, local_opt); - nlopt_destroy(local_opt); local_opt = NULL; - } - - mxFree(tmp); - return opt; + user_function_data* d = ((user_function_data*)d_)->dpre; + d->plhs[0] = d->plhs[1] = NULL; + memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double)); + memcpy(mxGetPr(d->prhs[d->xrhs + 1]), v, n * sizeof(double)); + + CHECK0(0 == mexCallMATLAB(1, d->plhs, d->nrhs, d->prhs, d->f), + "error calling user function"); + + CHECK0(mxIsDouble(d->plhs[0]) && !mxIsComplex(d->plhs[0]) + && (mxGetM(d->plhs[0]) == 1 || mxGetN(d->plhs[0]) == 1) + && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == n, + "vpre vector from user function is the wrong size"); + memcpy(vpre, mxGetPr(d->plhs[0]), n * sizeof(double)); + mxDestroyArray(d->plhs[0]); + d->neval++; + if (d->verbose) + mexPrintf("nlopt_optimize precond eval #%d\n", d->neval); } -#define CHECK(cond, msg) if (!(cond)) { mxFree(dh); mxFree(dfc); nlopt_destroy(opt); mexErrMsgTxt(msg); } +#define CHECK1(cond, msg) \ + if (!(cond)) { \ + mxFree(tmp); \ + nlopt_destroy(opt); \ + nlopt_destroy(local_opt); \ + mexWarnMsgTxt(msg); \ + return NULL; \ + }; -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) +nlopt_opt make_opt(const mxArray* opts, unsigned n) { - unsigned n; - double *x, *x0, opt_f; - nlopt_result ret; - mxArray *x_mx, *mx; - user_function_data d, dpre, *dfc = NULL, *dh = NULL; - nlopt_opt opt = NULL; - - CHECK(nrhs == 2 && nlhs <= 3, "wrong number of arguments"); - - /* options = prhs[0] */ - CHECK(mxIsStruct(prhs[0]), "opt must be a struct"); - - /* x0 = prhs[1] */ - CHECK(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) - && (mxGetM(prhs[1]) == 1 || mxGetN(prhs[1]) == 1), - "x must be real row or column vector"); - n = mxGetM(prhs[1]) * mxGetN(prhs[1]), - x0 = mxGetPr(prhs[1]); - - CHECK(opt = make_opt(prhs[0], n), "error initializing nlopt options"); - - d.neval = 0; - d.verbose = (int) struct_val_default(prhs[0], "verbose", 0); - d.opt = opt; - - /* function f = prhs[1] */ - mx = struct_funcval(prhs[0], "min_objective"); - if (!mx) mx = struct_funcval(prhs[0], "max_objective"); - CHECK(mx, "either opt.min_objective or opt.max_objective must exist"); - if (mxIsChar(mx)) { - CHECK(mxGetString(mx, d.f, FLEN) == 0, - "error reading function name string (too long?)"); - d.nrhs = 1; - d.xrhs = 0; - } - else { - d.prhs[0] = mx; - strcpy(d.f, "feval"); - d.nrhs = 2; - d.xrhs = 1; - } - d.prhs[d.xrhs] = mxCreateDoubleMatrix(1, n, mxREAL); - - if ((mx = struct_funcval(prhs[0], "pre"))) { - CHECK(mxIsChar(mx) || mxIsFunctionHandle(mx), - "pre must contain function handles or function names"); - if (mxIsChar(mx)) { - CHECK(mxGetString(mx, dpre.f, FLEN) == 0, - "error reading function name string (too long?)"); - dpre.nrhs = 2; - dpre.xrhs = 0; - } - else { - dpre.prhs[0] = mx; - strcpy(dpre.f, "feval"); - dpre.nrhs = 3; - dpre.xrhs = 1; - } - dpre.verbose = d.verbose > 2; - dpre.opt = opt; - dpre.neval = 0; - dpre.prhs[dpre.xrhs] = d.prhs[d.xrhs]; - dpre.prhs[d.xrhs+1] = mxCreateDoubleMatrix(1, n, mxREAL); - d.dpre = &dpre; - - if (struct_funcval(prhs[0], "min_objective")) - nlopt_set_precond_min_objective(opt, user_function,user_pre,&d); - else - nlopt_set_precond_max_objective(opt, user_function,user_pre,&d); - } - else { - dpre.nrhs = 0; - if (struct_funcval(prhs[0], "min_objective")) - nlopt_set_min_objective(opt, user_function, &d); - else - nlopt_set_max_objective(opt, user_function, &d); - } - - if ((mx = mxGetField(prhs[0], 0, "fc"))) { - int j, m; - double *fc_tol; - - CHECK(mxIsCell(mx), "fc must be a Cell array"); - m = mxGetM(mx) * mxGetN(mx);; - dfc = (user_function_data *) mxCalloc(m, sizeof(user_function_data)); - fc_tol = struct_arrval(prhs[0], "fc_tol", m, NULL); - - for (j = 0; j < m; ++j) { - mxArray *fc = mxGetCell(mx, j); - CHECK(mxIsChar(fc) || mxIsFunctionHandle(fc), - "fc must contain function handles or function names"); - if (mxIsChar(fc)) { - CHECK(mxGetString(fc, dfc[j].f, FLEN) == 0, - "error reading function name string (too long?)"); - dfc[j].nrhs = 1; - dfc[j].xrhs = 0; - } - else { - dfc[j].prhs[0] = fc; - strcpy(dfc[j].f, "feval"); - dfc[j].nrhs = 2; - dfc[j].xrhs = 1; - } - dfc[j].verbose = d.verbose > 1; - dfc[j].opt = opt; - dfc[j].neval = 0; - dfc[j].prhs[dfc[j].xrhs] = d.prhs[d.xrhs]; - CHECK(nlopt_add_inequality_constraint(opt, user_function, - dfc + j, - fc_tol ? fc_tol[j] : 0) - > 0, "nlopt error adding inequality constraint"); - } - } - - - if ((mx = mxGetField(prhs[0], 0, "h"))) { - int j, m; - double *h_tol; - - CHECK(mxIsCell(mx), "h must be a Cell array"); - m = mxGetM(mx) * mxGetN(mx);; - dh = (user_function_data *) mxCalloc(m, sizeof(user_function_data)); - h_tol = struct_arrval(prhs[0], "h_tol", m, NULL); - - for (j = 0; j < m; ++j) { - mxArray *h = mxGetCell(mx, j); - CHECK(mxIsChar(h) || mxIsFunctionHandle(h), - "h must contain function handles or function names"); - if (mxIsChar(h)) { - CHECK(mxGetString(h, dh[j].f, FLEN) == 0, - "error reading function name string (too long?)"); - dh[j].nrhs = 1; - dh[j].xrhs = 0; - } - else { - dh[j].prhs[0] = h; - strcpy(dh[j].f, "feval"); - dh[j].nrhs = 2; - dh[j].xrhs = 1; - } - dh[j].verbose = d.verbose > 1; - dh[j].opt = opt; - dh[j].neval = 0; - dh[j].prhs[dh[j].xrhs] = d.prhs[d.xrhs]; - CHECK(nlopt_add_equality_constraint(opt, user_function, - dh + j, - h_tol ? h_tol[j] : 0) - > 0, "nlopt error adding equality constraint"); - } - } - - - x_mx = mxCreateDoubleMatrix(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL); - x = mxGetPr(x_mx); - memcpy(x, x0, sizeof(double) * n); - - ret = nlopt_optimize(opt, x, &opt_f); - - mxFree(dh); - mxFree(dfc); - mxDestroyArray(d.prhs[d.xrhs]); - if (dpre.nrhs > 0) mxDestroyArray(dpre.prhs[d.xrhs+1]); - nlopt_destroy(opt); - - plhs[0] = x_mx; - if (nlhs > 1) { - plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); - *(mxGetPr(plhs[1])) = opt_f; - } - if (nlhs > 2) { - plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL); - *(mxGetPr(plhs[2])) = (int) ret; - } + nlopt_opt opt = NULL, local_opt = NULL; + nlopt_algorithm algorithm; + double* tmp = NULL; + unsigned i; + + algorithm = (nlopt_algorithm) + struct_val_default(opts, "algorithm", NLOPT_NUM_ALGORITHMS); + CHECK1(((int)algorithm) >= 0 && algorithm < NLOPT_NUM_ALGORITHMS, + "invalid opt.algorithm"); + + tmp = (double*)mxCalloc(n, sizeof(double)); + opt = nlopt_create(algorithm, n); + CHECK1(opt, "nlopt: out of memory"); + + nlopt_set_lower_bounds(opt, struct_arrval(opts, "lower_bounds", n, fill(tmp, n, -HUGE_VAL))); + nlopt_set_upper_bounds(opt, struct_arrval(opts, "upper_bounds", n, fill(tmp, n, +HUGE_VAL))); + + nlopt_set_stopval(opt, struct_val_default(opts, "stopval", -HUGE_VAL)); + nlopt_set_ftol_rel(opt, struct_val_default(opts, "ftol_rel", 0.0)); + nlopt_set_ftol_abs(opt, struct_val_default(opts, "ftol_abs", 0.0)); + nlopt_set_xtol_rel(opt, struct_val_default(opts, "xtol_rel", 0.0)); + nlopt_set_xtol_abs(opt, struct_arrval(opts, "xtol_abs", n, fill(tmp, n, 0.0))); + nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0.0) < 0 ? 0 : struct_val_default(opts, "maxeval", 0.0)); + nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0)); + + nlopt_set_population(opt, struct_val_default(opts, "population", 0)); + nlopt_set_vector_storage(opt, struct_val_default(opts, "vector_storage", 0)); + + if (struct_arrval(opts, "initial_step", n, NULL)) + nlopt_set_initial_step(opt, + struct_arrval(opts, "initial_step", n, NULL)); + + if (mxGetField(opts, 0, "local_optimizer")) { + const mxArray* local_opts = mxGetField(opts, 0, "local_optimizer"); + CHECK1(mxIsStruct(local_opts), + "opt.local_optimizer must be a structure"); + CHECK1(local_opt = make_opt(local_opts, n), + "error initializing local optimizer"); + nlopt_set_local_optimizer(opt, local_opt); + nlopt_destroy(local_opt); + local_opt = NULL; + } + + mxFree(tmp); + return opt; } + +#define CHECK(cond, msg) \ + if (!(cond)) { \ + mxFree(dh); \ + mxFree(dfc); \ + nlopt_destroy(opt); \ + mexErrMsgTxt(msg); \ + } + +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + unsigned n; + double *x, *x0, opt_f; + nlopt_result ret; + mxArray *x_mx, *mx; + user_function_data d, dpre, dmfc, dmh, *dfc = NULL, *dh = NULL; + nlopt_opt opt = NULL; + + CHECK(nrhs == 2 && nlhs <= 3, "wrong number of arguments"); + + /* options = prhs[0] */ + CHECK(mxIsStruct(prhs[0]), "opt must be a struct"); + + /* x0 = prhs[1] */ + CHECK(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) + && (mxGetM(prhs[1]) == 1 || mxGetN(prhs[1]) == 1), + "x must be real row or column vector"); + n = mxGetM(prhs[1]) * mxGetN(prhs[1]), + x0 = mxGetPr(prhs[1]); + + CHECK(opt = make_opt(prhs[0], n), "error initializing nlopt options"); + + d.neval = 0; + d.verbose = (int)struct_val_default(prhs[0], "verbose", 0); + d.opt = opt; + + /* function f = prhs[1] */ + mx = struct_funcval(prhs[0], "min_objective"); + if (!mx) + mx = struct_funcval(prhs[0], "max_objective"); + CHECK(mx, "either opt.min_objective or opt.max_objective must exist"); + if (mxIsChar(mx)) { + CHECK(mxGetString(mx, d.f, FLEN) == 0, + "error reading function name string (too long?)"); + d.nrhs = 1; + d.xrhs = 0; + } else { + d.prhs[0] = mx; + strcpy(d.f, "feval"); + d.nrhs = 2; + d.xrhs = 1; + } + d.prhs[d.xrhs] = mxCreateDoubleMatrix(1, n, mxREAL); + + if ((mx = struct_funcval(prhs[0], "pre"))) { + CHECK(mxIsChar(mx) || mxIsFunctionHandle(mx), + "pre must contain function handles or function names"); + if (mxIsChar(mx)) { + CHECK(mxGetString(mx, dpre.f, FLEN) == 0, + "error reading function name string (too long?)"); + dpre.nrhs = 2; + dpre.xrhs = 0; + } else { + dpre.prhs[0] = mx; + strcpy(dpre.f, "feval"); + dpre.nrhs = 3; + dpre.xrhs = 1; + } + dpre.verbose = d.verbose > 2; + dpre.opt = opt; + dpre.neval = 0; + dpre.prhs[dpre.xrhs] = d.prhs[d.xrhs]; + dpre.prhs[d.xrhs + 1] = mxCreateDoubleMatrix(1, n, mxREAL); + d.dpre = &dpre; + + if (struct_funcval(prhs[0], "min_objective")) + nlopt_set_precond_min_objective(opt, user_function, user_pre, &d); + else + nlopt_set_precond_max_objective(opt, user_function, user_pre, &d); + } else { + dpre.nrhs = 0; + if (struct_funcval(prhs[0], "min_objective")) + nlopt_set_min_objective(opt, user_function, &d); + else + nlopt_set_max_objective(opt, user_function, &d); + } + + if ((mx = mxGetField(prhs[0], 0, "mfc"))) { + int m; + double* mfc_tol; + + m = (int)struct_val_default(prhs[0], "mfc_count", -1); + mfc_tol = struct_arrval(prhs[0], "mfc_tol", m, NULL); + + CHECK(mxIsChar(mx) || mxIsFunctionHandle(mx), + "mfc must contain function handles or function names"); + if (mxIsChar(mx)) { + CHECK(mxGetString(mx, dmfc.f, FLEN) == 0, + "error reading function name string (too long?)"); + dmfc.nrhs = 1; + dmfc.xrhs = 0; + } else { + dmfc.prhs[0] = mx; + strcpy(dmfc.f, "feval"); + dmfc.nrhs = 2; + dmfc.xrhs = 1; + } + dmfc.verbose = d.verbose > 1; + dmfc.opt = opt; + dmfc.neval = 0; + dmfc.prhs[dmfc.xrhs] = d.prhs[d.xrhs]; + + CHECK(nlopt_add_inequality_mconstraint(opt, m, + user_mfunction, &dmfc, mfc_tol ? mfc_tol : 0) + > 0, + "nlopt error adding multiple inequality constraints"); + } + + if ((mx = mxGetField(prhs[0], 0, "mh"))) { + int m; + double* mh_tol; + + m = (int)struct_val_default(prhs[0], "mh_count", -1); + mh_tol = struct_arrval(prhs[0], "mh_tol", m, NULL); + + CHECK(mxIsChar(mx) || mxIsFunctionHandle(mx), + "mh must contain function handles or function names"); + if (mxIsChar(mx)) { + CHECK(mxGetString(mx, dmh.f, FLEN) == 0, + "error reading function name string (too long?)"); + dmh.nrhs = 1; + dmh.xrhs = 0; + } else { + dmh.prhs[0] = mx; + strcpy(dmh.f, "feval"); + dmh.nrhs = 2; + dmh.xrhs = 1; + } + dmh.verbose = d.verbose > 1; + dmh.opt = opt; + dmh.neval = 0; + dmh.prhs[dmh.xrhs] = d.prhs[d.xrhs]; + + CHECK(nlopt_add_equality_mconstraint(opt, m, + user_mfunction, &dmh, mh_tol ? mh_tol : 0) + > 0, + "nlopt error adding multiple equality constraints"); + } + + if ((mx = mxGetField(prhs[0], 0, "fc"))) { + int j, m; + double* fc_tol; + + CHECK(mxIsCell(mx), "fc must be a Cell array"); + m = mxGetM(mx) * mxGetN(mx); + dfc = (user_function_data*)mxCalloc(m, sizeof(user_function_data)); + fc_tol = struct_arrval(prhs[0], "fc_tol", m, NULL); + + for (j = 0; j < m; ++j) { + mxArray* fc = mxGetCell(mx, j); + CHECK(mxIsChar(fc) || mxIsFunctionHandle(fc), + "fc must contain function handles or function names"); + if (mxIsChar(fc)) { + CHECK(mxGetString(fc, dfc[j].f, FLEN) == 0, + "error reading function name string (too long?)"); + dfc[j].nrhs = 1; + dfc[j].xrhs = 0; + } else { + dfc[j].prhs[0] = fc; + strcpy(dfc[j].f, "feval"); + dfc[j].nrhs = 2; + dfc[j].xrhs = 1; + } + dfc[j].verbose = d.verbose > 1; + dfc[j].opt = opt; + dfc[j].neval = 0; + dfc[j].prhs[dfc[j].xrhs] = d.prhs[d.xrhs]; + CHECK(nlopt_add_inequality_constraint(opt, user_function, + dfc + j, + fc_tol ? fc_tol[j] : 0) + > 0, + "nlopt error adding inequality constraint"); + } + } + + if ((mx = mxGetField(prhs[0], 0, "h"))) { + int j, m; + double* h_tol; + + CHECK(mxIsCell(mx), "h must be a Cell array"); + m = mxGetM(mx) * mxGetN(mx); + dh = (user_function_data*)mxCalloc(m, sizeof(user_function_data)); + h_tol = struct_arrval(prhs[0], "h_tol", m, NULL); + + for (j = 0; j < m; ++j) { + mxArray* h = mxGetCell(mx, j); + CHECK(mxIsChar(h) || mxIsFunctionHandle(h), + "h must contain function handles or function names"); + if (mxIsChar(h)) { + CHECK(mxGetString(h, dh[j].f, FLEN) == 0, + "error reading function name string (too long?)"); + dh[j].nrhs = 1; + dh[j].xrhs = 0; + } else { + dh[j].prhs[0] = h; + strcpy(dh[j].f, "feval"); + dh[j].nrhs = 2; + dh[j].xrhs = 1; + } + dh[j].verbose = d.verbose > 1; + dh[j].opt = opt; + dh[j].neval = 0; + dh[j].prhs[dh[j].xrhs] = d.prhs[d.xrhs]; + CHECK(nlopt_add_equality_constraint(opt, user_function, + dh + j, + h_tol ? h_tol[j] : 0) + > 0, + "nlopt error adding equality constraint"); + } + } + + x_mx = mxCreateDoubleMatrix(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL); + x = mxGetPr(x_mx); + memcpy(x, x0, sizeof(double) * n); + + ret = nlopt_optimize(opt, x, &opt_f); + + mxFree(dh); + mxFree(dfc); + mxDestroyArray(d.prhs[d.xrhs]); + if (dpre.nrhs > 0) + mxDestroyArray(dpre.prhs[d.xrhs + 1]); + nlopt_destroy(opt); + + plhs[0] = x_mx; + if (nlhs > 1) { + plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); + *(mxGetPr(plhs[1])) = opt_f; + } + if (nlhs > 2) { + plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL); + *(mxGetPr(plhs[2])) = (int)ret; + } +} \ No newline at end of file From 30fd78f1603641bec530120c6c167737dc37d7eb Mon Sep 17 00:00:00 2001 From: Nicola Piccinelli Date: Thu, 23 Jan 2020 11:58:11 +0100 Subject: [PATCH 2/3] isNaN check performed on all elements of result vector --- src/octave/nlopt_optimize-mex.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/octave/nlopt_optimize-mex.c b/src/octave/nlopt_optimize-mex.c index 2b094e37..45f4bd23 100644 --- a/src/octave/nlopt_optimize-mex.c +++ b/src/octave/nlopt_optimize-mex.c @@ -161,7 +161,7 @@ static void user_mfunction(unsigned m, double* result, unsigned n, const double* mexPrintf("%g\n", result[n - 1]); } - for (size_t i = 0; i < n - 1; i++) { + for (size_t i = 0; i < n; i++) { if (mxIsNaN(result[i])) { nlopt_force_stop(d->opt); break; @@ -492,4 +492,4 @@ void mexFunction(int nlhs, mxArray* plhs[], plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL); *(mxGetPr(plhs[2])) = (int)ret; } -} \ No newline at end of file +} From c8e51d9ed12681992ab232c653584117208f3d2a Mon Sep 17 00:00:00 2001 From: Nicola Piccinelli Date: Mon, 27 Jan 2020 12:25:47 +0100 Subject: [PATCH 3/3] Fixed wrong gradient matrix assignment --- src/octave/nlopt_optimize-mex.c | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/octave/nlopt_optimize-mex.c b/src/octave/nlopt_optimize-mex.c index 45f4bd23..8fb6708b 100644 --- a/src/octave/nlopt_optimize-mex.c +++ b/src/octave/nlopt_optimize-mex.c @@ -144,10 +144,14 @@ static void user_mfunction(unsigned m, double* result, unsigned n, const double* if (gradient) { CHECK0(mxIsDouble(d->plhs[1]) && !mxIsComplex(d->plhs[1]) - && (mxGetM(d->plhs[1]) == m || mxGetN(d->plhs[1]) == m) - && mxGetM(d->plhs[1]) * mxGetN(d->plhs[1]) == m * n, - "gradient vector from user mfunction is the wrong size"); - memcpy(gradient, mxGetPr(d->plhs[1]), m * n * sizeof(double)); + && (mxGetM(d->plhs[1]) == m && mxGetN(d->plhs[1]) == n), + "gradient vector from user mfunction is the wrong size (mxn)"); + double* ptr = mxGetPr(d->plhs[1]); + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < m; i++) { + gradient[(i * n) + j] = ptr[i + (j * m)]; + } + } mxDestroyArray(d->plhs[1]); }