function [soln, res, alpha] = RPTRS (G, d, x_bar)

% ------------------------------------------------------------------------
% Parameterized Trust Region Subproblem Regularization Algorithm
%
% solves the problem min ||Gx - d||
%
% INPUT:
%   G     -- operator matrix
%   d     -- rhs (observed data)
%   x_bar -- true solution (optional), this parameter is used for
%            accuracy calculations
%
% OUTPUT:
%   x     -- regularized solution
%   res   -- norm of the residual
%   alpha -- Tikhonov regularization parameter
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

global A a gamma delta bidiag_p bidiag_q

% ------------------------------------------------------------------------
% Initialization
% ------------------------------------------------------------------------

% x_bar is optional
if nargin < 3, x_bar = ones(size(d,1),1); end

% in fact we do not need A matrix explicitly, this is only needed to work
% around a bug in the implementation of eigs()
A  = G'*G; 
a  = G'*d;
dd = d'*d;

% configuration options
lslope_tol1 = 1;

% compute the largest singular value of G
time = cputime; sigmaLA = svds(G,1); time = cputime - time;

% initial interval [low,up) for t and lambda
t_low   = 0;
t_up    = dd;
l_low   = -sigmaLA^2;
l_up    = 0;
itcount = 1;                    % iterations counter
phist   = [];                   % points history
pkappaU = 1e10;                 % previous curvature
pkappaL = 1e10;
nx_bar  = norm(x_bar);

% largest singular value computation time
disp([' ']);
disp(['   +------------------------------------------------------------------------------------+']);
disp([sprintf('   | initial time for sigmaLA = %7g %48c|', time, ' ')]);

% bidiagonalize matrix G
time = cputime; [gamma, delta, bidiag_p, bidiag_q] = lbidiag2 (G, d); time = cputime - time;
disp([sprintf('   | initial time for bidiagG = %7g %48c|', time, ' ')]);

% solve for the initial point: lambda -> t
lambda = l_low;
[t,x,k,time] = l2t(lambda);

% info header
disp(['   +----+--------------+--------------+--------------+---------+---------+--------------+']);
disp(['   | ## |      norm(x) |   norm(Gx-d) | accuracy [%] |    time |       t |       lambda |']);
disp(['   +----+--------------+--------------+--------------+---------+---------+--------------+']);

% ------------------------------------------------------------------------
% Loop
% ------------------------------------------------------------------------
l_damp =  1e-11;
while lambda < l_up - l_damp
    
  % calculate the slope of L-curve and d(lambda)/dt = y(1)
  eps2   = x'*x;                % norm of the solution squared
  r2     = k+dd;                % norm of the residual squared
  lslope = lambda*eps2/r2;
  dldt   = 1/(1+eps2);

  % save current point
  pt     = [t lambda sqrt(eps2) sqrt(r2) lslope dldt];
  phist  = [pt ; phist];

  % update lower bounds on t and lambda
  t_low  = pt(1);
  l_low  = pt(2);

  % relative accuracy
  acc    = norm(x - x_bar)/nx_bar;

  % plot point on the L-curve
  hold on; loglog(pt(3),pt(4),'rx'); hold off;
  
  % information
  disp([sprintf('   | %2d | %12.4e | %12.4e | %12.2f | %7g | %7g | %12.4e |', ...
       itcount, pt(3), pt(4), acc*100, time, pt(1), pt(2))]);
  disp(['   +----+--------------+--------------+--------------+---------+---------+--------------+']);

  % calculate curvature
  [kappaL, kappaU, time] = curvature(G, d, pt(3), r2, pt(2));
  disp([sprintf('   | curvature  = (%12.4e,%12.4e)  time = %7g %26c|', kappaL, kappaU, time, ' ')]);

  % --- termination
  stop = false;
  done = false;
  while ~stop
    if kappaL > pkappaU
      % return previous solution
      res   = phist(2,4);
      alpha = sqrt(-phist(2,2));
      done  = true;
      break;
    end
    
    if kappaU < pkappaL, pkappaL = kappaL; pkappaU = kappaU; stop = true;
    else
      time = cputime;
      [gamma, delta, bidiag_p, bidiag_q] = lbidiag2 (G, d, gamma, delta, bidiag_p, bidiag_q, length(gamma)+1);
      [kappaL, kappaU, time0] = curvature(G, d, pt(3), r2, pt(2));
      time = cputime - time;
      
      disp([sprintf('   | curvature  = (%12.4e,%12.4e)  time = %7g %26c|', kappaL, kappaU, time, ' ')]);

      % recalculate previous kappa bounds
      [pkappaL, pkappaU, time] = curvature(G, d, phist(2,3), phist(2,4)^2, phist(2,2));
      disp([sprintf('   | previous   = (%12.4e,%12.4e)  time = %7g %26c|', pkappaL, pkappaU, time, ' ')]);
    end
  end

  % true solution should sit between two last points
  if done, break; end

  disp([sprintf('   | slope      = %12.4e %57c|', pt(5),  ' ')]);

  % possible pause before next iteration
  % keyboard

  % target epsilon is the current one  
  tgt_eps = pt(3);
      
  % do triangle interpolation on k(t) to obtain the estimate for t
  e2p1 = tgt_eps^2 + 1;
    
  % we choose the point t=dd,k(t)=-dd for the second pivot, as we
  % know this point is always there
  k1    = e2p1*pt(2)-pt(1);
  l1    = e2p1*pt(6)-1;
  k2    = -dd;
  l2    =  -1;                     % slope ~ -1 there, since y(1)->0
  [t,u] = tinterpl (pt(1), dd, k1, k2, l1, l2);

  % display target epsilon
  disp([sprintf('   | target eps = %12.4e %57c|', tgt_eps, ' ')]);

  % information about projected point
  disp([sprintf('   | target t   = %12.4e %57c|', t, ' ')]);
  disp(['   +----+--------------+--------------+--------------+---------+---------+--------------+']);
  
  % now tgt_t is the projected point, check bounds, however both
  % MUST be good
  if t <= t_low || t >= t_up
    disp('>> PANIC: this should never happen!');
    return;
  end

  % save previous solution
  soln = x;

  % need previous lambda as an initial guess
  if itcount > 1, [lambda,x,k,time,y] = t2l(t,lambda,y);
  else,           [lambda,x,k,time,y] = t2l(t,lambda); end

  itcount = itcount + 1;
end

% info footer
disp(['   +------------------------------------------------------------------------------------+']);

% ------------------------------------------------------------------------
% Refine between the last two points
% ------------------------------------------------------------------------
ptR     = phist(1,:);           % right  point
ptC     = phist(2,:);           % center point
ptL     = phist(3,:);           % left   point

% largest known curvature so far
crvC    = (pkappaL + pkappaU)/2;

disp([sprintf('   | curvature [largest] = %12.4e             %36c|', crvC, ' ')]);
disp(['   +-------------------+--------------+--------------+---------+---------+--------------+']);

while ~stop

  done     = false;
  interval = -1;

  while ~done
    % left/right intervals
    if interval < 0, lambda = (ptC(2)+ptL(2))/2;
    else             lambda = (ptC(2)+ptR(2))/2; end

    [t,x,k,time] = l2t(lambda);

    eps2   = x'*x;              % norm of the solution squared
    r2     = k+dd;              % norm of the residual squared

    % relative accuracy
    acc    = norm(x - x_bar)/nx_bar;

    % plot point on the L-curve
    hold on; loglog(sqrt(eps2),sqrt(r2),'co'); hold off;

    % information
    disp([sprintf('   |      %12.4e | %12.4e | %12.2f | %7g | %7g | %12.4e |', ...
         sqrt(eps2), sqrt(r2), acc*100, time, t, lambda)]);
    disp(['   +-------------------+--------------+--------------+---------+---------+--------------+']);

    % calculate curvature
    [kappaL, kappaU, time] = curvature(G, d, sqrt(eps2), r2, lambda);
    disp([sprintf('   | curvature  = (%12.4e,%12.4e)  time = %7g %26c|', kappaL, kappaU, time, ' ')]);
    disp(['   +-------------------+--------------+--------------+---------+---------+--------------+']);

    % curvature
    crv = (kappaL+kappaU)/2;

    if crv < crvC
        % solution
        soln  = x;
        res   = sqrt(r2);
        alpha = sqrt(-lambda);

        done = true;
        stop = true;
    end

    if interval < 0, ptL = [t lambda sqrt(eps2) sqrt(r2)];
    else             ptR = [t lambda sqrt(eps2) sqrt(r2)]; end

    if interval < 0, interval = 1;
    else             done  = true; end
  end
end

% info footer
disp([sprintf('   | norm(x) = %12.4e, alpha = %12.4e, accuracy = %5.2f %20c|', norm(soln), alpha, acc*100, ' ')]);
disp(['   +------------------------------------------------------------------------------------+']);

% ------------------------------------------------------------------------
% --- END of main function
% ------------------------------------------------------------------------
    
% ------------------------------------------------------------------------
% --- conversion lambda -> t
% ------------------------------------------------------------------------
function [t,x,k,time] = l2t(l)
global A a

time = cputime; x = (A-l*speye(size(A)))\a; time = cputime - time;
t = l + a'*x;
k = (x'*x+1)*l - t;

% ------------------------------------------------------------------------
% --- conversion t -> lambda
% ------------------------------------------------------------------------
function [l,x,k,time,y] = t2l(t, lambda, y)
global A a

% options for eigs
opts.tol   = 1e6*eps;
opts.issym = 1;
opts.disp  = 0;

% starting eigenvector (seems that this doesn't help much)
if nargin > 2, opts.v0 = y; end

% construct matrix explicitly or eigs is doing crazy things
% note: this is confirmed to be a bug
time = cputime;
D = [t -a' ; -a A]; [y,l,inf] = eigs(D, 1, lambda, opts);
time = cputime-time;

if inf > 0
  disp('>> WARNING: EIGS did not converge!');
end

% normalize the sign of the eigenvector
if y(1) < 0, y = -y; end

eps2 = (1 - y(1)^2)/y(1)^2;
k    = (eps2+1)*l - t;
x    = y(2:end)/y(1);

% ------------------------------------------------------------------------
% --- triangle interpolation
% ------------------------------------------------------------------------
function [t,y] = tinterpl(t1, t2, k1, k2, slope1, slope2)

tt = pinv([-slope1 1 ; -slope2 1])*[k1 - t1 * slope1 ; k2 - t2 * slope2];
t  = tt(1);
y  = tt(2);

% ------------------------------------------------------------------------
% --- compute a'(A-l*I)^-3a: naive
% ------------------------------------------------------------------------
function [v,time] = est3naive (l)
global A a

time  = cputime;
v     = (A-l*speye(size(A)))\a;
w     = (A-l*speye(size(A)))\v;
v     = v'*w;
time  = cputime - time;

% ------------------------------------------------------------------------
% --- compute a'(A-l*I)^-3a: robust
% ------------------------------------------------------------------------
function [Gp,Rp,time] = est3robust (G, d, l)
global gamma delta bidiag_p bidiag_q

time    = cputime;
[Gp,Rp] = estgr(gamma, delta, G, d, -l, -3);
time    = cputime - time;

% ------------------------------------------------------------------------
% --- curvature
% ------------------------------------------------------------------------
function [kappaL,kappaU,time] = curvature(G, d, e, m, l)

time   = cputime;
[Gp,Rp,time0] = est3robust(G,d,l);

const0 = (e^2)*m;
const1 = const0*((e^4)*(l^2)+(m^2))^(-3/2);
const2 = 2*(e^2)*(l^2)-2*m*l;

kappaL = const1*(const2-const0/Gp);
kappaU = const1*(const2-const0/Rp);
time   = cputime - time;
