% L1 regularization for compressive sensing:
% Solving the model
%        min ||Lx||_1, st. Ax = b
% using CVX

clear, close all
rng(2)

m = 400;
n = 2*m;
r = 4;

k = ceil(m/r);
A = randn(m,n);
fprintf('\n[m n k] = [%i %i %i]\n',m,n,k)

snr = 50;
noise = randn(m,1);
noise = 10^(-snr/20)/norm(noise)*noise;

range = 50;
p = randperm(n);
xo = zeros(n,1);
xo(p(1:k))= range*randn(k,1);

D = cell(3,1);
e = ones(n,1);
D{1} = spdiags(e,0,n,n);
D{2} = spdiags([-e e],-1:0,n-1,n);
eqs = ' ================= ';
warning off

tic
x0 = xo; figure
for p = 1:2
    b = A*x0;
    b = b + norm(b)*noise;
    L = D{p};
    t0 = tic;
    
    % Using CVX to solve: min ||Lx||_1, st. Ax = b
    %---------------------------------------------
    % write your code here
    %---------------------------------------------
    
    t1 = toc(t0);
    relerr = norm(x-x0)/norm(x0);
    if p==1, fprintf(['\n' eqs ' CVX ' eqs '\n']); end
    fprintf('p = %1i: rel_err: %e  time: %f\n',p,relerr,t1)
    subplot(210+p); plot(1:n,x0,'b-',1:n,x,'r.');
    axis auto; shg
    title('CVX Results')
    x0 = cumsum(x0);
end

fprintf('\n----------------------------------------\n')
toc; disp(datestr(now))
fprintf('----------------------------------------\n\n')
