From: stevenj Date: Wed, 12 May 2010 04:00:09 +0000 (-0400) Subject: stab at Matlab nlopt_optimize interface X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=0ddaa97e5a9322558955a237058eaa5cb9185b78;p=nlopt.git stab at Matlab nlopt_optimize interface darcs-hash:20100512040009-c8de0-edd99ecc6904eca7c0fbb4d42a263b0ddb8cd00e.gz --- diff --git a/octave/Makefile.am b/octave/Makefile.am index c03c40c..b2e5cf7 100644 --- a/octave/Makefile.am +++ b/octave/Makefile.am @@ -16,14 +16,6 @@ oct_DATA = nlopt_optimize.oct m_DATA = $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m nlopt_optimize.m endif -nlopt_minimize_constrained.oct: nlopt_minimize_constrained-oct.cc nlopt_minimize_constrained_usage.h $(top_builddir)/.libs/libnlopt@NLOPT_SUFFIX@.la - $(MKOCTFILE) -o $@ $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(srcdir)/nlopt_minimize_constrained-oct.cc $(LDFLAGS) -L$(top_builddir)/.libs -lnlopt@NLOPT_SUFFIX@ - -nlopt_minimize_constrained_usage.h: $(srcdir)/nlopt_minimize_constrained.m - echo "#define NLOPT_MINIMIZE_CONSTRAINED_USAGE \\" > $@ - sed 's/\"/\\"/g' $(srcdir)/nlopt_minimize_constrained.m | sed 's,^% ,\",;s,^%,\",;s,$$,\\n\" \\,' >> $@ - echo "" >> $@ - nlopt_optimize.oct: nlopt_optimize-oct.cc nlopt_optimize_usage.h $(top_builddir)/.libs/libnlopt@NLOPT_SUFFIX@.la $(MKOCTFILE) -o $@ $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(srcdir)/nlopt_optimize-oct.cc $(LDFLAGS) -L$(top_builddir)/.libs -lnlopt@NLOPT_SUFFIX@ @@ -35,14 +27,14 @@ nlopt_optimize_usage.h: $(srcdir)/nlopt_optimize.m ####################################################################### mexdir = $(MEX_INSTALL_DIR) if WITH_MATLAB -mex_DATA = $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m +mex_DATA = nlopt_optimize.$(MEXSUFF) $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m nlopt_optimize.m endif -nlopt_minimize_constrained.$(MEXSUFF): nlopt_minimize_constrained-mex.c $(top_builddir)/.libs/libnlopt@NLOPT_SUFFIX@.la - $(MEX) -output nlopt_minimize_constrained -O $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(srcdir)/nlopt_minimize_constrained-mex.c $(LDFLAGS) -L$(top_builddir)/.libs -lnlopt@NLOPT_SUFFIX@ +nlopt_optimize.$(MEXSUFF): nlopt_optimize-mex.c $(top_builddir)/.libs/libnlopt@NLOPT_SUFFIX@.la + $(MEX) -output nlopt_optimize -O $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(srcdir)/nlopt_optimize-mex.c $(LDFLAGS) -L$(top_builddir)/.libs -lnlopt@NLOPT_SUFFIX@ ####################################################################### -EXTRA_DIST = nlopt_minimize_constrained-oct.cc nlopt_optimize-oct.cc nlopt_minimize_constrained-mex.c $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m nlopt_optimize.m +EXTRA_DIST = nlopt_optimize-oct.cc nlopt_optimize-mex.c $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m nlopt_optimize.m -CLEANFILES = nlopt_minimize_constrained.oct nlopt_minimize_constrained_usage.h nlopt_optimize.oct nlopt_optimize_usage.h nlopt_minimize_constrained.$(MEXSUFF) nlopt_minimize_constrained-oct.o nlopt_optimize-oct.o +CLEANFILES = nlopt_optimize.oct nlopt_optimize.h nlopt_optimize.$(MEXSUFF) nlopt_optimize-oct.o diff --git a/octave/nlopt_minimize_constrained-mex.c b/octave/nlopt_minimize_constrained-mex.c deleted file mode 100644 index 504739f..0000000 --- a/octave/nlopt_minimize_constrained-mex.c +++ /dev/null @@ -1,253 +0,0 @@ -/* Copyright (c) 2007-2010 Massachusetts Institute of Technology - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE - * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION - * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* Matlab MEX interface to NLopt, and in particular to nlopt_minimize_constrained */ - -#include -#include -#include -#include -#include - -#include "nlopt.h" - -#define xSTRIZE(x) #x -#define STRIZE(x) xSTRIZE(x) -#define CHECK(cond, msg) if (!(cond)) mexErrMsgTxt(msg); - -static double struct_val_default(const mxArray *s, const char *name, double dflt) -{ - mxArray *val = mxGetField(s, 0, name); - if (val) { - CHECK(mxIsNumeric(val) && !mxIsComplex(val) - && mxGetM(val) * mxGetN(val) == 1, - "stop fields, other than xtol_abs, must be real scalars"); - return mxGetScalar(val); - } - return dflt; -} - -#define FLEN 128 /* max length of user function name */ -#define MAXRHS 128 /* max nrhs for user function */ -typedef struct { - char f[FLEN]; - mxArray *plhs[2]; - mxArray *prhs[MAXRHS]; - int xrhs, nrhs; - int verbose, neval; -} user_function_data; - -static double user_function(int 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)); - - CHECK(0 == mexCallMATLAB(gradient ? 2 : 1, d->plhs, - d->nrhs, d->prhs, d->f), - "error calling user function"); - - CHECK(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) { - CHECK(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_minimize_constrained eval #%d: %g\n", d->neval, f); - return f; -} - -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) -{ - nlopt_algorithm algorithm; - int n, m, i, j; - double *lb, *ub, *x, *x0; - double minf_max, ftol_rel, ftol_abs, xtol_rel, *xtol_abs, maxtime; - int maxeval; - nlopt_result ret; - mxArray *x_mx; - double minf = HUGE_VAL; - user_function_data d, *dc; - - CHECK(nrhs == 9 && nlhs <= 3, "wrong number of arguments"); - - /* algorithm = prhs[0] */ - CHECK(mxIsNumeric(prhs[0]) && !mxIsComplex(prhs[0]) - && mxGetM(prhs[0]) * mxGetN(prhs[0]) == 1, - "algorithm must be real (integer) scalar"); - algorithm = (nlopt_algorithm) (mxGetScalar(prhs[0]) + 0.5); - CHECK(algorithm >= 0 && algorithm < NLOPT_NUM_ALGORITHMS, - "unknown algorithm"); - - /* function f = prhs[1] */ - CHECK(mxIsChar(prhs[1]) || mxIsFunctionHandle(prhs[1]), - "f must be a function handle or function name"); - if (mxIsChar(prhs[1])) { - CHECK(mxGetString(prhs[1], d.f, FLEN) == 0, - "error reading function name string (too long?)"); - d.nrhs = 1; - d.xrhs = 0; - } - else { - d.prhs[0] = (mxArray *) prhs[1]; - strcpy(d.f, "feval"); - d.nrhs = 2; - d.xrhs = 1; - } - - /* Cell f_data = prhs[2] */ - CHECK(mxIsCell(prhs[2]), "f_data must be a Cell array"); - CHECK(mxGetM(prhs[2]) * mxGetN(prhs[2]) + 1 <= MAXRHS, - "user function cannot have more than " STRIZE(MAXRHS) " arguments"); - d.nrhs += mxGetM(prhs[2]) * mxGetN(prhs[2]); - for (i = 0; i < d.nrhs - (1+d.xrhs); ++i) - d.prhs[(1+d.xrhs)+i] = mxGetCell(prhs[2], i); - - /* m = length(fc = prhs[3]) = length(fc_data = prhs[4]) */ - CHECK(mxIsCell(prhs[3]), "fc must be a Cell array"); - CHECK(mxIsCell(prhs[4]), "fc_data must be a Cell array"); - m = mxGetM(prhs[3]) * mxGetN(prhs[3]); - CHECK(m == mxGetM(prhs[4]) * mxGetN(prhs[4]), "fc and fc_data must have the same length"); - dc = (user_function_data *) malloc(sizeof(user_function_data) * m); - - for (j = 0; j < m; ++j) { - mxArray *fc, *fc_data; - - /* function fc = phrs[3] */ - fc = mxGetCell(prhs[3], j); - CHECK(mxIsChar(fc) || mxIsFunctionHandle(fc), - "fc must be Cell array of function handles or function names"); - if (mxIsChar(fc)) { - CHECK(mxGetString(fc, dc[j].f, FLEN) == 0, - "error reading function name string (too long?)"); - dc[j].nrhs = 1; - dc[j].xrhs = 0; - } - else { - dc[j].prhs[0] = fc; - strcpy(dc[j].f, "feval"); - dc[j].nrhs = 2; - dc[j].xrhs = 1; - } - - /* Cell fc_data = prhs[4] */ - fc_data = mxGetCell(prhs[4], j); - CHECK(mxIsCell(fc_data), "fc_data must be a Cell array of Cell arrays"); - CHECK(mxGetM(fc_data) * mxGetN(fc_data) + 1 <= MAXRHS, - "user function cannot have more than " STRIZE(MAXRHS) " arguments"); - dc[j].nrhs += mxGetM(fc_data) * mxGetN(fc_data); - for (i = 0; i < dc[j].nrhs - (1+dc[j].xrhs); ++i) - dc[j].prhs[(1+dc[j].xrhs)+i] = mxGetCell(fc_data, i); - } - - /* lb = prhs[5] */ - CHECK(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) - && (mxGetM(prhs[5]) == 1 || mxGetN(prhs[5]) == 1), - "lb must be real row or column vector"); - lb = mxGetPr(prhs[5]); - n = mxGetM(prhs[5]) * mxGetN(prhs[5]); - - /* ub = prhs[6] */ - CHECK(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) - && (mxGetM(prhs[6]) == 1 || mxGetN(prhs[6]) == 1) - && mxGetM(prhs[6]) * mxGetN(prhs[6]) == n, - "ub must be real row or column vector of same length as lb"); - ub = mxGetPr(prhs[6]); - - /* x0 = prhs[7] */ - CHECK(mxIsDouble(prhs[7]) && !mxIsComplex(prhs[7]) - && (mxGetM(prhs[7]) == 1 || mxGetN(prhs[7]) == 1) - && mxGetM(prhs[7]) * mxGetN(prhs[7]) == n, - "x must be real row or column vector of same length as lb"); - x0 = mxGetPr(prhs[7]); - - /* stopping criteria = prhs[8] */ - CHECK(mxIsStruct(prhs[8]), "stopping criteria must be a struct"); - minf_max = struct_val_default(prhs[8], "minf_max", -HUGE_VAL); - ftol_rel = struct_val_default(prhs[8], "ftol_rel", 0); - ftol_abs = struct_val_default(prhs[8], "ftol_abs", 0); - xtol_rel = struct_val_default(prhs[8], "xtol_rel", 0); - maxeval = (int) (struct_val_default(prhs[8], "maxeval", -1) + 0.5); - maxtime = struct_val_default(prhs[8], "maxtime", -1); - d.verbose = (int) struct_val_default(prhs[8], "verbose", 0); - d.neval = 0; - for (i = 0; i < m; ++i) { - dc[i].verbose = d.verbose > 1; - dc[i].neval = 0; - } - { - mxArray *val = mxGetField(prhs[8], 0, "xtol_abs"); - if (val) { - CHECK(mxIsNumeric(val) && !mxIsComplex(val) - && (mxGetM(val) == 1 || mxGetN(val) == 1) - && mxGetM(val) * mxGetN(val) == n, - "stop.xtol_abs must be real row/col vector of length n"); - xtol_abs = mxGetPr(val); - } - else - xtol_abs = NULL; - } - - - x_mx = mxCreateDoubleMatrix(1, n, mxREAL); - x = mxGetPr(x_mx); - memcpy(x, x0, sizeof(double) * n); - - d.prhs[d.xrhs] = mxCreateDoubleMatrix(1, n, mxREAL); - for (i = 0; i < m;++i) - dc[i].prhs[dc[i].xrhs] = d.prhs[d.xrhs]; - - ret = nlopt_minimize_constrained(algorithm, - n, - user_function, &d, - m, user_function, dc, - sizeof(user_function_data), - lb, ub, x, &minf, minf_max, - ftol_rel, ftol_abs, xtol_rel, xtol_abs, - maxeval, maxtime); - - mxDestroyArray(d.prhs[d.xrhs]); - free(dc); - - plhs[0] = x_mx; - if (nlhs > 1) { - plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); - *(mxGetPr(plhs[1])) = minf; - } - if (nlhs > 2) { - plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL); - *(mxGetPr(plhs[2])) = (int) ret; - } -} diff --git a/octave/nlopt_minimize_constrained-oct.cc b/octave/nlopt_minimize_constrained-oct.cc deleted file mode 100644 index 9d2dfcd..0000000 --- a/octave/nlopt_minimize_constrained-oct.cc +++ /dev/null @@ -1,190 +0,0 @@ -/* Copyright (c) 2007-2010 Massachusetts Institute of Technology - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE - * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION - * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -#include -#include -#include -#include -#include - -#include "nlopt.h" -#include "nlopt_minimize_constrained_usage.h" - -static double struct_val_default(Octave_map &m, const std::string& k, - double dflt) -{ - if (m.contains(k)) { - if (m.contents(k).length() == 1 && (m.contents(k))(0).is_real_scalar()) - return (m.contents(k))(0).double_value(); - } - return dflt; -} - -static Matrix struct_val_default(Octave_map &m, const std::string& k, - Matrix &dflt) -{ - if (m.contains(k)) { - if ((m.contents(k)).length() == 1) { - if ((m.contents(k))(0).is_real_scalar()) - return Matrix(1, dflt.length(), (m.contents(k))(0).double_value()); - else if ((m.contents(k))(0).is_real_matrix()) - return (m.contents(k))(0).matrix_value(); - } - } - return dflt; -} - -typedef struct { - octave_function *f; - Cell f_data; - int neval, verbose; -} user_function_data; - -static double user_function(int n, const double *x, - double *gradient, /* NULL if not needed */ - void *data_) -{ - user_function_data *data = (user_function_data *) data_; - octave_value_list args(1 + data->f_data.length(), 0); - Matrix xm(1,n); - for (int i = 0; i < n; ++i) - xm(i) = x[i]; - args(0) = xm; - for (int i = 0; i < data->f_data.length(); ++i) - args(1 + i) = data->f_data(i); - octave_value_list res = data->f->do_multi_index_op(gradient ? 2 : 1, args); - if (res.length() < (gradient ? 2 : 1)) - gripe_user_supplied_eval("nlopt_minimize_constrained"); - else if (!res(0).is_real_scalar() - || (gradient && !res(1).is_real_matrix() - && !(n == 1 && res(1).is_real_scalar()))) - gripe_user_returned_invalid("nlopt_minimize_constrained"); - else { - if (gradient) { - if (n == 1 && res(1).is_real_scalar()) - gradient[0] = res(1).double_value(); - else { - Matrix grad = res(1).matrix_value(); - for (int i = 0; i < n; ++i) - gradient[i] = grad(i); - } - } - data->neval++; - if (data->verbose) printf("nlopt_minimize_constrained eval #%d: %g\n", - data->neval, res(0).double_value()); - return res(0).double_value(); - } - return 0; -} - -#define CHECK(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); print_usage("nlopt_minimize_constrained"); return retval; } - -DEFUN_DLD(nlopt_minimize_constrained, args, nargout, NLOPT_MINIMIZE_CONSTRAINED_USAGE) -{ - octave_value_list retval; - double A; - - CHECK(args.length() == 9 && nargout <= 3, "wrong number of args"); - - CHECK(args(0).is_real_scalar(), "algorithm must be integer"); - nlopt_algorithm algorithm = nlopt_algorithm(args(0).int_value()); - - user_function_data d; - CHECK(args(1).is_function() || args(1).is_function_handle(), - "f must be a function handle"); - d.f = args(1).function_value(); - CHECK(args(2).is_cell(), "f_data must be cell array"); - d.f_data = args(2).cell_value(); - - CHECK(args(3).is_cell(), "fc must be cell array"); - CHECK(args(4).is_cell(), "fc_data must be cell array"); - int m = args(3).length(); - CHECK(m == args(4).length(), "fc and fc_data must have the same length"); - user_function_data *dc = new user_function_data[m+1]; - Cell fc = args(3).cell_value(); - Cell fc_data = args(4).cell_value(); - for (int i = 0; i < m; ++i) { - CHECK(fc(i).is_function() || fc(i).is_function_handle(), - "fc must be a cell array of function handles"); - dc[i].f = fc(i).function_value(); - CHECK(fc_data(i).is_cell(), "fc_data must be cell array of cell arrays"); - dc[i].f_data = fc_data(i).cell_value(); - } - - CHECK(args(5).is_real_matrix() || args(5).is_real_scalar(), - "lb must be real vector"); - Matrix lb = args(5).is_real_scalar() ? - Matrix(1, 1, args(5).double_value()) : args(5).matrix_value(); - int n = lb.length(); - - CHECK(args(6).is_real_matrix() || args(6).is_real_scalar(), - "ub must be real vector"); - Matrix ub = args(6).is_real_scalar() ? - Matrix(1, 1, args(6).double_value()) : args(6).matrix_value(); - CHECK(n == ub.length(), "lb and ub must have same length"); - - CHECK(args(7).is_real_matrix() || args(7).is_real_scalar(), - "x must be real vector"); - Matrix x = args(7).is_real_scalar() ? - Matrix(1, 1, args(7).double_value()) : args(7).matrix_value(); - CHECK(n == x.length(), "x and lb/ub must have same length"); - - CHECK(args(8).is_map(), "stop must be structure"); - Octave_map stop = args(8).map_value(); - double minf_max = struct_val_default(stop, "minf_max", -HUGE_VAL); - double ftol_rel = struct_val_default(stop, "ftol_rel", 0); - double ftol_abs = struct_val_default(stop, "ftol_abs", 0); - double xtol_rel = struct_val_default(stop, "xtol_rel", 0); - Matrix zeros(1, n, 0.0); - Matrix xtol_abs = struct_val_default(stop, "xtol_abs", zeros); - CHECK(n == xtol_abs.length(), "stop.xtol_abs must have same length as x"); - int maxeval = int(struct_val_default(stop, "maxeval", -1)); - double maxtime = struct_val_default(stop, "maxtime", -1); - - d.neval = 0; - d.verbose = (int) struct_val_default(stop, "verbose", 0); - for (int i = 0; i < m; ++i) { - dc[i].neval = 0; - dc[i].verbose = d.verbose > 1; - } - - double minf = HUGE_VAL; - nlopt_result ret = nlopt_minimize_constrained(algorithm, - n, user_function, &d, - m, user_function, dc, - sizeof(user_function_data), - lb.data(), ub.data(), - x.fortran_vec(), &minf, - minf_max, ftol_rel, ftol_abs, - xtol_rel, xtol_abs.data(), - maxeval, maxtime); - - retval(0) = x; - if (nargout > 1) - retval(1) = minf; - if (nargout > 2) - retval(2) = int(ret); - - delete[] dc; - - return retval; -} diff --git a/octave/nlopt_optimize-mex.c b/octave/nlopt_optimize-mex.c new file mode 100644 index 0000000..4956192 --- /dev/null +++ b/octave/nlopt_optimize-mex.c @@ -0,0 +1,317 @@ +/* Copyright (c) 2007-2010 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* Matlab MEX interface to NLopt, and in particular to nlopt_optimize */ + +#include +#include +#include +#include +#include + +#include "nlopt.h" + +#define CHECK0(cond, msg) if (!(cond)) mexErrMsgTxt(msg); + +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; +} + +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; +} + +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; +} + +static double *fill(double *arr, unsigned n, double val) +{ + unsigned i; + for (i = 0; i < n; ++i) arr[i] = val; + return arr; +} + +#define FLEN 128 /* max length of user function name */ +#define MAXRHS 2 /* max nrhs for user function */ +typedef struct { + char f[FLEN]; + mxArray *plhs[2]; + mxArray *prhs[MAXRHS]; + int xrhs, nrhs; + int verbose, neval; +} user_function_data; + +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); + return f; +} + +#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) +{ + 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, "lb", n, + fill(tmp, n, -HUGE_VAL))); + nlopt_set_upper_bounds(opt, struct_arrval(opts, "ub", 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_xtol_rel(opt, struct_val_default(opts, "maxeval", 0.0) < 0 ? + 0 : struct_val_default(opts, "maxeval", 0.0)); + nlopt_set_xtol_rel(opt, struct_val_default(opts, "maxtime", 0.0)); + + nlopt_set_population(opt, struct_val_default(opts, "population", 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, *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); + + /* 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 (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].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].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]); + 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; + } +} diff --git a/octave/nlopt_optimize-oct.cc b/octave/nlopt_optimize-oct.cc index 3ec290c..f95f7ec 100644 --- a/octave/nlopt_optimize-oct.cc +++ b/octave/nlopt_optimize-oct.cc @@ -29,10 +29,6 @@ #include "nlopt.h" #include "nlopt_optimize_usage.h" -static nlopt_func struct_val_func(Octave_map &m, const std::string& k) -{ -} - static int struct_val_default(Octave_map &m, const std::string& k, int dflt) { diff --git a/octave/nlopt_optimize.m b/octave/nlopt_optimize.m index 4230e22..6d0d85b 100644 --- a/octave/nlopt_optimize.m +++ b/octave/nlopt_optimize.m @@ -151,9 +151,9 @@ % optimization algorithm as a subroutine, typically for local optimization. % By default, they use MMA or COBYLA for gradient-based or derivative-free % searching, respectively. However, you can change this by specifying -% opt.local_opt: this is a structure with the same types of fields as opt +% opt.local_optimizer: this is a structure with the same types of fields as opt % (stopping criteria, algorithm, etcetera). The objective function -% and nonlinear constraint parameters of opt.local_opt are ignored. +% and nonlinear constraint parameters of opt.local_optimizer are ignored. % % INITIAL STEP SIZE: %