forge: main/optim/inst/leasqr.m annotate

annotate main/optim/inst/leasqr.m @ 9930:d30cfca46e8a octave-forge

optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses

author	carandraug
date	Fri, 30 Mar 2012 15:14:48 +0000
parents	88dc7961d556
children	fba8cdd5f9ad

rev	line source
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	1 %% Copyright (C) 1992-1994 Richard Shrager
9930 d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	2 %% Copyright (C) 1992-1994 Arthur Jutan <jutan@charon.engga.uwo.ca>
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	3 %% Copyright (C) 1992-1994 Ray Muzic <rfm2@ds2.uh.cwru.edu>
7951 b0fc599dcdb1 Update copyright notices of functions for 'subpackage' nonlin_residmin. i7tiol parents: 7684 diff changeset	4 %% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de>
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	5 %%
9930 d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	6 %% This program is free software; you can redistribute it and/or modify it under
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	7 %% the terms of the GNU General Public License as published by the Free Software
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	8 %% Foundation; either version 3 of the License, or (at your option) any later
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	9 %% version.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	10 %%
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	11 %% This program is distributed in the hope that it will be useful, but WITHOUT
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	12 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	13 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	14 %% details.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	15 %%
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	16 %% You should have received a copy of the GNU General Public License along with
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	17 %% this program; if not, see <http://www.gnu.org/licenses/>.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	18
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	19 %%function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	20 %% leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	21 %%
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	22 %% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	23 %%
9930 d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	24 %% Version 3.beta
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	25 %% Optional parameters are in braces {}.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	26 %% x = vector or matrix of independent variables.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	27 %% y = vector or matrix of observed values.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	28 %% wt = statistical weights (same dimensions as y). These should be
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	29 %% set to be proportional to (sqrt of var(y))^-1; (That is, the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	30 %% covariance matrix of the data is assumed to be proportional to
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	31 %% diagonal with diagonal equal to (wt.^2)^-1. The constant of
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	32 %% proportionality will be estimated.); default = ones( size (y)).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	33 %% pin = vec of initial parameters to be adjusted by leasqr.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	34 %% dp = fractional increment of p for numerical partial derivatives;
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	35 %% default = .001*ones(size(pin))
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	36 %% dp(j) > 0 means central differences on j-th parameter p(j).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	37 %% dp(j) < 0 means one-sided differences on j-th parameter p(j).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	38 %% dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	39 %% F = name of function in quotes or function handle; the function
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	40 %% shall be of the form y=f(x,p), with y, x, p of the form y, x, pin
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	41 %% as described above.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	42 %% dFdp = name of partial derivative function in quotes or function
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	43 %% handle; default is 'dfdp', a slow but general partial derivatives
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	44 %% function; the function shall be of the form
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	45 %% prt=dfdp(x,f,p,dp,F[,bounds]). For backwards compatibility, the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	46 %% function will only be called with an extra 'bounds' argument if the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	47 %% 'bounds' option is explicitely specified to leasqr (see dfdp.m).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	48 %% stol = scalar tolerance on fractional improvement in scalar sum of
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	49 %% squares = sum((wt.*(y-f))^2); default stol = .0001;
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	50 %% niter = scalar maximum number of iterations; default = 20;
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	51 %% options = structure, currently recognized fields are 'fract_prec',
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	52 %% 'max_fract_change', 'inequc', 'bounds', and 'equc'. For backwards
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	53 %% compatibility, 'options' can also be a matrix whose first and
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	54 %% second column contains the values of 'fract_prec' and
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	55 %% 'max_fract_change', respectively.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	56 %% Field 'options.fract_prec': column vector (same length as 'pin')
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	57 %% of desired fractional precisions in parameter estimates.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	58 %% Iterations are terminated if change in parameter vector (chg)
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	59 %% relative to current parameter estimate is less than their
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	60 %% corresponding elements in 'options.fract_prec' [ie. all (abs
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	61 %% (chg) < abs (options.fract_prec .* current_parm_est))] on two
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	62 %% consecutive iterations, default = zeros().
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	63 %% Field 'options.max_fract_change': column vector (same length as
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	64 %% 'pin) of maximum fractional step changes in parameter vector.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	65 %% Fractional change in elements of parameter vector is constrained to
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	66 %% be at most 'options.max_fract_change' between sucessive iterations.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	67 %% [ie. abs(chg(i))=abs(min([chg(i)
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	68 %% options.max_fract_change(i)*current param estimate])).], default =
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	69 %% Inf*ones().
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	70 %% Field 'options.inequc': cell-array containing up to four entries,
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	71 %% two entries for linear inequality constraints and/or one or two
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	72 %% entries for general inequality constraints. Initial parameters
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	73 %% must satisfy these constraints. Either linear or general
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	74 %% constraints may be the first entries, but the two entries for
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	75 %% linear constraints must be adjacent and, if two entries are given
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	76 %% for general constraints, they also must be adjacent. The two
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	77 %% entries for linear constraints are a matrix (say m) and a vector
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	78 %% (say v), specifying linear inequality constraints of the form
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	79 %% `m.' * parameters + v >= 0'. If the constraints are just bounds,
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	80 %% it is suggested to specify them in 'options.bounds' instead,
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	81 %% since then some sanity tests are performed, and since the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	82 %% function 'dfdp.m' is guarantied not to violate constraints during
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	83 %% determination of the numeric gradient only for those constraints
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	84 %% specified as 'bounds' (possibly with violations due to a certain
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	85 %% inaccuracy, however, except if no constraints except bounds are
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	86 %% specified). The first entry for general constraints must be a
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	87 %% differentiable vector valued function (say h), specifying general
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	88 %% inequality constraints of the form `h (p[, idx]) >= 0'; p is the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	89 %% column vector of optimized paraters and the optional argument idx
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	90 %% is a logical index. h has to return the values of all constraints
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	91 %% if idx is not given, and has to return only the indexed
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	92 %% constraints if idx is given (so computation of the other
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	93 %% constraints can be spared). If a second entry for general
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	94 %% constraints is given, it must be a function (say dh) which
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	95 %% returnes a matrix whos rows contain the gradients of the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	96 %% constraint function h with respect to the optimized parameters.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	97 %% It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	98 %% the column vector of optimized parameters, and idx is a logical
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	99 %% index --- only the rows indexed by idx must be returned (so
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	100 %% computation of the others can be spared). The other arguments of
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	101 %% dh are for the case that dh computes numerical gradients: vh is
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	102 %% the column vector of the current values of the constraint
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	103 %% function h, with idx already applied. h is a function h (p) to
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	104 %% compute the values of the constraints for parameters p, it will
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	105 %% return only the values indexed by idx. dp is a suggestion for
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	106 %% relative step width, having the same value as the argument 'dp'
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	107 %% of leasqr above. If bounds were specified to leasqr, they are
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	108 %% provided in the argument bounds of dh, to enable their
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	109 %% consideration in determination of numerical gradients. If dh is
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	110 %% not specified to leasqr, numerical gradients are computed in the
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	111 %% same way as with 'dfdp.m' (see above). If some constraints are
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	112 %% linear, they should be specified as linear constraints (or
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	113 %% bounds, if applicable) for reasons of performance, even if
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	114 %% general constraints are also specified.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	115 %% Field 'options.bounds': two-column-matrix, one row for each
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	116 %% parameter in 'pin'. Each row contains a minimal and maximal value
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	117 %% for each parameter. Default: [-Inf, Inf] in each row. If this
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	118 %% field is used with an existing user-side function for 'dFdp'
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	119 %% (see above) the functions interface might have to be changed.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	120 %% Field 'options.equc': equality constraints, specified the same
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	121 %% way as inequality constraints (see field 'options.inequc').
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	122 %% Initial parameters must satisfy these constraints.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	123 %% Note that there is possibly a certain inaccuracy in honoring
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	124 %% constraints, except if only bounds are specified.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	125 %% _Warning_: If constraints (or bounds) are set, returned guesses
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	126 %% of corp, covp, and Z are generally invalid, even if no constraints
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	127 %% are active for the final parameters. If equality constraints are
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	128 %% specified, corp, covp, and Z are not guessed at all.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	129 %% Field 'options.cpiv': Function for complementary pivot algorithm
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	130 %% for inequality constraints, default: cpiv_bard. No different
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	131 %% function is supplied.
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	132 %%
9930 d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	133 %% OUTPUT VARIABLES
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	134 %% f = column vector of values computed: f = F(x,p).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	135 %% p = column vector trial or final parameters. i.e, the solution.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	136 %% cvg = scalar: = 1 if convergence, = 0 otherwise.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	137 %% iter = scalar number of iterations used.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	138 %% corp = correlation matrix for parameters.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	139 %% covp = covariance matrix of the parameters.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	140 %% covr = diag(covariance matrix of the residuals).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	141 %% stdresid = standardized residuals.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	142 %% Z = matrix that defines confidence region (see comments in the source).
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	143 %% r2 = coefficient of multiple determination, intercept form.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	144 %%
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	145 %% Not suitable for non-real residuals.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	146 %%
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	147 %% References:
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	148 %% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
d30cfca46e8a optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses carandraug parents: 9928 diff changeset	149 %% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	150
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	151 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	152 leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	153
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	154 %% The following two blocks of comments are chiefly from the original
06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	155 %% version for Matlab. For later changes the logs of the Octave Forge
06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	156 %% svn repository should also be consulted.
06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	157
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	158 %% A modified version of Levenberg-Marquardt
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	159 %% Non-Linear Regression program previously submitted by R.Schrager.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	160 %% This version corrects an error in that version and also provides
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	161 %% an easier to use version with automatic numerical calculation of
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	162 %% the Jacobian Matrix. In addition, this version calculates statistics
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	163 %% such as correlation, etc....
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	164 %%
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	165 %% Version 3 Notes
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	166 %% Errors in the original version submitted by Shrager (now called
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	167 %% version 1) and the improved version of Jutan (now called version 2)
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	168 %% have been corrected.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	169 %% Additional features, statistical tests, and documentation have also been
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	170 %% included along with an example of usage. BEWARE: Some the the input and
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	171 %% output arguments were changed from the previous version.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	172 %%
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	173 %% Ray Muzic <rfm2@ds2.uh.cwru.edu>
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	174 %% Arthur Jutan <jutan@charon.engga.uwo.ca>
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	175
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	176 %% Richard I. Shrager (301)-496-1122
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	177 %% Modified by A.Jutan (519)-679-2111
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	178 %% Modified by Ray Muzic 14-Jul-1992
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	179 %% 1) add maxstep feature for limiting changes in parameter estimates
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	180 %% at each step.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	181 %% 2) remove forced columnization of x (x=x(:)) at beginning. x
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	182 %% could be a matrix with the ith row of containing values of
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	183 %% the independent variables at the ith observation.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	184 %% 3) add verbose option
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	185 %% 4) add optional return arguments covp, stdresid, chi2
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	186 %% 5) revise estimates of corp, stdev
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	187 %% Modified by Ray Muzic 11-Oct-1992
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	188 %% 1) revise estimate of Vy. remove chi2, add Z as return values
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	189 %% (later remark: the current code contains no variable Vy)
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	190 %% Modified by Ray Muzic 7-Jan-1994
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	191 %% 1) Replace ones(x) with a construct that is compatible with versions
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	192 %% newer and older than v 4.1.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	193 %% 2) Added global declaration of verbose (needed for newer than v4.x)
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	194 %% 3) Replace return value var, the variance of the residuals
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	195 %% with covr, the covariance matrix of the residuals.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	196 %% 4) Introduce options as 10th input argument. Include
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	197 %% convergence criteria and maxstep in it.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	198 %% 5) Correct calculation of xtx which affects coveraince estimate.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	199 %% 6) Eliminate stdev (estimate of standard deviation of
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	200 %% parameter estimates) from the return values. The covp is a
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	201 %% much more meaningful expression of precision because it
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	202 %% specifies a confidence region in contrast to a confidence
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	203 %% interval.. If needed, however, stdev may be calculated as
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	204 %% stdev=sqrt(diag(covp)).
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	205 %% 7) Change the order of the return values to a more logical order.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	206 %% 8) Change to more efficent algorithm of Bard for selecting epsL.
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	207 %% 9) Tighten up memory usage by making use of sparse matrices (if
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	208 %% MATLAB version >= 4.0) in computation of covp, corp, stdresid.
9928 88dc7961d556 maint: change text file encoding to UTF-8. No code change carandraug parents: 8717 diff changeset	209 %% Modified by Francesco Potortì
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	210 %% for use in Octave
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	211 %% Added linear inequality constraints with quadratic programming to
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	212 %% this file and special case bounds to this file and to dfdp.m
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	213 %% (24-Feb-2010) and later also general inequality constraints
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	214 %% (12-Apr-2010) (Reference: Bard, Y., 'An eclectic approach to
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	215 %% nonlinear programming', Proc. ANU Sem. Optimization, Canberra,
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	216 %% Austral. Nat. Univ.). Differences from the reference: adaption to
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	217 %% svd-based algorithm, linesearch or stepwidth adaptions to ensure
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	218 %% decrease in objective function omitted to rather start a new
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	219 %% overall cycle with a new epsL, some performance gains from linear
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	220 %% constraints even if general constraints are specified. Equality
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	221 %% constraints also implemented. Olaf Till
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	222 %% Now split into files leasqr.m and __lm_svd__.m.
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	223
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	224 %% needed for some anonymous functions
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	225 if (exist ('ifelse') ~= 5)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	226 ifelse = @ scalar_ifelse;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	227 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	228
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	229 __plot_cmds__ (); % flag persistent variables invalid
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	230
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	231 global verbose;
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	232
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	233 %% argument processing
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	234 %%
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	235
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	236 if (nargin > 8)
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	237 if (ischar (dFdp))
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	238 dfdp = str2func (dFdp);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	239 else
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	240 dfdp = dFdp;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	241 end
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	242 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	243
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	244 if (nargin <= 7) dp=.001(pin0+1); end %DT
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	245 if (nargin <= 6) wt = ones (size (y)); end % SMB modification
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	246 if (nargin <= 5) niter = []; end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	247 if (nargin == 4) stol=.0001; end
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	248 if (ischar (F)) F = str2func (F); end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	249 %%
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	250
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	251 if (any (size (y) ~= size (wt)))
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	252 error ('dimensions of observations and weights do not match');
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	253 end
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	254 wtl = wt(:);
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	255 pin=pin(:); dp=dp(:); %change all vectors to columns
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	256 [rows_y, cols_y] = size (y);
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	257 m = rows_y * cols_y; n=length(pin);
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	258 f_pin = F (x, pin);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	259 if (any (size (f_pin) ~= size (y)))
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	260 error ('dimensions of returned values of model function and of observations do not match');
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	261 end
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	262 f_pin = y - f_pin;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	263
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	264 dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, F);
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	265
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	266 %% processing of 'options'
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	267 pprec = zeros (n, 1);
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	268 maxstep = Inf * ones (n, 1);
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	269 have_gencstr = false; % no general constraints
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	270 have_genecstr = false; % no general equality constraints
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	271 n_gencstr = 0;
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	272 mc = zeros (n, 0);
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	273 vc = zeros (0, 1); rv = 0;
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	274 emc = zeros (n, 0);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	275 evc = zeros (0, 1); erv = 0;
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	276 bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	277 pin_cstr.inequ.lin_except_bounds = [];;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	278 pin_cstr.inequ.gen = [];;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	279 pin_cstr.equ.lin = [];;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	280 pin_cstr.equ.gen = [];;
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	281 dfdp_bounds = {};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	282 cpiv = @ cpiv_bard;
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	283 eq_idx = []; % numerical index for equality constraints in all
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	284 % constraints, later converted to
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	285 % logical index
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	286 if (nargin > 9)
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	287 if (ismatrix (options)) % backwards compatibility
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	288 tp = options;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	289 options = struct ('fract_prec', tp(:, 1));
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	290 if (columns (tp) > 1)
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	291 options.max_fract_change = tp(:, 2);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	292 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	293 end
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	294 if (isfield (options, 'cpiv') && ~isempty (options.cpiv))
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	295 %% As yet there is only one cpiv function distributed with leasqr,
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	296 %% but this may change; the algorithm of cpiv_bard is said to be
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	297 %% relatively fast, but may have disadvantages.
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	298 if (ischar (options.cpiv))
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	299 cpiv = str2func (options.cpiv);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	300 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	301 cpiv = options.cpiv;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	302 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	303 end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	304 if (isfield (options, 'fract_prec'))
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	305 pprec = options.fract_prec;
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	306 if (any (size (pprec) ~= [n, 1]))
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	307 error ('fractional precisions: wrong dimensions');
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	308 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	309 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	310 if (isfield (options, 'max_fract_change'))
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	311 maxstep = options.max_fract_change;
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	312 if (any (size (maxstep) ~= [n, 1]))
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	313 error ('maximum fractional step changes: wrong dimensions');
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	314 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	315 end
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	316 if (isfield (options, 'inequc'))
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	317 inequc = options.inequc;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	318 if (ismatrix (inequc{1}))
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	319 mc = inequc{1};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	320 vc = inequc{2};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	321 if (length (inequc) > 2)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	322 have_gencstr = true;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	323 f_gencstr = inequc{3};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	324 if (length (inequc) > 3)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	325 df_gencstr = inequc{4};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	326 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	327 df_gencstr = @ dcdp;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	328 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	329 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	330 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	331 lid = 0; % no linear constraints
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	332 have_gencstr = true;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	333 f_gencstr = inequc{1};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	334 if (length (inequc) > 1)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	335 if (ismatrix (inequc{2}))
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	336 lid = 2;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	337 df_gencstr = @ dcdp;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	338 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	339 df_gencstr = inequc{2};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	340 if (length (inequc) > 2)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	341 lid = 3;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	342 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	343 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	344 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	345 df_gencstr = @ dcdp;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	346 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	347 if (lid)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	348 mc = inequc{lid};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	349 vc = inequc{lid + 1};
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	350 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	351 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	352 if (have_gencstr)
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	353 if (ischar (f_gencstr))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	354 f_gencstr = str2func (f_gencstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	355 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	356 tp = f_gencstr (pin);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	357 n_gencstr = length (tp);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	358 f_gencstr = @ (p, idx) tf_gencstr (p, idx, f_gencstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	359 if (ischar (df_gencstr))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	360 df_gencstr = str2func (df_gencstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	361 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	362 if (strcmp (func2str (df_gencstr), 'dcdp'))
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	363 df_gencstr = @ (f, p, dp, idx, db) ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	364 df_gencstr (f(idx), p, dp, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	365 @ (tp) f_gencstr (tp, idx), db{:});
7411 2cf09f65806b - optim_problems.m: added unconstrained problems schittkowski_281, _289, i7tiol parents: 7248 diff changeset	366 else
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	367 df_gencstr = @ (f, p, dp, idx, db) ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	368 df_gencstr (f(idx), p, dp, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	369 @ (tp) f_gencstr (tp, idx), idx, db{:});
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	370 end
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	371 end
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	372 [rm, cm] = size (mc);
4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	373 [rv, cv] = size (vc);
4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	374 if (rm ~= n \|\| cm ~= rv \|\| cv ~= 1)
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	375 error ('linear inequality constraints: wrong dimensions');
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	376 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	377 pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	378 if (have_gencstr)
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	379 pin_cstr.inequ.gen = tp;
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	380 end
4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	381 end
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	382 if (isfield (options, 'equc'))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	383 equc = options.equc;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	384 if (ismatrix (equc{1}))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	385 emc = equc{1};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	386 evc = equc{2};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	387 if (length (equc) > 2)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	388 have_genecstr = true;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	389 f_genecstr = equc{3};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	390 if (length (equc) > 3)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	391 df_genecstr = equc{4};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	392 else
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	393 df_genecstr = @ dcdp;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	394 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	395 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	396 else
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	397 lid = 0; % no linear constraints
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	398 have_genecstr = true;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	399 f_genecstr = equc{1};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	400 if (length (equc) > 1)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	401 if (ismatrix (equc{2}))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	402 lid = 2;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	403 df_genecstr = @ dcdp;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	404 else
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	405 df_genecstr = equc{2};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	406 if (length (equc) > 2)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	407 lid = 3;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	408 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	409 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	410 else
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	411 df_genecstr = @ dcdp;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	412 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	413 if (lid)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	414 emc = equc{lid};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	415 evc = equc{lid + 1};
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	416 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	417 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	418 if (have_genecstr)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	419 if (ischar (f_genecstr))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	420 f_genecstr = str2func (f_genecstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	421 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	422 tp = f_genecstr (pin);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	423 n_genecstr = length (tp);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	424 f_genecstr = @ (p, idx) tf_gencstr (p, idx, f_genecstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	425 if (ischar (df_genecstr))
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	426 df_genecstr = str2func (df_genecstr);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	427 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	428 if (strcmp (func2str (df_genecstr), 'dcdp'))
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	429 df_genecstr = @ (f, p, dp, idx, db) ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	430 df_genecstr (f, p, dp, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	431 @ (tp) f_genecstr (tp, idx), db{:});
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	432 else
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	433 df_genecstr = @ (f, p, dp, idx, db) ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	434 df_genecstr (f, p, dp, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	435 @ (tp) f_genecstr (tp, idx), idx, db{:});
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	436 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	437 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	438 [erm, ecm] = size (emc);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	439 [erv, ecv] = size (evc);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	440 if (erm ~= n \|\| ecm ~= erv \|\| ecv ~= 1)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	441 error ('linear equality constraints: wrong dimensions');
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	442 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	443 pin_cstr.equ.lin = emc.' * pin + evc;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	444 if (have_genecstr)
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	445 pin_cstr.equ.gen = tp;
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	446 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	447 end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	448 if (isfield (options, 'bounds'))
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	449 bounds = options.bounds;
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	450 if (any (size (bounds) ~= [n, 2]))
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	451 error ('bounds: wrong dimensions');
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	452 end
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	453 idx = bounds(:, 1) > bounds(:, 2);
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	454 tp = bounds(idx, 2);
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	455 bounds(idx, 2) = bounds(idx, 1);
805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	456 bounds(idx, 1) = tp;
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	457 %% It is possible to take this decision here, since this frontend
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	458 %% is used only with one certain backend. The backend will check
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	459 %% this again; but it will not reference 'dp' in its message,
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	460 %% thats why the additional check here.
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	461 idx = bounds(:, 1) == bounds(:, 2);
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	462 if (any (idx))
7248 6f03b0449b74 optim_problems: added Schittkowski no. 327, 372, 373 (general constraints). i7tiol parents: 7240 diff changeset	463 warning ('leasqr:constraints', 'lower and upper bounds identical for some parameters, setting the respective elements of dp to zero');
6034 805204dbd13a Support bounds (from Olaf Till and Roessli Bertrand hauberg parents: 4737 diff changeset	464 dp(idx) = 0;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	465 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	466 %%
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	467 tp = eye (n);
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	468 lidx = ~isinf (bounds(:, 1));
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	469 uidx = ~isinf (bounds(:, 2));
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	470 mc = cat (2, mc, tp(:, lidx), - tp(:, uidx));
4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	471 vc = cat (1, vc, - bounds(lidx, 1), bounds(uidx, 2));
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	472 [rm, cm] = size (mc);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	473 [rv, cv] = size (vc);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	474 dfdp_bounds = {bounds};
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	475 dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	476 F, bounds);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	477 end
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	478 %% concatenate inequality and equality constraint functions, mc, and
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	479 %% vc; update eq_idx, rv, n_gencstr, have_gencstr
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	480 if (erv > 0)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	481 mc = cat (2, mc, emc);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	482 vc = cat (1, vc, evc);
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	483 eq_idx = rv + 1 : rv + erv;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	484 rv = rv + erv;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	485 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	486 if (have_genecstr)
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	487 eq_idx = cat (2, eq_idx, ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	488 rv + n_gencstr + 1 : rv + n_gencstr + n_genecstr);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	489 nidxi = 1 : n_gencstr;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	490 nidxe = n_gencstr + 1 : n_gencstr + n_genecstr;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	491 n_gencstr = n_gencstr + n_genecstr;
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	492 if (have_gencstr)
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	493 f_gencstr = @ (p, idx) cat (1, ...
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	494 f_gencstr (p, idx(nidxi)), ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	495 f_genecstr (p, idx(nidxe)));
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	496 df_gencstr = @ (f, p, dp, idx, db) ...
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	497 cat (1, ...
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	498 df_gencstr (f(nidxi), p, dp, idx(nidxi), db), ...
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	499 df_genecstr (f(nidxe), p, dp, idx(nidxe), db));
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	500 else
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	501 f_gencstr = f_genecstr;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	502 df_gencstr = df_genecstr;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	503 have_gencstr = true;
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	504 end
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	505 end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	506 end
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	507 if (have_gencstr)
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	508 nidxl = 1:rv;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	509 nidxh = rv+1:rv+n_gencstr;
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	510 f_cstr = @ (p, idx) ...
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	511 cat (1, mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), ...
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	512 f_gencstr (p, idx(nidxh)));
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	513 %% in the case of this interface, diffp is already zero at fixed;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	514 %% also in this special case, dfdp_bounds can be filled in directly
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	515 %% --- otherwise it would be a field of hook in the called function
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	516 df_cstr = @ (p, idx, dfdp_hook) ...
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	517 cat (1, mc(:, idx(nidxl)).', ...
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	518 df_gencstr (dfdp_hook.f(nidxh), p, dp, ...
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	519 idx(nidxh), ...
7411 2cf09f65806b - optim_problems.m: added unconstrained problems schittkowski_281, _289, i7tiol parents: 7248 diff changeset	520 dfdp_bounds));
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	521 else
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	522 f_cstr = @ (p, idx) mc(:, idx).' * p + vc(idx, 1);
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	523 df_cstr = @ (p, idx, dfdp_hook) mc(:, idx).';
7033 9ae2523ecbcd Fix: initial parameters should not be corrected to bounds if other constraints are given. i7tiol parents: 6992 diff changeset	524 end
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	525
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	526
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	527
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	528 %% in a general interface, check for all(fixed) here
6768 4b92ae35f351 leasqr.m i7tiol parents: 6722 diff changeset	529
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	530 %% passed constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	531 hook.mc = mc; % matrix of linear constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	532 hook.vc = vc; % vector of linear constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	533 hook.f_cstr = f_cstr; % function of all constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	534 hook.df_cstr = df_cstr; % function of derivatives of all constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	535 hook.n_gencstr = n_gencstr; % number of non-linear constraints
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	536 hook.eq_idx = false (size (vc, 1) + n_gencstr, 1);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	537 hook.eq_idx(eq_idx) = true; % logical index of equality constraints in
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	538 % all constraints
7953 546070338f68 Fix new bounds handling. Fix demo 2. i7tiol parents: 7951 diff changeset	539 hook.lbound = bounds(:, 1); % bounds, subset of linear inequality
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	540 % constraints in mc and vc
7953 546070338f68 Fix new bounds handling. Fix demo 2. i7tiol parents: 7951 diff changeset	541 hook.ubound = bounds(:, 2);
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	542
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	543 %% passed values of constraints for initial parameters
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	544 hook.pin_cstr = pin_cstr;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	545
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	546 %% passed derivative of model function
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	547 hook.dfdp = dFdp;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	548
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	549 %% passed function for complementary pivoting
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	550 hook.cpiv = cpiv;
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	551
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	552 %% passed value of residual function for initial parameters
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	553 hook.f_pin = f_pin;
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	554
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	555 %% passed options
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	556 hook.max_fract_change = maxstep;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	557 hook.fract_prec = pprec;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	558 hook.TolFun = stol;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	559 hook.MaxIter = niter;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	560 hook.weights = wt;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	561 hook.fixed = dp == 0;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	562 if (verbose)
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	563 hook.Display = 'iter';
8717 da32bc7df365 Move private functions into private/. Add corresponding workarounds for bug #31484 (Octave <= 3.2.3). i7tiol parents: 8282 diff changeset	564 __plot_cmds__ = @ __plot_cmds__; # for bug #31484 (Octave <= 3.2.4)
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	565 hook.plot_cmd = @ (f) __plot_cmds__ (x, y, y - f);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	566 else
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	567 hook.Display = 'off';
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	568 end
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	569
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	570 %% only preliminary, for testing
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	571 hook.testing = false;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	572 hook.new_s = false;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	573 if (exist ('options'))
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	574 if (isfield (options, 'testing'))
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	575 hook.testing = options.testing;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	576 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	577 if (isfield (options, 'new_s'))
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	578 hook.new_s = options.new_s;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	579 end
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	580 end
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	581
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	582 [p, resid, cvg, outp] = __lm_svd__ (@ (p) y - F (x, p), pin, hook);
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	583 f = y - resid;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	584 iter = outp.niter;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	585 cvg = cvg > 0;
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	586
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	587 if (~cvg) disp(' CONVERGENCE NOT ACHIEVED! '); end
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	588
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	589 if (~(verbose \|\| nargout > 4)) return; end
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	590
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	591 yl = y(:);
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	592 f = f(:);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	593 %% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	594 %% re-evaluate the Jacobian at optimal values
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	595 jac = dFdp (p, struct ('f', f));
d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	596 msk = ~hook.fixed;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	597 n = sum(msk); % reduce n to equal number of estimated parameters
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	598 jac = jac(:, msk); % use only fitted parameters
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	599
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	600 %% following section is Ray Muzic's estimate for covariance and correlation
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	601 %% assuming covariance of data is a diagonal matrix proportional to
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	602 %% diag(1/wt.^2).
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	603 %% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	604
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	605 tp = wtl.^2;
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	606 if (exist('sparse')) % save memory
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	607 Q = sparse (1:m, 1:m, 1 ./ tp);
06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	608 Qinv = sparse (1:m, 1:m, tp);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	609 else
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	610 Q = diag (ones (m, 1) ./ tp);
06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	611 Qinv = diag (tp);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	612 end
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	613 resid = resid(:); % un-weighted residuals
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	614 if (~isreal (resid)) error ('residuals are not real'); end
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	615 tp = resid.' * Qinv * resid;
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	616 covr = (tp / m) * Q; %covariance of residuals
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	617
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	618 %% Matlab compatibility and avoiding recomputation make the following
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	619 %% logic clumsy.
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	620 compute = 1;
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	621 if (m <= n \|\| any (eq_idx))
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	622 compute = 0;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	623 else
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	624 Qinv = ((m - n) / tp) * Qinv;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	625 %% simplified Eq. 7-5-13, Bard; cov of parm est, inverse; outer
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	626 %% parantheses contain inverse of guessed covariance matrix of data
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	627 covpinv = jac.' * Qinv * jac;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	628 if (exist ('rcond') && rcond (covpinv) <= eps)
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	629 compute = 0;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	630 elseif (rank (covpinv) < n)
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	631 %% above test is not equivalent to 'rcond' and may unnecessarily
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	632 %% reject some matrices
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	633 compute = 0;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	634 end
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	635 end
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	636 if (compute)
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	637 covp = inv (covpinv);
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	638 d=sqrt(diag(covp));
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	639 corp = covp ./ (d * d.');
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	640 else
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	641 covp = NA * ones (n);
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	642 corp = covp;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	643 end
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	644
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	645 if (exist('sparse'))
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	646 covr=spdiags(covr,0);
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	647 else
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	648 covr=diag(covr); % convert returned values to
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	649 % compact storage
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	650 end
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	651 covr = reshape (covr, rows_y, cols_y);
bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	652 stdresid = resid .* abs (wtl) / sqrt (tp / m); % equivalent to resid ./
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	653 % sqrt (covr)
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	654 stdresid = reshape (stdresid, rows_y, cols_y);
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	655
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	656 if (~(verbose \|\| nargout > 8)) return; end
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	657
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	658 if (m > n && ~any (eq_idx))
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	659 Z = ((m - n) / (n * resid.' * Qinv * resid)) * covpinv;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	660 else
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	661 Z = NA * ones (n);
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	662 end
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	663
24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	664 %%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	665 %%disp('Alternate estimate of cov. of param. est.')
6469 06fc77624665 main/optim/inst/leasqr.m: correct covariance matrices and standard residuals i7tiol parents: 6466 diff changeset	666 %%acovp=resid'Qinvresid/(m-n)inv(jac'Qinv*jac);
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	667
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	668 if (~(verbose \|\| nargout > 9)) return; end
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	669
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	670 %%Calculate R^2, intercept form
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	671 %%
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	672 tp = sumsq (yl - mean (yl));
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	673 if (tp > 0)
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	674 r2 = 1 - sumsq (resid) / tp;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	675 else
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	676 r2 = NA;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	677 end
2370 24d6a5cdedfe Changed the directory structure to match the package system hauberg parents: diff changeset	678
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	679 %% if someone has asked for it, let them have it
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	680 %%
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	681 if (verbose)
7684 d4f8cb9c7a93 Added nonlin_residmin and nonlin_curvefit with a common backend with leasqr. i7tiol parents: 7411 diff changeset	682 __plot_cmds__ (x, y, f);
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	683 disp(' Least Squares Estimates of Parameters')
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	684 disp(p.')
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	685 disp(' Correlation matrix of parameters estimated')
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	686 disp(corp)
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	687 disp(' Covariance matrix of Residuals' )
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	688 disp(covr)
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	689 disp(' Correlation Coefficient R^2')
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	690 disp(r2)
7110 bef6c2eea087 Removed restriction that observations must be a vector with length == rows (independents). i7tiol parents: 7059 diff changeset	691 fprintf(' 95%% conf region: F(0.05)(%.0f,%.0f)>= delta_pvec.%sZdelta_pvec\n', n, m - n, char (39)); % works with ' and '
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	692 Z
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	693 %% runs test according to Bard. p 201.
6476 11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	694 n1 = sum (resid > 0);
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	695 n2 = sum (resid < 0);
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	696 nrun=sum(abs(diff(resid > 0)))+1;
11209d9a9e63 main/optim/inst/leasqr.m: only compute additional data if requested and if reliably possible; reject non-real data; accept row-vector for x i7tiol parents: 6469 diff changeset	697 if ((n1 > 10) && (n2 > 10)) % sufficent data for test?
6466 24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	698 zed=(nrun-(2n1n2/(n1+n2)+1)+0.5)/(2n1n2(2n1*n2-n1-n2)...
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	699 /((n1+n2)^2*(n1+n2-1)));
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	700 if (zed < 0)
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	701 prob = erfc(-zed/sqrt(2))/2*100;
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	702 disp([num2str(prob),'% chance of fewer than ',num2str(nrun),' runs.']);
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	703 else
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	704 prob = erfc(zed/sqrt(2))/2*100;
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	705 disp([num2str(prob),'% chance of greater than ',num2str(nrun),' runs.']);
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	706 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	707 end
24e2fe3ba878 optim/inst/leasqr.m, optim/inst/dfdp.m: style cleanup, restored Matlab compatibility i7tiol parents: 6034 diff changeset	708 end
6778 424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	709
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	710 function ret = tf_gencstr (p, idx, f)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	711
7059 c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	712 %% necessary since user function f_gencstr might return [] or a row
c83d46a18b74 leasqr: equality constraints implemented, bugfix in general inequality constraints i7tiol parents: 7033 diff changeset	713 %% vector
6992 45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	714
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	715 ret = f (p, idx);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	716 if (isempty (ret))
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	717 ret = zeros (0, 1);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	718 elseif (size (ret, 2) > 1)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	719 ret = ret(:);
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	720 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	721
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	722 function fval = scalar_ifelse (cond, tval, fval)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	723
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	724 %% needed for some anonymous functions, builtin ifelse only available
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	725 %% in Octave > 3.2; we need only the scalar case here
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	726
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	727 if (cond)
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	728 fval = tval;
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	729 end
45a51449b113 - leasqr.m: i7tiol parents: 6875 diff changeset	730
6778 424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	731 %!demo
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	732 %! % Define functions
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	733 %! leasqrfunc = @(x, p) p(1) * exp (-p(2) * x);
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	734 %! leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)x), -p(1)x.exp(-p(2)x)];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	735 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	736 %! % generate test data
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	737 %! t = [1:10:100]';
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	738 %! p = [1; 0.1];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	739 %! data = leasqrfunc (t, p);
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	740 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	741 %! rnd = [0.352509; -0.040607; -1.867061; -1.561283; 1.473191; ...
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	742 %! 0.580767; 0.841805; 1.632203; -0.179254; 0.345208];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	743 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	744 %! % add noise
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	745 %! % wt1 = 1 /sqrt of variances of data
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	746 %! % 1 / wt1 = sqrt of var = standard deviation
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	747 %! wt1 = (1 + 0 * t) ./ sqrt (data);
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	748 %! data = data + 0.05 * rnd ./ wt1;
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	749 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	750 %! % Note by Thomas Walter <walter@pctc.chemie.uni-erlangen.de>:
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	751 %! %
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	752 %! % Using a step size of 1 to calculate the derivative is WRONG !!!!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	753 %! % See numerical mathbooks why.
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	754 %! % A derivative calculated from central differences need: s
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	755 %! % step = 0.001...1.0e-8
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	756 %! % And onesided derivative needs:
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	757 %! % step = 1.0e-5...1.0e-8 and may be still wrong
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	758 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	759 %! F = leasqrfunc;
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	760 %! dFdp = leasqrdfdp; % exact derivative
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	761 %! % dFdp = dfdp; % estimated derivative
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	762 %! dp = [0.001; 0.001];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	763 %! pin = [.8; .05];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	764 %! stol=0.001; niter=50;
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	765 %! minstep = [0.01; 0.01];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	766 %! maxstep = [0.8; 0.8];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	767 %! options = [minstep, maxstep];
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	768 %!
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	769 %! global verbose;
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	770 %! verbose = 1;
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	771 %! [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	772 %! leasqr (t, data, pin, F, stol, niter, wt1, dp, dFdp, options);
424555ba0b02 Move 'leasrdemo' into a demo block in 'leasqr' hauberg parents: 6768 diff changeset	773
6875 9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	774 %!demo
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	775 %! %% Example for linear inequality constraints.
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	776 %! %% model function:
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	777 %! F = @ (x, p) p(1) * exp (p(2) * x);
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	778 %! %% independents and dependents:
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	779 %! x = 1:5;
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	780 %! y = [1, 2, 4, 7, 14];
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	781 %! %% initial values:
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	782 %! init = [.25; .25];
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	783 %! %% other configuration (default values):
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	784 %! tolerance = .0001;
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	785 %! max_iterations = 20;
7953 546070338f68 Fix new bounds handling. Fix demo 2. i7tiol parents: 7951 diff changeset	786 %! weights = ones (1, 5);
6875 9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	787 %! dp = [.001; .001]; % bidirectional numeric gradient stepsize
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	788 %! dFdp = 'dfdp'; % function for gradient (numerical)
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	789 %!
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	790 %! %% linear constraints, A.' * parametervector + B >= 0
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	791 %! A = [1; -1]; B = 0; % p(1) >= p(2);
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	792 %! options.inequc = {A, B};
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	793 %!
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	794 %! %% start leasqr, be sure that 'verbose' is not set
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	795 %! global verbose; verbose = false;
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	796 %! [f, p, cvg, iter] = ...
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	797 %! leasqr (x, y, init, F, tolerance, max_iterations, ...
9264531021cf Changed version number and index file, include a constraints example into leasqr.m i7tiol parents: 6778 diff changeset	798 %! weights, dp, dFdp, options)

Mercurial > forge

annotate main/optim/inst/leasqr.m @ 9930:d30cfca46e8a octave-forge