% This code solves   
%  min mu*||H(z)||_*+1/2*||q.*(z-bz)||^2
% where H(z) is a block Hankel matrix formed from z,
% q is the weight matrix, and bx come from measurements: 
% z is m by n*(j+k-1) matrix.
% H(z) is m*j by n*k Hankel matrix. 
%
% The dual of this problem is
%  min   0.5*||nu||_2^2 + (q.*bz)'*nu
%  s.t.  H^*(y) + q.*nu = 0, sigma_max(y) <= mu
%
% This code applies PPG to the primal problem.
%
% Input:  
% q       - m by n*(j+k-1) matrix of zero and ones
% bz      - a m by n*(j+k-1) matrix (measured observations)
% m,n,j,k - parameters
% mu      - regularization parameter
% tol     - termination tolerance 
% freq    - frequency of checking termination criterion.
% maxiter - maximum number of iterations
%
% Output: 
% zbest   - m by n*(i+j-1) matrix, the approximate primal optimal solution
% fbest   - approximate primal optimal value
% dval    - approximate dual optimal value
% time    - cputime
% iter    - number of iterations
% dfeas   - relative dual feasibility at termination
%
% Written by Ting Kei Pong
% Last modified XXX

function [zbest,fbest,dval,time,iter,dfeas] = rcompletePPG(q,bz,m,n,j,k,mu,tol,freq,maxiter)

tstart = tic;
fprintf(' ********** PPG for system realization ***********\n')
% Set parameters
L = 1;
beta = 0.05/L; % seems to be faster when smaller
tau = beta*min(j,k);
qz = q.*bz;
tq = 1 - q;
gamma = 1 + 0.95*min(0.5,1/beta/L - 0.5);
gammabeta = gamma*beta;

% Initialization
z = zeros(m,n*(j+k-1));
y = zeros(m*j,n*k);
Hy2 = Hkxad(y,m,n,j,k);

svd_time = 0;
iter = 0;
fbest = inf;
zbest = z;

% mainloop
while iter <= maxiter
  
  % compute x = grad h(z)
  x = q.*z - qz;
  
  % update y
  tmp = y + Hkx(( z - beta*(x + Hy2))/tau,m,n,j,k);
  
  svd_time0 = clock;
  [R,S,T] = svd(tmp,'econ');
  svd_time = svd_time + etime(clock,svd_time0);  
  
  y = R*min(S,mu)*T';
  Hy2 = Hkxad(y,m,n,j,k);
  
  % update z
  z = z - gammabeta*(x + Hy2);
  
  iter = iter + 1;
  
  % compute function values and check for termination
  if mod(iter,freq) == 0
    
    Hz = Hkx(z,m,n,j,k);
    fval = 0.5*norm(q.*z-qz,'fro')^2 + mu*sum(svd(Hz));
    if fval < fbest
      fbest = fval;
      zbest = z;
    end
    
    dval = 0.5*norm(q.*Hy2,'fro')^2 - sum(sum(qz.*Hy2));
    dfeas = norm(tq.*Hy2,'fro')/max([1,norm(Hy2,'fro')]);
    gap = abs(fbest + dval)/max(fbest,1);
      
%     fprintf(' iter = %d, fval = %g, dval = %g, dfeas = %g\n', iter, fval, dval, dfeas)
  
    if gap < tol && dfeas < tol/5
      break
    end
    
  end
end
time = toc(tstart);
fprintf(' Termination: iter = %d, fval = %g, dval = %g, dfeas = %g, time = %g\n', iter, fval, dval, dfeas, time)
