1 /* Copyright (c) 2007-2010 Massachusetts Institute of Technology
3 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including
6 * without limitation the rights to use, copy, modify, merge, publish,
7 * distribute, sublicense, and/or sell copies of the Software, and to
8 * permit persons to whom the Software is furnished to do so, subject to
9 * the following conditions:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 /* Matlab MEX interface to NLopt, and in particular to nlopt_optimize */
33 #define CHECK0(cond, msg) if (!(cond)) mexErrMsgTxt(msg);
35 static double struct_val_default(const mxArray *s, const char *name, double dflt)
37 mxArray *val = mxGetField(s, 0, name);
39 CHECK0(mxIsNumeric(val) && !mxIsComplex(val)
40 && mxGetM(val) * mxGetN(val) == 1,
41 "opt fields, other than xtol_abs, must be real scalars");
42 return mxGetScalar(val);
47 static double *struct_arrval(const mxArray *s, const char *name, unsigned n,
50 mxArray *val = mxGetField(s, 0, name);
52 CHECK0(mxIsNumeric(val) && !mxIsComplex(val)
53 && mxGetM(val) * mxGetN(val) == n,
54 "opt vector field is not of length n");
60 static mxArray *struct_funcval(const mxArray *s, const char *name)
62 mxArray *val = mxGetField(s, 0, name);
64 CHECK0(mxIsChar(val) || mxIsFunctionHandle(val),
65 "opt function field is not a function handle/name");
71 static double *fill(double *arr, unsigned n, double val)
74 for (i = 0; i < n; ++i) arr[i] = val;
78 #define FLEN 128 /* max length of user function name */
79 #define MAXRHS 2 /* max nrhs for user function */
83 mxArray *prhs[MAXRHS];
88 static double user_function(unsigned n, const double *x,
89 double *gradient, /* NULL if not needed */
92 user_function_data *d = (user_function_data *) d_;
95 d->plhs[0] = d->plhs[1] = NULL;
96 memcpy(mxGetPr(d->prhs[d->xrhs]), x, n * sizeof(double));
98 CHECK0(0 == mexCallMATLAB(gradient ? 2 : 1, d->plhs,
99 d->nrhs, d->prhs, d->f),
100 "error calling user function");
102 CHECK0(mxIsNumeric(d->plhs[0]) && !mxIsComplex(d->plhs[0])
103 && mxGetM(d->plhs[0]) * mxGetN(d->plhs[0]) == 1,
104 "user function must return real scalar");
105 f = mxGetScalar(d->plhs[0]);
106 mxDestroyArray(d->plhs[0]);
108 CHECK0(mxIsDouble(d->plhs[1]) && !mxIsComplex(d->plhs[1])
109 && (mxGetM(d->plhs[1]) == 1 || mxGetN(d->plhs[1]) == 1)
110 && mxGetM(d->plhs[1]) * mxGetN(d->plhs[1]) == n,
111 "gradient vector from user function is the wrong size");
112 memcpy(gradient, mxGetPr(d->plhs[1]), n * sizeof(double));
113 mxDestroyArray(d->plhs[1]);
116 if (d->verbose) mexPrintf("nlopt_optimize eval #%d: %g\n", d->neval, f);
120 #define CHECK1(cond, msg) if (!(cond)) { mxFree(tmp); nlopt_destroy(opt); nlopt_destroy(local_opt); mexWarnMsgTxt(msg); return NULL; };
122 nlopt_opt make_opt(const mxArray *opts, unsigned n)
124 nlopt_opt opt = NULL, local_opt = NULL;
125 nlopt_algorithm algorithm;
129 algorithm = (nlopt_algorithm)
130 struct_val_default(opts, "algorithm", NLOPT_NUM_ALGORITHMS);
131 CHECK1(((int)algorithm) >= 0 && algorithm < NLOPT_NUM_ALGORITHMS,
132 "invalid opt.algorithm");
134 tmp = (double *) mxCalloc(n, sizeof(double));
135 opt = nlopt_create(algorithm, n);
136 CHECK1(opt, "nlopt: out of memory");
138 nlopt_set_lower_bounds(opt, struct_arrval(opts, "lower_bounds", n,
139 fill(tmp, n, -HUGE_VAL)));
140 nlopt_set_upper_bounds(opt, struct_arrval(opts, "upper_bounds", n,
141 fill(tmp, n, +HUGE_VAL)));
143 nlopt_set_stopval(opt, struct_val_default(opts, "stopval", -HUGE_VAL));
144 nlopt_set_ftol_rel(opt, struct_val_default(opts, "ftol_rel", 0.0));
145 nlopt_set_ftol_abs(opt, struct_val_default(opts, "ftol_abs", 0.0));
146 nlopt_set_xtol_rel(opt, struct_val_default(opts, "xtol_rel", 0.0));
147 nlopt_set_xtol_abs(opt, struct_arrval(opts, "xtol_abs", n,
149 nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0.0) < 0 ?
150 0 : struct_val_default(opts, "maxeval", 0.0));
151 nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0));
153 nlopt_set_population(opt, struct_val_default(opts, "population", 0));
155 if (struct_arrval(opts, "initial_step", n, NULL))
156 nlopt_set_initial_step(opt,
157 struct_arrval(opts, "initial_step", n, NULL));
159 if (mxGetField(opts, 0, "local_optimizer")) {
160 const mxArray *local_opts = mxGetField(opts, 0, "local_optimizer");
161 CHECK1(mxIsStruct(local_opts),
162 "opt.local_optimizer must be a structure");
163 CHECK1(local_opt = make_opt(local_opts, n),
164 "error initializing local optimizer");
165 nlopt_set_local_optimizer(opt, local_opt);
166 nlopt_destroy(local_opt); local_opt = NULL;
173 #define CHECK(cond, msg) if (!(cond)) { mxFree(dh); mxFree(dfc); nlopt_destroy(opt); mexErrMsgTxt(msg); }
175 void mexFunction(int nlhs, mxArray *plhs[],
176 int nrhs, const mxArray *prhs[])
179 double *x, *x0, opt_f;
182 user_function_data d, *dfc = NULL, *dh = NULL;
183 nlopt_opt opt = NULL;
185 CHECK(nrhs == 2 && nlhs <= 3, "wrong number of arguments");
187 /* options = prhs[0] */
188 CHECK(mxIsStruct(prhs[0]), "opt must be a struct");
191 CHECK(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1])
192 && (mxGetM(prhs[1]) == 1 || mxGetN(prhs[1]) == 1),
193 "x must be real row or column vector");
194 n = mxGetM(prhs[1]) * mxGetN(prhs[1]),
195 x0 = mxGetPr(prhs[1]);
197 CHECK(opt = make_opt(prhs[0], n), "error initializing nlopt options");
200 d.verbose = (int) struct_val_default(prhs[0], "verbose", 0);
202 /* function f = prhs[1] */
203 mx = struct_funcval(prhs[0], "min_objective");
204 if (!mx) mx = struct_funcval(prhs[0], "max_objective");
205 CHECK(mx, "either opt.min_objective or opt.max_objective must exist");
207 CHECK(mxGetString(mx, d.f, FLEN) == 0,
208 "error reading function name string (too long?)");
214 strcpy(d.f, "feval");
218 d.prhs[d.xrhs] = mxCreateDoubleMatrix(1, n, mxREAL);
219 if (struct_funcval(prhs[0], "min_objective"))
220 nlopt_set_min_objective(opt, user_function, &d);
222 nlopt_set_max_objective(opt, user_function, &d);
225 if ((mx = mxGetField(prhs[0], 0, "fc"))) {
229 CHECK(mxIsCell(mx), "fc must be a Cell array");
230 m = mxGetM(mx) * mxGetN(mx);;
231 dfc = (user_function_data *) mxCalloc(m, sizeof(user_function_data));
232 fc_tol = struct_arrval(prhs[0], "fc_tol", m, NULL);
234 for (j = 0; j < m; ++j) {
235 mxArray *fc = mxGetCell(mx, j);
236 CHECK(mxIsChar(fc) || mxIsFunctionHandle(fc),
237 "fc must contain function handles or function names");
239 CHECK(mxGetString(fc, dfc[j].f, FLEN) == 0,
240 "error reading function name string (too long?)");
246 strcpy(dfc[j].f, "feval");
250 dfc[j].verbose = d.verbose > 1;
252 dfc[j].prhs[dfc[j].xrhs] = d.prhs[d.xrhs];
253 CHECK(nlopt_add_inequality_constraint(opt, user_function,
255 fc_tol ? fc_tol[j] : 0)
256 > 0, "nlopt error adding inequality constraint");
261 if ((mx = mxGetField(prhs[0], 0, "h"))) {
265 CHECK(mxIsCell(mx), "h must be a Cell array");
266 m = mxGetM(mx) * mxGetN(mx);;
267 dh = (user_function_data *) mxCalloc(m, sizeof(user_function_data));
268 h_tol = struct_arrval(prhs[0], "h_tol", m, NULL);
270 for (j = 0; j < m; ++j) {
271 mxArray *h = mxGetCell(mx, j);
272 CHECK(mxIsChar(h) || mxIsFunctionHandle(h),
273 "h must contain function handles or function names");
275 CHECK(mxGetString(h, dh[j].f, FLEN) == 0,
276 "error reading function name string (too long?)");
282 strcpy(dh[j].f, "feval");
286 dh[j].verbose = d.verbose > 1;
288 dh[j].prhs[dh[j].xrhs] = d.prhs[d.xrhs];
289 CHECK(nlopt_add_equality_constraint(opt, user_function,
291 h_tol ? h_tol[j] : 0)
292 > 0, "nlopt error adding equality constraint");
297 x_mx = mxCreateDoubleMatrix(mxGetM(prhs[1]), mxGetN(prhs[1]), mxREAL);
299 memcpy(x, x0, sizeof(double) * n);
301 ret = nlopt_optimize(opt, x, &opt_f);
305 mxDestroyArray(d.prhs[d.xrhs]);
310 plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
311 *(mxGetPr(plhs[1])) = opt_f;
314 plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL);
315 *(mxGetPr(plhs[2])) = (int) ret;