% /home/apiola/opetus/projektit/odense2000/v2
% June 11, 2000
% Script File: ShowSDN
% Minimization by SD + Newton     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Default: 5 SD-steps followed by 10 Newton steps.
SDsteps   = 5 ; 
Nsteps=10;      % Newton steps    
% Add or remove comments according to needs.
%
format compact
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.
                    % Set the scale manually
                    
Nsteps=5;           % Change the above default
con_size = 60;      % Resolution of contourplot
xmin=-1; xmax=2;     % Contourplot axes limits
ymin=-0.5;ymax=2;     %               
%% Contours of f(x,y) and init point graphically
%
% Let's take the Rosenbrock banana function as our objective function.
% Define f in "array smart" way. 

% (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 (long) vektors.
% Siispa:

f=@(x,y) 100*(y-x.^2).^2 + (1-x).^2;  % Rosenbrock

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);
grid on
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'   % This is nice in many instances.
%xc=[0;2]  % If you want to give manually, just edit.
hold on
plot(xc(1),xc(2),'*')
% end of visualization of objf.
%% Use symbolic toolbox for gradient:
%{
>> syms x y
>> f(x,y)
ans =
(x - 1)^2 + 100*(- x^2 + y)^2
>> grf=gradient(f(x,y),[x,y])
grf =
 2*x - 400*x*(- x^2 + y) - 2
           - 200*x^2 + 200*y
>> %Check:
>> diff(f(x,y),x)
ans =
2*x - 400*x*(- x^2 + y) - 2
>> diff(f(x,y),y)
ans =
- 200*x^2 + 200*y
%}
% Thus gradient can be defined as:
objg=@(x) [2*x(1) - 400*x(1)*(- x(1)^2 + x(2)) - 2;...    
         - 200*x(1)^2 + 200*x(2)]
%%  
objf=@(x) f(x(1),x(2))
fc=objf(xc);
gc=objg(xc);
gc=gc/norm(gc);


%hold on
%plot(xvals(1,end),xvals(2,end),'o')
%% Trial step starting at xc in the direction of negative gradient.
%Tmax=5; % 1^st trial => minf = 0.4521
% t=0:0.1:Tmax;
t=linspace(0,Tmax);  % 2^nd trial Tmax =0.2525 minf = 0.4068
N=length(t);
xvals=zeros(2,N);
gvals=xvals;
fvals=zeros(1,N);
for k=1:N
    xvals(:,k)=xc-t(k)*gc;
    fvals(k)=objf(xvals(:,k));
end    
[minf,minind]=min(fvals)
Tmax=t(minind)
% Return above
%%
hold on
plot(xvals(1,:),xvals(2,:),'.-')
Tmax=t(minind)
%{

>> fvals(1:5)
ans =
    9.6619    4.6990    1.7526    0.8419    1.9859
>> xvals(:,1:5)
ans =
    0.0472    0.0598    0.0724    0.0849    0.0975
    0.2981    0.1989    0.0997    0.0005   -0.0987
%}
t(3:4)

%% 

%%

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)
