function xest = cspocs_l1(Phi, aPhi, l1val, y, niter)
% Syntax :
%  xest = cspocs_l1(Phi, aPhi, l1val, y, niter)
%
% Description : Testing Candes and Romberg POCS (alternate Projection Onto
% Convex Sets) for CS recovery, aka recovering of x from y = Phi*x when
% size(Phi,1) < size(Phi,2), using l1 criterion.
% 
% In :
% * Phi, aPhi : Measurement matrix and its reconstruction
% * l1val : the l1 norm of the intial signal x, i.e. norm(x,1)
% * niter : maximun number of alternate projections.
%
% Out :
% * xest : the recovered (or estimated) signal x
%
% Author of this mfile : L. Jacques, LTS2/EPFL, 2008.
%    
% Reference : 
%    Algo described in "Practical Signal Recovery from Random Projections",
%    Emmanuel Candès and Justin Romberg
% Example :
% >> N=128; K=20;
% >> x=[rand(1,K) zeros(1,N-K)]'; x=x(randperm(128));    
% >> m=floor(3.5*K);Phi=randn(m,N)/sqrt(m);aPhi=Phi';          
% >> y=Phi*x;
% >> figure; plot(x);
% >> nx=cspocs_l1(Phi,aPhi,norm(x,1),y,10000);      
% Stop at n=2157, score=9.987234e-11 ...
% Final score = 9.987234e-11 ...
% >> hold on; plot(nx,'ro');
% 
% Remark: This algo seems highly sensitive to the a priori l1 norm of x. To
% convince yourself of that, repeat the same experiment as above with
% norm(x,1)*0.99 and norm(x,1)*1.01 and observe the
% results. Recovery/NoRecovery transition points seems located around 3.2K
%
% This script/program is released under the Commons Creative Licence
% with Attribution Non-commercial Share Alike (by-nc-sa)
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% Short Disclaimer: this script is for educational purpose only.
% Longer Disclaimer see  http://igorcarron.googlepages.com/disclaimer
    
    [m,N] = size(Phi);
    
    pPhi = pinv(Phi);
    
    %% First projector (on the hyperplane Phi*x=y)
    P = @(u) u + pPhi*(y - Phi*u);
    
    %% Second projector on the l1 ball of radius l1val
    H = @(u, l1val) sft_th(u, l1toth(u,l1val));
    
    %% Intializations
    xest = pPhi*y;
    Tol = 1e-10;
    score = 0;
    
    for n = 1:niter;
        xest = P(xest);
        xest = H(xest, l1val);
        
        oldscore = score;
        score = norm(Phi*xest - y) / norm(xest);
        
        if (score < Tol)
            fprintf('Tolerance reached. Stop at n=%i\n',n);
            break
        end
    end

    fprintf('Final score: ||y - Phi*xest|| = %e\n', score);
    
    
function out = l1toth(x,l1val)
% (Internal function)
% Obtain the threshold level corresponding to the one such that the l1
% norm of the corresponding soft thresholded signal is lesser than l1val
    N = length(x);
    k = (1:N)';
    ax = abs(x);
    
    xs = sort(ax, 'descend');
    cxs = cumsum(xs);
    pos = (cxs - k.*xs) < l1val;
    out = xs(pos);
    out = out(end);
    
function out = sft_th(x, gamma)
% (Internal function)
% Soft Thresholding
    out = sign(x) .* (abs(x) - gamma) .* (abs(x) >= gamma);