function [fval,x,iter] = Alt_proj(A,b,k,tol,maxiter,freq)

% This computes the intersection of 
%
% ||u||_0 <= k, ||u||_infty <= 1e6 and Au = b, A has full row rank, m < n
%
% using alternating projection.

n = size(A,2);

% initialization
AAT = A*A';
x = zeros(n,1);
projx = x + A'*(AAT\(b - A*x));  % projection onto affine space

iter = 0;
while iter <= maxiter
  
  xold = x; 

  % projection onto sparsity
  tmp0 = max(min(projx,1e6),-1e6); % projection onto box
  tmp2 = projx.^2 - (tmp0 - projx).^2;
  [~,I] = sort(tmp2,'descend');
  x = zeros(n,1);
  x(I(1:k)) = tmp0(I(1:k));
  
  iter = iter + 1;  
  err = norm(x - xold);
  if err/max(norm(xold),1) < tol
    break
  end
  
  % projection onto affine space
  projx = x + A'*(AAT\(b - A*x)); 
  if mod(iter,freq) == 0
    fval = 0.5*norm(x - projx)^2; 
    fprintf(' iter %d fval %g err %g\n', iter, fval, err)
  end
  
end
projx = x + A'*(AAT\(b - A*x));
fval = 0.5*norm(x - projx)^2; 
fprintf(' Termination: ')
fprintf(' iter %d fval %g err %g\n', iter, fval, err)