% This code solves   
%  min mu*||H(x)||_*+1/2*||q.*(x-bx)||^2
% where H(x) is a block Hankel matrix formed from x,
% q is the weight matrix, and bx come from measurements: 
% x is m by n*(i+j-1) matrix.
% H(x) is m*i by n*j Hankel matrix. 
%
% The dual of this problem is
%  min   0.5*||gamma||_2^2 - (q.*bx)'*gamma
%  s.t.  H^*(Lambda) + q.*gamma = 0, sigma_max(Lambda) <= mu
%
% This code applies ADM to the dual problem, with proximal terms added to
% regularize the mapping H^*(x).
%
% Input:  
% q       - m by n*(i+j-1) matrix of zero and ones
% bx      - a m by n*(i+j-1) matrix (measured observations)
% m,n,i,j - parameters
% mu      - regularization parameter
% tol     - termination tolerance 
% init    - 1 initialize all variables at zero
%           2 initialize x at bx
% freq    - frequency of checking termination criterion.
% maxiter - maximum number of iterations
%
% Output: 
% x       - m by n*(i+j-1) 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
%
% Written by Maryam Fazel, Ting Kei Pong and Defeng Sun
% Last modified Dec 2, 2012

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

tic

% parameters

meps = 1e-8; % machine epsilon
beta = norm(bx,2)/(8*min(i,j)*mu);
taubeta = 1.61*beta;
sigma2 = 1/min(i,j); % sigma2*(||H*H||)<=1

% Initialization

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;
tq = 1 - q;

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

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

fbest = inf;
xbest = x;

% Loop

iter = 0;

while iter < maxiter
  
  xold = x;
  gammaold = gamma;
  lambdaold = lambda;
  
  gamma = (qx - q.*(x + beta*hL))/(1 + beta);
  
  lgx = q.*gamma + hL + x/beta;
  Hankeltime0 = clock;
  g = Hkx(sigma2*lgx,m,n,i,j);
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  lsg = lambda - 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);
  
  x = x + taubeta*(hL + q.*gamma);
  
  iter = iter + 1;
  
  if max([norm(x-xold,'fro'),norm(lambda-lambdaold,'fro'),norm(gamma-gammaold,'fro')]) < meps % early termination
    
    Hankeltime0 = clock;
    Hx = Hkx(x,m,n,i,j);
    Hankeltime = Hankeltime + etime(clock,Hankeltime0);
    
    fval = 1/2*norm(q.*x-qx,'fro')^2 + mu*sum(svd(Hx));
    dfeas = norm(tq.*hL,'fro')/max([1,norm(hL,'fro')]);
    dval = -0.5*norm(q.*hL,'fro')^2 - mytrace(qx,hL);
%    fprintf(' iter, fval, dval, dfea, time = %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas, toc)
    gap = abs(fval-dval)/max(abs(dval),1);
    fbest = fval;
    xbest = x;
    
    fprintf(' early termination \n')
%    fprintf(' iter, fval, dval, dfea, time = %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas, toc)
    break
    
  end 
  
  if mod(iter,freq) == 0
    
    Hankeltime0 = clock;
    Hx = Hkx(x,m,n,i,j);
    Hankeltime = Hankeltime + etime(clock,Hankeltime0);
    
    fval = 1/2*norm(q.*x-qx,'fro')^2 + mu*sum(svd(Hx));
    dfeas = norm(tq.*hL,'fro')/max([1,norm(hL,'fro')]);
    dval = -0.5*norm(q.*hL,'fro')^2 - mytrace(qx,hL);
    
    if fbest > fval
      fbest = fval;
      xbest = x;
    end
%    fprintf(' iter, fval, dval, dfea, time = %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas, toc)
    gap = abs(fbest-dval)/max(abs(dval),1);
    if max(gap,dfeas) < tol
      break
    end
    
  end
end

fval = fbest;
x = xbest;
t = toc;

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

fprintf(' #iter= %g  dobj= %4.3e  rel. dual feas= %g  pobj= %4.3e  \n rank(Ex)= %g  time= %g \n',iter,dval,dfeas,fval,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);
