NC/mp5/Project_5_Maggioni_Claudio/code_template.m

close all;
clear; clc;
addpath /Users/maggicl/Git/matlab2tikz/src/;

%% Load Default Img Data
load('blur_data/B.mat');
B=double(B);
load('blur_data/A.mat');
A=double(A);
ciao = A;

% Show Image
figure
im_l=min(min(B));
im_u=max(max(B));
imshow(B,[im_l,im_u])
title('Blured Image')

% Vectorize the image (row by row)
b=B';
b=b(:);

AT = A * A';
bt = A' * b;

IL = ichol(AT, struct('type', 'nofill', 'diagcomp', 0.01));
y = IL \ b;
x0 = IL' \ y;

[x, rvec] = myCG(AT, bt, x0, 200, 1e-6);
[x2, ~, ~, ~, rvec2] = pcg(AT, bt, 1e-6, 200, IL, IL');

figure;
semilogy(rvec);
hold on;
semilogy(rvec2);
hold off;
title('Residual norms over iteration (y is log)')
matlab2tikz('showInfo', false, '../res_log.tex');


X = zeros(250, 250);
for i = 0:249
    for j = 1:250
        X(i + 1, j) = x(i * 250 + j);
    end
end

figure;
im_l=min(min(X));
im_u=max(max(X));
imshow(X,[im_l,im_u])
title('Sharp Image (myCG)')
matlab2tikz('showInfo', false, '../img_my.tex');


X2 = zeros(250, 250);
for i = 0:249
    for j = 1:250
        X2(i + 1, j) = x2(i * 250 + j);
    end
end

figure;
im_l=min(min(X2));
im_u=max(max(X2));
imshow(X2,[im_l,im_u])
title('Sharp Image (rcg)')
matlab2tikz('showInfo', false, '../img_rcg.tex');



%% Validate Test values
load('test_data/A_test.mat');
load('test_data/x_test_exact.mat');
load('test_data/b_test.mat');

%res=||x^*-A^{-1}b||
res=x_test_exact-inv(A_test)*b_test;
norm(res);
%(Now do it with your CG and Matlab's PCG routine!!!)