% http://www.math.hut.fi/teaching/numsym/04/H/ShowSDN.m
% Script File: ShowSDN
% Modifioitu C. van Loanin: Scientific Computing using Matlab-kirjasta
% Minimization by SD + Newton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some Finnish instructions
%
% 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 **
  %        (aina ei näy toimivan)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Default: 5 SD-steps followed by 10 Newton steps.
SDsteps   = 5 ; Nsteps=10;    
% Add or remove comments according to needs.
%
close all           % Close graphics windows.
Lmax = 1;           % Max line step 
SDsteps   = 20;      % Total number of steps.
manualSteps = 2;    % Give user (yourself) control for initial max step.
                    % 2 manual steps might be enough but choose as you like.
                    % ei kasva, kun iteraatio etenee.  Tässä päästään asetta-
                    % maan tehtävän skaala manuaalisesti oikeaan luokkaan.

Nsteps=15;           % Change the above default
con_size = 60;      % Resolution of contourplot
xmin=-1; xmax=2;     % Contourplot axes limits
ymin=-0.5;ymax=2;     %               

%%%%%%%%%%%%%% 2-muuttujan fkt. korkeusk. ja alkup. graafisesti %%%%%%%%%%%%%%%%%%
% Let's take the Rosenbrock banana function as our objective function.
% Some Finnish :
% 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('100*(y-x.^2).^2 + (1-x).^2','x','y');  % Rosenbrock

% This is the standard Matlab-way of producing surface and contourplot
% datamatrix:
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,40);
C=contour(x,y,Z,0:20);

xlabel('x')
ylabel('y')
title('Enter Contour Labels (Optional). To quit, strike <return>.')

clabel(C,'manual') 
% Get ready for the iteration.

title('Enter Starting Point for Steepest Descent')
xc = ginput(1); xc=xc'   % This is nice in many instances (Not in Maple).
%xc=[0;2]  % If you want to give manually, just edit.

%%%%%%%%%%%%%% end of visualization of objf. %%%%%%%%%%%%%%%%%%%%%
%Määritellään objf ja objg inline-funktioina, edellinen on sama kuin yllä oleva
%f, paitsi rakennetaan ne vektorimuuttujan funktioiksi, kuten Maple-ws:ssä.
% Muistathan funktiot f ja fv Maplessa

%objf=inline('f(x(1),x(2))','x')   % Inline ei hallitse sisäkkäisiä 
                                  %  funktiokutsuja.
%Tehdään sittenkin  objf.m-tiedosto
%%% objf.m %%%%%%%%
function z=objf(X)
x=X(1); y=X(2);
z=100*(y-x.^2).^2 + (1-x).^2;
%%%%%%%%%%%%%%%%%%

% Testaus:
>> [f(1,2),objf([1;2])]
ans =
   100   100
% No nyt toimii!

function v=objg(X)
%  z=objg(X)
%  Syöte: n-sarakevektori X
%  Tulos: kohdefunktion arvon X:ssä, skalaari  
x=X(1);y=X(2);
v=[-400*(y-x^2)*x-2+2*x;  % Rosenbrock'n funktion gradeintti
   200*y-200*x^2];



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

>> objg(xc)
ans =
    0.2307
   -2.7669

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


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


for step = 1:SDsteps  
   if step <=manualSteps
     Lmaxtry=1
     while Lmaxtry > 0
        [xnew,fnew,gnew] = SDaskel(xc,fc,gc,Lmaxtry,0);
        title('If you want to change Lmax, continue with line search.')
        Lmaxtry=input('Press 0, if ok., give new Lmax if you go on ')
        if  Lmaxtry > 0
            Lmax=Lmaxtry
        end 
     end
     disp(['Lmax=',num2str(Lmax)])
     disp('Next search direction')
     disp('Paina ENTER, kun haluat jatkaa')
     pause
% If the minimum is at right endpoint or if it is increasing rapidly,
% change larger or smaller Lmax accordingly. (You will see in the figure.)

    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
SDtable = [xvals' fvals' gvals']

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

%%%%%%%%% Start Newton  %%%%%%%%%%%%%%%%%%%%%

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 
SDNtable = [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)
pause
figure
disp('Surface plot, enter')
surfl(x,y,Z)
