Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/GSP_Script_ISBI_2017.m
f85bf7ad46f1master

Last Change
Raw File

GSP_Script_ISBI_2017.m
View Options

%% 1. Loading of required data
addpath(genpath('Utilities'));
cd('..');
cd('Data/100307');
load('X');
load('W');
cd('..');
cd('..');
cd('Code');
% Defines the HRF
A.dt = 0.72;
A.name = 'hrf';
HRF = spm_get_bf(A);
PathToReg = '/Volumes/Lorena/GSP_HCP_Test_Subject/100307/Task/Events';
[time,Regressor] = GSP_Paradigm(PathToReg,10);
Regressor(Regressor == 5) = 1;
Regressor(Regressor == 6) = 2;
Regressor(Regressor == 7) = 3;
Regressor(Regressor == 8) = 4;
% Convolves each non-constant regressor with the HRF
U.u = Regressor';
U.name = {'Whatever'}; % U needs a name, but the string has no meaning
Regressor_real = spm_Volterra(U, HRF.bf);
%% 2. Analysis at the level of independent Slepian vectors
% To find the interesting nodes, I do a GLM beforehand
[SIGN,t] = GSP_MakeGLM({X},[ones(395,1),Regressor'],[0;1]);
% Indices of the regions to retain
idx_GLM = find(SIGN == 1);
% Bandwidth to use for the Slepian design; for now we go for something
% large enough
BW = 838;
% Number of regions
n_ROI = size(W,1);
% Computation of Slepian vectors
[Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,idx_GLM,{X});
PlotBrainGraph(zeros(n_ROI,n_ROI),abs(Slep(:,sidx)),Slep(:,sidx),CodeBook,0.5,max(abs(Slep(:,sidx))),...
max(abs(Slep(:,sidx))),2,1,'jet','jet',1,...
0.6,[]);
X_SLEP{1} = Slep' * X;
% Computation of Slepian vectors using a random set of nodes as a selection
% idx_random = randperm(n_ROI);
% [Slep_rand,Slep_coefs_rand,Xi_rand,Mu_rand] = GSP_ComputeSlepians(W,BW,idx_random(1:length(idx_GLM)),{X});
TB1 = figure;
set(gcf,'color','w');
imagesc(X_SLEP{1});
caxis([-100 100]);
% TB2 = figure;
% set(gcf,'color','w');
% imagesc(X_slep_rand{1});
% caxis([-100 100]);
% Options for the decomposition: we have a real symmetric matrix at
% hand
opts.issym=1;
opts.isreal=1;
opts.maxit=2500;
opts.disp=0;
% U contains the eigenvectors, while Lambda is a diagonal matrix
% containing the eigenvalues of the decomposition. They are sorted in
% increasing order of amplitude
[V,Lambda] = eigs(diag(sum(W,2))-W,n_ROI,'SA',opts);
[V_loc,Lambda_loc] = eigs(diag(sum(W(idx_GLM,idx_GLM),2))-W(idx_GLM,idx_GLM),length(idx_GLM),'SA',opts);
X_GFT{1} = V' * X;
X_GFT_loc{1} = V_loc' * X(idx_GLM,:);
TB3 = figure;
set(gcf,'color','w');
imagesc(X_GFT{1});
caxis([-100 100]);
% Number of null realisations to create
n_null = 100;
% Number of subjects (for now, 1)
n_subjects = 1;
% Generation of my null data
for s = 1:n_subjects
disp(num2str(s));
for n = 1:n_null
% Filter that we will use for null data generation
Phi_2 = zeros(n_ROI,n_ROI);
for i = 1:n_ROI
tmp = randn;
if tmp > 0
Phi_2(i,i) = 1;
else
Phi_2(i,i) = -1;
end
end
% Filter that we will use for null data generation (Slepians)
Phi_3 = zeros(BW,BW);
for i = 1:BW
tmp = randn;
if tmp > 0
Phi_3(i,i) = 1;
else
Phi_3(i,i) = -1;
end
end
XNull_GFT{s}(:,:,n) = V*Phi_2*V'*X;
XNull_SLEP{s}(:,:,n) = Slep*Phi_3*Slep'*X;
end
TNull_GFT{s} = prctile(XNull_GFT{s},[2.5 97.5],3);
TNull_SLEP{s} = prctile(XNull_SLEP{s},[2.5 97.5],3);
EXC_GFT{s} = zeros(n_ROI,395);
EXC_GFT{s}(X > TNull_GFT{1}(:,:,2)) = 1;
EXC_GFT{s}(X < TNull_GFT{1}(:,:,1)) = -1;
EXC_SLEP{s} = zeros(n_ROI,395);
EXC_SLEP{s}(X > TNull_SLEP{1}(:,:,2)) = 1;
EXC_SLEP{s}(X < TNull_SLEP{1}(:,:,1)) = -1;
end
TB4 = figure;
set(gcf,'color','w');
imagesc(EXC_GFT{1});
TB5 = figure;
set(gcf,'color','w');
imagesc(EXC_SLEP{1});
Slep = -1 * Slep;
% Cutoff on Xi to retain Slepians
epsilon = 1e-2;
% Computation of concentration
for t = 1:395
C(t) = ((Slep(:,diag(Xi)>epsilon)'*X(:,t)))' * Mu(diag(Xi)>epsilon)';
for n = 1:n_null
CN(t,n) = ((Slep(:,diag(Xi)>epsilon)'*squeeze(XNull_SLEP{1}(:,t,n))))' * Mu(diag(Xi)>epsilon)';
end
end
figure;
set(gcf,'color','w');
plot_ci(time',[prctile(CN,2.5,2),prctile(CN,97.5,2)],'PatchColor',...
[0.2 0.2 0.2], 'PatchAlpha', 0.9, 'MainLineWidth', 2);
hold on;
plot(time,C,'r','LineWidth',2);
plot(time,600*(Regressor_real-mean(Regressor_real)),'color',[1 0.8 0.2]);
plot(time,600*(Regressor-mean(Regressor)),'color',[0.8 0.4 0]);
xlabel('Time [s]');
ylabel('Concentration within visual system');
xlim([time(1) time(end)]);
% Going for trajectory plotting now
% Option 1: non linear dimensonality reduction using graphs
W1 = construct_knn_graph(X_GFT{1}',20,'euclidean');
W2 = construct_knn_graph(X_SLEP{1}',20,'euclidean');
W3 = construct_knn_graph(X_SLEP{1}(diag(Xi)>epsilon,:)',20,'euclidean');
W4 = construct_knn_graph(X_GFT_loc{1}',20,'euclidean');
[X1,Y1,Z1] = compute_non_linear_dim_reduction(W1);
[X2,Y2,Z2] = compute_non_linear_dim_reduction(W2);
[X3,Y3,Z3] = compute_non_linear_dim_reduction(W3);
[X4,Y4,Z4] = compute_non_linear_dim_reduction(W4);
CV = Regressor_real;
figure;
scatter3(X1,Y1,Z1,10,CV,'filled');
colormap(jet);
caxis([min(Regressor_real),max(Regressor_real)]);
CV = Regressor_real;
figure;
scatter3(X3,Y3,Z3,10,CV,'filled');
colormap(jet);
caxis([min(Regressor_real),max(Regressor_real)]);
CV = Regressor_real;
figure;
scatter3(X2,Y2,Z2,10,CV,'filled');
colormap(jet);
caxis([min(Regressor_real),max(Regressor_real)]);
CV = Regressor_real;
figure;
scatter3(X4,Y4,Z4,10,CV,'filled');
colormap(jet);
caxis([min(Regressor_real),max(Regressor_real)]);
% Option 2: PCA (linear reduction of dimensionality)
% Data has size N x D
Data = X_SLEP{1}(diag(Xi)>epsilon,:)';
% Demeaning
Mu_data = mean(Data,1);
Data = Data - repmat(Mu_data,[size(Data,1),1]);
% Computation of covariance matrix
C_Data = 1/size(Data,1) * (Data' * Data);
% Computation of its eigenvectors
[V_PCA,Lambda_PCA] = eigs(C_Data,size(C_Data,1),'LA',opts);
c4science · Help