1 ## Copyright (C) 2010, 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 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.
23 ## disabled PKG_ADD: __all_opts__ ("__nonlin_residmin__");
25 function [p, resid, cvg, outp] = \
26 __nonlin_residmin__ (f, pin, settings, hook)
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__;
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__;
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
45 cstep_default = 1e-20;
47 if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
48 p = optimset ("param_config", [], \
51 "f_inequc_pstruct", false, \
52 "f_equc_pstruct", false, \
54 "df_inequc_pstruct", false, \
55 "df_equc_pstruct", false, \
56 "dfdp_pstruct", false, \
60 "cpiv", @ cpiv_bard, \
61 "max_fract_change", [], \
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, \
72 "inequc_f_idx", false, \
73 "inequc_df_idx", false, \
74 "equc_f_idx", false, \
75 "equc_df_idx", false, \
77 "TolFun", stol_default, \
80 "Algorithm", "lm_svd_feasible", \
81 "plot_cmd", @ (f) 0, \
83 "lm_svd_feasible_alt_s", false);
87 assign = @ assign; # Is this faster in repeated calls?
90 error ("incorrect number of arguments");
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");
103 #### processing of settings and consistency checks
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", \
114 df_equc_pstruct = optimget (settings, "df_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");
131 optimget (settings, "complex_step_derivative_inequc", false);
132 do_cstep_equc = optimget (settings, "complex_step_derivative_equc", false);
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));
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);
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{:});
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{:});
165 if (! optimget (settings, "inequc_df_idx", false))
166 df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
169 df_inequc_pstruct = false;
171 if (user_df_genecstr)
172 if (! optimget (settings, "equc_df_idx", false))
173 df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
176 df_equc_pstruct = false;
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 || \
186 if (any_vector_conf || \
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");
197 pord = fieldnames (pin);
200 error ("given settings require specification of parameter order or initial parameters in the form of a structure");
204 if (pin_struct && ! all (isfield (pin, pord)))
205 error ("some initial parameters lacking");
207 if ((nnames = rows (unique (pord))) < rows (pord))
208 error ("duplicate parameter names in 'param_order'");
213 (@ size, fields2cell (pin, pord), "UniformOutput", false);
215 pdims = num2cell (ones (nnames, 2), 2);
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");
225 if (nnames != rows (pdims))
226 error ("lengths of 'param_order' and 'param_dims' not equal");
228 pnel = cellfun (@ prod, pdims);
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){:});
242 pnonscalar = false; # some less expensive interfaces later
243 prepidx = (1:nnames).';
245 psubidx = ones (nnames, 1);
248 pord = []; # spares checks for given but not needed
255 if (! isempty (pord) && np != sum (pnel))
256 error ("number of initial parameters not correct");
260 plabels = num2cell (num2cell ((1:np).'));
261 if (! isempty (pord))
262 plabels = cat (2, plabels, num2cell (epord), \
263 num2cell (num2cell (psubidx)));
266 ## some useful vectors
267 zerosvec = zeros (np, 1);
269 Infvec = Inf (np, 1);
270 falsevec = false (np, 1);
273 ## collect parameter-related configuration
274 if (! isempty (pconf))
275 ## use supplied configuration structure
277 ## parameter-related configuration is either allowed by a structure
280 error ("if param_config is given, its potential items must not \
281 be configured in another way");
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);
289 pconf = structcat (1, fields2cell (pconf, pord){:});
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
296 if (isfield (pconf, "lbound"))
297 idx = ! fieldempty (pconf, "lbound");
299 lbound (idx(prepidx), 1) = \
300 cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
301 {pconf(idx).lbound}.', \
302 cpnel(idx), "UniformOutput", false){:});
304 lbound(idx, 1) = cat (1, pconf.lbound);
309 if (isfield (pconf, "ubound"))
310 idx = ! fieldempty (pconf, "ubound");
312 ubound (idx(prepidx), 1) = \
313 cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
314 {pconf(idx).ubound}.', \
315 cpnel(idx), "UniformOutput", false){:});
317 ubound(idx, 1) = cat (1, pconf.ubound);
321 max_fract_change = fract_prec = NAvec;
323 if (isfield (pconf, "max_fract_change"))
324 idx = ! fieldempty (pconf, "max_fract_change");
326 max_fract_change(idx(prepidx)) = \
327 cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
328 {pconf(idx).max_fract_change}.', \
330 "UniformOutput", false){:});
332 max_fract_change(idx) = [pconf.max_fract_change];
336 if (isfield (pconf, "fract_prec"))
337 idx = ! fieldempty (pconf, "fract_prec");
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){:});
344 fract_prec(idx) = [pconf.fract_prec];
349 diffp(:) = diffp_default;
350 if (isfield (pconf, "diffp"))
351 idx = ! fieldempty (pconf, "diffp");
353 diffp(idx(prepidx)) = \
354 cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
355 {pconf(idx).diffp}.', cpnel(idx), \
356 "UniformOutput", false){:});
358 diffp(idx) = [pconf.diffp];
362 diff_onesided = fixed = falsevec;
364 if (isfield (pconf, "diff_onesided"))
365 idx = ! fieldempty (pconf, "diff_onesided");
367 diff_onesided(idx(prepidx)) = \
369 (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
370 {pconf(idx).diff_onesided}.', cpnel(idx), \
371 "UniformOutput", false){:}));
373 diff_onesided(idx) = logical ([pconf.diff_onesided]);
377 if (isfield (pconf, "fixed"))
378 idx = ! fieldempty (pconf, "fixed");
380 fixed(idx(prepidx)) = \
382 (cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
383 {pconf(idx).fixed}.', cpnel(idx), \
384 "UniformOutput", false){:}));
386 fixed(idx) = logical ([pconf.fixed]);
390 ## use supplied configuration vectors
392 if (isempty (lbound))
394 elseif (any (size (lbound) != sizevec))
395 error ("bounds: wrong dimensions");
398 if (isempty (ubound))
400 elseif (any (size (ubound) != sizevec))
401 error ("bounds: wrong dimensions");
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");
410 if (isempty (fract_prec))
412 elseif (any (size (fract_prec) != sizevec))
413 error ("fract_prec: wrong dimensions");
418 diffp(:) = diffp_default;
420 if (any (size (diffp) != sizevec))
421 error ("diffp: wrong dimensions");
423 diffp(isna (diffp)) = diffp_default;
426 if (isempty (diff_onesided))
427 diff_onesided = falsevec;
429 if (any (size (diff_onesided) != sizevec))
430 error ("diff_onesided: wrong dimensions")
432 diff_onesided(isna (diff_onesided)) = false;
433 diff_onesided = logical (diff_onesided);
439 if (any (size (fixed) != sizevec))
440 error ("fixed: wrong dimensions");
442 fixed(isna (fixed)) = false;
443 fixed = logical (fixed);
447 ## guaranty all (lbound <= ubound)
448 if (any (lbound > ubound))
449 error ("some lower bounds larger than upper bounds");
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
456 ## initial parameters
459 pin = cat (1, cellfun (@ (x, n) reshape (x, n, 1), \
460 fields2cell (pin, pord), cpnel, \
461 "UniformOutput", false){:});
463 pin = cat (1, fields2cell (pin, pord){:});
470 f = @ (p, varargin) \
472 (cellfun (@ reshape, mat2cell (p, ppartidx), \
473 pdims, "UniformOutput", false), \
474 pord, 1), varargin{:});
476 f = @ (p, varargin) \
477 f (cell2struct (num2cell (p), pord, 1), varargin{:});
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");
485 f = @ (p) f (p) - obs;
489 ## jacobian of model function
492 dfdp = @ (p, hook) jacobs (p, f, hook);
494 dfdp = @ (p, hook) __dfdp__ (p, f, hook);
503 (cellfun (@ reshape, mat2cell (p, ppartidx), \
504 pdims, "UniformOutput", false), \
511 (dfdp (cell2struct (num2cell (p), pord, 1), hook), \
516 ## function for general inequality constraints
517 if (f_inequc_pstruct)
519 f_genicstr = @ (p, varargin) \
520 f_genicstr (cell2struct \
521 (cellfun (@ reshape, mat2cell (p, ppartidx), \
522 pdims, "UniformOutput", false), \
523 pord, 1), varargin{:});
525 f_genicstr = @ (p, varargin) \
527 (cell2struct (num2cell (p), pord, 1), varargin{:});
532 possibly_pstruct_f_genicstr = f_genicstr;
534 ## jacobian of general inequality constraints
535 if (df_inequc_pstruct)
537 df_gencstr = @ (p, func, idx, hook) \
542 (cellfun (@ reshape, mat2cell (p, ppartidx), \
543 pdims, "UniformOutput", false), pord, 1), \
547 df_gencstr = @ (p, func, idx, hook) \
550 (df_gencstr (cell2struct (num2cell (p), pord, 1), \
556 ## function for general equality constraints
559 f_genecstr = @ (p, varargin) \
560 f_genecstr (cell2struct \
561 (cellfun (@ reshape, mat2cell (p, ppartidx), \
562 pdims, "UniformOutput", false), \
563 pord, 1), varargin{:});
565 f_genecstr = @ (p, varargin) \
567 (cell2struct (num2cell (p), pord, 1), varargin{:});
572 possibly_pstruct_f_genecstr = f_genecstr;
574 ## jacobian of general equality constraints
577 df_genecstr = @ (p, func, idx, hook) \
582 (cellfun (@ reshape, mat2cell (p, ppartidx), \
583 pdims, "UniformOutput", false), pord, 1), \
587 df_genecstr = @ (p, func, idx, hook) \
590 (df_genecstr (cell2struct (num2cell (p), pord, 1), \
596 ## linear inequality constraints
598 idx = isfield (mc, pord);
599 if (rows (fieldnames (mc)) > sum (idx))
600 error ("unknown fields in structure of linear inequality constraints");
603 mc = zeros (np, rows (vc));
604 mc(idx(prepidx), :) = cat (1, fields2cell (smc, pord(idx)){:});
607 ## linear equality constraints
609 idx = isfield (emc, pord);
610 if (rows (fieldnames (emc)) > sum (idx))
611 error ("unknown fields in structure of linear equality constraints");
614 emc = zeros (np, rows (evc));
615 emc(idx(prepidx), :) = cat (1, fields2cell (semc, pord(idx)){:});
618 ## parameter-related configuration for jacobian functions
619 if (dfdp_pstruct || df_inequc_pstruct || df_equc_pstruct)
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 \
637 (@ reshape, mat2cell (cat (1, x{:}), ppartidx), \
638 pdims, "UniformOutput", false), \
639 num2cell (plabels, 1), "UniformOutput", false){:}), \
642 s_orig_fixed = cell2struct \
643 (cellfun (@ reshape, mat2cell (fixed, ppartidx), \
644 pdims, "UniformOutput", false), pord, 1);
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);
655 #### some further values and checks
657 if (any (fixed & (pin < lbound | pin > ubound)))
658 warning ("some fixed parameters outside bounds");
661 ## dimensions of linear constraints
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");
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");
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");
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;
691 ## note number and initial values of general constraints
692 if (isempty (f_genicstr))
693 pin_cstr.inequ.gen = [];
696 n_genicstr = length (pin_cstr.inequ.gen = f_genicstr (pin));
698 if (isempty (f_genecstr))
699 pin_cstr.equ.gen = [];
702 n_genecstr = length (pin_cstr.equ.gen = f_genecstr (pin));
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);
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);
719 #### handle fixing of parameters
720 orig_lbound = lbound;
721 orig_ubound = ubound;
724 error ("no free parameters");
729 ## backend (returned values and initial parameters)
730 backend = @ (f, pin, hook) \
731 backend_wrapper (backend, fixed, f, pin, hook);
734 f = @ (p, varargin) f (assign (pin, nonfixed, p), varargin{:});
736 ## jacobian of model function
738 dfdp (assign (pin, nonfixed, p), hook)(:, nonfixed);
740 ## function for general inequality constraints
741 f_genicstr = @ (p, varargin) \
742 f_genicstr (assign (pin, nonfixed, p), varargin{:});
744 ## jacobian of general inequality constraints
745 df_gencstr = @ (p, func, idx, hook) \
746 df_gencstr (assign (pin, nonfixed, p), func, idx, hook) \
749 ## function for general equality constraints
750 f_genecstr = @ (p, varargin) \
751 f_genecstr (assign (pin, nonfixed, p), varargin{:});
753 ## jacobian of general equality constraints
754 df_genecstr = @ (p, func, idx, hook) \
755 df_genecstr (assign (pin, nonfixed, p), func, idx, hook) \
758 ## linear inequality constraints
759 vc += mc(fixed, :).' * (tp = pin(fixed));
760 mc = mc(nonfixed, :);
762 ## linear equality constraints
763 evc += emc(fixed, :).' * tp;
764 emc = emc(nonfixed, :);
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);
775 #### supplement constants to jacobian functions
777 ## jacobian of model function
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"}, \
790 dfdp (p, cell2fields \
791 ({diffp, diff_onesided, orig_lbound, orig_ubound, \
792 plabels, assign(orig_fixed, nonfixed, hook.fixed), \
794 {"diffp", "diff_onesided", "lbound", "ubound", \
795 "plabels", "fixed", "h"}, \
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"}, \
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"}, \
821 ## jacobian of general equality constraints
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"}, \
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"}, \
843 #### interfaces to constraints
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);
852 ## concatenate linear inequality and equality constraints
853 mc = cat (2, mc, emc);
854 vc = cat (1, vc, evc);
855 n_lincstr = rows (vc);
857 ## concatenate general inequality and equality constraints
860 nidxi = 1 : n_genicstr;
861 nidxe = n_genicstr + 1 : n_genicstr + n_genecstr;
862 f_gencstr = @ (p, idx, varargin) \
864 f_genicstr (p, idx(nidxi), varargin{:}), \
865 f_genecstr (p, idx(nidxe), varargin{:}));
866 df_gencstr = @ (p, idx, hook) \
868 df_gencstr (p, @ (p, varargin) \
869 possibly_pstruct_f_genicstr \
870 (p, idx(nidxi), varargin{:}), \
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{:}), \
878 setfield (hook, "f", \
879 hook.f(nidxe(idx(nidxe))))));
881 f_gencstr = f_genecstr;
882 df_gencstr = @ (p, idx, hook) \
885 possibly_pstruct_f_genecstr \
886 (p, idx, varargin{:}), \
888 setfield (hook, "f", hook.f(idx)));
891 f_gencstr = f_genicstr;
892 df_gencstr = @ (p, idx, hook) \
895 possibly_pstruct_f_genicstr (p, idx, varargin{:}), \
897 setfield (hook, "f", hook.f(idx)));
899 n_gencstr = n_genicstr + n_genecstr;
901 ## concatenate linear and general constraints, defining the final
902 ## function interfaces
905 nidxh = n_lincstr + 1 : n_lincstr + n_gencstr;
906 f_cstr = @ (p, idx, varargin) \
908 mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), \
909 f_gencstr (p, idx(nidxh), varargin{:}));
910 df_cstr = @ (p, idx, hook) \
912 mc(:, idx(nidxl)).', \
913 df_gencstr (p, idx(nidxh), \
914 setfield (hook, "f", \
917 f_cstr = @ (p, idx, varargin) mc(:, idx).' * p + vc(idx, 1);
918 df_cstr = @ (p, idx, hook) mc(:, idx).';
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;
928 #### prepare interface hook
930 ## passed constraints
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;
940 ## passed values of constraints for initial parameters
941 hook.pin_cstr = pin_cstr;
943 ## passed function for derivative of model function
946 ## passed function for complementary pivoting
947 ## hook.cpiv = cpiv; # set before
949 ## passed value of residual function for initial parameters
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;
962 [p, resid, cvg, outp] = backend (f, pin, hook);
967 (cellfun (@ reshape, mat2cell (p, ppartidx), \
968 pdims, "UniformOutput", false), \
971 p = cell2struct (num2cell (p), pord, 1);
977 function backend = map_matlab_algorithm_names (backend)
980 case "levenberg-marquardt"
981 backend = "lm_svd_feasible";
982 warning ("algorithm 'levenberg-marquardt' mapped to 'lm_svd_feasible'");
987 function backend = map_backend (backend)
990 case "lm_svd_feasible"
991 backend = "__lm_svd__";
993 error ("no backend implemented for algorithm '%s'", backend);
996 backend = str2func (backend);
1000 function [p, resid, cvg, outp] = backend_wrapper (backend, fixed, f, p, hook)
1002 [tp, resid, cvg, outp] = backend (f, p(! fixed), hook);
1008 function lval = assign (lval, lidx, rval)
1014 function ret = __optimget__ (s, name, default)
1016 if (isfield (s, name))
1026 function ret = apply_idx_if_given (ret, varargin)
1029 ret = ret(varargin{1});