On Wed, Aug 05, 2009 at 10:31:22AM +0200, Olaf Till wrote:
> On Tue, Aug 04, 2009 at 09:57:28PM +0200, Søren Hauberg wrote:
> > tir, 04 08 2009 kl. 13:02 +0200, skrev Olaf Till:
> > > A lot of
> > > cosmetic changes have been made to make leasqr.m better editable under
> > > emacs, so almost each line is changed and its not much use making
> > > diffs :-(.
> >
> > Then could you please start out by doing a diff where you do not change
> > coding style. Once we have the content in place, you can always change
> > the coding style.
>
> Ok, attached is an SVN diff inside main/optim/inst/.
As a result of a discussion on octave-help I send a new
patch. leasqr.m now only calls dfdp.m (or its replacement) with an
extra argument if the user specifies the new 'bounds' option. So users
who don't use the new option do not have to change their possible
dfdp-replacements while those who want to use the option are informed
by the help that they might have to change the dfdp-replacement.
Olaf
Index: leasqr.m
===================================================================
--- leasqr.m (revision 6079)
+++ leasqr.m (working copy)
@@ -42,21 +42,37 @@
% with y, x, p of the form y, x, pin as described above.
% dFdp = name of partial derivative function in quotes; default is "dfdp", a
% slow but general partial derivatives function; the function shall be
-% of the form prt=dfdp(x,f,p,dp,F) (see dfdp.m).
+% of the form prt=dfdp(x,f,p,dp,F[,bounds]). For backwards
+% compatibility, the function will only be called with an extra
+% 'bounds' argument if the 'bounds' option is explicitely specified
+% to leasqr (see dfdp.m).
% stol = scalar tolerance on fractional improvement in scalar sum of
% squares = sum((wt.*(y-f))^2); default stol = .0001;
% niter = scalar maximum number of iterations; default = 20;
-% options = matrix of n rows (same number of rows as pin) containing
-% column 1: desired fractional precision in parameter estimates.
-% Iterations are terminated if change in parameter vector (chg) on two
-% consecutive iterations is less than their corresponding elements
-% in options(:,1). [ie. all(abs(chg*current parm est) < options(:,1))
-% on two consecutive iterations.], default = zeros().
-% column 2: maximum fractional step change in parameter vector.
-% Fractional change in elements of parameter vector is constrained to be
-% at most options(:,2) between sucessive iterations.
-% [ie. abs(chg(i))=abs(min([chg(i) options(i,2)*current param estimate])).],
-% default = Inf*ones().
+% options = structure, currently recognized fields are 'fract_prec',
+% 'max_fract_change', and 'bounds'. For backwards compatibility,
+% 'options' can also be a matrix whose first and second column
+% contains the values of 'fract_prec' and 'max_fract_change',
+% respectively.
+% Field 'options.fract_prec': column vector (same length as 'pin') of
+% desired fractional precisions in parameter estimates. Iterations
+% are terminated if change in parameter vector (chg) on two
+% consecutive iterations is less than their corresponding elements in
+% 'options.fract_prec'. [ie. all(abs(chg*current parm est) <
+% options.fract_prec) on two consecutive iterations.], default =
+% zeros().
+% Field 'options.max_fract_change': column vector (same length as
+% 'pin) of maximum fractional step changes in parameter vector.
+% Fractional change in elements of parameter vector is constrained to
+% be at most 'options.max_fract_change' between sucessive iterations.
+% [ie. abs(chg(i))=abs(min([chg(i)
+% options.max_fract_change(i)*current param estimate])).], default =
+% Inf*ones().
+% Field 'options.bounds': two-column-matrix, one row for each
+% parameter in 'pin'. Each row contains a minimal and maximal value
+% for each parameter. Default: [-Inf, Inf] in each row. If this
+% field is used with an existing user-side function for 'dFdp'
+% (see above) the functions interface might have to be changed.
%
% OUTPUT VARIABLES
% f = column vector of values computed: f = F(x,p).
@@ -122,6 +138,7 @@
% MATLAB version >= 4.0) in computation of covp, corp, stdresid.
% Modified by Francesco Potortì
% for use in Octave
+% Added bounds feature to this file and to dfdp.m, Olaf Till 4-Aug-2009
%
% References:
% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
@@ -158,27 +175,66 @@
y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns
% check data vectors- same length?
-m=length(y); n=length(pin); p=pin;[m1,m2]=size(x);
+m=length(y); n=length(pin); [m1,m2]=size(x);
if m1~=m ,error('input(x)/output(y) data must have same number of rows ') ,end;
-if (nargin <= 9),
- options=[zeros(n,1), Inf*ones(n,1)];
- nor = n; noc = 2;
-else
- [nor, noc]=size(options);
- if (nor ~= n),
- error('options and parameter matrices must have same number of rows'),
- end;
- if (noc ~= 2),
- options=[options(:,1), Inf*ones(nor,1)];
- end;
-end;
-pprec=options(:,1);
-maxstep=options(:,2);
-%
+ %% processing of 'options'
+ pprec = zeros (n, 1);
+ maxstep = Inf * ones (n, 1);
+ bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
+ dfdp_cmd = "feval(dFdp,x,fbest,pprev,dp,F);"; % will possibly be redefined
+ if (nargin > 9)
+ if (ismatrix (options)) % backwards compatibility
+ tp = options;
+ options = struct ("fract_prec", tp(:, 1));
+ if (columns (tp) > 1)
+ options.max_fract_change = tp(:, 2);
+ endif
+ endif
+ if (isfield (options, "fract_prec"))
+ pprec = options.fract_prec;
+ if (rows (pprec) != n || columns (pprec) != 1)
+ error ("fractional precisions: wrong dimensions");
+ endif
+ endif
+ if (isfield (options, "max_fract_change"))
+ maxstep = options.max_fract_change;
+ if (rows (maxstep) != n || columns (maxstep) != 1)
+ error ("maximum fractional step changes: wrong dimensions");
+ endif
+ endif
+ if (isfield (options, "bounds"))
+ dfdp_cmd = "feval(dFdp,x,fbest,pprev,dp,F,bounds);";
+ bounds = options.bounds;
+ if (rows (bounds) != n || columns (bounds) != 2)
+ error ("bounds: wrong dimensions");
+ endif
+ idx = bounds(:, 1) > bounds(:, 2);
+ tp = bounds(idx, 2);
+ bounds(idx, 2) = bounds(idx, 1);
+ bounds(idx, 1) = tp;
+ if (any (idx = bounds(:, 1) == bounds(:, 2)))
+ warning ("lower and upper bounds identical for some parameters, setting the respective elements of 'dp' to zero");
+ dp(idx) = 0;
+ endif
+ if (any (idx = pin < bounds(:, 1)))
+ warning ("some initial parameters set to lower bound");
+ pin(idx) = bounds(idx, 1);
+ endif
+ if (any (idx = pin > bounds(:, 2)))
+ warning ("some initial parameters set to upper bound");
+ pin(idx) = bounds(idx, 2);
+ endif
+ endif
+ endif
+ if (all (dp == 0))
+ error ("no free parameters");
+ endif
+
% set up for iterations
%
+p = pin;
f=feval(F,x,p); fbest=f; pbest=p;
r=wt.*(y-f);
ss=r'*r;
@@ -193,7 +249,7 @@
%
for iter=1:niter,
pprev=pbest;
- prt=feval(dFdp,x,fbest,pprev,dp,F);
+ prt = eval (dfdp_cmd);
r=wt.*(y-fbest);
sprev=sbest;
sgoal=(1-stol)*sprev;
@@ -234,6 +290,20 @@
% ss=scalar sum of squares=sum((wt.*(y-f))^2).
if (any(abs(chg) > 0.1*aprec)),%--- % only worth evaluating function if
p=chg+pprev; % there is some non-miniscule change
+ %% apply bounds; preserving the direction of the attempted step
+ %% might lead to fixing _all_ parameters at their current value,
+ %% so decided not to do that and to simply reset parameters to
+ %% bounds
+ lidx = p < bounds(:, 1);
+ p(lidx) = bounds(lidx, 1);
+ uidx = p > bounds(:, 2);
+ p(uidx) = bounds(uidx, 2);
+ if (verbose && any (idx = lidx | uidx))
+ printf ("bounds apply for parameters number %s %i\n", \
+ sprintf ("%i, ", find (idx)(1:end - 1)), \
+ find (idx)(end));
+ endif
+ %%
f=feval(F,x,p);
r=wt.*(y-f);
ss=r'*r;
@@ -280,7 +350,7 @@
% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
% re-evaluate the Jacobian at optimal values
-jac=feval(dFdp,x,f,p,dp,F);
+jac = eval (dfdp_cmd);
msk = dp ~= 0;
n = sum(msk); % reduce n to equal number of estimated parameters
jac = jac(:, msk); % use only fitted parameters
Index: dfdp.m
===================================================================
--- dfdp.m (revision 6079)
+++ dfdp.m (working copy)
@@ -15,7 +15,8 @@
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.
-function prt=dfdp(x,f,p,dp,func)
+function prt=dfdp(x,f,p,dp,func,bounds)
+% function prt = dfdp (x, f, p, dp, func[, bounds])
% numerical partial derivatives (Jacobian) df/dp for use with leasqr
% --------INPUT VARIABLES---------
% x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
@@ -27,26 +28,62 @@
% dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
% func=string naming the function (.m) file
% e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
+% bounds=two-column-matrix of lower and upper bounds for parameters
+% If no 'bounds' options is specified to leasqr, it will call
+% dfdp without the 'bounds' argument.
%----------OUTPUT VARIABLES-------
% prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
%================================
m=size(x,1); if (m==1), m=size(x,2); end %# PAK: in case #cols > #rows
n=length(p); %dimensions
-ps=p; prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero
-for j=1:n
- del(j)=dp(j) .*p(j); %cal delx=fract(dp)*param value(p)
- if p(j)==0
- del(j)=dp(j); %if param=0 delx=fraction
- end
- p(j)=ps(j) + del(j);
- if del(j)~=0, f1=feval(func,x,p);
- if dp(j) < 0, prt(:,j)=(f1-f)./del(j);
- else
- p(j)=ps(j)- del(j);
- prt(:,j)=(f1-feval(func,x,p))./(2 .*del(j));
- end
- end
- p(j)=ps(j); %restore p(j)
-end
-return
+if (nargin < 6)
+ bounds = ones (n, 2);
+ bounds(:, 1) *= -Inf;
+ bounds(:, 2) *= Inf;
+endif
+prt=zeros(m,n); % initialise Jacobian to Zero
+del = dp .* p; %cal delx=fract(dp)*param value(p)
+idx = p == 0;
+del(idx) = dp(idx); %if param=0 delx=fraction
+idx = dp > 0;
+del(idx) = abs (del(idx)); % not for one-sided intervals, changed
+ % direction of intervals could change
+ % behavior of optimization without bounds
+min_del = min (abs (del), bounds(:, 2) - bounds(:, 1));
+ for j=1:n
+ ps = p;
+ if (dp(j)~=0)
+ if (dp(j) < 0)
+ ps(j) = p(j) + del(j);
+ if (ps(j) < bounds(j, 1) || ps(j) > bounds(j, 2))
+ t_del1 = max (bounds(j, 1) - p(j), - abs (del(j))); %
+ %non-positive
+ t_del2 = min (bounds(j, 2) - p(j), abs (del(j))); %
+ %non-negative
+ if (- t_del1 > t_del2)
+ del(j) = t_del1;
+ else
+ del(j) = t_del2;
+ endif
+ ps(j) = p(j) + del(j);
+ endif
+ prt(:, j) = (feval (func, x, ps) - f) / del(j);
+ else
+ if (p(j) - del(j) < bounds(j, 1))
+ tp = bounds(j, 1);
+ ps(j) = tp + min_del(j);
+ elseif (p(j) + del(j) > bounds(j, 2))
+ ps(j) = bounds(j, 2);
+ tp = ps(j) - min_del(j);
+ else
+ ps(j) = p(j) + del(j);
+ tp = p(j) - del(j);
+ min_del(j) = 2 * del(j);
+ endif
+ f1 = feval (func, x, ps);
+ ps(j) = tp;
+ prt(:, j) = (f1 - feval (func, x, ps)) / (min_del(j));
+ endif
+ endif
+ endfor
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
