% ------------------------------------------------------------------------
% 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;
[G, d_bar, x_bar] = blur(m, 5, 1);
[U,S,V] = svd(full(G));
e = normrnd(0,0.05,size(d_bar,1),1);

A = G'*G;
s = diag(S);

d = d_bar + e;
a = G'*d;

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

[svd_coeff_avg,foo]=movavg(svd_coeff,50,50);

figure;
h = semilogy((1:n),s,'-rx',(1:n),svd_coeff_avg,'-bo',(1:n),(svd_coeff_avg./s),'-g.');
hold on;
temp = svd_coeff_avg - s;
for i=n:-1:1, if temp(i)<0, break; end; end;
plot ([i+1,i+1],[min(svd_coeff),max(svd_coeff./s)],':k');
hold off;

xlabel ('index, i');

% l-curve
figure;
l_build (G, d, x_bar);

% around cross point
hold on;
for j=i-1:i+1,
  alpha = S(j,j);

  x = zeros(size(x_bar));
  for k=1:n, x = x + (S(k,k)^2/(S(k,k)^2+alpha^2))*(V(:,k)*U(:,k)'*d/S(k,k)); end
  r = G*x - d;
  plot(norm(x),norm(r),'rx');
end;
hold off;