1 /* Copyright (c) 2007-2014 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 #include <octave/oct.h>
24 #include <octave/oct-map.h>
25 #include <octave/ov.h>
26 #include <octave/parse.h>
31 #include "nlopt_optimize_usage.h"
33 #include <octave/version.h>
34 #if OCTAVE_MAJOR_VERSION < 3 || (OCTAVE_MAJOR_VERSION == 3 && OCTAVE_MINOR_VERSION < 8)
35 # define octave_map Octave_map
38 static int struct_val_default(octave_map &m, const std::string& k,
42 if (m.contents(k).numel() == 1 && (m.contents(k))(0).is_real_scalar())
43 return (m.contents(k))(0).int_value();
48 static double struct_val_default(octave_map &m, const std::string& k,
52 if (m.contents(k).numel() == 1 && (m.contents(k))(0).is_real_scalar())
53 return (m.contents(k))(0).double_value();
58 static Matrix struct_val_default(octave_map &m, const std::string& k,
62 if ((m.contents(k)).numel() == 1) {
63 if ((m.contents(k))(0).is_real_scalar())
64 return Matrix(1, dflt.numel(), (m.contents(k))(0).double_value());
65 else if ((m.contents(k))(0).is_real_matrix())
66 return (m.contents(k))(0).matrix_value();
78 static double user_function(unsigned n, const double *x,
79 double *gradient, /* NULL if not needed */
82 user_function_data *data = (user_function_data *) data_;
83 octave_value_list args(1, 0);
85 for (unsigned i = 0; i < n; ++i)
89 #if (OCTAVE_MAJOR_VERSION == 4 && OCTAVE_MINOR_VERSION > 2)
90 = octave::feval(data->f, args, gradient ? 2 : 1);
92 = data->f->do_multi_index_op(gradient ? 2 : 1, args);
94 if (res.length() < (gradient ? 2 : 1))
95 err_user_supplied_eval("nlopt_optimize");
96 else if (!res(0).is_real_scalar()
97 || (gradient && !res(1).is_real_matrix()
98 && !(n == 1 && res(1).is_real_scalar())))
99 err_user_returned_invalid("nlopt_optimize");
102 if (n == 1 && res(1).is_real_scalar())
103 gradient[0] = res(1).double_value();
105 Matrix grad = res(1).matrix_value();
106 for (unsigned i = 0; i < n; ++i)
107 gradient[i] = grad(i);
111 if (data->verbose) printf("nlopt_optimize eval #%d: %g\n",
112 data->neval, res(0).double_value());
113 double f = res(0).double_value();
114 if (f != f /* isnan(f) */) nlopt_force_stop(data->opt);
120 static double user_function1(unsigned n, const double *x,
121 double *gradient, /* NULL if not needed */
124 octave_function *f = (octave_function *) data_;
125 octave_value_list args(1, 0);
127 for (unsigned i = 0; i < n; ++i)
130 octave_value_list res
131 #if (OCTAVE_MAJOR_VERSION == 4 && OCTAVE_MINOR_VERSION > 2)
132 = octave::feval(f, args, gradient ? 2 : 1);
134 = f->do_multi_index_op(gradient ? 2 : 1, args);
136 if (res.length() < (gradient ? 2 : 1))
137 err_user_supplied_eval("nlopt_optimize");
138 else if (!res(0).is_real_scalar()
139 || (gradient && !res(1).is_real_matrix()
140 && !(n == 1 && res(1).is_real_scalar())))
141 err_user_returned_invalid("nlopt_optimize");
144 if (n == 1 && res(1).is_real_scalar())
145 gradient[0] = res(1).double_value();
147 Matrix grad = res(1).matrix_value();
148 for (unsigned i = 0; i < n; ++i)
149 gradient[i] = grad(i);
152 return res(0).double_value();
157 #define CHECK1(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); nlopt_destroy(opt); nlopt_destroy(local_opt); return NULL; }
159 nlopt_opt make_opt(octave_map &opts, int n)
161 nlopt_opt opt = NULL, local_opt = NULL;
163 nlopt_algorithm algorithm =
164 nlopt_algorithm(struct_val_default(opts, "algorithm",
165 NLOPT_NUM_ALGORITHMS));
166 CHECK1(((int)algorithm) >= 0 && algorithm < NLOPT_NUM_ALGORITHMS,
167 "invalid opt.algorithm");
169 opt = nlopt_create(algorithm, n);
170 CHECK1(opt, "nlopt: out of memory");
172 Matrix m_inf(1, n, -HUGE_VAL);
173 Matrix lb = struct_val_default(opts, "lower_bounds", m_inf);
174 CHECK1(n == lb.numel(), "wrong length of opt.lower_bounds");
175 CHECK1(nlopt_set_lower_bounds(opt, lb.data()) > 0, "nlopt: out of memory");
177 Matrix p_inf(1, n, +HUGE_VAL);
178 Matrix ub = struct_val_default(opts, "upper_bounds", p_inf);
179 CHECK1(n == ub.numel(), "wrong length of opt.upper_bounds");
180 CHECK1(nlopt_set_upper_bounds(opt, ub.data()) > 0, "nlopt: out of memory");
182 nlopt_set_stopval(opt, struct_val_default(opts, "stopval", -HUGE_VAL));
183 nlopt_set_ftol_rel(opt, struct_val_default(opts, "ftol_rel", 0.0));
184 nlopt_set_ftol_abs(opt, struct_val_default(opts, "ftol_abs", 0.0));
185 nlopt_set_xtol_rel(opt, struct_val_default(opts, "xtol_rel", 0.0));
188 Matrix zeros(1, n, 0.0);
189 Matrix xtol_abs = struct_val_default(opts, "xtol_abs", zeros);
190 CHECK1(n == xtol_abs.numel(), "stop.xtol_abs must have same length as x");
191 CHECK1(nlopt_set_xtol_abs(opt, xtol_abs.data())>0, "nlopt: out of memory");
194 nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0) < 0 ?
195 0 : struct_val_default(opts, "maxeval", 0));
196 nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0));
198 nlopt_set_population(opt, struct_val_default(opts, "population", 0));
199 nlopt_set_vector_storage(opt, struct_val_default(opts, "vector_storage", 0));
201 if (opts.contains("initial_step")) {
202 Matrix zeros(1, n, 0.0);
203 Matrix initial_step = struct_val_default(opts, "initial_step", zeros);
204 CHECK1(n == initial_step.numel(),
205 "stop.initial_step must have same length as x");
206 CHECK1(nlopt_set_initial_step(opt, initial_step.data()) > 0,
207 "nlopt: out of memory");
210 if (opts.contains("local_optimizer")) {
211 CHECK1(opts.contents("local_optimizer").numel() == 1
212 && (opts.contents("local_optimizer"))(0).is_map(),
213 "opt.local_optimizer must be a structure");
214 octave_map local_opts = (opts.contents("local_optimizer"))(0).map_value();
215 CHECK1((local_opt = make_opt(local_opts, n)),
216 "error initializing local optimizer");
217 nlopt_set_local_optimizer(opt, local_opt);
218 nlopt_destroy(local_opt); local_opt = NULL;
224 #define CHECK(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); nlopt_destroy(opt); return retval; }
226 DEFUN_DLD(nlopt_optimize, args, nargout, NLOPT_OPTIMIZE_USAGE)
228 octave_value_list retval;
230 nlopt_opt opt = NULL;
232 CHECK(args.length() == 2 && nargout <= 3, "wrong number of args");
234 CHECK(args(0).is_map(), "opt must be structure")
235 octave_map opts = args(0).map_value();
237 CHECK(args(1).is_real_matrix() || args(1).is_real_scalar(),
238 "x must be real vector");
239 Matrix x = args(1).is_real_scalar() ?
240 Matrix(1, 1, args(1).double_value()) : args(1).matrix_value();
243 CHECK((opt = make_opt(opts, n)), "error initializing nlopt options");
245 user_function_data d;
247 d.verbose = struct_val_default(opts, "verbose", 0);
249 if (opts.contains("min_objective")) {
250 CHECK(opts.contents("min_objective").numel() == 1
251 && (opts.contents("min_objective"))(0).is_function_handle(),
252 "opt.min_objective must be a function");
253 d.f = (opts.contents("min_objective"))(0).function_value();
254 nlopt_set_min_objective(opt, user_function, &d);
256 else if (opts.contains("max_objective")) {
257 CHECK(opts.contents("max_objective").numel() == 1
258 && (opts.contents("max_objective"))(0).is_function_handle(),
259 "opt.max_objective must be a function");
260 d.f = (opts.contents("max_objective"))(0).function_value();
261 nlopt_set_max_objective(opt, user_function, &d);
264 CHECK(0,"either opt.min_objective or opt.max_objective must exist");
267 if (opts.contains("fc") && opts.contents("fc").numel() == 1) {
268 CHECK((opts.contents("fc"))(0).is_cell(), "opt.fc must be cell array");
269 Cell fc = (opts.contents("fc"))(0).cell_value();
270 Matrix zeros(1, fc.numel(), 0.0);
271 Matrix fc_tol = struct_val_default(opts, "fc_tol", zeros);
272 CHECK(fc_tol.numel() == fc.numel(),
273 "opt.fc must have same length as opt.fc_tol");
274 for (int i = 0; i < fc.numel(); ++i) {
275 CHECK(fc(i).is_function() || fc(i).is_function_handle(),
276 "opt.fc must be a cell array of function handles");
277 CHECK(nlopt_add_inequality_constraint(opt, user_function1,
278 fc(i).function_value(),
280 "nlopt error adding inequality constraint");
284 if (opts.contains("h") && opts.contents("h").numel() == 1) {
285 CHECK((opts.contents("h"))(0).is_cell(), "opt.h must be cell array");
286 Cell h = (opts.contents("h"))(0).cell_value();
287 Matrix zeros(1, h.numel(), 0.0);
288 Matrix h_tol = struct_val_default(opts, "h_tol", zeros);
289 CHECK(h_tol.numel() == h.numel(),
290 "opt.h must have same length as opt.h_tol");
291 for (int i = 0; i < h.numel(); ++i) {
292 CHECK(h(i).is_function() || h(i).is_function_handle(),
293 "opt.h must be a cell array of function handles");
294 CHECK(nlopt_add_equality_constraint(opt, user_function1,
295 h(i).function_value(),
297 "nlopt error adding equality constraint");
303 nlopt_result ret = nlopt_optimize(opt, x.fortran_vec(), &opt_f);
309 retval(2) = int(ret);