--- /dev/null
+## Copyright (C) 2009 Levente Torok <TorokLev@gmail.com>
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## -*- 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<s, # New point is in ]a,s[
+
+ if gt > 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
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff