octave_packages/optim-1.2.0/nonlin_min.m

   1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin})
  18 ## @deftypefnx {Function File} {[@var{p}, @var{objf}, @var{cvg}, @var{outp}] =} nonlin_min (@var{f}, @var{pin}, @var{settings})
  19 ##
  20 ## Frontend for constrained nonlinear minimization of a scalar objective
  21 ## function. The functions supplied by the user have a minimal
  22 ## interface; any additionally needed constants can be supplied by
  23 ## wrapping the user functions into anonymous functions.
  24 ##
  25 ## The following description applies to usage with vector-based
  26 ## parameter handling. Differences in usage for structure-based
  27 ## parameter handling will be explained in a separate section below.
  28 ##
  29 ## @var{f}: objective function. It gets a column vector of real
  30 ## parameters as argument. In gradient determination, this function may
  31 ## be called with an informational second argument, whose content
  32 ## depends on the function for gradient determination.
  33 ##
  34 ## @var{pin}: real column vector of initial parameters.
  35 ##
  36 ## @var{settings}: structure whose fields stand for optional settings
  37 ## referred to below. The fields can be set by @code{optimset()} with
  38 ## Octave versions 3.3.55 or greater; with older Octave versions, the
  39 ## fields must be set directly as structure-fields in the correct case.
  40 ##
  41 ## The returned values are the column vector of final parameters
  42 ## @var{p}, the final value of the objective function @var{objf}, an
  43 ## integer @var{cvg} indicating if and how optimization succeeded or
  44 ## failed, and a structure @var{outp} with additional information,
  45 ## curently with only one field: @var{niter}, the number of iterations.
  46 ## @var{cvg} is greater than zero for success and less than or equal to
  47 ## zero for failure; its possible values depend on the used backend and
  48 ## currently can be @code{0} (maximum number of iterations exceeded),
  49 ## @code{1} (fixed number of iterations completed, e.g. in stochastic
  50 ## optimizers), @code{2} (parameter change less than specified precision
  51 ## in two consecutive iterations), @code{3} (improvement in objective
  52 ## function less than specified), or @code{-4} (algorithm got stuck).
  53 ##
  54 ## @var{settings}:
  55 ##
  56 ## @code{Algorithm}: String specifying the backend. Currently available
  57 ## are @code{"lm_feasible"} (default) and @code{"siman"}. They are
  58 ## described in separate sections below.
  59 ##
  60 ## @code{objf_grad}: Function computing the gradient of the objective
  61 ## function with respect to the parameters, assuming residuals are
  62 ## reshaped to a vector. Default: finite differences. Will be called
  63 ## with the column vector of parameters and an informational structure
  64 ## as arguments. The structure has the fields @code{f}: value of
  65 ## objective function for current parameters, @code{fixed}: logical
  66 ## vector indicating which parameters are not optimized, so these
  67 ## partial derivatives need not be computed and can be set to zero,
  68 ## @code{diffp}, @code{diff_onesided}, @code{lbound}, @code{ubound}:
  69 ## identical to the user settings of this name, @code{plabels}:
  70 ## 1-dimensional cell-array of column-cell-arrays, each column with
  71 ## labels for all parameters, the first column contains the numerical
  72 ## indices of the parameters. The default gradient function will call
  73 ## the objective function with the second argument set with fields
  74 ## @code{f}: as the @code{f} passed to the gradient function,
  75 ## @code{plabels}: cell-array of 1x1 cell-arrays with the entries of the
  76 ## column-cell-arrays of @code{plabels} as passed to the jacobian
  77 ## function corresponding to current parameter, @code{side}: @code{0}
  78 ## for one-sided interval, @code{1} or @code{2}, respectively, for the
  79 ## sides of a two-sided interval, and @code{parallel}: logical scalar
  80 ## indicating parallel computation of partial derivatives.
  81 ##
  82 ## @code{objf_hessian}: Function computing the Hessian of the objective
  83 ## function with respect to the parameters. The default is backend
  84 ## specific. Will be called with the column vector of parameters as
  85 ## argument.
  86 ##
  87 ## @code{diffp}: column vector of fractional intervals (doubled for
  88 ## central intervals) supposed to be used by gradient functions
  89 ## performing finite differencing. Default: @code{.001 * ones (size
  90 ## (parameters))}. The default gradient function will use these as
  91 ## absolute intervals for parameters with value zero.
  92 ##
  93 ## @code{diff_onesided}: logical column vector indicating that one-sided
  94 ## intervals should be used by gradient functions performing finite
  95 ## differencing. Default: @code{false (size (parameters))}.
  96 ##
  97 ## @code{complex_step_derivative_objf},
  98 ## @code{complex_step_derivative_inequc},
  99 ## @code{complex_step_derivative_equc}: logical scalars, default: false.
 100 ## Estimate gradient of objective function, general inequality
 101 ## constraints, and general equality constraints, respectively, with
 102 ## complex step derivative approximation. Use only if you know that your
 103 ## objective function, function of general inequality constraints, or
 104 ## function of general equality constraints, respectively, is suitable
 105 ## for this. No user function for the respective gradient must be
 106 ## specified.
 107 ##
 108 ## @code{cstep}: scalar step size for complex step derivative
 109 ## approximation. Default: 1e-20.
 110 ##
 111 ## @code{fixed}: logical column vector indicating which parameters
 112 ## should not be optimized, but kept to their inital value. Fixing is
 113 ## done independently of the backend, but the backend may choose to fix
 114 ## additional parameters under certain conditions.
 115 ##
 116 ## @code{lbound}, @code{ubound}: column vectors of lower and upper
 117 ## bounds for parameters. Default: @code{-Inf} and @code{+Inf},
 118 ## respectively. The bounds are non-strict, i.e. parameters are allowed
 119 ## to be exactly equal to a bound. The default gradient function will
 120 ## respect bounds (but no further inequality constraints) in finite
 121 ## differencing.
 122 ##
 123 ## @code{inequc}: Further inequality constraints. Cell-array containing
 124 ## up to four entries, two entries for linear inequality constraints
 125 ## and/or one or two entries for general inequality constraints. Either
 126 ## linear or general constraints may be the first entries, but the two
 127 ## entries for linear constraints must be adjacent and, if two entries
 128 ## are given for general constraints, they also must be adjacent. The
 129 ## two entries for linear constraints are a matrix (say @code{m}) and a
 130 ## vector (say @code{v}), specifying linear inequality constraints of
 131 ## the form @code{m.' * parameters + v >= 0}. The first entry for
 132 ## general constraints must be a differentiable vector valued function
 133 ## (say @code{h}), specifying general inequality constraints of the form
 134 ## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
 135 ## paraters and the optional argument @code{idx} is a logical index.
 136 ## @code{h} has to return the values of all constraints if @code{idx} is
 137 ## not given. It may choose to return only the indexed constraints if
 138 ## @code{idx} is given (so computation of the other constraints can be
 139 ## spared); in this case, the additional setting @code{inequc_f_idx} has
 140 ## to be set to @code{true}. In gradient determination, this function
 141 ## may be called with an informational third argument, whose content
 142 ## depends on the function for gradient determination. If a second entry
 143 ## for general inequality constraints is given, it must be a function
 144 ## computing the jacobian of the constraints with respect to the
 145 ## parameters. For this function, the description of @code{dfdp} above
 146 ## applies, with 2 exceptions: 1) it is called with 3 arguments since it
 147 ## has an additional argument @code{idx} --- a logical index --- at
 148 ## second position, indicating which rows of the jacobian must be
 149 ## returned (if the function chooses to return only indexed rows, the
 150 ## additional setting @code{inequc_df_idx} has to be set to
 151 ## @code{true}). 2) the default jacobian function calls @code{h} with 3
 152 ## arguments, since the argument @code{idx} is also supplied. Note that
 153 ## specifying linear constraints as general constraints will generally
 154 ## waste performance, even if further, non-linear, general constraints
 155 ## are also specified.
 156 ##
 157 ## @code{equc}: Equality constraints. Specified the same way as
 158 ## inequality constraints (see @code{inequc}). The respective additional
 159 ## settings are named @code{equc_f_idx} and @code{equc_df_idx}.
 160 ##
 161 ## @code{cpiv}: Function for complementary pivoting, usable in
 162 ## algorithms for constraints. Default: @ cpiv_bard. Only the default
 163 ## function is supplied with the package.
 164 ##
 165 ## @code{TolFun}: Minimum fractional improvement in objective function
 166 ## in an iteration (abortion criterium). Default: .0001.
 167 ##
 168 ## @code{MaxIter}: Maximum number of iterations (abortion criterium).
 169 ## Default: backend-specific.
 170 ##
 171 ## @code{fract_prec}: Column Vector, minimum fractional change of
 172 ## parameters in an iteration (abortion criterium if violated in two
 173 ## consecutive iterations). Default: backend-specific.
 174 ##
 175 ## @code{max_fract_change}: Column Vector, enforced maximum fractional
 176 ## change in parameters in an iteration. Default: backend-specific.
 177 ##
 178 ## @code{Display}: String indicating the degree of verbosity. Default:
 179 ## @code{"off"}. Possible values are currently @code{"off"} (no
 180 ## messages) and @code{"iter"} (some messages after each iteration).
 181 ## Support of this setting and its exact interpretation are
 182 ## backend-specific.
 183 ##
 184 ## @code{debug}: Logical scalar, default: @code{false}. Will be passed
 185 ## to the backend, which might print debugging information if true.
 186 ##
 187 ## Structure-based parameter handling
 188 ##
 189 ## The setting @code{param_order} is a cell-array with names of the
 190 ## optimized parameters. If not given, and initial parameters are a
 191 ## structure, all parameters in the structure are optimized. If initial
 192 ## parameters are a structure, it is an error if @code{param_order} is
 193 ## not given and there are any non-structure-based configuration items
 194 ## or functions.
 195 ##
 196 ## The initial parameters @var{pin} can be given as a structure
 197 ## containing at least all fields named in @code{param_order}. In this
 198 ## case the returned parameters @var{p} will also be a structure.
 199 ##
 200 ## Each user-supplied function can be called with the argument
 201 ## containing the current parameters being a structure instead of a
 202 ## column vector. For this, a corresponding setting must be set to
 203 ## @code{true}: @code{objf_pstruct} (objective function),
 204 ## @code{objf_grad_pstruct} (gradient of objective function),
 205 ## @code{objf_hessian_pstruct} (hessian of objective function),
 206 ## @code{f_inequc_pstruct} (general inequality constraints),
 207 ## @code{df_inequc_pstruct} (jacobian of general inequality
 208 ## constraints), @code{f_equc_pstruct} (general equality constraints),
 209 ## and @code{df_equc_pstruct} (jacobian of general equality
 210 ## constraints). If a gradient (jacobian) function is configured in such
 211 ## a way, it must return the entries (columns) of the gradient
 212 ## (jacobian) as fields of a structure under the respective parameter
 213 ## names. If the hessian function is configured in such a way, it must
 214 ## return a structure (say @code{h}) with fields e.g. as
 215 ## @code{h.a.b = value} for @code{value} being the 2nd partial derivative
 216 ## with respect to @code{a} and @code{b}. There is no need to also
 217 ## specify the field @code{h.b.a} in this example.
 218 ##
 219 ## Similarly, for specifying linear constraints, instead of the matrix
 220 ## (called @code{m} above), a structure containing the rows of the
 221 ## matrix in fields under the respective parameter names can be given.
 222 ## In this case, rows containing only zeros need not be given.
 223 ##
 224 ## The vector-based settings @code{lbound}, @code{ubound},
 225 ## @code{fixed}, @code{diffp}, @code{diff_onesided}, @code{fract_prec},
 226 ## and @code{max_fract_change} can be replaced by the setting
 227 ## @code{param_config}. It is a structure that can contain fields named
 228 ## in @code{param_order}. For each such field, there may be subfields
 229 ## with the same names as the above vector-based settings, but
 230 ## containing a scalar value for the respective parameter. If
 231 ## @code{param_config} is specified, none of the above
 232 ## vector/matrix-based settings may be used.
 233 ##
 234 ## Additionally, named parameters are allowed to be non-scalar real
 235 ## arrays. In this case, their dimensions are given by the setting
 236 ## @code{param_dims}, a cell-array of dimension vectors, each containing
 237 ## at least two dimensions; if not given, dimensions are taken from the
 238 ## initial parameters, if these are given in a structure. Any
 239 ## vector-based settings or not structure-based linear constraints then
 240 ## must correspond to an order of parameters with all parameters
 241 ## reshaped to vectors and concatenated in the user-given order of
 242 ## parameter names. Structure-based settings or structure-based initial
 243 ## parameters must contain arrays with dimensions reshapable to those of
 244 ## the respective parameters.
 245 ##
 246 ## Description of backends
 247 ##
 248 ## "lm_feasible"
 249 ##
 250 ## A Levenberg/Marquardt-like optimizer, attempting to honour
 251 ## constraints throughout the course of optimization. This means that
 252 ## the initial parameters must not violate constraints (to find an
 253 ## initial feasible set of parameters, e.g. Octaves @code{sqp} can be
 254 ## used, by specifying an objective function which is constant or which
 255 ## returns the quadratic distance to the initial values). If the
 256 ## constraints need only be honoured in the result of the optimization,
 257 ## Octaves @code{sqp} may be preferable. The Hessian is either supplied
 258 ## by the user or is approximated by the BFGS algorithm.
 259 ##
 260 ## Returned value @var{cvg} will be @code{2} or @code{3} for success and
 261 ## @code{0} or @code{-4} for failure (see above for meaning).
 262 ##
 263 ## Backend-specific defaults are: @code{MaxIter}: 20, @code{fract_prec}:
 264 ## @code{zeros (size (parameters))}, @code{max_fract_change}: @code{Inf}
 265 ## for all parameters.
 266 ##
 267 ## Interpretation of @code{Display}: if set to @code{"iter"}, currently
 268 ## only information on applying @code{max_fract_change} is printed.
 269 ##
 270 ## "siman"
 271 ##
 272 ## A simulated annealing (stochastic) optimizer, changing all parameters
 273 ## at once in a single step, so being suitable for non-bound
 274 ## constraints.
 275 ##
 276 ## No gradient or hessian of the objective function is used. The
 277 ## settings @code{MaxIter}, @code{fract_prec}, @code{TolFun}, and
 278 ## @code{max_fract_change} are not honoured.
 279 ##
 280 ## Accepts the additional settings @code{T_init} (initial temperature,
 281 ## default 0.01), @code{T_min} (final temperature, default 1.0e-5),
 282 ## @code{mu_T} (factor of temperature decrease, default 1.005),
 283 ## @code{iters_fixed_T} (iterations within one temperature step, default
 284 ## 10), @code{max_rand_step} (column vector or structure-based
 285 ## configuration of maximum random steps for each parameter, default
 286 ## 0.005 * @var{pin}), @code{stoch_regain_constr} (if @code{true},
 287 ## regain constraints after a random step, otherwise take new random
 288 ## value until constraints are met, default false), @code{trace_steps}
 289 ## (set field @code{trace} of @var{outp} with a matrix with a row for
 290 ## each step, first column iteration number, second column repeat number
 291 ## within iteration, third column value of objective function, rest
 292 ## columns parameter values, default false), and @code{siman_log} (set
 293 ## field @code{log} of @var{outp} with a matrix with a row for each
 294 ## iteration, first column temperature, second column value of objective
 295 ## function, rest columns numbers of tries with decrease, no decrease
 296 ## but accepted, and no decrease and rejected.
 297 ##
 298 ## Steps with increase @code{diff} of objective function are accepted if
 299 ## @code{rand (1) < exp (- diff / T)}, where @code{T} is the temperature
 300 ## of the current iteration.
 301 ##
 302 ## If regaining of constraints failed, optimization will be aborted and
 303 ## returned value of @var{cvg} will be @code{0}. Otherwise, @var{cvg}
 304 ## will be @code{1}.
 305 ##
 306 ## Interpretation of @code{Display}: if set to @code{"iter"}, an
 307 ## informational line is printed after each iteration.
 308 ##
 309 ## @end deftypefn
 310
 311 ## disabled PKG_ADD: __all_opts__ ("nonlin_min");
 312
 313 function [p, objf, cvg, outp] = nonlin_min (f, pin, settings)
 314
 315   if (compare_versions (version (), "3.3.55", "<"))
 316     ## optimset mechanism was fixed for option names with underscores
 317     ## sometime in 3.3.54+, if I remember right
 318     optimget = @ __optimget__;
 319   endif
 320
 321   if (compare_versions (version (), "3.2.4", "<="))
 322     ## For bug #31484; but Octave 3.6... shows bug #36288 due to this
 323     ## workaround. Octave 3.7... seems to be all right.
 324     __dfdp__ = @ __dfdp__;
 325   endif
 326
 327   ## some scalar defaults; some defaults are backend specific, so
 328   ## lacking elements in respective constructed vectors will be set to
 329   ## NA here in the frontend
 330   diffp_default = .001;
 331   stol_default = .0001;
 332   cstep_default = 1e-20;
 333
 334   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
 335     p = optimset ("param_config", [], \
 336                   "param_order", [], \
 337                   "param_dims", [], \
 338                   "f_inequc_pstruct", false, \
 339                   "f_equc_pstruct", false, \
 340                   "objf_pstruct", false, \
 341                   "df_inequc_pstruct", false, \
 342                   "df_equc_pstruct", false, \
 343                   "objf_grad_pstruct", false, \
 344                   "objf_hessian_pstruct", false, \
 345                   "lbound", [], \
 346                   "ubound", [], \
 347                   "objf_grad", [], \
 348                   "objf_hessian", [], \
 349                   "cpiv", @ cpiv_bard, \
 350                   "max_fract_change", [], \
 351                   "fract_prec", [], \
 352                   "diffp", [], \
 353                   "diff_onesided", [], \
 354                   "complex_step_derivative_objf", false, \
 355                   "complex_step_derivative_inequc", false, \
 356                   "complex_step_derivative_equc", false, \
 357                   "cstep", cstep_default, \
 358                   "fixed", [], \
 359                   "inequc", [], \
 360                   "equc", [], \
 361                   "inequc_f_idx", false, \
 362                   "inequc_df_idx", false, \
 363                   "equc_f_idx", false, \
 364                   "equc_df_idx", false, \
 365                   "TolFun", stol_default, \
 366                   "MaxIter", [], \
 367                   "Display", "off", \
 368                   "Algorithm", "lm_feasible", \
 369                   "T_init", .01, \
 370                   "T_min", 1.0e-5, \
 371                   "mu_T", 1.005, \
 372                   "iters_fixed_T", 10, \
 373                   "max_rand_step", [], \
 374                   "stoch_regain_constr", false, \
 375                   "trace_steps", false, \
 376                   "siman_log", false, \
 377                   "debug", false);
 378     return;
 379   endif
 380
 381   if (nargin < 2 || nargin > 3)
 382     print_usage ();
 383   endif
 384
 385   if (nargin == 2)
 386     settings = struct ();
 387   endif
 388
 389   if (ischar (f))
 390     f = str2func (f);
 391   endif
 392
 393   if (! (pin_struct = isstruct (pin)))
 394     if (! isvector (pin) || columns (pin) > 1)
 395       error ("initial parameters must be either a structure or a column vector");
 396     endif
 397   endif
 398
 399   #### processing of settings and consistency checks
 400
 401   pconf = optimget (settings, "param_config");
 402   pord = optimget (settings, "param_order");
 403   pdims = optimget (settings, "param_dims");
 404   f_inequc_pstruct = optimget (settings, "f_inequc_pstruct", false);
 405   f_equc_pstruct = optimget (settings, "f_equc_pstruct", false);
 406   f_pstruct = optimget (settings, "objf_pstruct", false);
 407   dfdp_pstruct = optimget (settings, "objf_grad_pstruct", f_pstruct);
 408   hessian_pstruct = optimget (settings, "objf_hessian_pstruct", f_pstruct);
 409   df_inequc_pstruct = optimget (settings, "df_inequc_pstruct", \
 410                                 f_inequc_pstruct);
 411   df_equc_pstruct = optimget (settings, "df_equc_pstruct", \
 412                               f_equc_pstruct);
 413   lbound = optimget (settings, "lbound");
 414   ubound = optimget (settings, "ubound");
 415   dfdp = optimget (settings, "objf_grad");
 416   if (ischar (dfdp)) dfdp = str2func (dfdp); endif
 417   hessian = optimget (settings, "objf_hessian");
 418   max_fract_change = optimget (settings, "max_fract_change");
 419   fract_prec = optimget (settings, "fract_prec");
 420   diffp = optimget (settings, "diffp");
 421   diff_onesided = optimget (settings, "diff_onesided");
 422   fixed = optimget (settings, "fixed");
 423   do_cstep = optimget (settings, "complex_step_derivative_objf", false);
 424   cstep = optimget (settings, "cstep", cstep_default);
 425   if (do_cstep && ! isempty (dfdp))
 426     error ("both 'complex_step_derivative_objf' and 'objf_grad' are set");
 427   endif
 428   do_cstep_inequc = \
 429       optimget (settings, "complex_step_derivative_inequc", false);
 430   do_cstep_equc = optimget (settings, "complex_step_derivative_equc", \
 431                             false);
 432   max_rand_step = optimget (settings, "max_rand_step");
 433
 434   any_vector_conf = ! (isempty (lbound) && isempty (ubound) && \
 435                        isempty (max_fract_change) && \
 436                        isempty (fract_prec) && isempty (diffp) && \
 437                        isempty (diff_onesided) && isempty (fixed) && \
 438                        isempty (max_rand_step));
 439
 440   ## collect constraints
 441   [mc, vc, f_genicstr, df_gencstr, user_df_gencstr] = \
 442       __collect_constraints__ (optimget (settings, "inequc"), \
 443                                do_cstep_inequc, "inequality constraints");
 444   [emc, evc, f_genecstr, df_genecstr, user_df_genecstr] = \
 445       __collect_constraints__ (optimget (settings, "equc"), \
 446                                do_cstep_equc, "equality constraints");
 447   mc_struct = isstruct (mc);
 448   emc_struct = isstruct (emc);
 449
 450   ## correct "_pstruct" settings if functions are not supplied, handle
 451   ## constraint functions not honoring indices
 452   if (isempty (dfdp)) dfdp_pstruct = false; endif
 453   if (isempty (hessian)) hessian_pstruct = false; endif
 454   if (isempty (f_genicstr))
 455     f_inequc_pstruct = false;
 456   elseif (! optimget (settings, "inequc_f_idx", false))
 457     f_genicstr = @ (p, varargin) apply_idx_if_given \
 458         (f_genicstr (p, varargin{:}), varargin{:});
 459   endif
 460   if (isempty (f_genecstr))
 461     f_equc_pstruct = false;
 462   elseif (! optimget (settings, "equc_f_idx", false))
 463     f_genecstr = @ (p, varargin) apply_idx_if_given \
 464         (f_genecstr (p, varargin{:}), varargin{:});
 465   endif
 466   if (user_df_gencstr)
 467     if (! optimget (settings, "inequc_df_idx", false))
 468       df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
 469     endif
 470   else
 471     df_inequc_pstruct = false;
 472   endif
 473   if (user_df_genecstr)
 474     if (! optimget (settings, "equc_df_idx", false))
 475       df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
 476     endif
 477   else
 478     df_equc_pstruct = false;
 479   endif
 480
 481   ## some settings require a parameter order
 482   if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
 483       f_equc_pstruct || f_pstruct || dfdp_pstruct || \
 484       hessian_pstruct || df_inequc_pstruct || df_equc_pstruct || \
 485       mc_struct || emc_struct)
 486     if (isempty (pord))
 487       if (pin_struct)
 488         if (any_vector_conf || \
 489             ! (f_pstruct && \
 490                (f_inequc_pstruct || isempty (f_genicstr)) && \
 491                (f_equc_pstruct || isempty (f_genecstr)) && \
 492                (dfdp_pstruct || isempty (dfdp)) && \
 493                (hessian_pstruct || isempty (hessian)) && \
 494                (df_inequc_pstruct || ! user_df_gencstr) && \
 495                (df_equc_pstruct || ! user_df_genecstr) && \
 496                (mc_struct || isempty (mc)) && \
 497                (emc_struct || isempty (emc))))
 498           error ("no parameter order specified and constructing a parameter order from the structure of initial parameters can not be done since not all configuration or given functions are structure based");
 499         else
 500           pord = fieldnames (pin);
 501         endif
 502       else
 503         error ("given settings require specification of parameter order or initial parameters in the form of a structure");
 504       endif
 505     endif
 506     pord = pord(:);
 507     if (pin_struct && ! all (isfield (pin, pord)))
 508       error ("some initial parameters lacking");
 509     endif
 510     if ((nnames = rows (unique (pord))) < rows (pord))
 511       error ("duplicate parameter names in 'param_order'");
 512     endif
 513     if (isempty (pdims))
 514       if (pin_struct)
 515         pdims = cellfun \
 516             (@ size, fields2cell (pin, pord), "UniformOutput", false);
 517       else
 518         pdims = num2cell (ones (nnames, 2), 2);
 519       endif
 520     else
 521       pdims = pdims(:);
 522       if (pin_struct && \
 523           ! all (cellfun (@ (x, y) prod (size (x)) == prod (y), \
 524                           struct2cell (pin), pdims)))
 525         error ("given param_dims and dimensions of initial parameters do not match");
 526       endif
 527     endif
 528     if (nnames != rows (pdims))
 529       error ("lengths of 'param_order' and 'param_dims' not equal");
 530     endif
 531     pnel = cellfun (@ prod, pdims);
 532     ppartidx = pnel;
 533     if (any (pnel > 1))
 534       pnonscalar = true;
 535       cpnel = num2cell (pnel);
 536       prepidx = cat (1, cellfun \
 537                      (@ (x, n) x(ones (1, n), 1), \
 538                       num2cell ((1:nnames).'), cpnel, \
 539                       "UniformOutput", false){:});
 540       epord = pord(prepidx, 1);
 541       psubidx = cat (1, cellfun \
 542                      (@ (n) (1:n).', cpnel, \
 543                       "UniformOutput", false){:});
 544     else
 545       pnonscalar = false; # some less expensive interfaces later
 546       prepidx = (1:nnames).';
 547       epord = pord;
 548       psubidx = ones (nnames, 1);
 549     endif
 550   else
 551     pord = []; # spares checks for given but not needed
 552   endif
 553
 554   if (pin_struct)
 555     np = sum (pnel);
 556   else
 557     np = length (pin);
 558     if (! isempty (pord) && np != sum (pnel))
 559       error ("number of initial parameters not correct");
 560     endif
 561   endif
 562
 563   plabels = num2cell (num2cell ((1:np).'));
 564   if (! isempty (pord))
 565     plabels = cat (2, plabels, num2cell (epord), \
 566                    num2cell (num2cell (psubidx)));
 567   endif
 568
 569   ## some useful vectors
 570   zerosvec = zeros (np, 1);
 571   NAvec = NA (np, 1);
 572   Infvec = Inf (np, 1);
 573   falsevec = false (np, 1);
 574   sizevec = [np, 1];
 575
 576   ## collect parameter-related configuration
 577   if (! isempty (pconf))
 578     ## use supplied configuration structure
 579
 580     ## parameter-related configuration is either allowed by a structure
 581     ## or by vectors
 582     if (any_vector_conf)
 583       error ("if param_config is given, its potential items must not \
 584           be configured in another way");
 585     endif
 586
 587     ## supplement parameter names lacking in param_config
 588     nidx = ! isfield (pconf, pord);
 589     pconf = cell2fields ({struct()}(ones (1, sum (nidx))), \
 590                          pord(nidx), 2, pconf);
 591
 592     pconf = structcat (1, fields2cell (pconf, pord){:});
 593
 594     ## in the following, use reshape with explicit dimensions (instead
 595     ## of x(:)) so that errors are thrown if a configuration item has
 596     ## incorrect number of elements
 597
 598     lbound = - Infvec;
 599     if (isfield (pconf, "lbound"))
 600       idx = ! fieldempty (pconf, "lbound");
 601       if (pnonscalar)
 602         lbound (idx(prepidx), 1) = \
 603             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 604                              {pconf(idx).lbound}.', \
 605                              cpnel(idx), "UniformOutput", false){:});
 606       else
 607         lbound(idx, 1) = cat (1, pconf.lbound);
 608       endif
 609     endif
 610
 611     ubound = Infvec;
 612     if (isfield (pconf, "ubound"))
 613       idx = ! fieldempty (pconf, "ubound");
 614       if (pnonscalar)
 615         ubound (idx(prepidx), 1) = \
 616             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 617                              {pconf(idx).ubound}.', \
 618                              cpnel(idx), "UniformOutput", false){:});
 619       else
 620         ubound(idx, 1) = cat (1, pconf.ubound);
 621       endif
 622     endif
 623
 624     max_fract_change = fract_prec = NAvec;
 625
 626     if (isfield (pconf, "max_fract_change"))
 627       idx = ! fieldempty (pconf, "max_fract_change");
 628       if (pnonscalar)
 629         max_fract_change(idx(prepidx)) = \
 630             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 631                              {pconf(idx).max_fract_change}.', \
 632                              cpnel(idx), \
 633                              "UniformOutput", false){:});
 634       else
 635         max_fract_change(idx) = [pconf.max_fract_change];
 636       endif
 637     endif
 638
 639     if (isfield (pconf, "fract_prec"))
 640       idx = ! fieldempty (pconf, "fract_prec");
 641       if (pnonscalar)
 642         fract_prec(idx(prepidx)) = \
 643             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 644                              {pconf(idx).fract_prec}.', cpnel(idx), \
 645                              "UniformOutput", false){:});
 646       else
 647         fract_prec(idx) = [pconf.fract_prec];
 648       endif
 649     endif
 650
 651     diffp = zerosvec;
 652     diffp(:) = diffp_default;
 653     if (isfield (pconf, "diffp"))
 654       idx = ! fieldempty (pconf, "diffp");
 655       if (pnonscalar)
 656         diffp(idx(prepidx)) = \
 657             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 658                              {pconf(idx).diffp}.', cpnel(idx), \
 659                              "UniformOutput", false){:});
 660       else
 661         diffp(idx) = [pconf.diffp];
 662       endif
 663     endif
 664
 665     diff_onesided = fixed = falsevec;
 666
 667     if (isfield (pconf, "diff_onesided"))
 668       idx = ! fieldempty (pconf, "diff_onesided");
 669       if (pnonscalar)
 670         diff_onesided(idx(prepidx)) = \
 671             logical \
 672             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 673                               {pconf(idx).diff_onesided}.', cpnel(idx), \
 674                              "UniformOutput", false){:}));
 675       else
 676         diff_onesided(idx) = logical ([pconf.diff_onesided]);
 677       endif
 678     endif
 679
 680     if (isfield (pconf, "fixed"))
 681       idx = ! fieldempty (pconf, "fixed");
 682       if (pnonscalar)
 683         fixed(idx(prepidx)) = \
 684             logical \
 685             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 686                               {pconf(idx).fixed}.', cpnel(idx), \
 687                              "UniformOutput", false){:}));
 688       else
 689         fixed(idx) = logical ([pconf.fixed]);
 690       endif
 691     endif
 692
 693     max_rand_step = NAvec;
 694
 695     if (isfield (pconf, "max_rand_step"))
 696       idx = ! fieldempty (pconf, "max_rand_step");
 697       if (pnonscalar)
 698         max_rand_step(idx(prepidx)) = \
 699             logical \
 700             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 701                               {pconf(idx).max_rand_step}.',
 702                               cpnel(idx), \
 703                               "UniformOutput", false){:}));
 704       else
 705         max_rand_step(idx) = logical ([pconf.max_rand_step]);
 706       endif
 707     endif
 708
 709   else
 710     ## use supplied configuration vectors
 711
 712     if (isempty (lbound))
 713       lbound = - Infvec;
 714     elseif (any (size (lbound) != sizevec))
 715       error ("bounds: wrong dimensions");
 716     endif
 717
 718     if (isempty (ubound))
 719       ubound = Infvec;
 720     elseif (any (size (ubound) != sizevec))
 721       error ("bounds: wrong dimensions");
 722     endif
 723
 724     if (isempty (max_fract_change))
 725       max_fract_change = NAvec;
 726     elseif (any (size (max_fract_change) != sizevec))
 727       error ("max_fract_change: wrong dimensions");
 728     endif
 729
 730     if (isempty (fract_prec))
 731       fract_prec = NAvec;
 732     elseif (any (size (fract_prec) != sizevec))
 733       error ("fract_prec: wrong dimensions");
 734     endif
 735
 736     if (isempty (diffp))
 737       diffp = zerosvec;
 738       diffp(:) = diffp_default;
 739     else
 740       if (any (size (diffp) != sizevec))
 741         error ("diffp: wrong dimensions");
 742       endif
 743       diffp(isna (diffp)) = diffp_default;
 744     endif
 745
 746     if (isempty (diff_onesided))
 747       diff_onesided = falsevec;
 748     else
 749       if (any (size (diff_onesided) != sizevec))
 750         error ("diff_onesided: wrong dimensions")
 751       endif
 752       diff_onesided(isna (diff_onesided)) = false;
 753       diff_onesided = logical (diff_onesided);
 754     endif
 755
 756     if (isempty (fixed))
 757       fixed = falsevec;
 758     else
 759       if (any (size (fixed) != sizevec))
 760         error ("fixed: wrong dimensions");
 761       endif
 762       fixed(isna (fixed)) = false;
 763       fixed = logical (fixed);
 764     endif
 765
 766     if (isempty (max_rand_step))
 767       max_rand_step = NAvec;
 768     elseif (any (size (max_rand_step) != sizevec))
 769       error ("max_rand_step: wrong dimensions");
 770     endif
 771
 772   endif
 773
 774   ## guaranty all (lbound <= ubound)
 775   if (any (lbound > ubound))
 776     error ("some lower bounds larger than upper bounds");
 777   endif
 778
 779   #### consider whether initial parameters and functions are based on
 780   #### parameter structures or parameter vectors; wrappers for call to
 781   #### default function for jacobians
 782
 783   ## initial parameters
 784   if (pin_struct)
 785     if (pnonscalar)
 786       pin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 787                              fields2cell (pin, pord), cpnel, \
 788                              "UniformOutput", false){:});
 789     else
 790       pin = cat (1, fields2cell (pin, pord){:});
 791     endif
 792   endif
 793
 794   ## objective function
 795   if (f_pstruct)
 796     if (pnonscalar)
 797       f = @ (p, varargin) \
 798           f (cell2struct \
 799              (cellfun (@ reshape, mat2cell (p, ppartidx), \
 800                        pdims, "UniformOutput", false), \
 801               pord, 1), varargin{:});
 802     else
 803       f = @ (p, varargin) \
 804           f (cell2struct (num2cell (p), pord, 1), varargin{:});
 805     endif
 806   endif
 807   f_pin = f (pin);
 808
 809   ## gradient of objective function
 810   if (isempty (dfdp))
 811     if (do_cstep)
 812       dfdp = @ (p, hook) jacobs (p, f, hook);
 813     else
 814       dfdp = @ (p, hook) __dfdp__ (p, f, hook);
 815     endif
 816   endif
 817   if (dfdp_pstruct)
 818     if (pnonscalar)
 819       dfdp = @ (p, hook) \
 820           cat (2, \
 821                fields2cell \
 822                (dfdp (cell2struct \
 823                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 824                                 pdims, "UniformOutput", false), \
 825                        pord, 1), hook), \
 826                 pord){:});
 827     else
 828       dfdp = @ (p, hook) \
 829           cat (2, \
 830                fields2cell \
 831                (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
 832                 pord){:});
 833     endif
 834   endif
 835
 836   ## hessian of objective function
 837   if (hessian_pstruct)
 838     if (pnonscalar)
 839       hessian = @ (p) \
 840           hessian_struct2mat \
 841           (hessian (cell2struct \
 842                     (cellfun (@ reshape, mat2cell (p, ppartidx), \
 843                               pdims, "UniformOutput", false), \
 844                      pord, 1)), pord);
 845     else
 846       hessian = @ (p) \
 847           hessian_struct2mat \
 848           (hessian (cell2struct (num2cell (p), pord, 1)), pord);
 849     endif
 850   endif
 851
 852   ## function for general inequality constraints
 853   if (f_inequc_pstruct)
 854     if (pnonscalar)
 855       f_genicstr = @ (p, varargin) \
 856           f_genicstr (cell2struct \
 857                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 858                                 pdims, "UniformOutput", false), \
 859                        pord, 1), varargin{:});
 860     else
 861       f_genicstr = @ (p, varargin) \
 862           f_genicstr \
 863           (cell2struct (num2cell (p), pord, 1), varargin{:});
 864     endif
 865   endif
 866
 867   ## note this stage
 868   possibly_pstruct_f_genicstr = f_genicstr;
 869
 870   ## jacobian of general inequality constraints
 871   if (df_inequc_pstruct)
 872     if (pnonscalar)
 873       df_gencstr = @ (p, func, idx, hook) \
 874           cat (2, \
 875                fields2cell \
 876                (df_gencstr \
 877                 (cell2struct \
 878                  (cellfun (@ reshape, mat2cell (p, ppartidx), \
 879                            pdims, "UniformOutput", false), pord, 1), \
 880                  func, idx, hook), \
 881                 pord){:});
 882     else
 883       df_gencstr = @ (p, func, idx, hook) \
 884           cat (2, \
 885                fields2cell \
 886                (df_gencstr (cell2struct (num2cell (p), pord, 1), \
 887                             func, idx, hook), \
 888                 pord){:});
 889     endif
 890   endif
 891
 892   ## function for general equality constraints
 893   if (f_equc_pstruct)
 894     if (pnonscalar)
 895       f_genecstr = @ (p, varargin) \
 896           f_genecstr (cell2struct \
 897                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 898                                 pdims, "UniformOutput", false), \
 899                        pord, 1), varargin{:});
 900     else
 901       f_genecstr = @ (p, varargin) \
 902           f_genecstr \
 903           (cell2struct (num2cell (p), pord, 1), varargin{:});
 904     endif
 905   endif
 906
 907   ## note this stage
 908   possibly_pstruct_f_genecstr = f_genecstr;
 909
 910   ## jacobian of general equality constraints
 911   if (df_equc_pstruct)
 912     if (pnonscalar)
 913       df_genecstr = @ (p, func, idx, hook) \
 914           cat (2, \
 915                fields2cell \
 916                (df_genecstr \
 917                 (cell2struct \
 918                  (cellfun (@ reshape, mat2cell (p, ppartidx), \
 919                            pdims, "UniformOutput", false), pord, 1), \
 920                  func, idx, hook), \
 921                 pord){:});
 922     else
 923       df_genecstr = @ (p, func, idx, hook) \
 924           cat (2, \
 925                fields2cell \
 926                (df_genecstr (cell2struct (num2cell (p), pord, 1), \
 927                              func, idx, hook), \
 928                 pord){:});
 929     endif
 930   endif
 931
 932   ## linear inequality constraints
 933   if (mc_struct)
 934     idx = isfield (mc, pord);
 935     if (rows (fieldnames (mc)) > sum (idx))
 936       error ("unknown fields in structure of linear inequality constraints");
 937     endif
 938     smc = mc;
 939     mc = zeros (np, rows (vc));
 940     mc(idx(prepidx), :) = cat (1, fields2cell (smc, pord(idx)){:});
 941   endif
 942
 943   ## linear equality constraints
 944   if (emc_struct)
 945     idx = isfield (emc, pord);
 946     if (rows (fieldnames (emc)) > sum (idx))
 947       error ("unknown fields in structure of linear equality constraints");
 948     endif
 949     semc = emc;
 950     emc = zeros (np, rows (evc));
 951     emc(idx(prepidx), :) = cat (1, fields2cell (semc, pord(idx)){:});
 952   endif
 953
 954   ## parameter-related configuration for jacobian functions
 955   if (dfdp_pstruct || df_inequc_pstruct || df_equc_pstruct)
 956     if(pnonscalar)
 957       s_diffp = cell2struct \
 958           (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
 959                     pdims, "UniformOutput", false), pord, 1);
 960       s_diff_onesided = cell2struct \
 961           (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
 962                     pdims, "UniformOutput", false), pord, 1);
 963       s_orig_lbound = cell2struct \
 964           (cellfun (@ reshape, mat2cell (lbound, ppartidx), \
 965                     pdims, "UniformOutput", false), pord, 1);
 966       s_orig_ubound = cell2struct \
 967           (cellfun (@ reshape, mat2cell (ubound, ppartidx), \
 968                     pdims, "UniformOutput", false), pord, 1);
 969       s_plabels = cell2struct \
 970           (num2cell \
 971            (cat (2, cellfun \
 972                  (@ (x) cellfun \
 973                   (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
 974                    pdims, "UniformOutput", false), \
 975                   num2cell (plabels, 1), "UniformOutput", false){:}), \
 976             2), \
 977            pord, 1);
 978       s_orig_fixed = cell2struct \
 979           (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
 980                     pdims, "UniformOutput", false), pord, 1);
 981     else
 982       s_diffp = cell2struct (num2cell (diffp), pord, 1);
 983       s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
 984       s_orig_lbound = cell2struct (num2cell (lbound), pord, 1);
 985       s_orig_ubound = cell2struct (num2cell (ubound), pord, 1);
 986       s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
 987       s_orig_fixed = cell2struct (num2cell (fixed), pord, 1);
 988     endif
 989   endif
 990
 991   #### some further values and checks
 992
 993   if (any (fixed & (pin < lbound | pin > ubound)))
 994     warning ("some fixed parameters outside bounds");
 995   endif
 996
 997   ## dimensions of linear constraints
 998   if (isempty (mc))
 999     mc = zeros (np, 0);
1000     vc = zeros (0, 1);
1001   endif
1002   if (isempty (emc))
1003     emc = zeros (np, 0);
1004     evc = zeros (0, 1);
1005   endif
1006   [rm, cm] = size (mc);
1007   [rv, cv] = size (vc);
1008   if (rm != np || cm != rv || cv != 1)
1009     error ("linear inequality constraints: wrong dimensions");
1010   endif
1011   [erm, ecm] = size (emc);
1012   [erv, ecv] = size (evc);
1013   if (erm != np || ecm != erv || ecv != 1)
1014     error ("linear equality constraints: wrong dimensions");
1015   endif
1016
1017   ## note initial values of linear constraits
1018   pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
1019   pin_cstr.equ.lin = emc.' * pin + evc;
1020
1021   ## note number and initial values of general constraints
1022   if (isempty (f_genicstr))
1023     pin_cstr.inequ.gen = [];
1024     n_genicstr = 0;
1025   else
1026     n_genicstr = length (pin_cstr.inequ.gen = f_genicstr (pin));
1027   endif
1028   if (isempty (f_genecstr))
1029     pin_cstr.equ.gen = [];
1030     n_genecstr = 0;
1031   else
1032     n_genecstr = length (pin_cstr.equ.gen = f_genecstr (pin));
1033   endif
1034
1035   #### collect remaining settings
1036   hook.TolFun = optimget (settings, "TolFun", stol_default);
1037   hook.MaxIter = optimget (settings, "MaxIter");
1038   if (ischar (hook.cpiv = optimget (settings, "cpiv", @ cpiv_bard)))
1039     hook.cpiv = str2func (hook.cpiv);
1040   endif
1041   hook.Display = optimget (settings, "Display", "off");
1042   hook.testing = optimget (settings, "debug", false);
1043   hook.siman.T_init = optimget (settings, "T_init", .01);
1044   hook.siman.T_min = optimget (settings, "T_min", 1.0e-5);
1045   hook.siman.mu_T = optimget (settings, "mu_T", 1.005);
1046   hook.siman.iters_fixed_T = optimget (settings, "iters_fixed_T", 10);
1047   hook.stoch_regain_constr = \
1048       optimget (settings, "stoch_regain_constr", false);
1049   hook.trace_steps = \
1050       optimget (settings, "trace_steps", false);
1051   hook.siman_log = \
1052       optimget (settings, "siman_log", false);
1053   backend = optimget (settings, "Algorithm", "lm_feasible");
1054   backend = map_matlab_algorithm_names (backend);
1055   backend = map_backend (backend);
1056
1057   #### handle fixing of parameters
1058   orig_lbound = lbound;
1059   orig_ubound = ubound;
1060   orig_fixed = fixed;
1061   if (all (fixed))
1062     error ("no free parameters");
1063   endif
1064
1065   nonfixed = ! fixed;
1066   if (any (fixed))
1067     ## backend (returned values and initial parameters)
1068     backend = @ (f, pin, hook) \
1069         backend_wrapper (backend, fixed, f, pin, hook);
1070
1071     ## objective function
1072     f = @ (p, varargin) f (assign (pin, nonfixed, p), varargin{:});
1073
1074     ## gradient of objective function
1075     dfdp = @ (p, hook) \
1076         dfdp (assign (pin, nonfixed, p), hook)(nonfixed);
1077
1078     ## hessian of objective function
1079     if (! isempty (hessian))
1080       hessian = @ (p) \
1081           hessian (assign (pin, nonfixed, p))(nonfixed, nonfixed);
1082     endif
1083
1084     ## function for general inequality constraints
1085     f_genicstr = @ (p, varargin) \
1086         f_genicstr (assign (pin, nonfixed, p), varargin{:});
1087
1088     ## jacobian of general inequality constraints
1089     df_gencstr = @ (p, func, idx, hook) \
1090         df_gencstr (assign (pin, nonfixed, p), func, idx, hook) \
1091         (:, nonfixed);
1092
1093     ## function for general equality constraints
1094     f_genecstr = @ (p, varargin) \
1095         f_genecstr (assign (pin, nonfixed, p), varargin{:});
1096
1097     ## jacobian of general equality constraints
1098     df_genecstr = @ (p, func, idx, hook) \
1099         df_genecstr (assign (pin, nonfixed, p), func, idx, hook) \
1100         (:, nonfixed);
1101
1102     ## linear inequality constraints
1103     vc += mc(fixed, :).' * (tp = pin(fixed));
1104     mc = mc(nonfixed, :);
1105
1106     ## linear equality constraints
1107     evc += emc(fixed, :).' * tp;
1108     emc = emc(nonfixed, :);
1109
1110     ## _last_ of all, vectors of parameter-related configuration,
1111     ## including "fixed" itself
1112     lbound = lbound(nonfixed, :);
1113     ubound = ubound(nonfixed, :);
1114     max_fract_change = max_fract_change(nonfixed);
1115     fract_prec = fract_prec(nonfixed);
1116     max_rand_step = max_rand_step(nonfixed);
1117     fixed = fixed(nonfixed);
1118   endif
1119
1120   #### supplement constants to jacobian functions
1121
1122   ## gradient of objective function
1123   if (dfdp_pstruct)
1124     dfdp = @ (p, hook) \
1125         dfdp (p, cell2fields \
1126               ({s_diffp, s_diff_onesided, s_orig_lbound, \
1127                 s_orig_ubound, s_plabels, \
1128                 cell2fields(num2cell(hook.fixed), pord(nonfixed), \
1129                             1, s_orig_fixed), cstep}, \
1130                {"diffp", "diff_onesided", "lbound", "ubound", \
1131                 "plabels", "fixed", "h"}, \
1132                2, hook));
1133   else
1134     dfdp = @ (p, hook) \
1135         dfdp (p, cell2fields \
1136               ({diffp, diff_onesided, orig_lbound, orig_ubound, \
1137                 plabels, assign(orig_fixed, nonfixed, hook.fixed), \
1138                 cstep}, \
1139                {"diffp", "diff_onesided", "lbound", "ubound", \
1140                 "plabels", "fixed", "h"}, \
1141                2, hook));
1142   endif
1143
1144   ## jacobian of general inequality constraints
1145   if (df_inequc_pstruct)
1146     df_gencstr = @ (p, func, idx, hook) \
1147         df_gencstr (p, func, idx, cell2fields \
1148                     ({s_diffp, s_diff_onesided, s_orig_lbound, \
1149                       s_orig_ubound, s_plabels, \
1150                       cell2fields(num2cell(hook.fixed), pord(nonfixed), \
1151                                   1, s_orig_fixed), cstep}, \
1152                      {"diffp", "diff_onesided", "lbound", "ubound", \
1153                       "plabels", "fixed", "h"}, \
1154                      2, hook));
1155   else
1156     df_gencstr = @ (p, func, idx, hook) \
1157         df_gencstr (p, func, idx, cell2fields \
1158                     ({diffp, diff_onesided, orig_lbound, \
1159                       orig_ubound, plabels, \
1160                       assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
1161                      {"diffp", "diff_onesided", "lbound", "ubound", \
1162                       "plabels", "fixed", "h"}, \
1163                      2, hook));
1164   endif
1165
1166   ## jacobian of general equality constraints
1167   if (df_equc_pstruct)
1168     df_genecstr = @ (p, func, idx, hook) \
1169         df_genecstr (p, func, idx, cell2fields \
1170                      ({s_diffp, s_diff_onesided, s_orig_lbound, \
1171                        s_orig_ubound, s_plabels, \
1172                        cell2fields(num2cell(hook.fixed), pord(nonfixed), \
1173                                    1, s_orig_fixed), cstep}, \
1174                       {"diffp", "diff_onesided", "lbound", "ubound", \
1175                        "plabels", "fixed", "h"}, \
1176                       2, hook));
1177   else
1178     df_genecstr = @ (p, func, idx, hook) \
1179         df_genecstr (p, func, idx, cell2fields \
1180                      ({diffp, diff_onesided, orig_lbound, \
1181                        orig_ubound, plabels, \
1182                        assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
1183                       {"diffp", "diff_onesided", "lbound", "ubound", \
1184                        "plabels", "fixed", "h"}, \
1185                       2, hook));
1186   endif
1187
1188   #### interfaces to constraints
1189
1190   ## include bounds into linear inequality constraints
1191   tp = eye (sum (nonfixed));
1192   lidx = lbound != - Inf;
1193   uidx = ubound != Inf;
1194   mc = cat (2, tp(:, lidx), - tp(:, uidx), mc);
1195   vc = cat (1, - lbound(lidx, 1), ubound(uidx, 1), vc);
1196
1197   ## concatenate linear inequality and equality constraints
1198   mc = cat (2, mc, emc);
1199   vc = cat (1, vc, evc);
1200   n_lincstr = rows (vc);
1201
1202   ## concatenate general inequality and equality constraints
1203   if (n_genecstr > 0)
1204     if (n_genicstr > 0)
1205       nidxi = 1 : n_genicstr;
1206       nidxe = n_genicstr + 1 : n_genicstr + n_genecstr;
1207       f_gencstr = @ (p, idx, varargin) \
1208           cat (1, \
1209                f_genicstr (p, idx(nidxi), varargin{:}), \
1210                f_genecstr (p, idx(nidxe), varargin{:}));
1211       df_gencstr = @ (p, idx, hook) \
1212           cat (1, \
1213                df_gencstr (p, @ (p, varargin) \
1214                            possibly_pstruct_f_genicstr \
1215                            (p, idx(nidxi), varargin{:}), \
1216                            idx(nidxi), \
1217                            setfield (hook, "f", \
1218                                      hook.f(nidxi(idx(nidxi))))), \
1219                df_genecstr (p, @ (p, varargin) \
1220                             possibly_pstruct_f_genecstr \
1221                             (p, idx(nidxe), varargin{:}), \
1222                             idx(nidxe), \
1223                             setfield (hook, "f", \
1224                                       hook.f(nidxe(idx(nidxe))))));
1225     else
1226       f_gencstr = f_genecstr;
1227       df_gencstr = @ (p, idx, hook) \
1228           df_genecstr (p, \
1229                        @ (p, varargin) \
1230                        possibly_pstruct_f_genecstr \
1231                        (p, idx, varargin{:}), \
1232                        idx, \
1233                        setfield (hook, "f", hook.f(idx)));
1234     endif
1235   else
1236     f_gencstr = f_genicstr;
1237     df_gencstr = @ (p, idx, hook) \
1238         df_gencstr (p, \
1239                     @ (p, varargin) \
1240                     possibly_pstruct_f_genicstr (p, idx, varargin{:}), \
1241                     idx, \
1242                     setfield (hook, "f", hook.f(idx)));
1243   endif
1244   n_gencstr = n_genicstr + n_genecstr;
1245
1246   ## concatenate linear and general constraints, defining the final
1247   ## function interfaces
1248   if (n_gencstr > 0)
1249     nidxl = 1:n_lincstr;
1250     nidxh = n_lincstr + 1 : n_lincstr + n_gencstr;
1251     f_cstr = @ (p, idx, varargin) \
1252         cat (1, \
1253              mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), \
1254              f_gencstr (p, idx(nidxh), varargin{:}));
1255     df_cstr = @ (p, idx, hook) \
1256         cat (1, \
1257              mc(:, idx(nidxl)).', \
1258              df_gencstr (p, idx(nidxh), \
1259                          setfield (hook, "f", \
1260                                    hook.f(nidxh))));
1261   else
1262     f_cstr = @ (p, idx, varargin) mc(:, idx).' * p + vc(idx, 1);
1263     df_cstr = @ (p, idx, hook) mc(:, idx).';
1264   endif
1265
1266   ## define eq_idx (logical index of equality constraints within all
1267   ## concatenated constraints
1268   eq_idx = false (n_lincstr + n_gencstr, 1);
1269   eq_idx(n_lincstr + 1 - rows (evc) : n_lincstr) = true;
1270   n_cstr = n_lincstr + n_gencstr;
1271   eq_idx(n_cstr + 1 - n_genecstr : n_cstr) = true;
1272
1273   #### prepare interface hook
1274
1275   ## passed constraints
1276   hook.mc = mc;
1277   hook.vc = vc;
1278   hook.f_cstr = f_cstr;
1279   hook.df_cstr = df_cstr;
1280   hook.n_gencstr = n_gencstr;
1281   hook.eq_idx = eq_idx;
1282   hook.lbound = lbound;
1283   hook.ubound = ubound;
1284
1285   ## passed values of constraints for initial parameters
1286   hook.pin_cstr = pin_cstr;
1287
1288   ## passed function for gradient of objective function
1289   hook.dfdp = dfdp;
1290
1291   ## passed function for hessian of objective function
1292   hook.hessian = hessian;
1293
1294   ## passed function for complementary pivoting
1295   ## hook.cpiv = cpiv; # set before
1296
1297   ## passed value of objective function for initial parameters
1298   hook.f_pin = f_pin;
1299
1300   ## passed options
1301   hook.max_fract_change = max_fract_change;
1302   hook.fract_prec = fract_prec;
1303   ## hook.TolFun = ; # set before
1304   ## hook.MaxIter = ; # set before
1305   hook.fixed = fixed;
1306   hook.max_rand_step = max_rand_step;
1307
1308   #### call backend
1309
1310   [p, objf, cvg, outp] = backend (f, pin, hook);
1311
1312   if (pin_struct)
1313     if (pnonscalar)
1314       p = cell2struct \
1315           (cellfun (@ reshape, mat2cell (p, ppartidx), \
1316                     pdims, "UniformOutput", false), \
1317            pord, 1);
1318     else
1319       p = cell2struct (num2cell (p), pord, 1);
1320     endif
1321   endif
1322
1323 endfunction
1324
1325 function backend = map_matlab_algorithm_names (backend)
1326
1327   ## nothing done here at the moment
1328
1329 endfunction
1330
1331 function backend = map_backend (backend)
1332
1333   switch (backend)
1334       ##    case "sqp_infeasible"
1335       ##      backend = "__sqp__";
1336       ##    case "sqp"
1337       ##      backend = "__sqp__";
1338     case "lm_feasible"
1339       backend = "__lm_feasible__";
1340     case "siman"
1341       backend = "__siman__";
1342     otherwise
1343       error ("no backend implemented for algorithm '%s'", backend);
1344   endswitch
1345
1346   backend = str2func (backend);
1347
1348 endfunction
1349
1350 function [p, resid, cvg, outp] = backend_wrapper (backend, fixed, f, p, hook)
1351
1352   [tp, resid, cvg, outp] = backend (f, p(! fixed), hook);
1353
1354   p(! fixed) = tp;
1355
1356 endfunction
1357
1358 function lval = assign (lval, lidx, rval)
1359
1360   lval(lidx) = rval;
1361
1362 endfunction
1363
1364 function m = hessian_struct2mat (s, pord)
1365
1366   m = cell2mat (fields2cell \
1367                 (structcat (1, NA, fields2cell (s, pord){:}), pord));
1368
1369   idx = isna (m);
1370
1371   m(idx) = (m.')(idx);
1372
1373 endfunction
1374
1375 function ret = __optimget__ (s, name, default)
1376
1377   if (isfield (s, name))
1378     ret = s.(name);
1379   elseif (nargin > 2)
1380     ret = default;
1381   else
1382     ret = [];
1383   endif
1384
1385 endfunction
1386
1387 function ret = apply_idx_if_given  (ret, varargin)
1388
1389   if (nargin > 1)
1390     ret = ret(varargin{1});
1391   endif
1392
1393 endfunction
1394
1395 %!demo
1396 %! ## Example for default optimization (Levenberg/Marquardt with
1397 %! ## BFGS), one non-linear equality constraint. Constrained optimum is
1398 %! ## at p = [0; 1].
1399 %! objective_function = @ (p) p(1)^2 + p(2)^2;
1400 %! pin = [-2; 5];
1401 %! constraint_function = @ (p) p(1)^2 + 1 - p(2);
1402 %! [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset ("equc", {constraint_function}))
1403
1404 %!demo
1405 %! ## Example for simulated annealing, two parameters, "trace_steps"
1406 %! ## is true;
1407 %! t_init = .2;
1408 %! t_min = .002;
1409 %! mu_t = 1.002;
1410 %! iters_fixed_t = 10;
1411 %! init_p = [2; 2];
1412 %! max_rand_step = [.2; .2];
1413 %! [p, objf, cvg, outp] = nonlin_min (@ (p) (p(1)/10)^2 + (p(2)/10)^2 + .1 * (-cos(4*p(1)) - cos(4*p(2))), init_p, optimset ("algorithm", "siman", "max_rand_step", max_rand_step, "t_init", t_init, "T_min", t_min, "mu_t", mu_t, "iters_fixed_T", iters_fixed_t, "trace_steps", true));
1414 %! p
1415 %! objf
1416 %! x = (outp.trace(:, 1) - 1) * iters_fixed_t + outp.trace(:, 2);
1417 %! x(1) = 0;
1418 %! plot (x, cat (2, outp.trace(:, 3:end), t_init ./ (mu_t .^ outp.trace(:, 1))))
1419 %! legend ({"objective function value", "p(1)", "p(2)", "Temperature"})
1420 %! xlabel ("subiteration")