% This code solves   
%  min mu*||H(x)*U||_*+1/2*||x-bx||^2
% where H(x) is a block Hankel matrix formed from x, and bx and U come from measurements: 
% U is a matrix such that U'U=I.
% x is p by N+1 matrix.
% U is N-r+1 by q matrix.
% H(x) is p(r+1) by N-r+1 Hankel matrix. 
%
% The Fenchel dual of this problem is
%  min   0.5*||E^*(Y)||_2^2 + bx'*E^*(Y)
%  s.t.  sigma_max(Y) <= mu,
% where E(x)=H(x)*U.
%
% This code applies PPA to the Fenchel dual problem.
%
% Input:  
% U       - N-r+1 by q matrix such that U'*U=I
% bx      - a p by N+1 real matrix (measured observations)
% r       - parameter
% mu      - regularization parameter
% tol     - termination tolerance 
% init    - 1. initialize at zero; 
%           2. initial point generated by DAGP with tol=0.005) 
% maxiter - maximum number of iterations
%
% Output: 
% x         - p by N+1 real matrix, the approximate primal optimal solution
% fval      - approximate primal optimal value
% dval      - approximate dual optimal value
% t         - cputime
% iter      - number of iterations
% gap       - duality gap
% iteration - iteration used in initialization by DAGP
% runtime   - runtime for DAGP
%
% Written by Maryam Fazel, Weimin Miao, Ting Kei Pong and Defeng Sun
% Last modified Nov 27, 2012


function [x,fval,dval,t,iter,gap2,iteration,runtime] = sysHankeldPPA(U,bx,r,mu,tol,init,maxiter)


fprintf(' Start dPPA\n')

tic

meps = 1e-8; % machine epsilon

% Initialization

Ut = U';
q = size(U,2);
p = size(bx,1);
N = size(bx,2) - 1;
m = (r+1)*p;
R0 = min(r+1,N-r+1);

if (init == 1) 
  x = zeros(p,N+1);
  lambda = zeros(m,q);
  iteration = 0;
  runtime = 0;
elseif (init == 2)
  AGPtol = 0.005;
  AGPiter = 2000;
  [x,lambda,iteration,runtime] = sysHankeldPPAinit(U,bx,r,mu,AGPtol,AGPiter);
end

Hankeltime = 0;
HankelAdjtime = 0;
svd_time = 0;

inmaxiter = 1000; % max iter for subproblem
index = 0; % the index for checking the criterion to update c

c = 1; % c represent lambda in paper, the initial value % The number of iteration increases as c decreases and the cost time in each step increases as c increases.
cmax = 1e8;
rho = 2; % parameter to update lambda

intolmin = tol*1e1; 
intol = intolmin*40;  % The initalized inner tolerance

Hankeltime0 = clock;
Hx = Hrx(x,r)*U;
Hankeltime = Hankeltime + etime(clock,Hankeltime0);

fbest = inf;
dbest = -inf;
xbest = x;

for iter = 1:maxiter
  if (iter == 1)
    L = 2*c*R0 + 1;
    s = 1.95/L;
  end
  alpha = 0.3; % parameter for line search

  GZ = lambda - c*Hx;
  
  svd_time0 = clock;
  [R,S,T] = svd(GZ,'econ');
  svd_time = svd_time + etime(clock,svd_time0);
      
  P = R*max(S - mu,0)*T';
          
  GZP = GZ - P;
  HankelAdjtime0 = clock;
  HaGZP = HradjX(GZP*Ut,N,r,p); % apply the adjoint of Hankel
  HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0); 

  diffx = x - bx;
  grad = HaGZP - diffx; 
  indval = mytrace(Hx,lambda) - c/2*norm(Hx,'fro')^2 + norm(P,'fro')^2/(2*c) - norm(diffx,'fro')^2/2;
         
  for initer = 1:inmaxiter
    indvalold = indval;
    for ininiter = 1:inmaxiter
      % fprintf('-----line search starts-----\n')
      xmid = x + s*grad; 
              
      Hankeltime0 = clock;
      Hx = Hrx(xmid,r)*U;
      Hankeltime = Hankeltime + etime(clock,Hankeltime0);
      
      GZ = lambda - c*Hx;
      
      svd_time0 = clock;
      [R,S,T] = svd(GZ,'econ');
      svd_time = svd_time + etime(clock,svd_time0);
      
      P = R*max(S - mu,0)*T';      
             
      diffx = xmid - bx;
      indval = mytrace(Hx,lambda) - c/2*norm(Hx,'fro')^2 + norm(P,'fro')^2/(2*c) - norm(diffx,'fro')^2/2;
%       fprintf(' s, ininiter, val, valold, armijoval1, armijoval = %g, %g, %g, %g, %g, %g\n', norm(s,'fro'), ininiter, indval, indvalold, indval-indvalold, 1e-4*s*norm(grad,'fro')^2)
%       fprintf('!') % show the iteration for line search
      
      if (indval - indvalold >= 1e-4*s*norm(grad,'fro')^2)
        % fprintf(' -----line search ends------\n')
        break
      else
        s = s*alpha;
      end
    end
    
    GZP = GZ - P;
    HankelAdjtime0 = clock;
    HaGZP = HradjX(GZP*Ut,N,r,p); % apply the adjoint of Hankel
    HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0); 
    
    grad = HaGZP - diffx;  
    if (norm(grad,'fro') < intol) 
%       fprintf(' initer, grad, indiff = %g, %g, %g, %g \n', initer, norm(grad,'fro'), indval-indvalold)
%       fprintf(' *******inner problem is solved!*******\n')
      x = xmid;
      break
    end
    if (ininiter == 1)
      s = s*1.11;  
    end
    % fprintf(' ++++++ s=%g, grad=%g, indval=%g, x-xmid=%g\n ',s,norm(grad,'fro'), indval, norm(x-xmid,'fro')')
    x = xmid;
    % fprintf('.') % show the iteration for subproblem
  end
  %fprintf('*')

  lambda = GZP;
  HaU = HaGZP;
  
  if iter > 1
    if max([norm(x-xold,'fro'), norm(lambda-lambdaold,'fro')]) < meps % early termination 
      fval = norm(diffx,'fro')^2/2 + mu*sum(svd(Hx,'econ'));
      
      x2 = HaU + bx;
      Hankeltime0 = clock;
      Hx2 = Hrx(x2,r)*U;
      Hankeltime = Hankeltime + etime(clock,Hankeltime0);
    
      nrm = 0.5*norm(HaU,'fro')^2;
      fval2 = nrm + mu*sum(svd(Hx2)); 
      dval = -(nrm + mytrace(HaU,bx));
      gap2 = abs(min(fval,fval2)-dval)/max(abs(dval),1);
      fprintf(' early termination \n')
      fprintf(' iter, fval, dval, time = %g, %g, %g, %g, %g\n', iter, fval, dval, toc)
      fbest = min(fval,fval2);
      dbest = dval;
      if fval > fval2
        xbest = x2;
      else
        xbest = x;
      end
      fprintf('.')
      break
    end
  end

  fval = norm(diffx,'fro')^2/2 + mu*sum(svd(Hx,'econ'));
  
  x2 = HaU + bx;
  Hankeltime0 = clock;
  Hx2 = Hrx(x2,r)*U;
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  nrm = 0.5*norm(HaU,'fro')^2;
  fval2 = nrm + mu*sum(svd(Hx2)); 
  dval = -(nrm + mytrace(HaU,bx));
  
  if dbest < dval
    dbest = dval;
  end
  if fbest > fval
    fbest = fval;
    xbest = x;
  end
  if fbest > fval2
    fbest = fval2;
    xbest = x2;
  end
  gap = abs(fval-dval)/max(abs(dval),1);
  gap2 = abs(fbest-dbest)/max(abs(dbest),1);
  % fprintf('\n c, intol, fval, dval, gap = %g, %g, %g, %g, %g\n', c, intol, fval, dval, gap)
  if gap2 < tol
    break
  end
% The following are the detailed skills to update the paramter c and intol, in order to guarantee the convergence.    
  if gap > 0.001
    alpha = 0.95;
  else
    alpha = 0.98;
  end
  if (iter > 1)
    relgapold = gap/gapold;
    if (intol > intolmin)
      relgap_benchmark = gap/gap_benchmark;
      if (relgap_benchmark < 1)
        crec = c;
        xrec = x;
        Hxrec = Hx;
        lambdarec = lambda;
        gaprec = gap;
      end
      if (relgap_benchmark > alpha^(index+1))
        index = index + 1;
        if (index == 3)
          c = crec;
          x = xrec;
          Hx = Hxrec;
          lambda = lambdarec;
          gap = gaprec;
          gap_benchmark = gaprec;
          index = 0;
          intol = max(intol/2, intolmin);
        else
          if (relgapold > alpha)
            c = rho*c;
          end
        end
      else
        index = 0;
        crec = c;
        xrec = x;
        Hxrec = Hx;
        lambdarec = lambda;
        gaprec = gap;
        gap_benchmark = gap;
      end
    else
      if (relgapold > alpha)
        c = rho*c;
      end
    end
  else
    index = 0;
    crec = c;
    xrec = x;
    Hxrec = Hx;
    lambdarec = lambda;
    gaprec = gap;
    gap_benchmark = gap;
  end
  if (c > cmax)
    fprintf(' early termination since c exceeds cmax!')
    break
  end
  gapold = gap;
  xold = x;
  lambdaold = lambda;
end

x = xbest;
fval = fbest;
dval = dbest;

t = toc;
if (init == 1)
    fprintf('\n Dual PPA terminates, initial =%g\n', init);
elseif (init == 2)
    fprintf('\n Dual PPA terminates, initial =%g, AGPtol=%g\n', init, AGPtol);
elseif (init == 3)
    fprintf('\n Dual PPA terminates, initial =%g, AGPiter=%g\n', init, AGPiter);
end
fprintf('#iter= %g, time= %g\n dobj= %g, pobj= %g\n',iter, t, dval, fval);
fprintf(' Total time for computing the the Hankel operator & its adjoint: %g  %g\n', Hankeltime, HankelAdjtime);
fprintf(' Total time for svd:  %g\n', svd_time);



    


