function [fval,z,x,y,iter] = PR_proj_hp(A,b,k,tol,maxiter,beta,eta)

% This computes the intersection of
%
% ||u||_0 <= k, ||u||_infty <= 1e6 and Au = b, A has full row rank, m < n
%
% Using PR splitting to minimize 0.5*d_C^2, where C is the affine space,
% subject to the sparsity and boundedness constraints
%
% This uses heuristics and prints out the gamma used

n = size(A,2);

% initialization
AAT = A*A';
x = zeros(n,1);
y = x;
z = x;
gamma0 = (beta - 2)/(beta + 1)^2; % has convergence guarantee when gamma < gamma0
gamma = eta/beta; % heuristic

iter = 0;
while iter <= maxiter

  xold = x;
  yold = y;
  zold = z;

  % update y: proj onto Ax = b
  tmp = x/(1 + beta*gamma);
  tmp2 = tmp + A'*(AAT\(b - A*tmp)); % projecting tmp onto A*x = b
  y = (x + gamma*tmp2)/(1 + (1 + beta)*gamma);

  % update z
  tmp = (2*y - x)/(1 - beta*gamma);
  tmp0 = max(min(tmp,1e6),-1e6);
  tmp2 = tmp.^2 - (tmp0 - tmp).^2;
  [~,I] = sort(tmp2,'descend');
  z = zeros(n,1);
  z(I(1:k)) = tmp0(I(1:k));

  x = x + 2*(z - y);

  iter = iter + 1;
  err = max([norm(x - xold), norm(y - yold), norm(z - zold)]);
  if err/max([norm(xold),norm(yold),norm(zold),1]) < tol
    break
  end

  if (norm(y - yold) >= 1e3/iter || norm(y) > 1e10) && gamma > gamma0 % decrease gamma to guarantee convergence
    gamma = max(gamma/2,gamma0*0.9999);
  end

%   if mod(iter,freq) == 0
%     projz = z + A'*(AAT\(b - A*z)); % Just for checking fval, not necessary for convergence
%     fval = 0.5*norm(z - projz)^2;
%     projy = y + A'*(AAT\(b - A*y)); % Just for checking D_gamma, not necessary for convergence
%     Dval = 0.5*norm(y - projy)^2 - 3*norm(y - z)^2/2/gamma + (x - y)'*(z - y)/gamma; % essentially decreasing
%     fprintf(' iter %d Dval %g fval %g err %g normx %g gamma %g\n', iter, Dval, fval, err, norm(x), gamma)
%   end

end
projz = z + A'*(AAT\(b - A*z)); % Just for checking fval, not necessary for convergence
fval = 0.5*norm(z - projz)^2;
fprintf(' Termination: ')
fprintf(' iter %d fval %g err %g\n', iter, fval, err)
