% This code uses the ADM approach to solve the Hankel matrix completion
% problem:
% min mu*||H(x)||_*+1/(2)*||q.*(x-bx)||^2
% where H(x) is the Hankel matrix formed from x, bx comes
% measurements and q is a vector with entries zero or 1.
% x and q are m by n(i+j-1) vectors.
% H(x) is mi by nj matrix.
% This does not use proximal term in the x update.
%
% Inputs:
% q     = m by n(i+j-1) matrix of 0 and 1
% bx    = m by n(i+j-1) matrix
% mu    = regularization parameter (1/6)
% tol   = termination tolerance
% init  = 1. zero initialization 2. least square initialization
% freq  = frequency of checking termination criterion and displaying
%         function values.
%
% Outputs:
% x     = recovered vector
% fval  = optimal value
% dval  = dual optimal value
% t     = cpu time
% iter  = number of iterations
%
% function called:
% mytrace, Hkx, Hkxad

function [x,fval,dval,t,iter,gap] = rcompleteADM2(q,bx,m,n,i,j,mu,tol,init,freq,maxiter)

tic

% parameters

meps = 1e-8; % machine epsilon
beta = min(i,j)*min(mu,1)/(2*norm(bx,2)); % Too small => slow improvement; too large => overshooting

% Initialization

if init == 1
  x = zeros(m,n*(i+j-1));
  X = zeros(m*i,n*j);
else
  x = q.*bx;
  X = -Hkx(x,m,n,i,j);
end
y = zeros(m,n*(i+j-1));
lambda = zeros(m*i,n*j);
gamma = y;
qx = q.*bx;

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

Hankeltime0 = clock;
Hx = Hkx(x,m,n,i,j);
Hankeltime = Hankeltime + etime(clock,Hankeltime0);

HankelAdjtime0 = clock;
hL = Hkxad(lambda,m,n,i,j);
HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);

% Loop

iter = 0;

while iter < maxiter
  
  yold = y;
  xold = x;
  Xold = X;
  lambdaold = lambda;
  gammaold = gamma;
  
  y = beta/(1+beta)*(gamma/beta-q.*x+qx);
  svd_time0 = clock;
  [R,S,T] = svd(-Hx+lambda/beta,'econ');
  svd_time = svd_time + etime(clock,svd_time0);
  smb = S-mu/beta;
  X = R*max(smb,0)*T'; % matrix shrinkage
  
  HankelAdjtime0 = clock;
  hxL = Hkxad(X,m,n,i,j);
  HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
  
  x = invmap((hL+q.*gamma)/beta - hxL - q.*y + qx,q,m,n,i,j);
  
  Hankeltime0 = clock;
  Hx = Hkx(x,m,n,i,j);
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  lambda = lambda - beta*(Hx+X);
  gamma = gamma - beta*(q.*x+y-qx);  
  
  HankelAdjtime0 = clock;
  hL = Hkxad(lambda,m,n,i,j);
  HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
  
  iter = iter + 1;
  
  if max([norm(x-xold,'fro'),norm(lambda-lambdaold,'fro'),norm(X-Xold,'fro'),norm(y-yold,'fro'),norm(gamma-gammaold,'fro')]) < meps % early termination
  
    fval = 1/2*norm(q.*x-qx,'fro')^2 + mu*sum(svd(Hx));
    dfeas1 = norm(hL+q.*gamma,'fro')/max([1,norm(hL,'fro'),norm(q.*gamma,'fro')]);
    dfeas2 = max(0,max(svd(lambda))-mu)/max(1,norm(lambda,'fro'));
    dval = -0.5*norm(gamma,'fro')^2 + mytrace(qx,gamma);
    gap = abs(fval-dval)/max(abs(fval),1);
    
    fprintf(' early termination \n')
%    fprintf(' iter, fval, dval, dfea1, dfea2, time = %g, %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas1, dfeas2, toc)
    break
    
  end 
  
  if mod(iter,freq) == 0
    
    fval = 1/2*norm(q.*x-qx,'fro')^2 + mu*sum(svd(Hx));
    dfeas1 = norm(hL+q.*gamma,'fro')/max([1,norm(hL,'fro'),norm(q.*gamma,'fro')]);
    dfeas2 = max(0,max(svd(lambda))-mu)/max(1,norm(lambda,'fro'));
    dval = -0.5*norm(gamma,'fro')^2 + mytrace(qx,gamma);
    
%     fprintf(' iter, fval, dval, dfea1, dfea2, time = %g, %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas1, dfeas2, toc)
    gap = abs(fval-dval)/max(abs(fval),1);
    if max([gap,dfeas1,dfeas2]) < tol
      break
    end
    
  end
end

t = toc;

rk = rank(Hx); 

fprintf('\n ADM terminates, init = %g \n', init);

fprintf(' #iter= %g  dobj= %4.3e  rel. dual feas= %g, %g  pobj= %4.3e  \n rank(Ex)= %g  time= %g \n',iter,dval,dfeas1,dfeas2,fval,rk,t);

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