function [x,y]=NajbrziSpust(f,a1,b1,a2,b2,x0,y0,epsilon)
% NAJBRZISPUST - Funkcija graficki prikazuje nalazenje minimuma funkcionala 
% metodom najbrzeg spusta
%
% Argumenti funkcije su :
%   f - string koji sadrzi opis funkcionala ciji minimum trazimo 
%   a1,b1,a2,b2 - granice intervala koji se graficki prikazuje,
%                 pozeljno je da se u ovom intervalu nalazi lokalni minimum
%   x0,y0 - pocetna iteracija
%   epsilon - zeljena tacnost.  Ukoliko se ne navede, podrazumeva se eps
%
% Funkcija vraca:
%   x,y - koordinate minimuma funkcionala

% Autor : Filip Maric
% Verzija : 1.0
% Datum : 30.3.2003.

if (nargin<8)
    epsilon=eps;
end

% Provera korektnosti argumenata
if nargin<7
    error(['Upotreba : ' mfilename '(f,a1,b1,a2,b2,x0,y0)']);
end

if (a1>=b1) | max(size(a1))>1 | max(size(b1))>1 | ~isnumeric(a1) | ~isnumeric(b1)
    error('Interval [a1,b1] je nesipravno zadat');
end

if (a2>=b2) | max(size(a2))>1 | max(size(b2))>1 | ~isnumeric(a2) | ~isnumeric(b2)
    error('Interval [a2,b2] je nesipravno zadat');
end

if (max(size(x0))>1) | ~isnumeric(x0) | (max(size(y0))>1) | ~isnumeric(y0)
    error('Pocetna iteracija nije ispravno zadana');
end

if ~isstr(f)
    error('Prvi argument mora da bude string koji sadrzi funkciju nezavisne promenjive x')'
end


%Gradimo simbolicku funkciju na osnovu stringa f
f=sym(f);

% Matrice xg i yg ce da opisuju ekvidistantnu mrezu cvorova u kojima
% cemo da crtamo vrednost funkcionala
[xg,yg]=meshgrid(linspace(a1,b1,20),linspace(a2,b2,20));
% Izracunavamo vrednost funkcionala u tackama mreze
zg=subs(subs(f,{'x','y'},{xg,yg}));
% Crtamo funkcional
mesh(xg,yg,zg);

alpha 0.5
hold on

% Izracunavamo parcijalne izvode funkcionala prvog reda
fx=diff(f,'x');
fy=diff(f,'y');
% Izracunavamo parcijalne izvode funkcionala drugog reda
fxx=diff(fx,'x');
fxy=diff(fx,'y');
fyx=diff(fy,'x');
fyy=diff(fy,'y');

% Izracunavamo vrednost funkcionala u pocetnoj tacki
z=subs(subs(f,{'x','y'},{x0,y0}));
% Iscrtavamo pocetnu tacku
plot3(x0,y0,z,'o');

% Postavljamo vrednost tekuce tacke na pocetnu
x=x0;
y=y0;

% Ispisujemo naslov grafika
title(['x=' num2str(x) ', y=' num2str(y)]);
% Cekamo da se pritisne taster
pause


while(1)
    % Izracunavamo vektor gradijenta u tekucoj tacki
    d=double(-[subs(fx,{'x','y'},{x,y}) subs(fy,{'x','y'},{x,y})]);

    %Iscrtavamo vektor gradijenta
    CrtajVektor([x y z],[d -1],'k');

    % Izracunavamo vrednost Hessiana u tekucoj tacki
    H=subs([subs(fxx,{'x','y'},{x,y}) subs(fxy,{'x','y'},{x,y});...
            subs(fyx,{'x','y'},{x,y}) subs(fyy,{'x','y'},{x,y})]);

    % Ukoliko tekuca tacka nije stacionarna izracunavamo koeficijenat lambda
    % koji odredjuje "koliko dugo" je potrebno spustati se u pravcu d 
    if (~isequal(d,[0 0]))
        l=d*d'/(d*H*d');
    else
        % Ukoliko je tacka stacionarna, nasli smo lokalni minimum
        break;
    end
    
    % Izracuanvamo vrednost nove tacke
    x1=x+l*d(1);
    y1=y+l*d(2);
    
    % Izracunavamo vrednost funkcionala u toj novoj tacki
    z=subs(subs(f,{'x','y'},{x1,y1}));
    % Iscrtavamo novu tacku na funkcionalu
    plot3(x1,y1,z,'o');
  
    % Krivom po funkcionalu spajamo staru i novu tacku
    xl=linspace(x,x1,20);
    yl=linspace(y,y1,20);
    zl=subs(subs(f,{'x','y'},{xl,yl}));
    plot3(xl,yl,zl,'r');
    
    % Ispisujemo naslov grafika
    title(['x=' num2str(x1) ', y=' num2str(y1)]);

    % Ukoliko je postignuta zeljena tacnost prekidamo iterativni postupak
    if (abs(x1-x)<epsilon & abs(y1-y)<epsilon)
        break;    
    end
    
    % Postavljamo tekucu tacku na upravno nadjenu novu
    x=x1;
    y=y1;    
    
    % Cekamo da korisnik pritisne taster
    pause
end

function CrtajVektor(A,B,boja)
nb=norm(B);
if nb~=0
    B=B/nb;
end
plot3([A(1) A(1)+B(1)],[A(2) A(2)+B(2)],[A(3) A(3)+B(3)],boja);