%% 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
close all

% Set to non-zero to generate LaTeX for graphs 
enable_m2tikz = 1;
if enable_m2tikz
    addpath /home/claudio/git/matlab2tikz/src
else
    matlab2tikz = @(a,b,c) 0;
end
    
syms x y;
f = (1 - x)^2 + 100 * (y - x^2)^2;
global fl
fl = matlabFunction(f);

%% 1.3 - Newton and GD solutions and energy plot 

[x1, xs1, gs1] = Newton(f, [0;0], 50000, 1e-6, true);
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;
[x2, xs2, gs2] = Newton(f, [0;0], 50000, 1e-6, false);
plot(xs2(1, :), xs2(2, :), 'Marker', '.');
fprintf("Newton: it=%d\n", size(xs2, 2)-1);

[x3, xs3, gs3] = GD(f, [0;0], 50000, 1e-6, true);
plot(xs3(1, :), xs3(2, :), 'Marker', '.');
fprintf("GD backtracking: it=%d\n", size(xs3, 2)-1);

[x4, xs4, gs4] = GD(f, [0;0], 50000, 1e-6, false);
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;
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([5e-19, 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');

figure;

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

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