From b9f61d603b314ef72620777f7fc28dddebb99e43 Mon Sep 17 00:00:00 2001 From: stevenj Date: Mon, 5 Apr 2010 22:50:26 -0400 Subject: [PATCH] initial stab at new-style Octave api darcs-hash:20100406025026-c8de0-653797e7b13d5e004082e273187d950c8db7e436.gz --- octave/Makefile.am | 16 +- octave/nlopt_minimize_constrained-oct.cc | 2 +- octave/nlopt_optimize-oct.cc | 297 +++++++++++++++++++++++ octave/nlopt_optimize.m | 16 ++ 4 files changed, 326 insertions(+), 5 deletions(-) create mode 100644 octave/nlopt_optimize-oct.cc create mode 100644 octave/nlopt_optimize.m diff --git a/octave/Makefile.am b/octave/Makefile.am index 1ef2b27..53708dd 100644 --- a/octave/Makefile.am +++ b/octave/Makefile.am @@ -12,8 +12,8 @@ dummy_SOURCES = dummy.c octdir = $(OCT_INSTALL_DIR) mdir = $(M_INSTALL_DIR) if WITH_OCTAVE -oct_DATA = nlopt_minimize_constrained.oct -m_DATA = $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m +oct_DATA = nlopt_minimize_constrained.oct 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 @@ -24,6 +24,14 @@ nlopt_minimize_constrained_usage.h: $(srcdir)/nlopt_minimize_constrained.m 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@ + +nlopt_optimize_usage.h: $(srcdir)/nlopt_optimize.m + echo "#define NLOPT_OPTIMIZE_USAGE \\" > $@ + sed 's/\"/\\"/g' $(srcdir)/nlopt_optimize.m | sed 's,^% ,\",;s,^%,\",;s,$$,\\n\" \\,' >> $@ + echo "" >> $@ + ####################################################################### mexdir = $(MEX_INSTALL_DIR) if WITH_MATLAB @@ -35,6 +43,6 @@ nlopt_minimize_constrained.$(MEXSUFF): nlopt_minimize_constrained-mex.c $(top_bu ####################################################################### -EXTRA_DIST = nlopt_minimize_constrained-oct.cc nlopt_minimize_constrained-mex.c $(MFILES) nlopt_minimize.m nlopt_minimize_constrained.m +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 -CLEANFILES = nlopt_minimize_constrained.oct nlopt_minimize_constrained_usage.h nlopt_minimize_constrained.$(MEXSUFF) nlopt_minimize_constrained-oct.o +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 diff --git a/octave/nlopt_minimize_constrained-oct.cc b/octave/nlopt_minimize_constrained-oct.cc index 74a7176..9d2dfcd 100644 --- a/octave/nlopt_minimize_constrained-oct.cc +++ b/octave/nlopt_minimize_constrained-oct.cc @@ -105,7 +105,7 @@ DEFUN_DLD(nlopt_minimize_constrained, args, nargout, NLOPT_MINIMIZE_CONSTRAINED_ CHECK(args.length() == 9 && nargout <= 3, "wrong number of args"); - CHECK(args(0).is_real_scalar(), "n must be real scalar"); + CHECK(args(0).is_real_scalar(), "algorithm must be integer"); nlopt_algorithm algorithm = nlopt_algorithm(args(0).int_value()); user_function_data d; diff --git a/octave/nlopt_optimize-oct.cc b/octave/nlopt_optimize-oct.cc new file mode 100644 index 0000000..5a393cf --- /dev/null +++ b/octave/nlopt_optimize-oct.cc @@ -0,0 +1,297 @@ +/* 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_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) +{ + if (m.contains(k)) { + if (m.contents(k).length() == 1 && (m.contents(k))(0).is_real_scalar()) + return (m.contents(k))(0).int_value(); + } + return dflt; +} + +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; + int neval, verbose; +} user_function_data; + +static double user_function(unsigned 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, 0); + Matrix xm(1,n); + for (unsigned i = 0; i < n; ++i) + xm(i) = x[i]; + args(0) = xm; + 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_optimize"); + 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_optimize"); + 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 (unsigned i = 0; i < n; ++i) + gradient[i] = grad(i); + } + } + data->neval++; + if (data->verbose) printf("nlopt_optimize eval #%d: %g\n", + data->neval, res(0).double_value()); + return res(0).double_value(); + } + return 0; +} + +static double user_function1(unsigned n, const double *x, + double *gradient, /* NULL if not needed */ + void *data_) +{ + octave_function *f = (octave_function *) data_; + octave_value_list args(1, 0); + Matrix xm(1,n); + for (unsigned i = 0; i < n; ++i) + xm(i) = x[i]; + args(0) = xm; + octave_value_list res = f->do_multi_index_op(gradient ? 2 : 1, args); + if (res.length() < (gradient ? 2 : 1)) + gripe_user_supplied_eval("nlopt_optimize"); + 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_optimize"); + 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 (unsigned i = 0; i < n; ++i) + gradient[i] = grad(i); + } + } + return res(0).double_value(); + } + return 0; +} + +#define CHECK1(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); nlopt_destroy(opt); nlopt_destroy(local_opt); return NULL; } + +nlopt_opt make_opt(Octave_map &opts, int n) +{ + nlopt_opt opt = NULL, local_opt = NULL; + + nlopt_algorithm algorithm = + nlopt_algorithm(struct_val_default(opts, "algorithm", + NLOPT_NUM_ALGORITHMS)); + CHECK1(((int)algorithm) >= 0 && algorithm < NLOPT_NUM_ALGORITHMS, + "invalid opt.algorithm"); + + opt = nlopt_create(algorithm, n); + CHECK1(opt, "nlopt: out of memory"); + + Matrix m_inf(1, n, -HUGE_VAL); + Matrix lb = struct_val_default(opts, "lb", m_inf); + CHECK1(n == lb.length(), "wrong length of opt.lb"); + CHECK1(nlopt_set_lower_bounds(opt, lb.data()) > 0, "nlopt: out of memory"); + + Matrix p_inf(1, n, +HUGE_VAL); + Matrix ub = struct_val_default(opts, "ub", p_inf); + CHECK1(n == ub.length(), "wrong length of opt.ub"); + CHECK1(nlopt_set_upper_bounds(opt, ub.data()) > 0, "nlopt: out of memory"); + + 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)); + + { + Matrix zeros(1, n, 0.0); + Matrix xtol_abs = struct_val_default(opts, "xtol_abs", zeros); + CHECK1(n == xtol_abs.length(), "stop.xtol_abs must have same length as x"); + CHECK1(nlopt_set_xtol_abs(opt, xtol_abs.data())>0, "nlopt: out of memory"); + } + + nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0) < 0 ? + 0 : struct_val_default(opts, "maxeval", 0)); + nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0)); + + nlopt_set_population(opt, struct_val_default(opts, "population", 0)); + + if (opts.contains("initial_step")) { + Matrix zeros(1, n, 0.0); + Matrix initial_step = struct_val_default(opts, "initial_step", zeros); + CHECK1(n == initial_step.length(), + "stop.initial_step must have same length as x"); + CHECK1(nlopt_set_initial_step(opt, initial_step.data()) > 0, + "nlopt: out of memory"); + } + + if (opts.contains("local_optimizer")) { + CHECK1(opts.contents("local_optimizer").length() == 1 + && (opts.contents("local_optimizer"))(0).is_map(), + "opt.local_optimizer must be a structure"); + Octave_map local_opts = (opts.contents("local_optimizer"))(0).map_value(); + 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; + } + + return opt; +} + +#define CHECK(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); print_usage("nlopt_optimize"); nlopt_destroy(opt); return retval; } + +DEFUN_DLD(nlopt_optimize, args, nargout, NLOPT_OPTIMIZE_USAGE) +{ + octave_value_list retval; + double A; + nlopt_opt opt = NULL; + + CHECK(args.length() == 2 && nargout <= 3, "wrong number of args"); + + CHECK(args(0).is_map(), "opt must be structure") + Octave_map opts = args(0).map_value(); + + CHECK(args(1).is_real_matrix() || args(1).is_real_scalar(), + "x must be real vector"); + Matrix x = args(1).is_real_scalar() ? + Matrix(1, 1, args(1).double_value()) : args(1).matrix_value(); + int n = x.length(); + + CHECK((opt = make_opt(opts, n)), "error initializing nlopt options"); + + user_function_data d; + d.neval = 0; + d.verbose = struct_val_default(opts, "verbose", 0); + if (opts.contains("min_objective")) { + CHECK(opts.contents("min_objective").length() == 1 + && (opts.contents("min_objective"))(0).is_function_handle(), + "opt.min_objective must be a function"); + d.f = (opts.contents("min_objective"))(0).function_value(); + nlopt_set_min_objective(opt, user_function, &d); + } + else if (opts.contains("max_objective")) { + CHECK(opts.contents("max_objective").length() == 1 + && (opts.contents("max_objective"))(0).is_function_handle(), + "opt.max_objective must be a function"); + d.f = (opts.contents("max_objective"))(0).function_value(); + nlopt_set_max_objective(opt, user_function, &d); + } + else { + CHECK(0,"either opt.min_objective or opt.max_objective must exist"); + } + + if (opts.contains("fc") && opts.contents("fc").length() == 1) { + CHECK((opts.contents("fc"))(0).is_cell(), "opt.fc must be cell array"); + Cell fc = (opts.contents("fc"))(0).cell_value(); + Matrix zeros(1, fc.length(), 0.0); + Matrix fc_tol = struct_val_default(opts, "fc_tol", zeros); + CHECK(fc_tol.length() == fc.length(), + "opt.fc must have same length as opt.fc_tol"); + for (int i = 0; i < fc.length(); ++i) { + CHECK(fc(i).is_function() || fc(i).is_function_handle(), + "opt.fc must be a cell array of function handles"); + CHECK(nlopt_add_inequality_constraint(opt, user_function1, + fc(i).function_value(), + fc_tol(i)) > 0, + "nlopt error adding inequality constraint"); + } + } + + if (opts.contains("h") && opts.contents("h").length() == 1) { + CHECK((opts.contents("h"))(0).is_cell(), "opt.h must be cell array"); + Cell h = (opts.contents("h"))(0).cell_value(); + Matrix zeros(1, h.length(), 0.0); + Matrix h_tol = struct_val_default(opts, "h_tol", zeros); + CHECK(h_tol.length() == h.length(), + "opt.h must have same length as opt.h_tol"); + for (int i = 0; i < h.length(); ++i) { + CHECK(h(i).is_function() || h(i).is_function_handle(), + "opt.h must be a cell array of function handles"); + CHECK(nlopt_add_equality_constraint(opt, user_function1, + h(i).function_value(), + h_tol(i)) > 0, + "nlopt error adding equality constraint"); + } + } + + + double opt_f; + nlopt_result ret = nlopt_optimize(opt, x.fortran_vec(), &opt_f); + + retval(0) = x; + if (nargout > 1) + retval(1) = opt_f; + if (nargout > 2) + retval(2) = int(ret); + + nlopt_destroy(opt); + + return retval; +} diff --git a/octave/nlopt_optimize.m b/octave/nlopt_optimize.m new file mode 100644 index 0000000..c2e486f --- /dev/null +++ b/octave/nlopt_optimize.m @@ -0,0 +1,16 @@ +% Usage: [xopt, fopt, retcode] = nlopt_optimize(opt, xinit) +% +% Optimizes (minimizes or maximizes) a nonlinear function under +% nonlinear constraints from the starting guess xinit, where the +% objective, constraints, stopping criteria, and other options are +% specified in the structure opt described below. A variety of local +% and global optimization algorithms can be used, as specified by the +% opt.algorithm parameter described below. Returns the optimum +% function value fopt, the location xopt of the optimum, and a +% return code retcode described below (> 0 on success). +% +% This function is a front-end for the external routine nlopt_optimize +% in the free NLopt nonlinear-optimization library, which is a wrapper +% around a number of free/open-source optimization subroutines. More +% details can be found on the NLopt web page (ab-initio.mit.edu/nlopt) +% and also under 'man nlopt_minimize' on Unix. -- 2.30.2