octave_packages/optim-1.2.0/private/__residmin_stat__.m

   1 ## Copyright (C) 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 residmin_stat --- see there --- and
  17 ## others. Calling __residmin_stat__ indirectly hides the argument
  18 ## "hook", usable by wrappers, from users. Currently, hook can contain
  19 ## the field "observations". Since much uf the interface code is taken
  20 ## from __nonlin_residmin__, it may be that not everything is ideal for
  21 ## the present case; but I think it's allright to leave it so.
  22 ##
  23 ## Some general considerations while making this function:
  24 ##
  25 ## Different Functions for optimization statistics should be made for
  26 ## mere objective function optimization (to be made yet) and
  27 ## residual-derived minimization (this function), since there are
  28 ## different computing aspects. Don't put the contained functionality
  29 ## (statistics) into the respective optimization functions (or
  30 ## backends), since different optimization algorithms can share a way to
  31 ## compute statistics (e.g. even stochastic optimizers can mimize
  32 ## (weighted) squares of residuals). Also, don't use the same frontend
  33 ## for optimization and statistics, since the differences in the
  34 ## interface for both uses may be confusing otherwise, also the optimset
  35 ## options only partially overlap.
  36
  37 ## disabled PKG_ADD: __all_opts__ ("__residmin_stat__");
  38
  39 function ret = __residmin_stat__ (f, pfin, settings, hook)
  40
  41   if (compare_versions (version (), "3.3.55", "<"))
  42     ## optimset mechanism was fixed for option names with underscores
  43     ## sometime in 3.3.54+, if I remember right
  44     optimget = @ __optimget__;
  45   endif
  46
  47   ## scalar defaults
  48   diffp_default = .001;
  49   cstep_default = 1e-20;
  50
  51   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
  52     ret = optimset ("param_config", [], \
  53                     "param_order", [], \
  54                     "param_dims", [], \
  55                     "f_pstruct", false, \
  56                     "dfdp_pstruct", false, \
  57                     "dfdp", [], \
  58                     "diffp", [], \
  59                     "diff_onesided", [], \
  60                     "complex_step_derivative", false, \
  61                     "cstep", cstep_default, \
  62                     "fixed", [], \
  63                     "weights", [], \
  64                     "residuals", [], \
  65                     "covd", [], \
  66                     "objf", [], \ # no default, e.g. "wls"
  67                     "ret_dfdp", false, \
  68                     "ret_covd", false, \
  69                     "ret_covp", false, \
  70                     "ret_corp", false);
  71     return;
  72   endif
  73
  74   assign = @ assign; # Is this faster in repeated calls?
  75
  76   if (nargin != 4)
  77     error ("incorrect number of arguments");
  78   endif
  79
  80   if (ischar (f))
  81     f = str2func (f);
  82   endif
  83
  84   if (! (p_struct = isstruct (pfin)))
  85     if (! isempty (pfin) && (! isvector (pfin) || columns (pfin) > 1))
  86       error ("parameters must be either a structure or a column vector");
  87     endif
  88   endif
  89
  90   #### processing of settings and consistency checks
  91
  92   pconf = optimget (settings, "param_config");
  93   pord = optimget (settings, "param_order");
  94   pdims = optimget (settings, "param_dims");
  95   f_pstruct = optimget (settings, "f_pstruct", false);
  96   dfdp_pstruct = optimget (settings, "dfdp_pstruct", f_pstruct);
  97   dfdp = optimget (settings, "dfdp");
  98   if (ischar (dfdp)) dfdp = str2func (dfdp); endif
  99   if (isstruct (dfdp)) dfdp_pstruct = true; endif
 100   diffp = optimget (settings, "diffp");
 101   diff_onesided = optimget (settings, "diff_onesided");
 102   fixed = optimget (settings, "fixed");
 103   residuals = optimget (settings, "residuals");
 104   do_cstep = optimget (settings, "complex_step_derivative", false);
 105   cstep = optimget (settings, "cstep", cstep_default);
 106   if (do_cstep && ! isempty (dfdp))
 107     error ("both 'complex_step_derivative' and 'dfdp' are set");
 108   endif
 109
 110   any_vector_conf = ! (isempty (diffp) && isempty (diff_onesided) && \
 111                        isempty (fixed));
 112
 113   ## correct "_pstruct" settings if functions are not supplied
 114   if (isempty (dfdp)) dfdp_pstruct = false; endif
 115   if (isempty (f)) f_pstruct = false; endif
 116
 117   ## some settings require a parameter order
 118   if (p_struct || ! isempty (pconf) || f_pstruct || dfdp_pstruct)
 119     if (isempty (pord))
 120       if (p_struct)
 121         if (any_vector_conf || \
 122             ! ((f_pstruct || isempty (f)) && \
 123                (dfdp_pstruct || isempty (dfdp))))
 124           error ("no parameter order specified and constructing a parameter order from the structure of parameters can not be done since not all configuration or given functions are structure based");
 125         else
 126           pord = fieldnames (pfin);
 127         endif
 128       else
 129         error ("given settings require specification of parameter order or parameters in the form of a structure");
 130       endif
 131     endif
 132     pord = pord(:);
 133     if (p_struct && ! all (isfield (pfin, pord)))
 134       error ("some parameters lacking");
 135     endif
 136     if ((nnames = rows (unique (pord))) < rows (pord))
 137       error ("duplicate parameter names in 'param_order'");
 138     endif
 139     if (isempty (pdims))
 140       if (p_struct)
 141         pdims = cellfun \
 142             (@ size, fields2cell (pfin, pord), "UniformOutput", false);
 143       else
 144         pdims = num2cell (ones (nnames, 2), 2);
 145       endif
 146     else
 147       pdims = pdims(:);
 148       if (p_struct && \
 149           ! all (cellfun (@ (x, y) prod (size (x)) == prod (y), \
 150                           struct2cell (pfin), pdims)))
 151         error ("given param_dims and dimensions of parameters do not match");
 152       endif
 153     endif
 154     if (nnames != rows (pdims))
 155       error ("lengths of 'param_order' and 'param_dims' not equal");
 156     endif
 157     pnel = cellfun (@ prod, pdims);
 158     ppartidx = pnel;
 159     if (any (pnel > 1))
 160       pnonscalar = true;
 161       cpnel = num2cell (pnel);
 162       prepidx = cat (1, cellfun \
 163                      (@ (x, n) x(ones (1, n), 1), \
 164                       num2cell ((1:nnames).'), cpnel, \
 165                       "UniformOutput", false){:});
 166       epord = pord(prepidx, 1);
 167       psubidx = cat (1, cellfun \
 168                      (@ (n) (1:n).', cpnel, \
 169                       "UniformOutput", false){:});
 170     else
 171       pnonscalar = false; # some less expensive interfaces later
 172       prepidx = (1:nnames).';
 173       epord = pord;
 174       psubidx = ones (nnames, 1);
 175     endif
 176   else
 177     pord = []; # spares checks for given but not needed
 178   endif
 179
 180   if (p_struct)
 181     np = sum (pnel);
 182   else
 183     np = length (pfin);
 184     if (! isempty (pord) && np != sum (pnel))
 185       error ("number of initial parameters not correct");
 186     endif
 187   endif
 188   if (ismatrix (dfdp) && ! ischar (dfdp) && ! isempty (dfdp) && \
 189       np == 0)
 190     np = columns (dfdp);
 191   endif
 192
 193   plabels = num2cell (num2cell ((1:np).'));
 194   if (! isempty (pord))
 195     plabels = cat (2, plabels, num2cell (epord), \
 196                    num2cell (num2cell (psubidx)));
 197   endif
 198
 199   ## some useful vectors
 200   zerosvec = zeros (np, 1);
 201   falsevec = false (np, 1);
 202   sizevec = [np, 1];
 203
 204   ## collect parameter-related configuration
 205   if (! isempty (pconf))
 206     ## use supplied configuration structure
 207
 208     ## parameter-related configuration is either allowed by a structure
 209     ## or by vectors
 210     if (any_vector_conf)
 211       error ("if param_config is given, its potential items must not \
 212           be configured in another way");
 213     endif
 214
 215     ## supplement parameter names lacking in param_config
 216     nidx = ! isfield (pconf, pord);
 217     pconf = cell2fields ({struct()}(ones (1, sum (nidx))), \
 218                          pord(nidx), 2, pconf);
 219
 220     pconf = structcat (1, fields2cell (pconf, pord){:});
 221
 222     ## in the following, use reshape with explicit dimensions (instead
 223     ## of x(:)) so that errors are thrown if a configuration item has
 224     ## incorrect number of elements
 225
 226     diffp = zerosvec;
 227     diffp(:) = diffp_default;
 228     if (isfield (pconf, "diffp"))
 229       idx = ! fieldempty (pconf, "diffp");
 230       if (pnonscalar)
 231         diffp(idx(prepidx)) = \
 232             cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 233                              {pconf(idx).diffp}.', cpnel(idx), \
 234                              "UniformOutput", false){:});
 235       else
 236         diffp(idx) = [pconf.diffp];
 237       endif
 238     endif
 239
 240     diff_onesided = fixed = falsevec;
 241
 242     if (isfield (pconf, "diff_onesided"))
 243       idx = ! fieldempty (pconf, "diff_onesided");
 244       if (pnonscalar)
 245         diff_onesided(idx(prepidx)) = \
 246             logical \
 247             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 248                               {pconf(idx).diff_onesided}.', cpnel(idx), \
 249                              "UniformOutput", false){:}));
 250       else
 251         diff_onesided(idx) = logical ([pconf.diff_onesided]);
 252       endif
 253     endif
 254
 255     if (isfield (pconf, "fixed"))
 256       idx = ! fieldempty (pconf, "fixed");
 257       if (pnonscalar)
 258         fixed(idx(prepidx)) = \
 259             logical \
 260             (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 261                               {pconf(idx).fixed}.', cpnel(idx), \
 262                              "UniformOutput", false){:}));
 263       else
 264         fixed(idx) = logical ([pconf.fixed]);
 265       endif
 266     endif
 267   else
 268     ## use supplied configuration vectors
 269
 270     if (isempty (diffp))
 271       diffp = zerosvec;
 272       diffp(:) = diffp_default;
 273     else
 274       if (any (size (diffp) != sizevec))
 275         error ("diffp: wrong dimensions");
 276       endif
 277       diffp(isna (diffp)) = diffp_default;
 278     endif
 279
 280     if (isempty (diff_onesided))
 281       diff_onesided = falsevec;
 282     else
 283       if (any (size (diff_onesided) != sizevec))
 284         error ("diff_onesided: wrong dimensions")
 285       endif
 286       diff_onesided(isna (diff_onesided)) = false;
 287       diff_onesided = logical (diff_onesided);
 288     endif
 289
 290     if (isempty (fixed))
 291       fixed = falsevec;
 292     else
 293       if (any (size (fixed) != sizevec))
 294         error ("fixed: wrong dimensions");
 295       endif
 296       fixed(isna (fixed)) = false;
 297       fixed = logical (fixed);
 298     endif
 299   endif
 300
 301   #### consider whether parameters and functions are based on parameter
 302   #### structures or parameter vectors; wrappers for call to default
 303   #### function for jacobians
 304
 305   ## parameters
 306   if (p_struct)
 307     if (pnonscalar)
 308       pfin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
 309                               fields2cell (pfin, pord), cpnel, \
 310                               "UniformOutput", false){:});
 311     else
 312       pfin = cat (1, fields2cell (pfin, pord){:});
 313     endif
 314   endif
 315
 316   ## model function
 317   if (f_pstruct)
 318     if (pnonscalar)
 319       f = @ (p, varargin) \
 320           f (cell2struct \
 321              (cellfun (@ reshape, mat2cell (p, ppartidx), \
 322                        pdims, "UniformOutput", false), \
 323               pord, 1), varargin{:});
 324     else
 325       f = @ (p, varargin) \
 326           f (cell2struct (num2cell (p), pord, 1), varargin{:});
 327     endif
 328   endif
 329   if (isempty (residuals))
 330     if (isempty (f))
 331       error ("neither model function nor residuals given");
 332     endif
 333     residuals = f (pfin);
 334   endif
 335   if (isfield (hook, "observations"))
 336     if (any (size (residuals) != size (obs = hook.observations)))
 337       error ("dimensions of observations and values of model function must match");
 338     endif
 339     f = @ (p) f (p) - obs;
 340     residuals -= obs;
 341   endif
 342
 343   ## jacobian of model function
 344   if (isempty (dfdp))
 345     if (! isempty (f))
 346       if (do_cstep)
 347         dfdp = @ (p, hook) jacobs (p, f, hook);
 348       else
 349         __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
 350         dfdp = @ (p, hook) __dfdp__ (p, f, hook);
 351       endif
 352     endif
 353   elseif (! isa (dfdp, "function_handle"))
 354     if (ismatrix (dfdp))
 355       if (any (size (dfdp) != [prod(size(residuals)), np]))
 356         error ("jacobian has wrong size");
 357       endif
 358     elseif (! dfdp_pstruct)
 359       error ("jacobian has wrong type");
 360     endif
 361     dfdp = @ (varargin) dfdp; # simply make a function returning it
 362   endif
 363   if (dfdp_pstruct)
 364     if (pnonscalar)
 365       dfdp = @ (p, hook) \
 366           cat (2, \
 367                fields2cell \
 368                (dfdp (cell2struct \
 369                       (cellfun (@ reshape, mat2cell (p, ppartidx), \
 370                                 pdims, "UniformOutput", false), \
 371                        pord, 1), hook), \
 372                 pord){:});
 373     else
 374       dfdp = @ (p, hook) \
 375           cat (2, \
 376                fields2cell \
 377                (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
 378                 pord){:});
 379     endif
 380   endif
 381
 382   ## parameter-related configuration for jacobian function
 383   if (dfdp_pstruct)
 384     if(pnonscalar)
 385       s_diffp = cell2struct \
 386           (cellfun (@ reshape, mat2cell (diffp, ppartidx), \
 387                     pdims, "UniformOutput", false), pord, 1);
 388       s_diff_onesided = cell2struct \
 389           (cellfun (@ reshape, mat2cell (diff_onesided, ppartidx), \
 390                     pdims, "UniformOutput", false), pord, 1);
 391       s_plabels = cell2struct \
 392           (num2cell \
 393            (cat (2, cellfun \
 394                  (@ (x) cellfun \
 395                   (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
 396                    pdims, "UniformOutput", false), \
 397                   num2cell (plabels, 1), "UniformOutput", false){:}), \
 398             2), \
 399            pord, 1);
 400       s_orig_fixed = cell2struct \
 401           (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
 402                     pdims, "UniformOutput", false), pord, 1);
 403     else
 404       s_diffp = cell2struct (num2cell (diffp), pord, 1);
 405       s_diff_onesided = cell2struct (num2cell (diff_onesided), pord, 1);
 406       s_plabels = cell2struct (num2cell (plabels, 2), pord, 1);
 407       s_fixed = cell2struct (num2cell (fixed), pord, 1);
 408     endif
 409   endif
 410
 411   #### further values and checks
 412
 413   ## check weights dimensions
 414   weights = optimget (settings, "weights", ones (size (residuals)));
 415   if (any (size (weights) != size (residuals)))
 416     error ("dimension of weights and residuals must match");
 417   endif
 418
 419
 420   #### collect remaining settings
 421   need_dfdp = false;
 422   covd = optimget (settings, "covd");
 423   need_objf_label = false;
 424   if ((ret_dfdp = optimget (settings, "ret_dfdp", false)))
 425     need_dfdp = true;
 426   endif
 427   if ((ret_covd = optimget (settings, "ret_covd", false)))
 428     need_objf_label = true;
 429     if (np == 0)
 430       error ("number of parameters must be known for 'covd', specify either parameters or a jacobian matrix");
 431     endif
 432   endif
 433   if ((ret_covp = optimget (settings, "ret_covp", false)))
 434     need_objf_label = true;
 435     need_dfdp = true;
 436   endif
 437   if ((ret_corp = optimget (settings, "ret_corp", false)))
 438     need_objf_label = true;
 439     need_dfdp = true;
 440   endif
 441   if (need_objf_label)
 442     if (isempty (objf = optimget (settings, "objf")))
 443       error ("label of objective function must be specified");
 444     else
 445       funs = map_objf (objf);
 446     endif
 447   else
 448     funs = struct ();
 449   endif
 450   if (isempty (dfdp) && need_dfdp)
 451     error ("jacobian required and default function for jacobian requires a model function");
 452   endif
 453
 454   ####
 455
 456   ## Everything which is computed is stored in a hook structure which is
 457   ## passed to and returned by every backend function. This hook is not
 458   ## identical to the returned structure, since some more results could
 459   ## be computed by the way.
 460
 461   #### handle fixing of parameters
 462
 463   orig_p = pfin;
 464
 465   if (all (fixed) && ! isempty (fixed))
 466     error ("no free parameters");
 467   endif
 468
 469   ## The policy should be that everything which is computed is left as
 470   ## it is up to the end --- since other computations might need it in
 471   ## this form --- and supplemented with values corresponding to fixed
 472   ## parameters (mostly NA, probably) not until then.
 473
 474   nonfixed = ! fixed;
 475
 476   np_after_fixing = sum (nonfixed);
 477
 478   if (any (fixed))
 479
 480     if (! isempty (pfin))
 481       pfin = pfin(nonfixed);
 482     endif
 483
 484     ## model function
 485     f = @ (p, varargin) f (assign (pfin, nonfixed, p), varargin{:});
 486
 487     ## jacobian of model function
 488     if (! isempty (dfdp))
 489       dfdp = @ (p, hook) \
 490           dfdp (assign (pfin, nonfixed, p), hook)(:, nonfixed);
 491     endif
 492
 493   endif
 494
 495   #### supplement constants to jacobian function
 496   if (dfdp_pstruct)
 497     dfdp = @ (p, hook) \
 498         dfdp (p, cell2fields \
 499               ({s_diffp, s_diff_onesided, s_plabels, s_fixed, cstep}, \
 500                {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
 501                2, hook));
 502   else
 503     if (! isempty (dfdp))
 504       dfdp = @ (p, hook) \
 505           dfdp (p, cell2fields \
 506                 ({diffp, diff_onesided, plabels, fixed, cstep}, \
 507                  {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
 508                  2, hook));
 509     endif
 510   endif
 511
 512   #### prepare interface hook
 513
 514   ## passed final parameters of an optimization
 515   hook.pfin = pfin;
 516
 517   ## passed function for derivative of model function
 518   hook.dfdp = dfdp;
 519
 520   ## passed function for complementary pivoting
 521   ## hook.cpiv = cpiv; # set before
 522
 523   ## passed value of residual function for initial parameters
 524   hook.residuals = residuals;
 525
 526   ## passed weights
 527   hook.weights = weights;
 528
 529   ## passed dimensions
 530   hook.np = np_after_fixing;
 531   hook.nm = prod (size (residuals));
 532
 533   ## passed statistics functions
 534   hook.funs = funs;
 535
 536   ## passed covariance matrix of data (if given by user)
 537   if (! isempty (covd))
 538     covd_dims = size (covd);
 539     if (length (covd_dims) != 2 || any (covd_dims != hook.nm))
 540       error ("wrong dimensions of covariance matrix of data");
 541     endif
 542     hook.covd = covd;
 543   endif
 544
 545   #### do the actual work
 546
 547   if (ret_dfdp)
 548     hook.jac = hook.dfdp (hook.pfin, hook);
 549   endif
 550
 551   if (ret_covd)
 552     hook = funs.covd (hook);
 553   endif
 554
 555   if (ret_covp || ret_corp)
 556     hook = funs.covp_corp (hook);
 557   endif
 558
 559   #### convert (consider fixing ...) and return results
 560
 561   ret = struct ();
 562
 563   if (ret_dfdp)
 564     ret.dfdp = zeros (hook.nm, np);
 565     ret.dfdp(:, nonfixed) = hook.jac;
 566   endif
 567
 568   if (ret_covd)
 569     ret.covd = hook.covd;
 570   endif
 571
 572   if (ret_covp)
 573     if (any (fixed))
 574       ret.covp = NA (np);
 575       ret.covp(nonfixed, nonfixed) = hook.covp;
 576     else
 577       ret.covp = hook.covp;
 578     endif
 579   endif
 580
 581   if (ret_corp)
 582     if (any (fixed))
 583       ret.corp = NA (np);
 584       ret.corp(nonfixed, nonfixed) = hook.corp;
 585     else
 586       ret.corp = hook.corp;
 587     endif
 588   endif
 589
 590 endfunction
 591
 592 function funs = map_objf (objf)
 593
 594   switch (objf)
 595     case "wls" # weighted least squares
 596       funs.covd = str2func ("__covd_wls__");
 597       funs.covp_corp = str2func ("__covp_corp_wls__");
 598     otherwise
 599       error ("no statistics implemented for objective function '%s'", \
 600              objf);
 601   endswitch
 602
 603 endfunction
 604
 605 function lval = assign (lval, lidx, rval)
 606
 607   lval(lidx) = rval;
 608
 609 endfunction
 610
 611 function ret = __optimget__ (s, name, default)
 612
 613   if (isfield (s, name))
 614     ret = s.(name);
 615   elseif (nargin > 2)
 616     ret = default;
 617   else
 618     ret = [];
 619   endif
 620
 621 endfunction