Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/GSP_Script_ISBI_2017.m
9f1d070fc71amaster

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 = 300;
% Number of regions
n_ROI = size(W,1);
IDX_OI = idx_GLM(idx_GLM > 820);
vec_idx = zeros(n_ROI,1);
vec_idx(idx_GLM(idx_GLM > 750)) = 1;
% PlotBrainGraph(zeros(n_ROI,n_ROI),abs(vec_idx),vec_idx,CodeBook,0.5,max(vec_idx),...
% max(abs(vec_idx)),2,1,'jet','jet',1,...
% 0.6,[]);
% Computation of Slepian vectors
[Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,IDX_OI,{X});
%
% PlotBrainGraph(zeros(n_ROI,n_ROI),abs(vec_idx),vec_idx,CodeBook,0.5,max(vec_idx),...
% max(abs(vec_idx)),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_OI,IDX_OI),2))-W(IDX_OI,IDX_OI),length(IDX_OI),'SA',opts);
X_GFT{1} = V' * X;
X_GFT_loc{1} = V_loc' * X(IDX_OI,:);
% 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;
%
% TEST = 1;
%
% % Computation of concentration
% for t = 1:395
%
% C(t) = ((Slep'*X(:,t)))' * Mu';
%
% for n = 1:n_null
% CN(t,n) = ((Slep'*squeeze(XNull_SLEP{1}(:,t,n))))' * Mu';
% 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
f1 = 2;
n_freq = sum(diag(Xi)>epsilon);
n_freq2 = sum(Mu>epsilon);
% Option 1: non linear dimensonality reduction using graphs
W1 = construct_knn_graph(X_GFT{1}(1:n_freq2-(f1-1),:)',20,'euclidean');
%W2 = construct_knn_graph(X_SLEP{1}(1:end,:)',20,'euclidean');
idx1 = ones(size(X_SLEP{1},1),1);
idx1(Mu <= epsilon) = 0;
idx1(1:f1-1) = 0;
W3 = construct_knn_graph(X_SLEP{1}(f1:sum(Mu>epsilon),:)',20,'euclidean');
W4 = construct_knn_graph(X_GFT_loc{1}(1:n_freq2-(f1-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(X2,Y2,Z2,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(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