% ------------------------------------------------------------------------
% 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;
n = 32;
[G, d_bar, x_bar] = shaw(n);
[U,S,V] = svd(G);

e = 1e-3*rand(size(d_bar));

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

figure;
h = semilogy((1:n),s,'-rx',(1:n),svd_coeff,'-bo',(1:n),(svd_coeff./s),'-g.',[1,n],[eps,eps],'--k');
hold on;
temp = svd_coeff - 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');
legend (h, '\sigma_i', 'U_{:i}^T d_{\rm true}', '|U_{:i}^T d_{\rm true} / \sigma_i|', 'machine \epsilon');

% 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;