src/octave/nlopt_optimize-oct.cc

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   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:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  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.
  21  */
  22
  23 #include <octave/oct.h>
  24 #include <octave/oct-map.h>
  25 #include <octave/ov.h>
  26 #include <octave/parse.h>
  27 #include <math.h>
  28 #include <stdio.h>
  29
  30 #include "nlopt.h"
  31 #include "nlopt_optimize_usage.h"
  32
  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
  36 #endif
  37
  38 static int struct_val_default(octave_map &m, const std::string& k,
  39                                  int dflt)
  40 {
  41   if (m.contains(k)) {
  42     if (m.contents(k).length() == 1 && (m.contents(k))(0).is_real_scalar())
  43       return (m.contents(k))(0).int_value();
  44   }
  45   return dflt;
  46 }
  47
  48 static double struct_val_default(octave_map &m, const std::string& k,
  49                                  double dflt)
  50 {
  51   if (m.contains(k)) {
  52     if (m.contents(k).length() == 1 && (m.contents(k))(0).is_real_scalar())
  53       return (m.contents(k))(0).double_value();
  54   }
  55   return dflt;
  56 }
  57
  58 static Matrix struct_val_default(octave_map &m, const std::string& k,
  59                                  Matrix &dflt)
  60 {
  61   if (m.contains(k)) {
  62     if ((m.contents(k)).length() == 1) {
  63       if ((m.contents(k))(0).is_real_scalar())
  64         return Matrix(1, dflt.length(), (m.contents(k))(0).double_value());
  65       else if ((m.contents(k))(0).is_real_matrix())
  66         return (m.contents(k))(0).matrix_value();
  67     }
  68   }
  69   return dflt;
  70 }
  71
  72 typedef struct {
  73   octave_function *f;
  74   int neval, verbose;
  75   nlopt_opt opt;
  76 } user_function_data;
  77
  78 static double user_function(unsigned n, const double *x,
  79                             double *gradient, /* NULL if not needed */
  80                             void *data_)
  81 {
  82   user_function_data *data = (user_function_data *) data_;
  83   octave_value_list args(1, 0);
  84   Matrix xm(1,n);
  85   for (unsigned i = 0; i < n; ++i)
  86     xm(i) = x[i];
  87   args(0) = xm;
  88   octave_value_list res
  89 #if (OCTAVE_MAJOR_VERSION == 4 && OCTAVE_MINOR_VERSION > 2)
  90     = octave::feval(data->f, args, gradient ? 2 : 1);
  91 #else
  92     = data->f->do_multi_index_op(gradient ? 2 : 1, args);
  93 #endif
  94   if (res.length() < (gradient ? 2 : 1))
  95     gripe_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     gripe_user_returned_invalid("nlopt_optimize");
 100   else {
 101     if (gradient) {
 102       if (n == 1 && res(1).is_real_scalar())
 103         gradient[0] = res(1).double_value();
 104       else {
 105         Matrix grad = res(1).matrix_value();
 106         for (unsigned i = 0; i < n; ++i)
 107           gradient[i] = grad(i);
 108       }
 109     }
 110     data->neval++;
 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);
 115     return f;
 116   }
 117   return 0;
 118 }
 119
 120 static double user_function1(unsigned n, const double *x,
 121                             double *gradient, /* NULL if not needed */
 122                             void *data_)
 123 {
 124   octave_function *f = (octave_function *) data_;
 125   octave_value_list args(1, 0);
 126   Matrix xm(1,n);
 127   for (unsigned i = 0; i < n; ++i)
 128     xm(i) = x[i];
 129   args(0) = xm;
 130   octave_value_list res
 131 #if (OCTAVE_MAJOR_VERSION == 4 && OCTAVE_MINOR_VERSION > 2)
 132     = octave::feval(f, args, gradient ? 2 : 1);
 133 #else
 134     = f->do_multi_index_op(gradient ? 2 : 1, args);
 135 #endif
 136   if (res.length() < (gradient ? 2 : 1))
 137     gripe_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     gripe_user_returned_invalid("nlopt_optimize");
 142   else {
 143     if (gradient) {
 144       if (n == 1 && res(1).is_real_scalar())
 145         gradient[0] = res(1).double_value();
 146       else {
 147         Matrix grad = res(1).matrix_value();
 148         for (unsigned i = 0; i < n; ++i)
 149           gradient[i] = grad(i);
 150       }
 151     }
 152     return res(0).double_value();
 153   }
 154   return 0;
 155 }
 156
 157 #define CHECK1(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); nlopt_destroy(opt); nlopt_destroy(local_opt); return NULL; }
 158
 159 nlopt_opt make_opt(octave_map &opts, int n)
 160 {
 161   nlopt_opt opt = NULL, local_opt = NULL;
 162
 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");
 168
 169   opt = nlopt_create(algorithm, n);
 170   CHECK1(opt, "nlopt: out of memory");
 171
 172   Matrix m_inf(1, n, -HUGE_VAL);
 173   Matrix lb = struct_val_default(opts, "lower_bounds", m_inf);
 174   CHECK1(n == lb.length(), "wrong length of opt.lower_bounds");
 175   CHECK1(nlopt_set_lower_bounds(opt, lb.data()) > 0, "nlopt: out of memory");
 176
 177   Matrix p_inf(1, n, +HUGE_VAL);
 178   Matrix ub = struct_val_default(opts, "upper_bounds", p_inf);
 179   CHECK1(n == ub.length(), "wrong length of opt.upper_bounds");
 180   CHECK1(nlopt_set_upper_bounds(opt, ub.data()) > 0, "nlopt: out of memory");
 181
 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));
 186
 187   {
 188     Matrix zeros(1, n, 0.0);
 189     Matrix xtol_abs = struct_val_default(opts, "xtol_abs", zeros);
 190     CHECK1(n == xtol_abs.length(), "stop.xtol_abs must have same length as x");
 191     CHECK1(nlopt_set_xtol_abs(opt, xtol_abs.data())>0, "nlopt: out of memory");
 192   }
 193
 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));
 197
 198   nlopt_set_population(opt, struct_val_default(opts, "population", 0));
 199   nlopt_set_vector_storage(opt, struct_val_default(opts, "vector_storage", 0));
 200
 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.length(),
 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");
 208   }
 209
 210   if (opts.contains("local_optimizer")) {
 211     CHECK1(opts.contents("local_optimizer").length() == 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;
 219   }
 220
 221   return opt;
 222 }
 223
 224 #define CHECK(cond, msg) if (!(cond)) { fprintf(stderr, msg "\n\n"); nlopt_destroy(opt); return retval; }
 225
 226 DEFUN_DLD(nlopt_optimize, args, nargout, NLOPT_OPTIMIZE_USAGE)
 227 {
 228   octave_value_list retval;
 229   double A;
 230   nlopt_opt opt = NULL;
 231
 232   CHECK(args.length() == 2 && nargout <= 3, "wrong number of args");
 233
 234   CHECK(args(0).is_map(), "opt must be structure")
 235   octave_map opts = args(0).map_value();
 236
 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();
 241   int n = x.length();
 242
 243   CHECK((opt = make_opt(opts, n)), "error initializing nlopt options");
 244
 245   user_function_data d;
 246   d.neval = 0;
 247   d.verbose = struct_val_default(opts, "verbose", 0);
 248   d.opt = opt;
 249   if (opts.contains("min_objective")) {
 250     CHECK(opts.contents("min_objective").length() == 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);
 255   }
 256   else if (opts.contains("max_objective")) {
 257     CHECK(opts.contents("max_objective").length() == 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);
 262   }
 263   else {
 264     CHECK(0,"either opt.min_objective or opt.max_objective must exist");
 265   }
 266
 267   if (opts.contains("fc") && opts.contents("fc").length() == 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.length(), 0.0);
 271     Matrix fc_tol = struct_val_default(opts, "fc_tol", zeros);
 272     CHECK(fc_tol.length() == fc.length(),
 273           "opt.fc must have same length as opt.fc_tol");
 274     for (int i = 0; i < fc.length(); ++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(),
 279                                             fc_tol(i)) > 0,
 280             "nlopt error adding inequality constraint");
 281     }
 282   }
 283
 284   if (opts.contains("h") && opts.contents("h").length() == 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.length(), 0.0);
 288     Matrix h_tol = struct_val_default(opts, "h_tol", zeros);
 289     CHECK(h_tol.length() == h.length(),
 290           "opt.h must have same length as opt.h_tol");
 291     for (int i = 0; i < h.length(); ++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(),
 296                                             h_tol(i)) > 0,
 297             "nlopt error adding equality constraint");
 298     }
 299   }
 300
 301
 302   double opt_f;
 303   nlopt_result ret = nlopt_optimize(opt, x.fortran_vec(), &opt_f);
 304
 305   retval(0) = x;
 306   if (nargout > 1)
 307     retval(1) = opt_f;
 308   if (nargout > 2)
 309     retval(2) = int(ret);
 310
 311   nlopt_destroy(opt);
 312
 313   return retval;
 314 }