%% Homework 3 - Optimization Methods
% Author: Claudio Maggioni
%
% Note: exercises are not in the right order due to matlab constraints of
% functions inside of scripts.

clc
clear

% Set to non-zero to generate LaTeX for graphs 
enable_m2tikz = 0;
if enable_m2tikz
    addpath /home/claudio/git/matlab2tikz/src
else
    matlab2tikz = @(a,b,c) 0;
end

%% 1.1 - Rosenbrock function definition and surf plot
close all
syms x y;
f = (1 - x)^2 + 100 * (y - x^2)^2;
global fl
fl = matlabFunction(f);

xs = -0.2:0.01:1.2;
Z = zeros(size(xs,2), size(xs,2));

for x = 1:size(xs, 2)
    for y = 1:size(xs, 2)
        Z(x,y) = fl(xs(x), xs(y));
    end
end

surf(xs, xs, Z, 'EdgeColor', 'none', 'FaceAlpha', 0.4);
hold on
plot3([0,1],[0,1],[fl(0,0),fl(1,1)],'-r.');
plot3([1],[1],[fl(1,1)],'-k.');
hold off

figure;


%% 1.2 - Minimizing the Rosenbrock function

[x1, xs1, gs1] = Newton(f, [0;0], 50000, 1e-6, true);
[x2, xs2, gs2] = Newton(f, [0;0], 50000, 1e-6, false);
[x3, xs3, gs3] = GD(f, [0;0], 50000, 1e-6, true);
[x4, xs4, gs4] = GD(f, [0;0], 50000, 1e-6, false);

%% 1.3 - Newton and GD solutions and energy plot 

plot(xs1(1, :), xs1(2, :), 'Marker', '.');
fprintf("Newton backtracking: it=%d\n", size(xs1, 2)-1);
xlim([-0.01 1.01]) 
ylim([-0.01 1.01])
hold on;
plot(xs2(1, :), xs2(2, :), 'Marker', '.');
fprintf("Newton: it=%d\n", size(xs2, 2)-1);

plot(xs3(1, :), xs3(2, :), 'Marker', '.');
fprintf("GD backtracking: it=%d\n", size(xs3, 2)-1);

plot(xs4(1, :), xs4(2, :), 'Marker', '.');
fprintf("GD: it=%d\n", size(xs4, 2)-1);
hold off;
legend('Newton + backtracking', 'Newton', 'GD + backtracking', 'GD (alpha=1)')
sgtitle("Iterations of Newton and Gradient descent methods over 2D energy landscape");
matlab2tikz('showInfo', false, './ex1-3.tex');

figure;
plot(xs4(1, :), xs4(2, :), 'Marker', '.');
legend('GD (alpha=1)')
sgtitle("Iterations of Newton and Gradient descent methods over 2D energy landscape");
matlab2tikz('showInfo', false, './ex1-3-large.tex');

%% 2.3 - BGFS solution and energy plot 

figure;
[x5, xs5, gs5] = BGFS(f, [0;0], eye(2), 50000, 1e-6);
xlim([-0.01 1.01]) 
ylim([-0.01 1.01])
plot(xs5(1, :), xs5(2, :), 'Marker', '.');
fprintf("BGFS backtracking: it=%d\n", size(xs5, 2)-1);

sgtitle("Iterations of BGFS method over 2D energy landscape");
matlab2tikz('showInfo', false, './ex2-3.tex');

%% 1.4 - Newton and GD gradient norm log

figure;
semilogy(0:size(xs1, 2)-1, gs1, 'Marker', '.');
ylim([5e-10, 1e12]);
xlim([-1, 30]);
hold on
semilogy(0:size(xs2, 2)-1, gs2, 'Marker', '.');
semilogy(0:size(xs3, 2)-1, gs3, 'Marker', '.');
semilogy(0:size(xs4, 2)-1, gs4, 'Marker', '.');
hold off

legend('Newton + backtracking', 'Newton', 'GD + backtracking', 'GD (alpha=1)')
sgtitle("Gradient norm w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-grad.tex');


figure;
plot(0:size(xs1, 2)-1, gs1, 'Marker', '.');
ylim([5e-10, 401]);
xlim([-1, 30]);
hold on
plot(0:size(xs2, 2)-1, gs2, 'Marker', '.');
plot(0:size(xs3, 2)-1, gs3, 'Marker', '.');
plot(0:size(xs4, 2)-1, gs4, 'Marker', '.');
hold off

legend('Newton + backtracking', 'Newton', 'GD + backtracking', 'GD (alpha=1)')
sgtitle("Gradient norm w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-grad.tex');

figure;
semilogy(0:size(xs3, 2)-1, gs3, 'Marker', '.');
ylim([1e-7, 1e10]);
hold on
semilogy(0:size(xs4, 2)-1, gs4, 'Marker', '.');
hold off

legend('GD + backtracking', 'GD (alpha=1)')
sgtitle("Gradient norm w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-grad-large.tex');

figure;
ys1 = funvalues(xs1);
ys2 = funvalues(xs2);
ys3 = funvalues(xs3);
ys4 = funvalues(xs4);
ys5 = funvalues(xs5);
semilogy(0:size(xs1, 2)-1, ys1, 'Marker', '.');
ylim([0, 1e12]);
xlim([-1, 30]);
hold on
semilogy(0:size(xs2, 2)-1, ys2, 'Marker', '.');
semilogy(0:size(xs3, 2)-1, ys3, 'Marker', '.');
semilogy(0:size(xs4, 2)-1, ys4, 'Marker', '.');
hold off

legend('Newton + backtracking', 'Newton', 'GD + backtracking', 'GD (alpha=1)')
sgtitle("Objective function value w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-ys.tex');

plot(0:size(xs1, 2)-1, ys1, 'Marker', '.');
ylim([0, 101]);
xlim([-1, 30]);
hold on
plot(0:size(xs2, 2)-1, ys2, 'Marker', '.');
plot(0:size(xs3, 2)-1, ys3, 'Marker', '.');
plot(0:size(xs4, 2)-1, ys4, 'Marker', '.');
hold off

legend('Newton + backtracking', 'Newton', 'GD + backtracking', 'GD (alpha=1)')
sgtitle("Objective function value w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-ys.tex');

figure;

semilogy(0:size(xs3, 2)-1, ys3, 'Marker', '.');
hold on
semilogy(0:size(xs4, 2)-1, ys4, 'Marker', '.');
hold off

legend('GD + backtracking', 'GD (alpha=1)')
sgtitle("Objective function value w.r.t. iteration number for Newton and GD methods");
matlab2tikz('showInfo', false, './ex1-4-ys-large.tex');

%% 2.4 - BGFS gradient norms plot

figure;
semilogy(0:size(xs5, 2)-1, gs5, 'Marker', '.');

sgtitle("Gradient norm w.r.t. iteration number for BGFS method");
matlab2tikz('showInfo', false, './ex2-4-grad.tex');

%% 2.4 - BGFS objective values plot

figure;
semilogy(0:size(xs5, 2)-1, ys5, 'Marker', '.');

sgtitle("Objective function value w.r.t. iteration number for BGFS methods");
matlab2tikz('showInfo', false, './ex2-4-ys.tex');

function ys = funvalues(xs)
    ys = zeros(1, size(xs, 2));
    global fl
    for i = 1:size(xs, 2)
        ys(i) = fl(xs(1,i), xs(2,i));
    end
end