[x,y,fcount]=globalmin(f,a,b,x0,m,e,t,varargin)

function [x,y,fcount]=globalmin(f,a,b,x0,m,e,t,varargin)
% globalmin finds the global minimum of f(x) with a<=x<=b in R
%
%  calling sequence :
%  [x,y,fcount]=globalmin(f,a,b,x0,m,e,t,varargin)
% 
% Input
% f               is the function name (handle) to be minimized
%                 y=f(x,...)
% a,b           bounds constraints
% x0            an initial guess for x with a<=x0<=b
% m             a value such that f''(x)<m,  and a<=x<=b
% e              precision computation of y=f(x,....)
% t               absolute tolerance on the minimum value of f(x,...)
% varargin  list of supplementary parameters for f(x,....)=f(x,varargin)
%                as usually within matlab
%
% Output
% x          global optimum position
% y          f(x) minimum value
% fcount  f(x) calls number
%
% Remark1 :  smaller is the difference between max(f'') and m
%                   smaller is fcount see Brent's book
%
% Remark2 :  if the user can set e such that |fl[f(x*(1+-eps))]-f(x)|<=e
%                   then it is guaranteed that the global minimum y
%                   satisfies [min(f)-e] <= y <= [min(f)+e] + t
%                  for more details see Brent's book
%
% Examples :
%           [x,y,fcount]=globalmin(@f,a,b,x0,m,e,t);  
%                              with some function y=f(x)
%
%           [x,y,fcount]=globalmin(@f,a,b,x0,m,e,t,P1,P2,...);  
%                              with some function y=f(x,P1,P2,...)
%
%           [x,y,fcount]=globalmin(@(x) f(x,P1,P2,...),a,b,x0,m,e,t);  
%                              with some function y=f(x,P1,P2,...)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Author A. BARRAUD from R. BRENT 
%                                  "Algorithms for minimization without
%                                   derivatives" Dover Publication 2002.
%  June 2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%initialization
a0=b;x=a0;a2=a;fcount=2;
y0=f(b,varargin{:});yb=y0;
y2=f(a,varargin{:});y=y2;
if y0>=y,x=a;else y=y0;end
if m<=0 || a>=b,return;end%trivial case
m2=0.5*(1+16*eps)*m;
sc=x0;
if sc<=a || sc>=b,sc=0.5*(a+b);end
y1=f(sc,varargin{:});k=3;d0=a2-sc;h=0.81818181818;fcount=fcount+1;
if y1<y,x=sc;y=y1;end
while 1% main loop
    d1=a2-a0;d2=sc-a0;z2=b-a2;z0=y2-y1;z1=y2-y0;
    r=d1*d1*z0-d0*d0*z1;p=r;
    qs=2*(d0*z1-d1*z0);q=qs;
    loop40=1;test=k>100000 & y<y2;
    while loop40
        if ~test
            %try to find a lower value of f using quadratic interpolation
            if q*(r*yb-y2)+z2*q*((y2-y)+t) < z2*m2*r*(z2*q-r)
                a3=a2+r/q;y3=f(a3,varargin{:});fcount=fcount+1;
                if y3<y,x=a3;y=y3;end
            end
        end
        %with a probability about 0.1 do a random search 
        %for a lower value of f
        k=fix(mod(1611*k,1048576));q=1;r=k*(b-a)*1e-5;test=loop40==0;
        if r<z2,loop40=1;else loop40=0;end
    end
    %prepare to step as far as possible
    r=m2*d0*d1*d2;s=sqrt(((y2-y)+t)/m2);
    h=0.5*(1+h);p=h*(p+2*r*s);q=r+0.5*qs;
    r=-0.5*(d0+(z0+2.01*e)/(d0*m2));
    if r>=s && d0>=0,r=a2+r;else r=a2+s;end
    %it is safe to step to r, but we may try to step further
    if p*q<=0,a3=r;else a3=a2+p/q;end
    loop90=1;
    while loop90
        if a3<r,a3=r;end
        if a3<b,y3=f(a3,varargin{:});fcount=fcount+1;else a3=b;y3=yb;end
        if y3<y,x=a3;y=y3;end
        d0=a3-a2;
        if a3>r
            %inspect the parabolic lower bound on f in [a2  a3]
            p=2*(y2-y3)/(m*d0);
            if~(abs(p)>=(1+9*eps)*d0)
                if~(0.5*m2*(d0*d0+p*p)<=(y2-y)+(y3-y)+2*t)
                    %halve the step and try again
                    a3=0.5*(a2+a3);h=0.9*h;loop90=1;
                else
                    loop90=0;
                end
            else
                loop90=0;
            end
        else
            loop90=0;
        end
    end
    if a3<b
        %update for the next step
        a0=sc;sc=a2;a2=a3;y0=y1;y1=y2;y2=y3;
    else
        break;end
end% main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%