brent0

(tmp_for_tar/deriv_min.doc/brent0.m)

Function Synopsis

[xmin,fx,niter] = brent0(dp,func,p,a,c,b)

Help text

        [x0,v,niter] = brent0(dx,func,x,[a,b,c]|[bsize])

 Returned values : 

 x0    : MxN  the found minimum of 'func'.
 v     : 1x1  the value of func( x0, all_va_args )
 niter : 1x1  the number of function evaluations. 

 Arguments :

 a,b,c,dx : Search will take place on the line segment 

         { args(nmin) + dx*w | min(a,b,c) <= w <= max(a,b,c) }

 func  : the name of the function to minimize

 args  : the arguments of func. The function will be minimized over
         the nmin'th argument. 

 Variables in the function :

 'ftol', ftol  : stop minimizing when abs(a-b) < ftol. This is the
                 default stopping criterion, with ftol==1e-6. Don't
                 take ftol<abs(eps).

 'maxeval',e   : The maximum number of function evaluations. A
                 negative value means 'no limit'. Default is 50000. 

 'bsize',s     : Find the bracketing range, starting from a 2*s wide
                 interval. 

 'report',n    : Show a short report every (approximately) n
                 evaluations. 0 (default) means never.

 'verbose'     : Produce a little message at end

 'debug',debug : Pass a struct 'debug'. At some points in the
                 function, some internal variables are checked
                 against the corresponding key of debug. If they
                 mismatch, a message is printed and a call to
                 keyboard is made. Recognized keys are :
                   M,N,x,ftol,func
                 See the code for more info.

 Last modified: October 2000

Listing of function file tmp_for_tar/deriv_min.doc/brent0.m

##        [x0,v,niter] = brent0(dx,func,x,[a,b,c]|[bsize])
## 
## Returned values : 
##
## x0    : MxN  the found minimum of 'func'.
## v     : 1x1  the value of func( x0, all_va_args )
## niter : 1x1  the number of function evaluations. 
##
## Arguments :
## 
## a,b,c,dx : Search will take place on the line segment 
##
##         { args(nmin) + dx*w | min(a,b,c) <= w <= max(a,b,c) }
##
## func  : the name of the function to minimize
## 
## args  : the arguments of func. The function will be minimized over
##         the nmin'th argument. 
## 
##
##
## Variables in the function :
## 
## 'ftol', ftol  : stop minimizing when abs(a-b) < ftol. This is the
##                 default stopping criterion, with ftol==1e-6. Don't
##                 take ftol<abs(eps).
## 
## 'maxeval',e   : The maximum number of function evaluations. A
##                 negative value means 'no limit'. Default is 50000. 
## 
## 'bsize',s     : Find the bracketing range, starting from a 2*s wide
##                 interval. 
##  
## 'report',n    : Show a short report every (approximately) n
##                 evaluations. 0 (default) means never.
## 
## 'verbose'     : Produce a little message at end
## 
## 'debug',debug : Pass a struct 'debug'. At some points in the
##                 function, some internal variables are checked
##                 against the corresponding key of debug. If they
##                 mismatch, a message is printed and a call to
##                 keyboard is made. Recognized keys are :
##                   M,N,x,ftol,func
##                 See the code for more info.
## 


## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>
## Last modified: October 2000

function [xmin,fx,niter] = brent0(dp,func,p,a,c,b)
# Set the defaults (not all) #########################################
ftol = sqrt(eps) ;			# maximum distance
maxeval = 5000 ;		# maximum number of evals
report = 0 ;			# Never report
verbose = 0 ;			# Not even at end
				# 
debug.dummy = 0 ;		# Dummy structure element
niter = 0 ;
bsize = 1 ;
tiny = 1e-10 ;		# too small a step
cgold = 0.381966 ;




if ! isstr( func ) ,
  printf("Argument 'func' is not a string\n");
  keyboard ;
end

				# Check for zero dp
if all(dp==0),
  xmin = 0 ;
  niter = 0;
  fx = feval(func,p) ;
end

if exist("b") != 1,
  ## printf("brent0 : finding bracket\n") ;
  if exist("a") != 1, bsize = 1 ;
  else                bsize = a ;
  end

  [a,c,b] = bracket_minimum0(func,bsize,p,dp) ;
  ## printf("brent0 : found bracket %8.3g %8.3g, %8.3g\n",a,b,c) ;
end



# ####################################################################

# ####################################################################
# Okay, after all that administrative stuff, we can start the        #
# optimization algorithm, at last.                                   # 
# ####################################################################

ee = 0 ;

nextprint = report ;

ax = a ; bx = b ; cx = c ;
if ax < cx , a=ax ; else a=cx; end
if ax > cx , b=ax ; else b=cx; end 


x = v = w = bx ;
fx = fv = fw = feval(func,p+x*dp) ;

while niter++ <= maxeval || maxeval < 0, 

  xm = 0.5*(a+b) ;
  
  t1 = ftol*abs(x)+tiny ;
  t2 = 2*t1 ;
  ## printf("xm=%-5d t1=%-5d t2=%-5d\n",xm,t1,t2);
  
  if abs( x-xm )  <= t2-0.5*(b-a) ,
    xmin = x ;

    if verbose ,
      printf("brent0:niter=%-5d   xm=%-8.3g fx=%-8.3g\n>>> Return\n",...
	     niter,xm,fx) ;
    end
    return 
  end

  if verbose ,
    printf("niter=%-5d   x=%-8.3g xm=%-8.3g fx=%-8.3g\n",...
	   niter,x,xm,fx) ;
  end
    

  if abs( ee ) > t1 ,
    rr = (x-w)*(fx-fv) ;
    qq = (x-v)*(fx-fw) ;
    pp = (x-v)*qq-(x-w)*rr ;
    qq = 2*(qq-rr) ;
    if qq>0, pp = -pp ; end
    qq = abs(qq) ;
    etmp = ee ;
    ee = dd ;
    if abs(pp) >= abs( 0.5*qq*etmp ) || pp < qq*(a-x) || pp>=qq*(b-x),
      if x >= xm, ee = a-x ; else ee = b-x ; end 
      dd = cgold * ee ;
    else
      dd = pp/qq ;
      u = x+dd ;
      if u-a < t2 || b-u<t2,
	## ((b) > 0.0 ? fabs(a) : -fabs(a))
	if xm-x>0; dd = t1 ; else dd = -t1 ; end
      end
    end
  else
    if x >= xm, ee = a-x ; else ee = b-x ; end 
    dd = cgold * ee ;
  end
  if abs(dd) > t1, 
    u = x+dd ; 
  else 
    if dd>=0 , u = x+t1; else u=x-t1 ; end
  end
  fu = feval(func,p+u*dp) ;
  
  if fu <= fx ,
    if u >= x , a=x ; else b=x ; end
    v = w ; w = x ; x = u ;
    fv = fw ; fw = fx ; fx = fu ;
  else
    if u < x , a=u; else b=u ; end
    if fu <= fw || w == x ,
      v = w; w=u; fv=fw; fw=fu ;
    elseif fu <= fv || v==x || v==w ,
      v=u ;
      fv=fu ;
    end
  end
end

if verbose
  printf("brent0:Too many iterations") ;
end
## keyboard