X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Foptim-1.2.0%2Fbrent_line_min.m;fp=octave_packages%2Foptim-1.2.0%2Fbrent_line_min.m;h=1c881e191fae2dec8939c1db413c4e2a74bd4fa3;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/optim-1.2.0/brent_line_min.m b/octave_packages/optim-1.2.0/brent_line_min.m new file mode 100644 index 0000000..1c881e1 --- /dev/null +++ b/octave_packages/optim-1.2.0/brent_line_min.m @@ -0,0 +1,233 @@ +## Copyright (C) 2009 Levente Torok +## +## 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 . + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{s},@var{v},@var{n}]} brent_line_min ( @var{f},@var{df},@var{args},@var{ctl} ) +## Line minimization of f along df +## +## Finds minimum of f on line @math{ x0 + dx*w | a < w < b } by +## bracketing. a and b are passed through argument ctl. +## +## @subheading Arguments +## @itemize @bullet +## @item @var{f} : string : Name of function. Must return a real value +## @item @var{args} : cell : Arguments passed to f or RxC : f's only argument. x0 must be at @var{args}@{ @var{ctl}(2) @} +## @item @var{ctl} : 5 : (optional) Control variables, described below. +## @end itemize +## +## @subheading Returned values +## @itemize @bullet +## @item @var{s} : 1 : Minimum is at x0 + s*dx +## @item @var{v} : 1 : Value of f at x0 + s*dx +## @item @var{nev} : 1 : Number of function evaluations +## @end itemize +## +## @subheading Control Variables +## @itemize @bullet +## @item @var{ctl}(1) : Upper bound for error on s Default=sqrt(eps) +## @item @var{ctl}(2) : Position of minimized argument in args Default= 1 +## @item @var{ctl}(3) : Maximum number of function evaluations Default= inf +## @item @var{ctl}(4) : a Default=-inf +## @item @var{ctl}(5) : b Default= inf +## @end itemize +## +## Default values will be used if ctl is not passed or if nan values are +## given. +## @end deftypefn + +function [s,gs,nev] = brent_line_min( f,dx,args,ctl ) + +verbose = 0; + +seps = sqrt (eps); + + # Default control variables +tol = 10*eps; # sqrt (eps); +narg = 1; +maxev = inf; +a = -inf; +b = inf; + +if nargin >= 4, # Read arguments + if !isnan (ctl (1)), tol = ctl(1); end + if length (ctl)>=2 && !isnan (ctl(2)), narg = ctl(2); end + if length (ctl)>=3 && !isnan (ctl(3)), maxev = ctl(3); end + if length (ctl)>=4 && !isnan (ctl(4)), a = ctl(4); end + if length (ctl)>=5 && !isnan (ctl(5)), b = ctl(5); end + +end # Otherwise, use defaults, def'd above + +if a>b, tmp=a; a=b; b=tmp; end + +if narg > length (args), + printf ("brent_line_min : narg==%i > length (args)==%i",\ + narg, length (args)); + keyboard +end + + +if ! iscell (args), + args = {args}; +endif + +x = args{ narg }; + +[R,C] = size (x); +N = R*C; # Size of argument + +gs0 = gs = feval (f, args); +nev = 1; + # Initial value +s = 0; + +if all (dx==0), return; end # Trivial case + + # If any of the bounds is infinite, find + # finite bounds that bracket minimum +if !isfinite (a) || !isfinite (b), + if !isfinite (a) && !isfinite (b), + [a,b, ga,gb, n] = __bracket_min (f, dx, narg, args); + elseif !isfinite (a), + [a,b, ga,gb, n] = __semi_bracket (f, -dx, -b, narg, args); + tmp = a; a = -b; b = -tmp; + tmp = ga; ga = gb; gb = tmp; + else + [a,b, ga,gb, n] = __semi_bracket (f, dx, a, narg, args); + end + nev += n; +else + args{narg} = x+a*dx; ga = feval( f, args ); + args{narg} = x+b*dx; gb = feval( f, args ); + nev += 2; +end # End of finding bracket for minimum + +if a > b, # Check assumptions + printf ("brent_line_min : a > b\n"); + keyboard +end + +s = 0.5*(a+b); +args{narg} = x+ s*dx; gs = feval( f, args ); +nev++; + +if verbose, + printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\ + a,s,b,ga,gs,gb); +end + +maxerr = 2*tol; + +while ( b-a > maxerr ) && nev < maxev, + + if gs > ga && gs > gb, # Check assumptions + printf ("brent_line_min : gs > ga && gs > gb\n"); + keyboard + end + + if ga == gb && gb == gs, break end + + # Don't trust poly_2_ex for glued points + # (see test_poly_2_ex). + + ## if min (b-s, s-a) > 10*seps, + + # If s is not glued to a or b and does not + # look linear + ## mydet = sum (l([2 3 1]).*f([3 1 2])-l([3 1 2]).*f([2 3 1])) + mydet = sum ([s b a].*[gb ga gs] - [b a s].*[gs gb ga]); + if min (b-s, s-a) > 10*seps && abs (mydet) > 10*seps && \ + (t = poly_2_ex ([a,s,b], [ga, gs, gb])) < b && t > a, + + # t has already been set + + ## if t>=b || t<=a, + ## printf ("brent_line_min : t is not in ]a,b[\n"); + ## keyboard + ## end + + # Otherwise, reduce the biggest of the + # intervals, but not too much + elseif s-a > b-s, + t = max (0.5*(a+s), s-100*seps); + else + t = min (0.5*(s+b), s+100*seps); + end + + if abs (t-s) < 0.51*maxerr, + #sayif (verbose, "ungluing t and s\n"); + t = s + (1 - 2*(s-a > b-s))*0.49*maxerr ; + end + + if a > s || s > b, # Check assumptions + printf ("brent_line_min : a > s || s > b\n"); + keyboard + end + + xt = args; + args{narg} = x+t*dx; + gt = feval( f, args ); + nev++; + + if verbose, + printf ("t = %.3e, gt = %.3e\n",t,gt); + end + + if t ga + seps, # Check assumptions + printf ("brent_line_min : gt > ga\n"); + keyboard + end + + if gt < gs, + b = s; gb = gs; + s = t; gs = gt; + else + a = t; ga = gt; + end + else # New point is in ]s,b[ + if gt > gb + seps, # Check assumptions + printf ("brent_line_min : gt > gb\n"); + keyboard + end + + if gt < gs, + a = s; ga = gs; + s = t; gs = gt; + else + b = t; gb = gt; + end + end + + if verbose, + printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\ + a,s,b,ga,gs,gb); + end + ## keyboard + ## [b-a, maxerr] +end + +s2 = 0.5*(a+b); +args{narg} = x + s2*dx; gs2 = feval (f, args ); +nev++; + +if gs2 < gs, + s = s2; gs = gs2; +end + +if gs > gs0, + printf ("brent_line_min : goes uphill by %8.3e\n",gs-gs0); + keyboard +end