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

% ------------------------------------------------------------------------
% Initialization
% ------------------------------------------------------------------------
clear all;
shaw_n   = 200;                 % default size of the Shaw problem
blur_n   = 25;                  % default size of the blurring problem
picture  = false;
fig_orig = false;
fig_data = false;
fig_trs  = false;
fig_cg   = false;
fig_lc   = false;

% ------------------------------------------------------------------------
% Loop
% ------------------------------------------------------------------------
while true

  % select problem to solve
  idx = menu('Select a test problem:','Shaw','Blurring','Quit');
  if idx == 3, break; end;

  % init random generator
  randn('seed',70957);

  if idx == 1                   % Shaw
    created = false;
    while true
      str_1 = sprintf('Size: %d', shaw_n);
      idxs  = menu('Options:','Generate',str_1,'->');

      if idxs == 3, break; end
      if idxs == 2, shaw_n = input('Size: '); end

      if idxs == 1
        disp(['>> Data is being generated...']);
        [G, d_bar, x_bar] = shaw(shaw_n);
        e = 1e-3*rand(size(d_bar));
        created = true;
        picture = false;
        disp(['>> ...done']);
        break;
      end
    end

    if ~created, break; end
  end

  if idx == 2                   % Blurring
    created = false;
    while true
      str_1 = sprintf('Size: %d', blur_n);
      idxb  = menu('Options:','Load','Generate',str_1,'->');

      if idxb == 4, break; end
      if idxb == 3, blur_n = input('Size: '); end

      if idxb == 2
        disp(['>> Data is being generated...']);
        [G, d_bar, x_bar] = blur(blur_n, 5, 1);
        e = normrnd(0,0.05,size(d_bar,1),1);
        created = true;
        picture = true;
        disp(['>> ...done']);
        break;
      end

      if idxb == 1
        name = input('Name: ','s');
        data = imread(name);
        disp(['>> Image is being loaded...']);
        [G, d_bar, x_bar] = blur(size(data,1), 5, 0.8);
        e = normrnd(0,0.05,size(d_bar,1),1);
        x_bar = double(data(:));
        d_bar = G*x_bar;
        picture = true;
        created = true;
        disp(['>> ...done']);
        break;
      end
    end

    if ~created, break; end
  end

  % rhs with error
  d = d_bar + e;

  % for CGLS
  tol.G = 0;
  tol.f = 1;
  tol.d = norm(e)/norm(d);
  xini  = zeros(size(d,1),1);
  soln_cg  = false;
  soln_trs = false;

  % processing
  while true
    if picture, idxp = menu('Options:','Build L-curve','CGLS','TRS','Images','->');
    else,       idxp = menu('Options:','Build L-curve','CGLS','TRS','->'); end

    if idxp == 4
      if ~picture, break; end

      if ~fig_orig,                        fig_orig = figure; end
      if ~fig_data,                        fig_data = figure; end
      if ~fig_cg  && size(soln_cg,1)  > 1, fig_cg   = figure; end
      if ~fig_trs && size(soln_trs,1) > 1, fig_trs  = figure; end

      % draw pictures
      figure(fig_orig); set(fig_orig, 'Name', 'true solution: x_bar'); img(x_bar);
      figure(fig_data); set(fig_data, 'Name', 'observed data: d');     img(d);

      if size(soln_cg,1) > 1
        figure(fig_cg); set(fig_cg, 'Name', 'CGLS solution: x_cg'); img(soln_cg);
      end

      if size(soln_trs,1) > 1
        figure(fig_trs); set(fig_trs, 'Name', 'RPTRS solution: x_trs'); img(soln_trs);
      end
    end

    if idxp == 5, break; end

    if idxp == 1,
      if ~fig_lc, fig_lc = figure; set(fig_lc, 'Name', 'L-curve'); end
      figure(fig_lc);
      l_build (G, d, x_bar);
    end

    if idxp == 2                % CGLS method
      if ~fig_lc, fig_lc = figure; set(fig_lc, 'Name', 'L-curve'); end
      figure(fig_lc);
      [soln_cg, phist] = CGLS(G, d, xini, x_bar, tol);
      disp(['>> |x| = ', num2str(phist(1,1)), ', |res| = ', num2str(phist(1,2)), ...
            ', accuracy = ', num2str(100*phist(1,3)), '%']);
    end
    
    if idxp == 3                % TRS method
      if ~fig_lc, fig_lc = figure; set(fig_lc, 'Name', 'L-curve'); end
      figure(fig_lc);
      [soln_trs, res, alpha] = RPTRS (G, d, x_bar);
      disp(['>> |x| = ', num2str(norm(soln_trs)), ', |res| = ', num2str(res), ...
            ', accuracy = ', num2str(100*norm(soln_trs - x_bar)/norm(x_bar)), '%']);
    end
  end

end
