octave_packages/optim-1.2.0/private/__nonlin_residmin__.m

   1 ## Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.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 ## Internal function, called by nonlin_residmin --- see there --- and
  17 ## others. Calling __nonlin_residmin__ indirectly hides the argument
  18 ## "hook", usable by wrappers, from users. Currently, hook can contain
  19 ## the field "observations", so that dimensions of observations and
  20 ## returned values of unchanged model function can be checked against
  21 ## each other exactly one time.
  22
  23 ## disabled PKG_ADD: __all_opts__ ("__nonlin_residmin__");
  24
  25 function [p, resid, cvg, outp] = \
  26       __nonlin_residmin__ (f, pin, settings, hook)
  27
  28   if (compare_versions (version (), "3.3.55", "<"))
  29     ## optimset mechanism was fixed for option names with underscores
  30     ## sometime in 3.3.54+, if I remember right
  31     optimget = @ __optimget__;
  32   endif
  33
  34   if (compare_versions (version (), "3.2.4", "<="))
  35     ## For bug #31484; but Octave 3.6... shows bug #36288 due to this
  36     ## workaround. Octave 3.7... seems to be all right.
  37     __dfdp__ = @ __dfdp__;
  38   endif
  39
  40   ## some scalar defaults; some defaults are backend specific, so
  41   ## lacking elements in respective constructed vectors will be set to
  42   ## NA here in the frontend
  43   diffp_default = .001;
  44   stol_default = .0001;
  45   cstep_default = 1e-20;
  46
  47   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
  48     p = optimset ("param_config", [], \
  49                   "param_order", [], \
  50                   "param_dims", [], \
  51                   "f_inequc_pstruct", false, \
  52                   "f_equc_pstruct", false, \
  53                   "f_pstruct", false, \
  54                   "df_inequc_pstruct", false, \
  55                   "df_equc_pstruct", false, \
  56                   "dfdp_pstruct", false, \
  57                   "lbound", [], \
  58                   "ubound", [], \
  59                   "dfdp", [], \
  60                   "cpiv", @ cpiv_bard, \
  61                   "max_fract_change", [], \
  62                   "fract_prec", [], \
  63                   "diffp", [], \
  64                   "diff_onesided", [], \
  65                   "complex_step_derivative_f", false, \
  66                   "complex_step_derivative_inequc", false, \
  67                   "complex_step_derivative_equc", false, \
  68                   "cstep", cstep_default, \
  69                   "fixed", [], \
  70                   "inequc", [], \
  71                   "equc", [], \
  72                   "inequc_f_idx", false, \
  73                   "inequc_df_idx", false, \
  74                   "equc_f_idx", false, \
  75                   "equc_df_idx", false, \
  76                   "weights", [], \
  77                   "TolFun", stol_default, \
  78                   "MaxIter", [], \
  79                   "Display", "off", \
  80                   "Algorithm", "lm_svd_feasible", \
  81                   "plot_cmd", @ (f) 0, \
  82                   "debug", false, \
  83                   "lm_svd_feasible_alt_s", false);
  84     return;
  85   endif
  86
  87   assign = @ assign; # Is this faster in repeated calls?
  88
  89   if (nargin != 4)
  90     error ("incorrect number of arguments");
  91   endif
  92
  93   if (ischar (f))
  94     f = str2func (f);
  95   endif
  96
  97   if (! (pin_struct = isstruct (pin)))
  98     if (! isvector (pin) || columns (pin) > 1)
  99       error ("initial parameters must be either a structure or a column vector");
 100     endif
 101   endif
 102
 103   #### processing of settings and consistency checks
 104
 105   pconf = optimget (settings, "param_config");
 106   pord = optimget (settings, "param_order");
 107   pdims = optimget (settings, "param_dims");
 108   f_inequc_pstruct = optimget (settings, "f_inequc_pstruct", false);
 109   f_equc_pstruct = optimget (settings, "f_equc_pstruct", false);
 110   f_pstruct = optimget (settings, "f_pstruct", false);
 111   dfdp_pstruct = optimget (settings, "dfdp_pstruct", f_pstruct);
 112   df_inequc_pstruct = optimget (settings, "df_inequc_pstruct", \
 113                                 f_inequc_pstruct);
 114   df_equc_pstruct = optimget (settings, "df_equc_pstruct", \
 115                               f_equc_pstruct);
 116   lbound = optimget (settings, "lbound");
 117   ubound = optimget (settings, "ubound");
 118   dfdp = optimget (settings, "dfdp");
 119   if (ischar (dfdp)) dfdp = str2func (dfdp); endif
 120   max_fract_change = optimget (settings, "max_fract_change");
 121   fract_prec = optimget (settings, "fract_prec");
 122   diffp = optimget (settings, "diffp");
 123   diff_onesided = optimget (settings, "diff_onesided");
 124   fixed = optimget (settings, "fixed");
 125   do_cstep = optimget (settings, "complex_step_derivative_f", false);
 126   cstep = optimget (settings, "cstep", cstep_default);
 127   if (do_cstep && ! isempty (dfdp))
 128     error ("both 'complex_step_derivative_f' and 'dfdp' are set");
 129   endif
 130   do_cstep_inequc = \
 131       optimget (settings, "complex_step_derivative_inequc", false);
 132   do_cstep_equc = optimget (settings, "complex_step_derivative_equc", false);
 133
 134   any_vector_conf = ! (isempty (lbound) && isempty (ubound) && \
 135                        isempty (max_fract_change) && \
 136                        isempty (fract_prec) && isempty (diffp) && \
 137                        isempty (diff_onesided) && isempty (fixed));
 138
 139   ## collect constraints
 140   [mc, vc, f_genicstr, df_gencstr, user_df_gencstr] = \
 141       __collect_constraints__ (optimget (settings, "inequc"), \
 142                                do_cstep_inequc, "inequality constraints");
 143   [emc, evc, f_genecstr, df_genecstr, user_df_genecstr] = \
 144       __collect_constraints__ (optimget (settings, "equc"), \
 145                                do_cstep_equc, "equality constraints");
 146   mc_struct = isstruct (mc);
 147   emc_struct = isstruct (emc);
 148
 149   ## correct "_pstruct" settings if functions are not supplied, handle
 150   ## constraint functions not honoring indices
 151   if (isempty (dfdp)) dfdp_pstruct = false; endif
 152   if (isempty (f_genicstr))
 153     f_inequc_pstruct = false;
 154   elseif (! optimget (settings, "inequc_f_idx", false))
 155     f_genicstr = @ (p, varargin) apply_idx_if_given \
 156         (f_genicstr (p, varargin{:}), varargin{:});
 157   endif
 158   if (isempty (f_genecstr))
 159     f_equc_pstruct = false;
 160   elseif (! optimget (settings, "equc_f_idx", false))
 161     f_genecstr = @ (p, varargin) apply_idx_if_given \
 162         (f_genecstr (p, varargin{:}), varargin{:});
 163   endif
 164   if (user_df_gencstr)
 165     if (! optimget (settings, "inequc_df_idx", false))
 166       df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
 167     endif
 168   else
 169     df_inequc_pstruct = false;
 170   endif
 171   if (user_df_genecstr)
 172     if (! optimget (settings, "equc_df_idx", false))
 173       df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
 174     endif
 175   else
 176     df_equc_pstruct = false;
 177   endif
 178
 179   ## some settings require a parameter order
 180   if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
 181       f_equc_pstruct || f_pstruct || dfdp_pstruct || \
 182       df_inequc_pstruct || df_equc_pstruct || mc_struct || \
 183       emc_struct)
 184     if (isempty (pord))
 185       if (pin_struct)
 186         if (any_vector_conf || \
 187             ! (f_pstruct && \
 188                (f_inequc_pstruct || isempty (f_genicstr)) && \
 189                (f_equc_pstruct || isempty (f_genecstr)) && \
 190                (dfdp_pstruct || isempty (dfdp)) && \
 191                (df_inequc_pstruct || ! user_df_gencstr) && \
 192                (df_equc_pstruct || ! user_df_genecstr) && \
 193                (mc_struct || isempty (mc)) && \
 194                (emc_struct || isempty (emc))))
 195           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");
 196         else
 197           pord = fieldnames (pin);
 198         endif
 199       else
 200         error ("given settings require specification of parameter order or initial parameters in the form of a structure");
 201       endif
 202     endif
 203     pord = pord(:);
 204     if (pin_struct && ! all (isfield (pin, pord)))
 205       error ("some initial parameters lacking");
 206     endif
 207     if ((nnames = rows (unique (pord))) < rows (pord))
 208       error ("duplicate parameter names in 'param_order'");
 209     endif
 210     if (isempty (pdims))
 211       if (pin_struct)
 212         pdims = cellfun \
 213             (@ size, fields2cell (pin, pord), "UniformOutput", false);
 214       else
 215         pdims = num2cell (ones (nnames, 2), 2);
 216       endif
 217     else
 218       pdims = pdims(:);
 219       if (pin_struct && \
 220           ! all (cellfun (@ (x, y) prod (size (x)) == prod (y), \
 221                           struct2cell (pin), pdims)))
 222         error ("given param_dims and dimensions of initial parameters do not match");
 223       endif
 224     endif
 225     if (nnames != rows (pdims))
 226       error ("lengths of 'param_order' and 'param_dims' not equal");
 227     endif
 228     pnel = cellfun (@ prod, pdims);
 229     ppartidx = pnel;
 230     if (any (pnel > 1))
 231       pnonscalar = true;
 232       cpnel = num2cell (pnel);
 233       prepidx = cat (1, cellfun \
 234                      (@ (x, n) x(ones (1, n), 1), \
 235                       num2cell ((1:nnames).'), cpnel, \
 236                       "UniformOutput", false){:});
 237       epord = pord(prepidx, 1);
 238       psubidx = cat (1, cellfun \
 239                      (@ (n) (1:n).', cpnel, \
 240                       "UniformOutput", false){:});
 241     else
 242       pnonscalar = false; # some less expensive interfaces later
 243       prepidx = (1:nnames).';
 244       epord = pord;
 245       psubidx = ones (nnames, 1);
 246     endif
 247   else
 248     pord = []; # spares checks for given but not needed
 249   endif
 250
 251   if (pin_struct)
 252     np = sum (pnel);
 253   else
 254     np = length (pin);
 255     if (! isempty (pord) && np != sum (pnel))
 256       error ("number of initial parameters not correct");
 257     endif
 258   endif
 259
 260   plabels = num2cell (num2cell ((1:np).'));
 261   if (! isempty (pord))
 262     plabels = cat (2, plabels, num2cell (epord), \
 263                    num2cell (num2cell (psubidx)));
 264   endif
 265
 266   ## some useful vectors
 267   zerosvec = zeros (np, 1);
 268   NAvec = NA (np, 1);
 269   Infvec = Inf (np, 1);
 270   falsevec = false (np, 1);
 271   sizevec = [np, 1];
 272
 273   ## collect parameter-related configuration
 274   if (! isempty (pconf))
 275     ## use supplied configuration structure
 276
 277     ## parameter-related configuration is either allowed by a structure
 278     ## or by vectors
 279     if (any_vector_conf)
 280       error ("if param_config is given, its potential items must not \
 281           be configured in another way");
 282     endif
 283
 284     ## supplement parameter names lacking in param_config
 285     nidx = ! isfield (pconf, pord);
 286     pconf = cell2fields ({struct()}(ones (1, sum (nidx))), \
 287                          pord(nidx), 2, pconf);
 288
 289     pconf = structcat (1, fields2cell (pconf, pord){:});
 290
 291     ## in the following, use reshape with explicit dimensions (instead
 292     ## of x(:)) so that errors are thrown if a configuration item has
 293     ## incorrect number of elements
 294
 295     lbound = - Infvec;
 296     if (isfield (pconf, "lbound"))
 297       idx = ! fieldempty (pconf, "lbound");
 298       if (pnonscalar)
 299         lbound (idx(prepidx), 1) = \
 300             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 301                              {pconf(idx).lbound}.', \
 302                              cpnel(idx), "UniformOutput", false){:});
 303       else
 304         lbound(idx, 1) = cat (1, pconf.lbound);
 305       endif
 306     endif
 307
 308     ubound = Infvec;
 309     if (isfield (pconf, "ubound"))
 310       idx = ! fieldempty (pconf, "ubound");
 311       if (pnonscalar)
 312         ubound (idx(prepidx), 1) = \
 313             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 314                              {pconf(idx).ubound}.', \
 315                              cpnel(idx), "UniformOutput", false){:});
 316       else
 317         ubound(idx, 1) = cat (1, pconf.ubound);
 318       endif
 319     endif
 320
 321     max_fract_change = fract_prec = NAvec;
 322
 323     if (isfield (pconf, "max_fract_change"))
 324       idx = ! fieldempty (pconf, "max_fract_change");
 325       if (pnonscalar)
 326         max_fract_change(idx(prepidx)) = \
 327             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 328                              {pconf(idx).max_fract_change}.', \
 329                              cpnel(idx), \
 330                              "UniformOutput", false){:});
 331       else
 332         max_fract_change(idx) = [pconf.max_fract_change];
 333       endif
 334     endif
 335
 336     if (isfield (pconf, "fract_prec"))
 337       idx = ! fieldempty (pconf, "fract_prec");
 338       if (pnonscalar)
 339         fract_prec(idx(prepidx)) = \
 340             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 341                              {pconf(idx).fract_prec}.', cpnel(idx), \
 342                              "UniformOutput", false){:});
 343       else
 344         fract_prec(idx) = [pconf.fract_prec];
 345       endif
 346     endif
 347
 348     diffp = zerosvec;
 349     diffp(:) = diffp_default;
 350     if (isfield (pconf, "diffp"))
 351       idx = ! fieldempty (pconf, "diffp");
 352       if (pnonscalar)
 353         diffp(idx(prepidx)) = \
 354             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 355                              {pconf(idx).diffp}.', cpnel(idx), \
 356                              "UniformOutput", false){:});
 357       else
 358         diffp(idx) = [pconf.diffp];
 359       endif
 360     endif
 361
 362     diff_onesided = fixed = falsevec;
 363
 364     if (isfield (pconf, "diff_onesided"))
 365       idx = ! fieldempty (pconf, "diff_onesided");
 366       if (pnonscalar)
 367         diff_onesided(idx(prepidx)) = \
 368             logical \
 369             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 370                               {pconf(idx).diff_onesided}.', cpnel(idx), \
 371                              "UniformOutput", false){:}));
 372       else
 373         diff_onesided(idx) = logical ([pconf.diff_onesided]);
 374       endif
 375     endif
 376
 377     if (isfield (pconf, "fixed"))
 378       idx = ! fieldempty (pconf, "fixed");
 379       if (pnonscalar)
 380         fixed(idx(prepidx)) = \
 381             logical \
 382             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 383                               {pconf(idx).fixed}.', cpnel(idx), \
 384                              "UniformOutput", false){:}));
 385       else
 386         fixed(idx) = logical ([pconf.fixed]);
 387       endif
 388     endif
 389   else
 390     ## use supplied configuration vectors
 391
 392     if (isempty (lbound))
 393       lbound = - Infvec;
 394     elseif (any (size (lbound) != sizevec))
 395       error ("bounds: wrong dimensions");
 396     endif
 397
 398     if (isempty (ubound))
 399       ubound = Infvec;
 400     elseif (any (size (ubound) != sizevec))
 401       error ("bounds: wrong dimensions");
 402     endif
 403
 404     if (isempty (max_fract_change))
 405       max_fract_change = NAvec;
 406     elseif (any (size (max_fract_change) != sizevec))
 407       error ("max_fract_change: wrong dimensions");
 408     endif
 409
 410     if (isempty (fract_prec))
 411       fract_prec = NAvec;
 412     elseif (any (size (fract_prec) != sizevec))
 413       error ("fract_prec: wrong dimensions");
 414     endif
 415
 416     if (isempty (diffp))
 417       diffp = zerosvec;
 418       diffp(:) = diffp_default;
 419     else
 420       if (any (size (diffp) != sizevec))
 421         error ("diffp: wrong dimensions");
 422       endif
 423       diffp(isna (diffp)) = diffp_default;
 424     endif
 425
 426     if (isempty (diff_onesided))
 427       diff_onesided = falsevec;
 428     else
 429       if (any (size (diff_onesided) != sizevec))
 430         error ("diff_onesided: wrong dimensions")
 431       endif
 432       diff_onesided(isna (diff_onesided)) = false;
 433       diff_onesided = logical (diff_onesided);
 434     endif
 435
 436     if (isempty (fixed))
 437       fixed = falsevec;
 438     else
 439       if (any (size (fixed) != sizevec))
 440         error ("fixed: wrong dimensions");
 441       endif
 442       fixed(isna (fixed)) = false;
 443       fixed = logical (fixed);
 444     endif
 445   endif
 446
 447   ## guaranty all (lbound <= ubound)
 448   if (any (lbound > ubound))
 449     error ("some lower bounds larger than upper bounds");
 450   endif
 451
 452   #### consider whether initial parameters and functions are based on
 453   #### parameter structures or parameter vectors; wrappers for call to
 454   #### default function for jacobians
 455
 456   ## initial parameters
 457   if (pin_struct)
 458     if (pnonscalar)
 459       pin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 460                              fields2cell (pin, pord), cpnel, \
 461                              "UniformOutput", false){:});
 462     else
 463       pin = cat (1, fields2cell (pin, pord){:});
 464     endif
 465   endif
 466
 467   ## model function
 468   if (f_pstruct)
 469     if (pnonscalar)
 470       f = @ (p, varargin) \
 471           f (cell2struct \
 472              (cellfun (@ reshape, mat2cell (p, ppartidx), \
 473                        pdims, "UniformOutput", false), \
 474               pord, 1), varargin{:});
 475     else
 476       f = @ (p, varargin) \
 477           f (cell2struct (num2cell (p), pord, 1), varargin{:});
 478     endif
 479   endif
 480   f_pin = f (pin);
 481   if (isfield (hook, "observations"))
 482     if (any (size (f_pin) != size (obs = hook.observations)))
 483       error ("dimensions of observations and values of model function must match");
 484     endif
 485     f = @ (p) f (p) - obs;
 486     f_pin -= obs;
 487   endif
 488
 489   ## jacobian of model function
 490   if (isempty (dfdp))
 491     if (do_cstep)
 492       dfdp = @ (p, hook) jacobs (p, f, hook);
 493     else
 494       dfdp = @ (p, hook) __dfdp__ (p, f, hook);
 495     endif
 496   endif
 497   if (dfdp_pstruct)
 498     if (pnonscalar)
 499       dfdp = @ (p, hook) \
 500           cat (2, \
 501                fields2cell \
 502                (dfdp (cell2struct \
 503                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 504                                 pdims, "UniformOutput", false), \
 505                        pord, 1), hook), \
 506                 pord){:});
 507     else
 508       dfdp = @ (p, hook) \
 509           cat (2, \
 510                fields2cell \
 511                (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
 512                 pord){:});
 513     endif
 514   endif
 515
 516   ## function for general inequality constraints
 517   if (f_inequc_pstruct)
 518     if (pnonscalar)
 519       f_genicstr = @ (p, varargin) \
 520           f_genicstr (cell2struct \
 521                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 522                                 pdims, "UniformOutput", false), \
 523                        pord, 1), varargin{:});
 524     else
 525       f_genicstr = @ (p, varargin) \
 526           f_genicstr \
 527           (cell2struct (num2cell (p), pord, 1), varargin{:});
 528     endif
 529   endif
 530
 531   ## note this stage
 532   possibly_pstruct_f_genicstr = f_genicstr;
 533
 534   ## jacobian of general inequality constraints
 535   if (df_inequc_pstruct)
 536     if (pnonscalar)
 537       df_gencstr = @ (p, func, idx, hook) \
 538           cat (2, \
 539                fields2cell \
 540                (df_gencstr \
 541                 (cell2struct \
 542                  (cellfun (@ reshape, mat2cell (p, ppartidx), \
 543                            pdims, "UniformOutput", false), pord, 1), \
 544                  func, idx, hook), \
 545                 pord){:});
 546     else
 547       df_gencstr = @ (p, func, idx, hook) \
 548           cat (2, \
 549                fields2cell \
 550                (df_gencstr (cell2struct (num2cell (p), pord, 1), \
 551                             func, idx, hook), \
 552                 pord){:});
 553     endif
 554   endif
 555
 556   ## function for general equality constraints
 557   if (f_equc_pstruct)
 558     if (pnonscalar)
 559       f_genecstr = @ (p, varargin) \
 560           f_genecstr (cell2struct \
 561                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 562                                 pdims, "UniformOutput", false), \
 563                        pord, 1), varargin{:});
 564     else
 565       f_genecstr = @ (p, varargin) \
 566           f_genecstr \
 567           (cell2struct (num2cell (p), pord, 1), varargin{:});
 568     endif
 569   endif
 570
 571   ## note this stage
 572   possibly_pstruct_f_genecstr = f_genecstr;
 573
 574   ## jacobian of general equality constraints
 575   if (df_equc_pstruct)
 576     if (pnonscalar)
 577       df_genecstr = @ (p, func, idx, hook) \
 578           cat (2, \
 579                fields2cell \
 580                (df_genecstr \
 581                 (cell2struct \
 582                  (cellfun (@ reshape, mat2cell (p, ppartidx), \
 583                            pdims, "UniformOutput", false), pord, 1), \
 584                  func, idx, hook), \
 585                 pord){:});
 586     else
 587       df_genecstr = @ (p, func, idx, hook) \
 588           cat (2, \
 589                fields2cell \
 590                (df_genecstr (cell2struct (num2cell (p), pord, 1), \
 591                              func, idx, hook), \
 592                 pord){:});
 593     endif
 594   endif
 595
 596   ## linear inequality constraints
 597   if (mc_struct)
 598     idx = isfield (mc, pord);
 599     if (rows (fieldnames (mc)) > sum (idx))
 600       error ("unknown fields in structure of linear inequality constraints");
 601     endif
 602     smc = mc;
 603     mc = zeros (np, rows (vc));
 604     mc(idx(prepidx), :) = cat (1, fields2cell (smc, pord(idx)){:});
 605   endif
 606
 607   ## linear equality constraints
 608   if (emc_struct)
 609     idx = isfield (emc, pord);
 610     if (rows (fieldnames (emc)) > sum (idx))
 611       error ("unknown fields in structure of linear equality constraints");
 612     endif
 613     semc = emc;
 614     emc = zeros (np, rows (evc));
 615     emc(idx(prepidx), :) = cat (1, fields2cell (semc, pord(idx)){:});
 616   endif
 617
 618   ## parameter-related configuration for jacobian functions
 619   if (dfdp_pstruct || df_inequc_pstruct || df_equc_pstruct)
 620     if(pnonscalar)
 621       s_diffp = cell2struct \
 622           (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
 623                     pdims, "UniformOutput", false), pord, 1);
 624       s_diff_onesided = cell2struct \
 625           (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
 626                     pdims, "UniformOutput", false), pord, 1);
 627       s_orig_lbound = cell2struct \
 628           (cellfun (@ reshape, mat2cell (lbound, ppartidx), \
 629                     pdims, "UniformOutput", false), pord, 1);
 630       s_orig_ubound = cell2struct \
 631           (cellfun (@ reshape, mat2cell (ubound, ppartidx), \
 632                     pdims, "UniformOutput", false), pord, 1);
 633       s_plabels = cell2struct \
 634           (num2cell \
 635            (cat (2, cellfun \
 636                  (@ (x) cellfun \
 637                   (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
 638                    pdims, "UniformOutput", false), \
 639                   num2cell (plabels, 1), "UniformOutput", false){:}), \
 640             2), \
 641            pord, 1);
 642       s_orig_fixed = cell2struct \
 643           (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
 644                     pdims, "UniformOutput", false), pord, 1);
 645     else
 646       s_diffp = cell2struct (num2cell (diffp), pord, 1);
 647       s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
 648       s_orig_lbound = cell2struct (num2cell (lbound), pord, 1);
 649       s_orig_ubound = cell2struct (num2cell (ubound), pord, 1);
 650       s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
 651       s_orig_fixed = cell2struct (num2cell (fixed), pord, 1);
 652     endif
 653   endif
 654
 655   #### some further values and checks
 656
 657   if (any (fixed & (pin < lbound | pin > ubound)))
 658     warning ("some fixed parameters outside bounds");
 659   endif
 660
 661   ## dimensions of linear constraints
 662   if (isempty (mc))
 663     mc = zeros (np, 0);
 664     vc = zeros (0, 1);
 665   endif
 666   if (isempty (emc))
 667     emc = zeros (np, 0);
 668     evc = zeros (0, 1);
 669   endif
 670   [rm, cm] = size (mc);
 671   [rv, cv] = size (vc);
 672   if (rm != np || cm != rv || cv != 1)
 673     error ("linear inequality constraints: wrong dimensions");
 674   endif
 675   [erm, ecm] = size (emc);
 676   [erv, ecv] = size (evc);
 677   if (erm != np || ecm != erv || ecv != 1)
 678     error ("linear equality constraints: wrong dimensions");
 679   endif
 680
 681   ## check weights dimensions
 682   weights = optimget (settings, "weights", ones (size (f_pin)));
 683   if (any (size (weights) != size (f_pin)))
 684     error ("dimension of weights and residuals must match");
 685   endif
 686
 687   ## note initial values of linear constraits
 688   pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
 689   pin_cstr.equ.lin = emc.' * pin + evc;
 690
 691   ## note number and initial values of general constraints
 692   if (isempty (f_genicstr))
 693     pin_cstr.inequ.gen = [];
 694     n_genicstr = 0;
 695   else
 696     n_genicstr = length (pin_cstr.inequ.gen = f_genicstr (pin));
 697   endif
 698   if (isempty (f_genecstr))
 699     pin_cstr.equ.gen = [];
 700     n_genecstr = 0;
 701   else
 702     n_genecstr = length (pin_cstr.equ.gen = f_genecstr (pin));
 703   endif
 704
 705   #### collect remaining settings
 706   hook.TolFun = optimget (settings, "TolFun", stol_default);
 707   hook.MaxIter = optimget (settings, "MaxIter");
 708   if (ischar (hook.cpiv = optimget (settings, "cpiv", @ cpiv_bard)))
 709     hook.cpiv = str2func (hook.cpiv);
 710   endif
 711   hook.Display = optimget (settings, "Display", "off");
 712   hook.plot_cmd = optimget (settings, "plot_cmd", @ (f) 0);
 713   hook.testing = optimget (settings, "debug", false);
 714   hook.new_s = optimget (settings, "lm_svd_feasible_alt_s", false);
 715   backend = optimget (settings, "Algorithm", "lm_svd_feasible");
 716   backend = map_matlab_algorithm_names (backend);
 717   backend = map_backend (backend);
 718
 719   #### handle fixing of parameters
 720   orig_lbound = lbound;
 721   orig_ubound = ubound;
 722   orig_fixed = fixed;
 723   if (all (fixed))
 724     error ("no free parameters");
 725   endif
 726
 727   nonfixed = ! fixed;
 728   if (any (fixed))
 729     ## backend (returned values and initial parameters)
 730     backend = @ (f, pin, hook) \
 731         backend_wrapper (backend, fixed, f, pin, hook);
 732
 733     ## model function
 734     f = @ (p, varargin) f (assign (pin, nonfixed, p), varargin{:});
 735
 736     ## jacobian of model function
 737     dfdp = @ (p, hook) \
 738         dfdp (assign (pin, nonfixed, p), hook)(:, nonfixed);
 739
 740     ## function for general inequality constraints
 741     f_genicstr = @ (p, varargin) \
 742         f_genicstr (assign (pin, nonfixed, p), varargin{:});
 743
 744     ## jacobian of general inequality constraints
 745     df_gencstr = @ (p, func, idx, hook) \
 746         df_gencstr (assign (pin, nonfixed, p), func, idx, hook) \
 747         (:, nonfixed);
 748
 749     ## function for general equality constraints
 750     f_genecstr = @ (p, varargin) \
 751         f_genecstr (assign (pin, nonfixed, p), varargin{:});
 752
 753     ## jacobian of general equality constraints
 754     df_genecstr = @ (p, func, idx, hook) \
 755         df_genecstr (assign (pin, nonfixed, p), func, idx, hook) \
 756         (:, nonfixed);
 757
 758     ## linear inequality constraints
 759     vc += mc(fixed, :).' * (tp = pin(fixed));
 760     mc = mc(nonfixed, :);
 761
 762     ## linear equality constraints
 763     evc += emc(fixed, :).' * tp;
 764     emc = emc(nonfixed, :);
 765
 766     ## _last_ of all, vectors of parameter-related configuration,
 767     ## including "fixed" itself
 768     lbound = lbound(nonfixed, :);
 769     ubound = ubound(nonfixed, :);
 770     max_fract_change = max_fract_change(nonfixed);
 771     fract_prec = fract_prec(nonfixed);
 772     fixed = fixed(nonfixed);
 773   endif
 774
 775   #### supplement constants to jacobian functions
 776
 777   ## jacobian of model function
 778   if (dfdp_pstruct)
 779     dfdp = @ (p, hook) \
 780         dfdp (p, cell2fields \
 781               ({s_diffp, s_diff_onesided, s_orig_lbound, \
 782                 s_orig_ubound, s_plabels, \
 783                 cell2fields(num2cell(hook.fixed), pord(nonfixed), \
 784                             1, s_orig_fixed), cstep}, \
 785                {"diffp", "diff_onesided", "lbound", "ubound", \
 786                 "plabels", "fixed", "h"}, \
 787                2, hook));
 788   else
 789     dfdp = @ (p, hook) \
 790         dfdp (p, cell2fields \
 791               ({diffp, diff_onesided, orig_lbound, orig_ubound, \
 792                 plabels, assign(orig_fixed, nonfixed, hook.fixed), \
 793                 cstep}, \
 794                {"diffp", "diff_onesided", "lbound", "ubound", \
 795                 "plabels", "fixed", "h"}, \
 796                2, hook));
 797   endif
 798
 799   ## jacobian of general inequality constraints
 800   if (df_inequc_pstruct)
 801     df_gencstr = @ (p, func, idx, hook) \
 802         df_gencstr (p, func, idx, cell2fields \
 803                     ({s_diffp, s_diff_onesided, s_orig_lbound, \
 804                       s_orig_ubound, s_plabels, \
 805                       cell2fields(num2cell(hook.fixed), pord(nonfixed), \
 806                                   1, s_orig_fixed), cstep}, \
 807                      {"diffp", "diff_onesided", "lbound", "ubound", \
 808                       "plabels", "fixed", "h"}, \
 809                      2, hook));
 810   else
 811     df_gencstr = @ (p, func, idx, hook) \
 812         df_gencstr (p, func, idx, cell2fields \
 813                     ({diffp, diff_onesided, orig_lbound, \
 814                       orig_ubound, plabels, \
 815                       assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
 816                      {"diffp", "diff_onesided", "lbound", "ubound", \
 817                       "plabels", "fixed", "h"}, \
 818                      2, hook));
 819   endif
 820
 821   ## jacobian of general equality constraints
 822   if (df_equc_pstruct)
 823     df_genecstr = @ (p, func, idx, hook) \
 824         df_genecstr (p, func, idx, cell2fields \
 825                      ({s_diffp, s_diff_onesided, s_orig_lbound, \
 826                        s_orig_ubound, s_plabels, \
 827                        cell2fields(num2cell(hook.fixed), pord(nonfixed), \
 828                                    1, s_orig_fixed), cstep}, \
 829                       {"diffp", "diff_onesided", "lbound", "ubound", \
 830                        "plabels", "fixed", "h"}, \
 831                       2, hook));
 832   else
 833     df_genecstr = @ (p, func, idx, hook) \
 834         df_genecstr (p, func, idx, cell2fields \
 835                      ({diffp, diff_onesided, orig_lbound, \
 836                        orig_ubound, plabels, \
 837                        assign(orig_fixed, nonfixed, hook.fixed), cstep}, \
 838                       {"diffp", "diff_onesided", "lbound", "ubound", \
 839                        "plabels", "fixed", "h"}, \
 840                       2, hook));
 841   endif
 842
 843   #### interfaces to constraints
 844
 845   ## include bounds into linear inequality constraints
 846   tp = eye (sum (nonfixed));
 847   lidx = lbound != - Inf;
 848   uidx = ubound != Inf;
 849   mc = cat (2, tp(:, lidx), - tp(:, uidx), mc);
 850   vc = cat (1, - lbound(lidx, 1), ubound(uidx, 1), vc);
 851
 852   ## concatenate linear inequality and equality constraints
 853   mc = cat (2, mc, emc);
 854   vc = cat (1, vc, evc);
 855   n_lincstr = rows (vc);
 856
 857   ## concatenate general inequality and equality constraints
 858   if (n_genecstr > 0)
 859     if (n_genicstr > 0)
 860       nidxi = 1 : n_genicstr;
 861       nidxe = n_genicstr + 1 : n_genicstr + n_genecstr;
 862       f_gencstr = @ (p, idx, varargin) \
 863           cat (1, \
 864                f_genicstr (p, idx(nidxi), varargin{:}), \
 865                f_genecstr (p, idx(nidxe), varargin{:}));
 866       df_gencstr = @ (p, idx, hook) \
 867           cat (1, \
 868                df_gencstr (p, @ (p, varargin) \
 869                            possibly_pstruct_f_genicstr \
 870                            (p, idx(nidxi), varargin{:}), \
 871                            idx(nidxi), \
 872                            setfield (hook, "f", \
 873                                      hook.f(nidxi(idx(nidxi))))), \
 874                df_genecstr (p, @ (p, varargin) \
 875                             possibly_pstruct_f_genecstr \
 876                             (p, idx(nidxe), varargin{:}), \
 877                             idx(nidxe), \
 878                             setfield (hook, "f", \
 879                                       hook.f(nidxe(idx(nidxe))))));
 880     else
 881       f_gencstr = f_genecstr;
 882       df_gencstr = @ (p, idx, hook) \
 883           df_genecstr (p, \
 884                        @ (p, varargin) \
 885                        possibly_pstruct_f_genecstr \
 886                        (p, idx, varargin{:}), \
 887                        idx, \
 888                        setfield (hook, "f", hook.f(idx)));
 889     endif
 890   else
 891     f_gencstr = f_genicstr;
 892     df_gencstr = @ (p, idx, hook) \
 893         df_gencstr (p, \
 894                     @ (p, varargin) \
 895                     possibly_pstruct_f_genicstr (p, idx, varargin{:}), \
 896                     idx, \
 897                     setfield (hook, "f", hook.f(idx)));
 898   endif
 899   n_gencstr = n_genicstr + n_genecstr;
 900
 901   ## concatenate linear and general constraints, defining the final
 902   ## function interfaces
 903   if (n_gencstr > 0)
 904     nidxl = 1:n_lincstr;
 905     nidxh = n_lincstr + 1 : n_lincstr + n_gencstr;
 906     f_cstr = @ (p, idx, varargin) \
 907         cat (1, \
 908              mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), \
 909              f_gencstr (p, idx(nidxh), varargin{:}));
 910     df_cstr = @ (p, idx, hook) \
 911         cat (1, \
 912              mc(:, idx(nidxl)).', \
 913              df_gencstr (p, idx(nidxh), \
 914                          setfield (hook, "f", \
 915                                    hook.f(nidxh))));
 916   else
 917     f_cstr = @ (p, idx, varargin) mc(:, idx).' * p + vc(idx, 1);
 918     df_cstr = @ (p, idx, hook) mc(:, idx).';
 919   endif
 920
 921   ## define eq_idx (logical index of equality constraints within all
 922   ## concatenated constraints
 923   eq_idx = false (n_lincstr + n_gencstr, 1);
 924   eq_idx(n_lincstr + 1 - rows (evc) : n_lincstr) = true;
 925   n_cstr = n_lincstr + n_gencstr;
 926   eq_idx(n_cstr + 1 - n_genecstr : n_cstr) = true;
 927
 928   #### prepare interface hook
 929
 930   ## passed constraints
 931   hook.mc = mc;
 932   hook.vc = vc;
 933   hook.f_cstr = f_cstr;
 934   hook.df_cstr = df_cstr;
 935   hook.n_gencstr = n_gencstr;
 936   hook.eq_idx = eq_idx;
 937   hook.lbound = lbound;
 938   hook.ubound = ubound;
 939
 940   ## passed values of constraints for initial parameters
 941   hook.pin_cstr = pin_cstr;
 942
 943   ## passed function for derivative of model function
 944   hook.dfdp = dfdp;
 945
 946   ## passed function for complementary pivoting
 947   ## hook.cpiv = cpiv; # set before
 948
 949   ## passed value of residual function for initial parameters
 950   hook.f_pin = f_pin;
 951
 952   ## passed options
 953   hook.max_fract_change = max_fract_change;
 954   hook.fract_prec = fract_prec;
 955   ## hook.TolFun = ; # set before
 956   ## hook.MaxIter = ; # set before
 957   hook.weights = weights;
 958   hook.fixed = fixed;
 959
 960   #### call backend
 961
 962   [p, resid, cvg, outp] = backend (f, pin, hook);
 963
 964   if (pin_struct)
 965     if (pnonscalar)
 966       p = cell2struct \
 967           (cellfun (@ reshape, mat2cell (p, ppartidx), \
 968                     pdims, "UniformOutput", false), \
 969            pord, 1);
 970     else
 971       p = cell2struct (num2cell (p), pord, 1);
 972     endif
 973   endif
 974
 975 endfunction
 976
 977 function backend = map_matlab_algorithm_names (backend)
 978
 979   switch (backend)
 980     case "levenberg-marquardt"
 981       backend = "lm_svd_feasible";
 982       warning ("algorithm 'levenberg-marquardt' mapped to 'lm_svd_feasible'");
 983   endswitch
 984
 985 endfunction
 986
 987 function backend = map_backend (backend)
 988
 989   switch (backend)
 990     case "lm_svd_feasible"
 991       backend = "__lm_svd__";
 992     otherwise
 993       error ("no backend implemented for algorithm '%s'", backend);
 994   endswitch
 995
 996   backend = str2func (backend);
 997
 998 endfunction
 999
1000 function [p, resid, cvg, outp] = backend_wrapper (backend, fixed, f, p, hook)
1001
1002   [tp, resid, cvg, outp] = backend (f, p(! fixed), hook);
1003
1004   p(! fixed) = tp;
1005
1006 endfunction
1007
1008 function lval = assign (lval, lidx, rval)
1009
1010   lval(lidx) = rval;
1011
1012 endfunction
1013
1014 function ret = __optimget__ (s, name, default)
1015
1016   if (isfield (s, name))
1017     ret = s.(name);
1018   elseif (nargin > 2)
1019     ret = default;
1020   else
1021     ret = [];
1022   endif
1023
1024 endfunction
1025
1026 function ret = apply_idx_if_given  (ret, varargin)
1027
1028   if (nargin > 1)
1029     ret = ret(varargin{1});
1030   endif
1031
1032 endfunction