function [xmin, fmin] = dbrent(fun,dfun, ax,bx,cx,tol)
%
% Brentova metoda s derivacemi pro hledani minima funkce.
% Podle rutiny DBRENT z Numerical Recipes in C.
%    
% Hleda Brentovou metodou minimum zadane funkce FUNC, ktere je jiz
%   dopredu ohraniceno hodnotami AX,BX,CX, kde BX lezi mezi AX a CX 
%   a funkce je v BX minimalni z dane trojice. 
% TOL zadava relativni presnost hledani minima,
% Vystup:
%   XMIN je nalezena poloha minima,
%   FMIN je hodnota funkce v tomto bodu.
%
% Pokud je na vstupu FUN inline funkce, predpokladame skalarni funkci     
% Pokud je na vstupu FUN = {FUNCND,PT0,DIR)}, provadime 1D minimalizaci 
%       na primce  PT0 + x*DIR
%

    %% prevadi fun na skalarni funkci
    function f = func(fun,x)
    % 
    % Pokud jsme v 1D, funguje klasicky: 
    %       FUNC(x) = FUN(x)
    % Pokud jsme ve vice dim, prevadi problem na 1D minimalizaci na primce:
    %       FUNC(x) = FUNCND( PT0 + x*DIR )
    %    
        if (iscell(fun))
            % funkce vektoru, pocatecni bod a smer
            funcnd = fun{1};
            pt0    = fun{2};
            dir    = fun{3};
            f = funcnd(pt0+dir*x);
        else
            % napriklad: fun=inline('sin(x)-0.5');
            f = fun(x);
        end;
        return;
    end

    %% prevadi dfun na skalar:
    function df = dfunc1d(fun,dfun,x)
    % 
    % Pokud jsme v 1D, funguje klasicky: 
    %       DFUNC1D(x) = DFUN(x)
    % Pokud jsme ve vice dim, vezme projekci gradientu do smeru:
    %       DFUNC1D(x) = scalprod( DIR , DFUN( PT0+x*DIR ) )
    %    
        if (iscell(fun))
            % funkce vektoru, pocatecni bod a smer
            funcnd = fun{1};
            pt0    = fun{2};
            dir    = fun{3};
            df = dot( dir , dfun(pt0+dir*x) );
       else
            % napriklad: dfun=inline('cos(x)');
            df = dfun(x);
        end;
        return;
    end
        

    %% Definice konstant
    ITMAX  = 100;
    ZEPS   = 1.0e-10;

    
    if (ax < cx) 
        a = ax;
    else 
        a = cx;
    end;
    if (ax > cx)
        b = ax;
    else 
        b = cx;
    end;
    v = bx;   w = bx;  x = bx;
    e = 0.0;  d = 0.0;
    fx = func(fun, x);           fv = fx;  fw = fx;
    dx = dfunc1d(fun,dfun, x);   dv = dx;  dw = dx;
    for iter=1:ITMAX
        xm   = 0.5*(a+b);
        tol1 = tol*abs(x) + ZEPS;
        tol2 = 2.0*tol1;
        if ( abs(x-xm) <= (tol2-0.5*(b-a)) ) 
            xmin = x;
            fmin = fx;
            return;
        end;
        if (abs(e) > tol1) 
            d1 = 2.0*(b-a);
            d2 = d1;
            if (dw ~= dx)
                d1 = (w-x)*dx/(dx-dw);
            end;
            if (dv ~= dx)
                d2 = (v-x)*dx/(dx-dv);
            end;
            u1 = x+d1;
            u2 = x+d2;
            ok1 = ( ((a-u1)*(u1-b)>0.0) && (dx*d1<=0.0) );
            ok2 = ( ((a-u2)*(u2-b)>0.0) && (dx*d2<=0.0) );
            olde = e;
            e    = d;
            if (ok1 || ok2) 
		if (ok1 && ok2)
                    if (abs(d1)<abs(d2))
                        d = d1;
                    else
                        d = d2;
                    end;
		elseif (ok1)
                    d = d1;
                else
                    d = d2;
                end;
		if (abs(d) <= abs(0.5*olde)) 
                    u = x+d;
                    if ( (u-a<tol2)||(b-u<tol2) )
                        d = abs(tol1)*sign(xm-x);
                    end;                    
                else 
                    if (dx>=0.0)
                        e = a-x;
                    else
                        e = b-x;
                    end;
                    d = 0.5*e;
                end;
            else 
                if (dx >= 0.0)
                    e = a-x;
                else
                    e = b-x;
                end;
                d = 0.5*e;
            end
        else 
            if (dx >= 0.0)
                e = a-x;
            else
                e = b-x;
            end;
            d = 0.5*e;
        end;
        if (abs(d)>=tol1)
            u  = x + d;
            fu = func(fun, u);
        else
            u = x + abs(tol1)*sign(d);
            fu = func(fun, u);
            if (fu > fx) 
                xmin = x;
                fmin = fx;
                return;
            end;
        end;
        du = dfunc1d(fun,dfun, u);
        if (fu <= fx) 
            if (u >= x) 
                a=x; 
            else 
                b=x;
            end;
            v = w;   fv = fw;   dv = dw;
            w = x;   fw = fx;   dw = dx;
            x = u;   fx = fu;   dx = du;
        else 
            if (u < x) 
                a=u; 
            else 
                b=u;
            end;
            if ( (fu<=fw) || (w==x) ) 
                v = w;   fv = fw;   dv = dw;
                w = u;   fw = fu;   dw = du;
            elseif ( (fu<fv) || (v==x) || (v==w) )
                v = u;   fv = fu;   dv = du;
            end;
        end;
    end
    disp(['Too many iterations in DBRENT']);
    xmin = x;
    fmin = fx;
    return;

end


