function x_hat = camp(y, A, lambda)
    [M,N] = size(A);
    assert(M == length(y) && isvector(y));
    
    A_H = A';
    x_hat = zeros(N,1); z = zeros(M,1);
    dx = zeros(N,1);
    dy = zeros(N,1);

    t = 0;
    x_hat_prev = ones(N,1);
    
    while t<= 50 && norm(x_hat_prev-x_hat) > 1e-5
        % update z. 
        z = y-A*x_hat;
        for a=1:M
            s = A_H(:,a)*z(a);
            z(a) = z(a) - (A(a,:)*(dx.*real(s)));
            z(a) = z(a) - (A(a,:)*(dy.*imag(s)));
        end
        sigma2 = norm(z)^2/M;

        % update x_hat. 
        v = A_H*z;
        x_hat_prev = x_hat;
        [x_hat, dx, dy] = shrink(x_hat+v, lambda, sigma2);
        t = t+1;
    end

end

function [x_hat, dx, dy] = shrink(r, lambda, sigma2)
    b = abs(r)>lambda*sqrt(sigma2);
    x_hat = (r-lambda*sqrt(sigma2)*exp(1j*angle(r))).*b;
    d = ((1-lambda*sqrt(sigma2)./(2*abs(r))).*b);
    dc = -(lambda*sqrt(sigma2)*(r.^2)./(2*abs(r).^3).*b);
    dx = d+dc;
    dy = 1j*(d-dc);
end

