% Statistical and Computational Methods in Brain Image Analysis
% Chapter 15. Sparse Networks
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu


%----------------------------------------
% Basis pursuit 

% set up paths

path(path, './l1magic-1.1');
path(path, './l1magic-1.1/Optimization');

% Lets solve b= Ax
% b and A are given. We solve for sparse x

% signal length
p = 100;
% number of spikes in the signal
T = 4;
% number of observations to make
n = 50;

% unknown ground truth: random +/- 1 signal
x = zeros(p,1);
q = randperm(p);
x(q(1:T)) = sign(randn(T,1));
figure; subplot(2,2,1); plot(x)
title('Unknown parameter')

% known measurement matrix is given
A = randn(n,p);

% known observation
b = A*x;
subplot(2,2,2); plot(b);
title('Observation')
% need to estimate s using LSE 
x_lse=inv(A'*A)*A'*b;
subplot(2,2,3); plot(x_lse)
title('Least squares estimation')

% Basis pursuit
% algorithm requires an initial guess. 
x0 = inv(A'*A)*A'*b;
x_bp = l1eq_pd(x0, A, [], b, 1e-3);

subplot(2,2,4);  plot(x_bp)
title('Basis pursuit')

print_pdf('CS1')

%--------------------------------
% LASSO

% Lets solve b= Ax + noise
% b and A are given. We solve for sparse x under noise

figure; subplot(2,2,1); plot(x)
title('Unknown parameter')

% noisy observations
sigma = 0.7;
e = sigma*randn(n,1);
subplot(2,2,2); plot(e);
title('Noise')


b1= A*x + e;
subplot(2,2,2);
plot(b,'--r')
hold on; plot(b1);
title('Observation')

% initial guess = min energy
x0 = A'*b1;

% solve basis pursuit
x_sp = l1eq_pd(x0, A, [], b1, 1e-3);

subplot(2,2,3);  plot(x_sp)
title('Basis pursuit')

epsilon = 5
x_lasso = l1qc_logbarrier(x0, A, [], b1, epsilon, 1e-3);
subplot(2,2,4); plot(x_lasso)
title('LASSO')

print_pdf('CS1')

%-----------------------------------
% g-LASSO
% read PDF for detail

n=10
p=100
J=normrnd(0, 1, n, p);
S=J'*J/n;
figure; subplot(2,2,1); imagesc(S); colorbar;
title('covariance for n=10, p=100')

[Theta W] = graphicalLasso(S, 5, 10, 0.1);
subplot(2,2,2); imagesc(W); colorbar;
title('g-Lasso with lambda=5')

[Theta W] = graphicalLasso(S, 0.5, 10, 0.1);
subplot(2,2,3); imagesc(W); colorbar;
title('g-Lasso with lambda=0.5')

[Theta W] = graphicalLasso(S, 0.2, 10, 0.1);
subplot(2,2,4); imagesc(W); colorbar;
title('g-Lasso with lambda=0.05')



