clear 
close all


a=-2; b=2; 

n = 50;

t = linspace(a,b,n); % Datenpunkte pro Raumdimension  (nur für Plot)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameter für das Gradientenverfahren 

K = 1000; % Maximale Anzahl von Iterationen 

% Startwerte für die Himmelblau-Funktion
x_0 = [ 0 ; 4];     % --> schnelle Konvergenz gegen [-2.8051 ; 3.1313]
x_0 = [ 0 ; -1];    % --> schnelle Konvergenz gegen [3.5844 ; -1.8481]
x_0 = [ 1.5 ; 2 ];  % --> schnelle Konvergenz gegen [3 ; 2]
x_0 = [ -1 ; 1];   % --> schnelle Konvergenz gegen [-2.8051 ; 3.1313]

% Startwerte für die Rosenbrock-Funktion: bis zu K=10000 Iterationen nötig
%   x_0 = [-1.5 3];    % --> sehr langsame Konvergenz gegen [1 1]
%   x_0 = [-0.5 4];     % --> extrem langsame Konvergenz gegen [1 1]  

% Parameter für die Amijo-Regel
  delta = 0.2; 
  beta = 0.5;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% % % Himmelblau-Funktion 
Funktion = @(x,y) (x.^2+y-11).^2 + (x+y.^2-7).^2;
Funktion_dx = @(x,y) 2*(x.^2+y-11).*2*x + 2*(x+y.^2-7);
Funktion_dy = @(x,y) 2*(x.^2+y-11) + 2*(x+y.^2-7).*2*y;


 


% % % Rosenbrock-Funktion
%  Funktion = @(x,y)    100*(y-x.^2).^2 + (x-1).^2;
%  Funktion_dx = @(x,y) 200*(y-x.^2).*(-2*x) + 2*(x-1);
%  Funktion_dy = @(x,y) 200*(y-x.^2);



% Alternative: Logarithmus der Rosenbrock-Funktion statt der Rosenbrock-Funktion minimieren
% Funktion = @(x,y)    log(100*(y-x.^2).^2 + (x-1).^2);
% Funktion_dx = @(x,y) (200*(y-x.^2).*(-2*x) + 2*(x-1)) / (100*(y-x.^2).^2 + (x-1).^2);
% Funktion_dy = @(x,y) (200*(y-x.^2)) / (100*(y-x.^2).^2 + (x-1).^2);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gradientenabstiegsverfahren %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1=zeros(K+1,1);
x2=zeros(K+1,1);
x1(1)=x_0(1); 
x2(1)=x_0(2);


% Speichern aller Schrittweiten, um zu sehen, wann die Iteration schnell bzw. langsam ist
sigma=zeros(K+1,1); 
Funktionswert=zeros(K+1,1);
Funktionswert(1)=Funktion(x_0(1),x_0(2));

Norm_Grad_F=1;
k=0;


while  Norm_Grad_F>1e-18 && k<K+1
    j=0;
    k=k+1;
    sigmahilf=1;
    
    
    gradient_x=Funktion_dx(x1(k),x2(k));     %1. Komponente des Gradienten an der Stelle x(k) 
    gradient_y=Funktion_dy(x1(k),x2(k));         %2. Komponente des Gradienten an der Stelle x(k)
    Norm_Grad_F = norm([gradient_x,gradient_y],2);


 while ( Funktion(x1(k),x2(k)) < Funktion(x1(k)-sigmahilf* gradient_x,x2(k)-sigmahilf*gradient_y) + delta *sigmahilf * Norm_Grad_F^2) ; 
        j = j+1;
        sigmahilf = beta^(j);
        
 end
    

    sigma(k+1)=sigmahilf;
    x1(k+1) = x1(k)-sigmahilf* gradient_x;
    x2(k+1) = x2(k) -sigmahilf* gradient_y;
    
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Funktionswert = Funktion(x1,x2);


display(['Es wurde ',num2str(k-1),' Iterationen zur Berechnung des Ergebnisses benötigt.'])
display(['Die Lösung ist: [',num2str(x1(k)),',',num2str(x2(k)),'].'])
    
    
h1=(max(x1)-min(x1))*0.2;
h2=(max(x2)-min(x2))*0.2;
a1=floor(min(x1)-h1);
a2=floor(min(x2)-h2);
b1=ceil(max(x1)+h1);
b2=ceil(max(x2)+h2);
t1 = linspace(a1,b1,n);
t2 = linspace(a2,b2,n);
 

% Durch die "plotverschoenerung" schweben die Punkte des Abstiegspfades etwas
% über dem eigentlichen Graph der Funktion, sind aber dadurch besser sichtbar
fh=zeros(n,n);
for i=1:n
    for j=1:n
        fh(i,j) = Funktion(t1(i),t2(j));
    end
end
plotverschoenerung=mean2(fh)/10;


figure;
plot3(x1(1:k),x2(1:k),Funktionswert(1:k)+plotverschoenerung,'ko--','LineWidth',2,'MarkerSize',8); hold on;
plot3(x1(1),x2(1),Funktionswert(1)+plotverschoenerung,'ko','LineWidth',2); 
plot3(x1(k),x2(k),Funktionswert(k)+plotverschoenerung,'go','LineWidth',2); 
surf(t1,t2,fh');
title('Abstiegspfad');


figure;    
plot(x1(1:k),x2(1:k),'ko-','LineWidth',1.5,'MarkerSize',6); hold on;
plot(x1(1),x2(1),'ko','LineWidth',1.5,'MarkerSize',6); 
plot(x1(k),x2(k),'go','LineWidth',1.5,'MarkerSize',6); 
contour(t1,t2,fh',80);   
axis equal
title('Höhenlinien'); 
% saveas(gcf,strcat('figureb2hoehe.fig'))    
    

    
    
figure;
plot(1:k,Funktionswert(1:k));  
title('Funktionswerte');
    
figure
plot(sigma(1:k),'.-');           
title('Schrittweiten');
saveas(gcf,strcat('figure1schritt.fig'))



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Zum Plotten der Funktion  auf [-2,2] x [-2,2]
f=zeros(n,n);
for i=1:n
    for j=1:n
        f(i,j) = Funktion(t(i),t(j));
    end
end

figure;
surf(t,t,f);
title('Plot der Funktion auf [-2,2] x [-2,2]');




% 
% figure;    
% plot(x1(1:k),x2(1:k),'ko-'); hold on;
% plot(x1(k),x2(k),'go'); hold on;
% contour(t,t,log(f'),50); 
% title('Höhenlinien von log(f)');


% figure;    
% contour(t,t,f',50); 
% title('Höhenlinien von log(f)');

