% Subroutine of sysHankeldPPA.m

function [x,Y,k,t] = sysHankeldPPAinit(U,bx,r,mu,tol,maxiter)

fprintf(' Calculate the initial point by DAGP \n')

tstart = tic;
meps = 1e-8;

% Initialization
Ut = U';

n = size(U,2);
p = size(bx,1);
N = size(bx,2) - 1;
m = (r+1)*p;

L = min(r+1,N-r+1);

freq = 10;

HankelTime = 0;
HankelAdjTime = 0;
svd_time = 0;

Y = zeros(m,n);
fbest = inf;
dbest = -inf;

Yold = Y;
% Compute E^*(Y)
HankelAdjTime0 = clock;
EadjY = HradjX(Y*Ut,N,r,p);
HankelAdjTime = HankelAdjTime + etime(clock,HankelAdjTime0);
  
% Compute E(x)
x = bx + EadjY;
HankelTime0 = clock;
Ex = Hrx(x,r)*U;
HankelTime = HankelTime + etime(clock,HankelTime0);
%   xavg = zeros(p,N+1);
theta = 1;
thetaold = theta;
k = 0;

while k < maxiter
  
  % Extrapolate Y in the direction Y-Yold
  Yxtrp = Y + (theta/thetaold-theta)*(Y - Yold);
  
  Exold = Ex;
  x = bx + EadjY;
      
  HankelTime0 = clock;
  Ex = Hrx(x,r)*U;
%   Exavg = Hrx(xavg,r)*U;
  HankelTime = HankelTime + etime(clock,HankelTime0);
    
  if (k > 0) && (mod(k,freq) == 0)
    nmEY = norm(EadjY,'fro')^2/2;
    pobj = nmEY + mu*sum(svd(Ex));
%     pavg = norm(xavg-bx,'fro')^2/2+mu*sum(svd(Exavg));
    dobj = -(nmEY + mytrace(bx,EadjY) );
    if fbest > pobj
      fbest = pobj;
      xbest = x;
    end
    if dbest < dobj
      dbest = dobj;
      Ybest = Y;
    end
%     fprintf(' iter= %g  dobj= %g  pobj= %g  time= %g \n',k,dobj,pavg,toc);
%     fprintf('.');
%     gap = abs(min(pavg,pobj)-dobj)/max(abs(dobj),1);
    gap = abs(fbest - dbest)/max(abs(dbest),1);
    if gap < tol
      break
    end
  end
    
  G = Ex + (theta/thetaold-theta)*(Ex - Exold);
 
  YxGL = Yxtrp - G/L;
  svd_time0 = clock;
  [R,D,S] = svd(YxGL,'econ');			%faster alternative?
  svd_time = svd_time + etime(clock,svd_time0);

  Yold = Y;
  Y = R*min(D,mu)*S';
  HankelAdjTime0 = clock;
  EadjY = HradjX(Y*Ut,N,r,p);
  HankelAdjTime = HankelAdjTime + etime(clock,HankelAdjTime0);
%   xavg = (1-theta)*xavg+theta*(bx+EadjYxtrp);
  thetaold = theta;
  theta = 2*theta/(sqrt(theta^2+4)+theta);
  k = k + 1;
    
  % early termination if change in iterate is small
  step = norm(Y-Yold,'fro');
  if step < meps
    Ybest = Y;
    xbest = x;
    fprintf(' termination due to small stepsize = %g \n',step)
    break
  end

end 

Y = Ybest;
x = xbest;
t = toc(tstart);

fprintf('#iter=%g, time=%g, \n ', k,t)