% 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 prox-type method in Tseng 2008 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] = rcompletemFBS_uv(q,bz,m,n,j,k,mu,tol,freq,maxiter)

tstart = tic;
fprintf(' ********** mFBS for system realization ***********\n')
% Set parameters
L = 0.5*(1 + sqrt(1 + 4*min(j,k)))/0.95;
qz = q.*bz;
tq = 1 - q;

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

barz = z;
bary = y;
svd_time = 0;
iter = 0;
fbest = inf;
zbest = z;

% mainloop
while iter <= maxiter
  
  % compute v
  gradv = Hkx(-z,m,n,j,k);
  tmp = y - gradv/L;
  svd_time0 = clock;
  [R,S,T] = svd(tmp,'econ');
  svd_time = svd_time + etime(clock,svd_time0);    
  v = R*min(S,mu)*T';
  
  % compute u
  gradu = q.*z - qz + Hkxad(y,m,n,j,k);
  u = z - gradu/L;
  
  % update z 
  z = u - (q.*u - qz + Hkxad(v,m,n,j,k) - gradu)/L;
  
  % update y
  y = v - (Hkx(-u,m,n,j,k) - gradv)/L;
  
  iter = iter + 1;  
  barz = u; % much slower with averaging
  bary = v;
  
  % compute function values and check for termination
  if mod(iter,freq) == 0
    
    Hz = Hkx(barz,m,n,j,k);
    fval = 0.5*norm(q.*barz-qz,'fro')^2 + mu*sum(svd(Hz));
    if fval < fbest
      fbest = fval;
      zbest = barz;
    end
    
    barHy2 = Hkxad(bary,m,n,j,k);
    dval = 0.5*norm(q.*barHy2,'fro')^2 - sum(sum(qz.*barHy2));
    dfeas = norm(tq.*barHy2,'fro')/max([1,norm(barHy2,'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)
