1 ## Copyright (C) 2011 Olaf Till <olaf.till@uni-jena.de>
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.
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.
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/>.
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.
23 ## Some general considerations while making this function:
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.
37 ## disabled PKG_ADD: __all_opts__ ("__residmin_stat__");
39 function ret = __residmin_stat__ (f, pfin, settings, hook)
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__;
49 cstep_default = 1e-20;
51 if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
52 ret = optimset ("param_config", [], \
56 "dfdp_pstruct", false, \
59 "diff_onesided", [], \
60 "complex_step_derivative", false, \
61 "cstep", cstep_default, \
66 "objf", [], \ # no default, e.g. "wls"
74 assign = @ assign; # Is this faster in repeated calls?
77 error ("incorrect number of arguments");
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");
90 #### processing of settings and consistency checks
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");
110 any_vector_conf = ! (isempty (diffp) && isempty (diff_onesided) && \
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
117 ## some settings require a parameter order
118 if (p_struct || ! isempty (pconf) || f_pstruct || dfdp_pstruct)
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");
126 pord = fieldnames (pfin);
129 error ("given settings require specification of parameter order or parameters in the form of a structure");
133 if (p_struct && ! all (isfield (pfin, pord)))
134 error ("some parameters lacking");
136 if ((nnames = rows (unique (pord))) < rows (pord))
137 error ("duplicate parameter names in 'param_order'");
142 (@ size, fields2cell (pfin, pord), "UniformOutput", false);
144 pdims = num2cell (ones (nnames, 2), 2);
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");
154 if (nnames != rows (pdims))
155 error ("lengths of 'param_order' and 'param_dims' not equal");
157 pnel = cellfun (@ prod, pdims);
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){:});
171 pnonscalar = false; # some less expensive interfaces later
172 prepidx = (1:nnames).';
174 psubidx = ones (nnames, 1);
177 pord = []; # spares checks for given but not needed
184 if (! isempty (pord) && np != sum (pnel))
185 error ("number of initial parameters not correct");
188 if (ismatrix (dfdp) && ! ischar (dfdp) && ! isempty (dfdp) && \
193 plabels = num2cell (num2cell ((1:np).'));
194 if (! isempty (pord))
195 plabels = cat (2, plabels, num2cell (epord), \
196 num2cell (num2cell (psubidx)));
199 ## some useful vectors
200 zerosvec = zeros (np, 1);
201 falsevec = false (np, 1);
204 ## collect parameter-related configuration
205 if (! isempty (pconf))
206 ## use supplied configuration structure
208 ## parameter-related configuration is either allowed by a structure
211 error ("if param_config is given, its potential items must not \
212 be configured in another way");
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);
220 pconf = structcat (1, fields2cell (pconf, pord){:});
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
227 diffp(:) = diffp_default;
228 if (isfield (pconf, "diffp"))
229 idx = ! fieldempty (pconf, "diffp");
231 diffp(idx(prepidx)) = \
232 cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
233 {pconf(idx).diffp}.', cpnel(idx), \
234 "UniformOutput", false){:});
236 diffp(idx) = [pconf.diffp];
240 diff_onesided = fixed = falsevec;
242 if (isfield (pconf, "diff_onesided"))
243 idx = ! fieldempty (pconf, "diff_onesided");
245 diff_onesided(idx(prepidx)) = \
247 (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
248 {pconf(idx).diff_onesided}.', cpnel(idx), \
249 "UniformOutput", false){:}));
251 diff_onesided(idx) = logical ([pconf.diff_onesided]);
255 if (isfield (pconf, "fixed"))
256 idx = ! fieldempty (pconf, "fixed");
258 fixed(idx(prepidx)) = \
260 (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
261 {pconf(idx).fixed}.', cpnel(idx), \
262 "UniformOutput", false){:}));
264 fixed(idx) = logical ([pconf.fixed]);
268 ## use supplied configuration vectors
272 diffp(:) = diffp_default;
274 if (any (size (diffp) != sizevec))
275 error ("diffp: wrong dimensions");
277 diffp(isna (diffp)) = diffp_default;
280 if (isempty (diff_onesided))
281 diff_onesided = falsevec;
283 if (any (size (diff_onesided) != sizevec))
284 error ("diff_onesided: wrong dimensions")
286 diff_onesided(isna (diff_onesided)) = false;
287 diff_onesided = logical (diff_onesided);
293 if (any (size (fixed) != sizevec))
294 error ("fixed: wrong dimensions");
296 fixed(isna (fixed)) = false;
297 fixed = logical (fixed);
301 #### consider whether parameters and functions are based on parameter
302 #### structures or parameter vectors; wrappers for call to default
303 #### function for jacobians
308 pfin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
309 fields2cell (pfin, pord), cpnel, \
310 "UniformOutput", false){:});
312 pfin = cat (1, fields2cell (pfin, pord){:});
319 f = @ (p, varargin) \
321 (cellfun (@ reshape, mat2cell (p, ppartidx), \
322 pdims, "UniformOutput", false), \
323 pord, 1), varargin{:});
325 f = @ (p, varargin) \
326 f (cell2struct (num2cell (p), pord, 1), varargin{:});
329 if (isempty (residuals))
331 error ("neither model function nor residuals given");
333 residuals = f (pfin);
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");
339 f = @ (p) f (p) - obs;
343 ## jacobian of model function
347 dfdp = @ (p, hook) jacobs (p, f, hook);
349 __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
350 dfdp = @ (p, hook) __dfdp__ (p, f, hook);
353 elseif (! isa (dfdp, "function_handle"))
355 if (any (size (dfdp) != [prod(size(residuals)), np]))
356 error ("jacobian has wrong size");
358 elseif (! dfdp_pstruct)
359 error ("jacobian has wrong type");
361 dfdp = @ (varargin) dfdp; # simply make a function returning it
369 (cellfun (@ reshape, mat2cell (p, ppartidx), \
370 pdims, "UniformOutput", false), \
377 (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
382 ## parameter-related configuration for jacobian function
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 \
395 (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
396 pdims, "UniformOutput", false), \
397 num2cell (plabels, 1), "UniformOutput", false){:}), \
400 s_orig_fixed = cell2struct \
401 (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
402 pdims, "UniformOutput", false), pord, 1);
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);
411 #### further values and checks
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");
420 #### collect remaining settings
422 covd = optimget (settings, "covd");
423 need_objf_label = false;
424 if ((ret_dfdp = optimget (settings, "ret_dfdp", false)))
427 if ((ret_covd = optimget (settings, "ret_covd", false)))
428 need_objf_label = true;
430 error ("number of parameters must be known for 'covd', specify either parameters or a jacobian matrix");
433 if ((ret_covp = optimget (settings, "ret_covp", false)))
434 need_objf_label = true;
437 if ((ret_corp = optimget (settings, "ret_corp", false)))
438 need_objf_label = true;
442 if (isempty (objf = optimget (settings, "objf")))
443 error ("label of objective function must be specified");
445 funs = map_objf (objf);
450 if (isempty (dfdp) && need_dfdp)
451 error ("jacobian required and default function for jacobian requires a model function");
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.
461 #### handle fixing of parameters
465 if (all (fixed) && ! isempty (fixed))
466 error ("no free parameters");
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.
476 np_after_fixing = sum (nonfixed);
480 if (! isempty (pfin))
481 pfin = pfin(nonfixed);
485 f = @ (p, varargin) f (assign (pfin, nonfixed, p), varargin{:});
487 ## jacobian of model function
488 if (! isempty (dfdp))
490 dfdp (assign (pfin, nonfixed, p), hook)(:, nonfixed);
495 #### supplement constants to jacobian function
498 dfdp (p, cell2fields \
499 ({s_diffp, s_diff_onesided, s_plabels, s_fixed, cstep}, \
500 {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
503 if (! isempty (dfdp))
505 dfdp (p, cell2fields \
506 ({diffp, diff_onesided, plabels, fixed, cstep}, \
507 {"diffp", "diff_onesided", "plabels", "fixed", "h"}, \
512 #### prepare interface hook
514 ## passed final parameters of an optimization
517 ## passed function for derivative of model function
520 ## passed function for complementary pivoting
521 ## hook.cpiv = cpiv; # set before
523 ## passed value of residual function for initial parameters
524 hook.residuals = residuals;
527 hook.weights = weights;
530 hook.np = np_after_fixing;
531 hook.nm = prod (size (residuals));
533 ## passed statistics functions
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");
545 #### do the actual work
548 hook.jac = hook.dfdp (hook.pfin, hook);
552 hook = funs.covd (hook);
555 if (ret_covp || ret_corp)
556 hook = funs.covp_corp (hook);
559 #### convert (consider fixing ...) and return results
564 ret.dfdp = zeros (hook.nm, np);
565 ret.dfdp(:, nonfixed) = hook.jac;
569 ret.covd = hook.covd;
575 ret.covp(nonfixed, nonfixed) = hook.covp;
577 ret.covp = hook.covp;
584 ret.corp(nonfixed, nonfixed) = hook.corp;
586 ret.corp = hook.corp;
592 function funs = map_objf (objf)
595 case "wls" # weighted least squares
596 funs.covd = str2func ("__covd_wls__");
597 funs.covp_corp = str2func ("__covp_corp_wls__");
599 error ("no statistics implemented for objective function '%s'", \
605 function lval = assign (lval, lidx, rval)
611 function ret = __optimget__ (s, name, default)
613 if (isfield (s, name))