% ------------------------------------------------------------------------
% Demonstration package for regularization algorithms.
%
% This requires Hansen package to be installed, more specifically blur
% and shaw functions should be present.
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

clear all;
close all;
m = 20; n = m^2; avgstep = 20;
[G, d_bar, x_bar] = blur(m, 5, 1);
[U,S,V] = svd(full(G));
s = diag(S);
e = normrnd(0,0.05,size(d_bar,1),1);
d = d_bar + e;

coeff_svd = zeros(size(x_bar));
for k=1:n, coeff_svd(k) = abs(U(:,k)'*d); end

coeff_e = zeros(size(x_bar));
for k=1:n, coeff_e(k) = abs(U(:,k)'*e); end

coeff_bar = zeros(size(x_bar));
for k=1:n, coeff_bar(k) = abs(U(:,k)'*d_bar); end

% -- Picard plots
f1 = figure; set(f1, 'Name', 'Picard plot for the unpertubed rhs');
[coeff_bar_avg,foo]=movavg(coeff_bar,avgstep,avgstep);
semilogy((1:n),s,'-rx'); hold on;
semilogy((1:n),coeff_bar,'-bo', 'MarkerSize', 5);
semilogy((1:n),(coeff_bar./s),'-g.');
semilogy((1:n),coeff_bar_avg,'-k','LineWidth',4); hold off;
xlabel ('index, i');
legend ('\sigma_i', '|U_{:i}^T d_{\rm true}|', '|U_{:i}^T d_{\rm true} / \sigma_i|', 'average of |U_{:i}^T d_{\rm true}|', 3);

f2 = figure; set(f2, 'Name', 'Picard plot for the noise vector');
[coeff_e_avg,foo]=movavg(coeff_e,avgstep,avgstep);
semilogy((1:n),s,'-rx'); hold on;
semilogy((1:n),coeff_e,'-bo', 'MarkerSize', 5);
semilogy((1:n),(coeff_e./s),'-g.');
semilogy((1:n),coeff_e_avg,'-k','LineWidth',4); hold off;
xlabel ('index, i');
legend ('\sigma_i', '|U_{:i}^T \eta|', '|U_{:i}^T \eta / \sigma_i|', 'average of |U_{:i}^T \eta|', 2);

f3 = figure; set(f3, 'Name', 'Picard plot for the perturbed rhs');
[coeff_svd_avg,foo]=movavg(coeff_svd,avgstep,avgstep);
semilogy((1:n),s,'-rx'); hold on;
semilogy((1:n),coeff_svd,'-bo', 'MarkerSize', 5);
semilogy((1:n),(coeff_svd./s),'-g.');
semilogy((1:n),coeff_svd_avg,'-k','LineWidth',4); hold off;
xlabel ('index, i');
legend ('\sigma_i', '|U_{:i}^T d|', '|U_{:i}^T d / \sigma_i|', 'average of |U_{:i}^T d|', 2);
