%-------------------------------------
% Lecture 11. 
% November 18, 2010.
% Compressed Sensing
% go to l1 magic directory /l1magic-1.1

%----------------------------------------
% Basis pursuit on underdetermined system

path(path, './Optimization');
path(path, './Data');

% 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;

% ground truth: random +/- 1 signal
x = zeros(p,1);
q = randperm(p);
x(q(1:T)) = sign(randn(T,1));
figure; plot(x)

% measurement matrix is given
A = randn(n,p);
figure; imagesc(A); colorbar


% This is what we observe
b = A*x;
figure; plot(b);

% LSE 
x_lse=inv(A'*A)*A'*b;
figure; plot(x_lse)


% algorithm requires an initial guess. 
% simply put anything here
x0 = A'*b;
%x0 = inv(A'*A)*A'*b;

% solve basis pursuit

x_sp = l1eq_pd(x0, A, [], b, 1e-3);
figure; plot(x_sp)

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

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


% noisy observations
sigma = 0.005;
e = sigma*randn(n,1);
b = A*x + e;
figure; plot(b);

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

epsilon = 1
x_sp = l1qc_logbarrier(x0, A, [], b, epsilon, 1e-3);
figure; plot(x_sp)


%----------------------------
% brain network application

% Bunch of file information. This is the data set used in
% kim.2011.ISBI.pdf

c=0;
p=1;

id=[
2	c	13	0
4	c	11	0
5	c	13	1
9	p	10	0
12	p	10	0
13	p	11	1
14	p	12	1
19	p	13	0
20	p	9	1
22	p	12	0
23	p	14	0
24	p	13	0
26	c	10	0
27	c	12	0
30	c	12	0
31	p	12	1
32	p	9	0
33	p	10	1
36	c	11	1
37	c	10	0
38	c	13	0
39	c	14	0
40	c	11	1
41	c	11	1
42	c	9	0
43	c	14	1
44	p	14	1
45	p	10	0
46	c	10	1
47	c	13	1
48	c	13	1
49	c	13	1
52	c	10	1
58	p	9	0
64	c	12	1
82	p	9	1
93	p	10	0
96	c	9	0
97	p	13	1
98	p	12	0
102	p	10	0
103	c	10	0
106	p	12	0
109	p	9	0
110	p	14	0
112	p	9	1
113	p	11	0
119	c	9	0
121	c	9	1
127	c	9	1
131	c	12	1
132	p	10	0
133	c	14	1
134	c	12	1
135	c	14	1
136	c	11	1
140	c	12	1
142	p	10	0
149	p	13	1
150	p	11	0
151	p	14	1
153	c	12	0
158	p	13	1
159	p	10	1
162	c	11	1]

% subject number:
n_subject=size(id,1)

% group variable:
group=id(:,2);

% subject numbers in each group
n1=sum(~group) % control = 33
n2=sum(group)   % neglected children =  32 


% variable length subject ids are handled as
for i=1:n_subject
    str=num2str(id(i,1));
    switch length(str)
        case 1 
            str=strcat('00',str);
        case 2 
            str=strcat('0',str);
    end;
    id_num(i,:)=str;
end;

% directory and file names associated with subject id
% this line has to be changed for different computers

directory='/Users/titipas/Desktop/jacobian/displacement/0';
str_x='_IR+acpc_corrected_finalfullinversewarpxvec.hdr';
str_y='_IR+acpc_corrected_finalfullinversewarpyvec.hdr';
str_z='_IR+acpc_corrected_finalfullinversewarpzvec.hdr';

fname_x=strcat(directory, id_num, str_x);
fname_y=strcat(directory, id_num, str_y);
fname_z=strcat(directory, id_num, str_z);


%--------------------------------------------------------
% Reading Jacobian determinant in white matter only


% Whitematter template has values 0 to 1 (white matter density).
% It has to be thresholded to obtain the binary white matter region
% Let's threshold at 0.8 and above and assign it to vol_re 


% white matter template ( vol_tem )
fname_template='/Users/titipas/Desktop/jacobian/template/templatewm.hdr';
file_name=strcat(fname_template)
image_info = analyze75info(file_name)
vol_tem = analyze75read(image_info);

surf_tem=isosurface(vol_tem);
surf_tem= reducepatch(surf_tem,0.1);
figure;figure_patch(surf_tem,[0.9 0.9 0.9],0.2);hold on;
axis off
%view([-90 90])
set(gcf,'Color','w');
hold on;
  
%------------------------------
% uniform node selection in white matter
%
% threshold white matter density at 0.8
[temp, vol_re]=threshold_image(vol_tem, 0.8);
figure_slices(vol_re, [100,100,80]) 

% identify all voxel coordinates in vol_re > 0.8

d=size(vol_re);
coord=[];
for k=1:d(3)
    [i,j] = find(vol_re(:,:,k) >0.8);
    temp=[i j ones(length(i),1)*k];
    coord=[coord ; temp];
end;

coord_sub=[];
% subsampling at 2000 gives 169 nodes
for i=1:2000:size(coord,1)
    coord_sub=[coord_sub ; coord(i,:)];
end;

%  node display
figure;figure_patch(surf_tem,[0.9 0.9 0.9],0.2);hold on;
plot3(coord_sub(:,2), coord_sub(:,1), coord_sub(:,3),'.r');


%-------------------------------
% Approximate Jacobian determinants over 169 nodes
% This has been computed for you.
%download NIP.2010/matlab/jacobian169nodes.mat

load jacobian169nodes.mat
% JJ  =  subject x nodes
% JJ is 65 x 169 matrix consisting of jacobian determinant at
% 169 uniformaly sampled nodes in white matter for 65 subjects

% Measurements are in X

X=JJ(~group,:);  % Let's do it for only control
n=n1;
p=size(JJ,2);

figure; imagesc(X); colorbar

%------------------
% positive correlation

corr= corrcoef(X);
[temp, corrplus]=threshold_image(corr, 0);
figure; imagesc(corrplus); colorbar



A = spalloc(n*p,p^2,n*(p-1)*p);
%Sparse matrix. Otherwise it will blow up your memory
%A=zeros(n*p,p^2);
for i=1:p
    temp = X;
    temp(:,i)=0;
    A((i-1)*n+1:i*n, (i-1)*p+1:i*p)=temp;
end;

y=[];
for i=1:p
    y=[y; X(:,i)];
end;

% partial correlation via LSE

%beta_lse = A'*inv(A*A')*y;

beta_lse = inv(A'*A)*A'*y;
r=y-A*beta_lse;
res=reshape(r,n,p);
pcorr= corrcoef(res);
[temp, pcorrplus]=threshold_image(pcorr, 0);
figure; imagesc(pcorrplus); colorbar
set(gcf,'Color','w') 

%print -r300 -dtiff partial_LSE

figure; hist(reshape(pcorr,p^2,1),100);

%--------------------------------
% Basis pursuit via L1-magic

% initial guess

x0 = A'*pinv(A*A')*y;
xp = l1eq_pd(x0, A, [], y, 1e-3);

beta_sp=reshape(xp,p,p);
figure;imagesc(beta_sp);

%-----------------------
% download l1_ls toolbox
%LASSO

beta=[];
for i=0.9:-0.1:0.1
          % regularization parameter
    lambda=i
    rel_tol = 0.01;     % relative target duality gap
    tic
    [x,status]=l1_ls(A,y,lambda,rel_tol);
    toc
    beta=[beta x];
end;

save LASSOresult_0.9to0.1.mat beta

for i=1:9
    beta_sp=reshape(beta(:,i),p,p);
    figure;imagesc(beta_sp);
end;





%------------------------
% partial correlation 
% residual vectors


for i=1:10
    beta_sp=beta(:,i);
    r=y-A*beta_sp;
    res=reshape(r,n,p);
    pcorr= corrcoef(res);
    [temp, pcorrplus]=threshold_image(pcorr, 0);
    figure; imagesc(pcorrplus); colorbar
    set(gcf,'Color','w')
end;



print -r300 -dtiff partial_sparse100

figure; hist(reshape(pcorr,p^2,1),100);

%---------------------------------------------
% graph construction using partial correalation

figure;figure_patch(surf_tem,[0.9 0.9 0.9],0.2);hold on;
plot3(coord_sub(:,2), coord_sub(:,1), coord_sub(:,3),'.r');

vertices(:,1)=coord_sub(:,2);
vertices(:,2)=coord_sub(:,1);
vertices(:,3)=coord_sub(:,3);

for i=1:10
    beta_sp=beta(:,i);
    r=y-A*beta_sp;
    res=reshape(r,n,p);
    pcorr= corrcoef(res);
    [temp, pcorrplus]=threshold_image(pcorr, 0);

    pcorrplus= pcorrplus- eye(p);
    %figure; imagesc(pcorrplus); colorbar

    adj=reshape(pcorrplus,p^2,1);
    ind=find(adj<0.4);
    adj(ind)=0;
    adj=reshape(adj,p,p);
    %figure; imagesc(adj); colorbar

    f=sum(pcorrplus);

    figure;figure_patch(surf_tem,[0.99 0.99 0.99],0.1);hold on;
    view([90 0])
    figure_graph_color(adj,vertices, f);
    print('-dtiff', '-r300',strcat('sparsegraph', num2str(i),'.tif'))
end;





