% /home/apiola/opetus/peruskurssi/v2-3/matlab/CV
% 31-May-2000
% Script File: ShowSDN
% Minimoidaan yhdistelmällä SD + Newton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Yleistä Matlab-skriptillä operoinnista:
%   - Skripti voidaan ymmärtää "pääohjelmaksi"
%   - Kun skriptin nimi kirjoitetaan Matlab-istunnossa, siinä olevat komennot
%     suorittuvat.
%   - Skriptin testaamiseksi voidaan lisäillä kommentteja. 
%     ** Komento break on erityisen kätevä, suoritus katkeaa siihen **
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Tässä on olennaisesti vain yhdistetty skriptit ShowSD ja ShowNmin
% Oletuksena: 5 SD-askelta, jonka jälkeen 10 Newton-askelta.
%     SDsteps   = 5 ; Nsteps=10;   alla   
% Lisää tai poista kommentteja "luovasti", muuttele parametreja järkevästi, ...
%
close all           % Sulje grafiikkaikkunat.
Lmax = 1;           % Max line step 
SDsteps   = 5;      % Total number of steps.
manualSteps = 2;    % Montako kuvaa piirretään 1-ul. minimoinnista
                    % Kaksi manuaalista askelta riittäköön yleensä Lmax
                    % ei kasva, kun iteraatio etenee.  Tässä päästään asetta-
                    % maan tehtävän skaala manuaalisesti oikeaan luokkaan.
Nsteps=10;
con_size = 20;      % Korkeusmatriisin resoluutio
xmin=-1;xmax=4;     % Korkeuskäyräpiirrosalueen
ymin=-5;ymax=5;     %               rajat

%%%%%%%%%%%%%% 2-muuttujan fkt. korkeusk. ja alkup. graafisesti %%%%%%%%%%%%%%%%%%
% Kohdefunktion korkeuskäyräpiirros.
% Määrittelemme f:n tässä vektoroidusti. Funktiota fungrad emme voi käyttää
% nätisti meshgrid:n kanssa, vaikka hajoittaisimme sen f- ja grad-osiin, koska
% haluamme yleisen, dimensiosta riippumattoman tavan. Mehgrid taas on raken-
% nettu 2 muuttfkt:lle f(x,y), jossa x ja y ovat (pitkiä) vektoreita.
% Siispä:

f=inline('(x-2).^4+(x-2*y).^2','x','y');

x=linspace(xmin,xmax,con_size);
y=linspace(ymin,ymax,con_size);
[X,Y]=meshgrid(x,y);
Z=f(X,Y);

C = contour(x,y,Z,20);
%C=contour(x,y,Z,0:20);
title('Enter Contour Labels (Optional). To quit, strike <return>.')
xlabel('x')
ylabel('y')
%clabel(C,'manual') 
% Get ready for the iteration.
title('Enter Starting Point for Steepest Descent')
xc = ginput(1); xc=xc'

%%%%%%%%%%%%%% 2-muuttujan fkt. visualisointi päättyy %%%%%%%%%%%%%%%%%%%%%

fc=objf(xc);
gc=objg(xc);

xvals = xc
fvals = fc;
gvals = norm(gc);

%hold on
%plot(xvals(1,end),xvals(2,end),'o')


title('')
% Steepest Descent with line search:

clc
disp('Step       objf            x(1)         x(2)       norm(grad)')
disp('----------------------------------------------------------------')
disp(sprintf('%2.0f     %10.6f     %10.6f   %10.6f      %8.3e ',0,fc,xc,norm(gc)))

for step = 1:SDsteps  
   if step <=manualSteps
     Lmaxtry=1
     while Lmaxtry > 0
        [xnew,fnew,gnew] = SDaskel(xc,fc,gc,Lmaxtry,0);
        title('Jos haluat muuttaa Lmax:ia, anna haun jatkua.')
        Lmaxtry=input('Paina 0, jos ok., anna uusi Lmax, jos haku jatkuu. ')
        if  Lmaxtry > 0
            Lmax=Lmaxtry
        end 
     end
     disp(['Lmax=',num2str(Lmax)])
     disp('Seuraava hakusuunta')
     disp('Paina ENTER, kun haluat jatkaa')
     pause
% Jos minimi on oikeassa päätepisteessä tai jos funktio kasvaa kovasti 
% välillä, valitse suurempi tai pienempi Lmax (kuvasta näet) ja kutsu
% uudestaan.

    else
   [xnew,fnew,gnew] = SDaskel(xc,fc,gc,Lmax,1);
    end

   xvals = [xvals xnew];       xc = xnew;
   fvals = [fvals fnew];       fc = fnew;
   gvals = [gvals norm(gnew)]; gc = gnew;

end
SDtaulukko=[xvals' fvals' gvals']

%%%%%%%%%% End SD %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% Newton alkaa %%%%%%%%%%%%%%%%%%%%%

for step = 1:Nsteps
    xnew=Nminstep(xc);
    xvals = [xvals xnew];xc = xnew;fc=objf(xnew);gc=norm(objg(xnew));
    fvals = [fvals fc];       
    gvals = [gvals gc]; 
end 
kokotaulukko=[xvals' fvals' gvals']

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Show the descent path on the contour plot:

figure(1)
hold on
plot(xvals(1,2:end),xvals(2,2:end),'bo') 
plot(xvals(1,:),xvals(2,:),'r',xvals(1,1),xvals(2,1),'kx')
hold off
%break
title(sprintf('xmin = (%8.3f,%8.3f)   norm(gmin)= %8.4e',xvals(1,end),xvals(2,end),gvals(end)))

% Plot the descent of objf and its gradient:

figure
subplot(2,1,1)
plot(fvals)
title('Value of objf')
xlabel('Iteration')
subplot(2,1,2)
semilogy(gvals)
title('Value of norm(objg).')
xlabel('Iteration')

figure(1)

