clear; close all; clc;

%% Exercise 8

Q = [6     0    -4     0
     0     6     0    -4
    -4     0     6     0
     0    -4     0     6];
c = [ 1 -1  2 -3 ]';
x0 = zeros(length(c),1);
tolerance = 1e-3;

%% Exercise 9

% data = [-5.0000  -96.2607 
%    -4.8000  -85.9893
%    -4.6000  -55.2451
%    -4.4000  -55.6153
%    -4.2000  -44.8827
%    -4.0000  -24.1306
%    -3.8000  -19.4970
%    -3.6000  -10.3972
%    -3.4000   -2.2633
%    -3.2000    0.2196
%    -3.0000    4.5852
%    -2.8000    7.1974
%    -2.6000    8.2207
%    -2.4000   16.0614
%    -2.2000   16.4224
%    -2.0000   17.5381
%    -1.8000   11.4895
%    -1.6000   14.1065
%    -1.4000   16.8220
%    -1.2000   16.1584
%    -1.0000   11.6846
%    -0.8000    5.9991
%    -0.6000    7.8277
%    -0.4000    2.8236
%    -0.2000    2.7129
%          0    1.1669
%     0.2000   -1.4223
%     0.4000   -3.8489
%     0.6000   -4.7101
%     0.8000   -8.1538
%     1.0000   -7.3364
%     1.2000  -13.6464
%     1.4000  -15.2607
%     1.6000  -14.8747
%     1.8000   -9.9271
%     2.0000  -10.5022
%     2.2000   -7.7297
%     2.4000  -11.7859
%     2.6000  -10.2662
%     2.8000   -7.1364
%     3.0000   -2.1166
%     3.2000    1.9428
%     3.4000    4.0905
%     3.6000   16.3151
%     3.8000   16.9854
%     4.0000   17.6418
%     4.2000   46.3117
%     4.4000   53.2609
%     4.6000   72.3538
%     4.8000   49.9166
%     5.0000   89.1652];
% 
% a = data(:,1) ;
% b = data(:,2) ;
% 
% A = [ ones(length(a),1) a a.^2 a.^3 ]; % Vandermonde matrix
% 
% Q = A'*A ;
% c = -A'*b ;
% x0 = zeros(length(c),1);
% tolerance = 1e-3;

%% method

fprintf('Conjugate Gradient method\n\n');
fprintf('iter \t f(x) \t\t ||grad f(x)||\n\n');

iter = 0;

% starting point
x = x0;

while true
    v = 0.5*x'*Q*x + c'*x;
    g = Q*x + c ;
    fprintf('%1.0f \t %1.6e \t %1.4e\n',iter,v,norm(g));
    
    % stopping criterion
    if norm(g) < tolerance
        break
    end
    
    %   search direction
    if iter == 0
        d = -g; 
    else
        beta = (norm(g)^2)/(norm(g_prev)^2);
        d = -g + beta*d_prev;
    end
    
    %   step size
    t = (norm(g)^2)/(d'*Q*d);
    
    %   new point
    iter = iter + 1;
    x = x + t*d;
    d_prev = d ;
    g_prev = g ; 
end

