X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=sidebyside;f=octave_packages%2Foptim-1.2.0%2Fd2_min.m;fp=octave_packages%2Foptim-1.2.0%2Fd2_min.m;h=befda03588c39db45691735d8322f8d5b5ded7ce;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/optim-1.2.0/d2_min.m b/octave_packages/optim-1.2.0/d2_min.m new file mode 100644 index 0000000..befda03 --- /dev/null +++ b/octave_packages/optim-1.2.0/d2_min.m @@ -0,0 +1,393 @@ +## Copyright (C) 2002 Etienne Grossmann +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + +## [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization +## +## Minimize f(x) using 1st and 2nd derivatives. Any function w/ second +## derivatives can be minimized, as in Newton. f(x) decreases at each +## iteration, as in Levenberg-Marquardt. This function is inspired from the +## Levenberg-Marquardt algorithm found in the book "Numerical Recipes". +## +## ARGUMENTS : +## f : string : Cost function's name +## +## d2f : string : Name of function returning the cost (1x1), its +## differential (1xN) and its second differential or it's +## pseudo-inverse (NxN) (see ctl(5) below) : +## +## [v,dv,d2v] = d2f (x). +## +## args : list : f and d2f's arguments. By default, minimize the 1st +## or matrix : argument. +## +## ctl : vector : Control arguments (see below) +## or struct +## +## code : string : code will be evaluated after each outer loop that +## produced some (any) improvement. Variables visible from +## "code" include "x", the best parameter found, "v" the +## best value and "args", the list of all arguments. All can +## be modified. This option can be used to re-parameterize +## the argument space during optimization +## +## CONTROL VARIABLE ctl : (optional). May be a struct or a vector of length +## ---------------------- 5 or less where NaNs are ignored. Default values +## are written . +## FIELD VECTOR +## NAME POS +## +## ftol, f N/A : Stop search when value doesn't improve, as tested by +## +## f > Deltaf/max(|f(x)|,1) +## +## where Deltaf is the decrease in f observed in the last +## iteration. <10*sqrt(eps)> +## +## utol, u N/A : Stop search when updates are small, as tested by +## +## u > max { dx(i)/max(|x(i)|,1) | i in 1..N } +## +## where dx is the change in the x that occured in the last +## iteration. +## +## dtol, d N/A : Stop search when derivative is small, as tested by +## +## d > norm (dv) +## +## crit, c ctl(1) : Set one stopping criterion, 'ftol' (c=1), 'utol' (c=2) +## or 'dtol' (c=3) to the value of by the 'tol' option. <1> +## +## tol, t ctl(2) : Threshold in termination test chosen by 'crit' <10*eps> +## +## narg, n ctl(3) : Position of the minimized argument in args <1> +## maxev,m ctl(4) : Maximum number of function evaluations +## maxout,m : Maximum number of outer loops +## id2f, i ctl(5) : 0 if d2f returns the 2nd derivatives, 1 if <0> +## it returns its pseudo-inverse. +## +## verbose, v N/A : Be more or less verbose (quiet=0) <0> + +function [xbest,vbest,nev,hbest,args] = d2_min (f,d2f,args,ctl,code) + +maxout = inf; +maxinner = 30 ; + +tcoeff = 0.5 ; # Discount on total weight +ncoeff = 0.5 ; # Discount on weight of newton +ocoeff = 1.5 ; # Factor for outwards searching + +report = 0 ; # Never report +verbose = 0 ; # Be quiet +prudent = 1 ; # Check coherence of d2f and f? + +niter = 0 ; + +crit = 0; # Default control variables +ftol = 10 * sqrt (eps); +dtol = eps; +utol = tol = nan; +narg = 1; +maxev = inf; +id2f = 0; + +if nargin >= 4 # Read arguments + if isnumeric (ctl) + if length (ctl)>=1 && !isnan (ctl(1)), crit = ctl(1); end + if length (ctl)>=2 && !isnan (ctl(2)), tol = ctl(2); end + if length (ctl)>=3 && !isnan (ctl(3)), narg = ctl(3); end + if length (ctl)>=4 && !isnan (ctl(4)), maxev = ctl(4); end + if length (ctl)>=5 && !isnan (ctl(5)), id2f = ctl(5); end + elseif isstruct (ctl) + if isfield (ctl, "crit") , crit = ctl.crit ; end + if isfield (ctl, "tol") , tol = ctl.tol ; end + if isfield (ctl, "narg") , narg = ctl.narg ; end + if isfield (ctl, "maxev") , maxev = ctl.maxev ; end + if isfield (ctl, "maxout") , maxout = ctl.maxout ; end + if isfield (ctl, "id2f") , id2f = ctl.id2f ; end + if isfield (ctl, "verbose"), verbose = ctl.verbose; end + if isfield (ctl, "code") , code = ctl.code ; end + else + error ("The 'ctl' argument should be either a vector or a struct"); + end +end + +if crit == 1, ftol = tol; +elseif crit == 2, utol = tol; +elseif crit == 3, dtol = tol; +elseif crit, error ("crit is %i. Should be 1,2 or 3.\n"); +end + + +if nargin < 5, code = "" ; end + +if iscell (args) # List of arguments + x = args{narg}; +else # Single argument + x = args; + args = {args}; +end + +############################## Checking ############################## +if narg > length (args) + error ("d2_min : narg==%i, length (args)==%i\n", + narg, length (args)); +end + +if tol <= 0 + printf ("d2_min : tol=%8.3g <= 0\n",tol) ; +end + +if !ischar (d2f) || !ischar (f) + printf ("d2_min : f and d2f should be strings!\n"); +end + +sz = size (x); N = prod (sz); + +v = feval (f, args{:}); +nev = [1,0]; + +if prudent && (! isnumeric (v) || isnan (v) || any (size (v)>1)) + error ("Function '%s' returns inadequate output", f); +end + +xbest = x = x(:); +vold = vbest = nan ; # Values of f +hbest = nan ; # Inv. Hessian + +if verbose + printf ( "d2_min : Initially, v=%8.3g\n",v); +end + +while niter <= maxout + niter += 1; + if nev(1) < maxev, break; end; + + [v,d,h] = feval (d2f, args{1:narg-1},reshape(x,sz),args{narg+1:end}); + nev(2)++; + + if prudent && niter <= 1 && \ + (! isnumeric (v) || isnan (v) || any (size (v)>1) || \ + ! isnumeric (d) || length (d(:)) != N || \ + ! isnumeric (h) || any (size (h) != N)) + error ("Function '%s' returns inadequate output", d2f); + end + + if ! id2f, h = pinv (h); end + d = d(:); + + if prudent + v2 = feval (f, args{1:narg-1},reshape(x,sz),args{narg+1:end}); + nev(1)++; + if abs(v2-v) > 0.001 * sqrt(eps) * max (abs(v2), 1) + printf ("d2_min : f and d2f disagree %8.3g\n",abs(v2-v)); + end + end + + xbest = x ; + if ! isnan (vbest) # Check that v ==vbest + if abs (vbest - v) > 1000*eps * max (vbest, 1) + printf ("d2_min : vbest changed at beginning of outer loop\n"); + end + end + vold = vbest = v ; + hbest = h ; + + if length (code), abest = args; end # Eventually stash all args + + if verbose || (report && rem(niter,max(report,1)) == 1) + printf ("d2_min : niter=%d, v=%8.3g\n",niter,v ); + end + + if norm (d) < dtol # Check for small derivative + if verbose || report + printf ("d2_min : Exiting because of low gradient\n"); + end + break; # Exit outer loop + end + + dnewton = -h*d ; # Newton step + # Heuristic for negative hessian + if dnewton'*d > 0, dnewton = -100*d; end + wn = 1 ; # Weight of Newton step + wt = 1 ; # Total weight + + ninner = 0; + done_inner = 0; # 0=not found. 1=Ready to quit inner. + + # ########################################## + while ninner < maxinner, # Inner loop ############################### + ninner += 1; + # Proposed step + dx = wt*(wn*dnewton - (1-wn)*d) ; + xnew = x+dx ; + + if verbose + printf (["Weight : total=%8.3g, newtons's=%8.3g vbest=%8.3g ",... + "Norm:Newton=%8.3g, deriv=%8.3g\n"],... + wt,wn,vbest,norm(wt*wn*dnewton),norm(wt*(1-wn)*d)); + end + if any(isnan(xnew)) + printf ("d2_min : Whoa!! any(isnan(xnew)) (1)\n"); + end + + vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:end}); + nev(1)++; + + if vnew= maxinner # There was a problem + if verbose + printf ( "d2_min : Too many inner loops (vnew=%8.3g)\n",vnew); + end + + # ########################################## + else # Look for improvement along dbest + wn = ocoeff ; + xnew = x+wn*dbest; + if any(isnan(xnew)), + printf ("d2_min : Whoa!! any(isnan(xnew)) (2)\n"); + end + vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:end}); + nev(1)++; + + while vnew < vbest, + vbest = vnew; # Stash best values + wbest = wn; + xbest = xnew; + wn = wn*ocoeff ; + xnew = x+wn*dbest; + vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:length(args)}); + if verbose + printf ( "Looking farther : v = %8.3g\n",vnew); + end + nev(1)++; + end + end # End of improving along dbest + # ########################################## + + if verbose || rem(niter,max(report,1)) == 1 + if vold, + if verbose + printf ("d2_min : Inner loop : vbest=%8.5g, vbest/vold=%8.5g\n",\ + vbest,vbest/vold); + end + else + if verbose + printf ( "d2_min : Inner loop : vbest=%8.5g, vold=0\n", vbest); + end + end + end + + if vbest < vold + ## "improvement found" + if prudent + tmpv = feval (f, args{1:narg-1},reshape(xbest,sz),args{2:end}); + nev(1)++; + + if abs (tmpv-vbest) > eps + printf ("d2_min : Whoa! Value at xbest changed by %g\n",\ + abs(tmpv-vbest)); + end + end + v = vbest; x = xbest; + if ! isempty (code) + if verbose + printf ("d2_min : Going to eval (\"%s\")\n",code); + end + + xstash = xbest; + astash = abest; + args = abest; # Here : added 2001/11/07. Is that right? + x = xbest; + eval (code, "printf (\"code fails\\n\");"); + xbest = x; + abest = args; + # Check whether eval (code) changes value + if prudent + tmpv = feval (f, args{1:narg-1},reshape(x,sz),args{2:end}); + nev(1)++; + if abs (tmpv-vbest) > max (min (100*eps,0.00001*abs(vbest)), eps) , + printf ("d2_min : Whoa! Value changes by %g after eval (code)\n",\ + abs (tmpv-vbest)); + end + end + end + end + + if ! isnan (ftol) && ftol > (vold-vbest)/max(vold,1), + if verbose || report + printf ("d2_min : Quitting, niter=%-3d v=%8.3g, ",niter,v); + if vold, printf ("v/vold=%8.3g \n",v/vold); + else printf ("vold =0 \n",v); + end + end + break ; # out of outer loop + end + if ! isnan (utol) && utol > max (abs (wbest*dbest)) / max(abs (xbest),1) + if verbose || report + printf ("d2_min : Quitting, niter=%-3d v=%8.3g, ",niter,v); + if vold, printf ("v/vold=%8.3g \n",v/vold); + else printf ("vold =0 \n",v); + end + end + break ; # out of outer loop + end +end # End of outer loop ################## + +xbest = reshape (xbest, sz); +if length (code) + args = abest; + args(narg) = xbest; +end + +if niter > maxout + if verbose + printf ( "d2_min : Outer loop lasts forever\n"); + end +end + + # One last check +if prudent + err = feval (f, args{1:narg-1},reshape(xbest,sz),args{2:end}); + nev(1)++; + + if abs (err-vbest) > eps, + printf ("d2_min : Whoa!! xbest does not eval to vbest\n"); + printf (" : %8.3e - %8.3e = %8.3e != 0\n",err,vbest,err-vbest); + end +end +