% This code uses the ADM approach to solve the dual of the following
% 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.
%
% 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

% Ting Kei Pong, May 30, 2010
% June 2, 2010: symbols changed to match with writeup.
% June 6, 2010: changed beta for faster convergence.
% June 17, 2010: added two different initializations.
% June 19, 2010: changed beta for faster convergence.

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

% Initialization

tic

% maxiter=1000;
meps=1e-8; % machine epsilon

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

%beta=norm(bx,2)/(2*min(i,j)*min(mu,1));

sigma1=0.95;  % sigma1*(||Q*Q||)<1
sigma2=0.95/min(i,j); % sigma2*(||H*H||)<1

if init==1
  x=zeros(m,n*(i+j-1));
else
  x=q.*bx;
end
lambda=zeros(m*i,n*j);
gamma=zeros(m,n*(i+j-1));
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
  
  xold=x;
  gammaold=gamma;
  lambdaold=lambda;
  
  G=q.*(x/beta+q.*gamma+hL);
  gamma=sigma1/(sigma1+beta)*(qx-beta*G+beta/sigma1*gamma);
  
  lg=q.*gamma+hL;
  
  Hankeltime0 = clock;
  Hkxlg=Hkx(lg,m,n,i,j);
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  g=Hx/beta+Hkxlg;
  
  lsg=lambda-sigma2*g;
  
  svd_time0 = clock;
  [R,S,T]=svd(lsg,'econ');
  svd_time = svd_time + etime(clock,svd_time0);
  
  lambda=R*min(S,mu)*T';
  
  HankelAdjtime0=clock;
  hL=Hkxad(lambda,m,n,i,j);
  HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
  
  G2=(hL+q.*gamma);
  x=x+beta*G2;
  
  Hankeltime0 = clock;
  Hx=Hkx(x,m,n,i,j);
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  iter=iter+1;
  
  if max([norm(x-xold,'fro'),norm(lambda-lambdaold,'fro'),norm(gamma-gammaold,'fro')])<meps % early termination
    fval=1/2*norm(q.*x-qx,'fro')^2+mu*sum(svd(Hx));
    dfeas=norm(hL+q.*gamma,'fro')/max([1,norm(hL,'fro'),norm(q.*gamma,'fro')]);
    dval=-0.5*norm(gamma,'fro')^2+mytrace(qx,gamma);
    gap=abs(fval-dval)/max(abs(dval),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));
    dfeas=norm(hL+q.*gamma,'fro')/max([1,norm(hL,'fro'),norm(q.*gamma,'fro')]);
    dval=-0.5*norm(gamma,'fro')^2+mytrace(qx,q.*gamma);
%    fprintf(' iter, fval, dval, dfea1, time = %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas, toc)
    gap=abs(fval-dval)/max(abs(dval),1);
    if max([gap,dfeas])<tol
      break
    end
  end
end
t=toc;
rk =rank(Hx); 
fprintf('\n DADM terminates, init = %g \n', init);
fprintf(' #iter= %g  dobj= %g  rel. dual feas= %g  pobj= %g  \n rank(Ex)= %g  time= %g \n',iter,dval,dfeas,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);