Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/GSP_Script_ISBI_2017_Simulated_TB.m
275d2d988655master

Last Change
Raw File

GSP_Script_ISBI_2017_Simulated_TB.m
View Options

%% 1. Generation of the graph
% Number of regions (to be close to HCP)
n_ROI = 500;
% According to Power et al. 2011 (Neuron) for the number of communities;
% else, we pick a small enough world density and reasonably high community
% numbers
param.Nc = 10;
param.min_comm = 20;
param.world_density = 1/n_ROI;
% Creates the graph, plots it and retrieves W and L
G = gsp_community(n_ROI,param);
gsp_plot_graph(G);
W = G.W;
L = G.L;
% We want to pick the largest community as our selected, ground truth set;
% for this, we determine the index of each node with spectral clustering,
% and pick the group with most entries
opts.issym=1;
opts.isreal=1;
opts.maxit=2500;
opts.disp=0;
[V,Lambda] = eigs(L,n_ROI,'SA',opts);
idx = kmeans(V(:,2:7),param.Nc);
n_elements = hist(idx,1:param.Nc);
tmp = find(n_elements == max(n_elements)); tmp = tmp(1);
clear V
clear Lambda
% S contains the indices of the nodes from our representation that are part
% of our ground truth community of interest
S = find(idx == tmp);
S_rest = find(idx ~= tmp);
n_S = length(S);
%% 2. Generation of the signals on the graph
% Number of time points
n_TP = 900;
% Mean of the time courses
Mu = zeros(n_ROI,1);
% Standard deviation of the noise
Std = 0.5*ones(n_ROI,1);
% Noise level of our time courses
X = repmat(Mu,1,n_TP);
% Paradigm at hand: we say we will have two paradigms, equivalently
% segregated within the nodes from the selected community and with two
% different timings
time = 1:1:n_TP;
Par1 = zeros(1,n_TP);
Par1(time > 100 & time < 200) = 1;
Par1(time > 300 & time < 400) = 1;
Par1(time > 500 & time < 600) = 1;
Par1(time > 700 & time < 800) = 1;
Par2 = zeros(1,n_TP);
Par2(time > 120 & time < 230) = 1;
%Par2(time > 300 & time < 400) = 1;
Par2(time > 520 & time < 630) = 1;
%Par2(time > 700 & time < 800) = 1;
TB1 = figure;
set(gcf,'color','w');
plot(time,Par1,'r','LineWidth',2);
hold on;
plot(time,Par2,'b','LineWidth',2);
xlabel('Time [s]');
ylabel('Paradigm');
legend('Type 1','Type 2');
type_vec = randperm(n_S);
S1 = S(type_vec(1:floor(n_S/2)));
S2 = S(type_vec(floor(n_S/2)+1:end));
X(S1,:) = X(S1,:) + repmat(Par1,length(S1),1);
X(S2,:) = X(S2,:) + repmat(Par2,length(S2),1);
X_noisy = X + repmat(Std,1,n_TP).*randn(n_ROI,n_TP);
TB2 = figure;
set(gcf,'color','w');
imagesc(X);
TB3 = figure;
set(gcf,'color','w');
imagesc(X_noisy);
% Error between X and X_noisy
MSE_start = GSP_ComputeMSE(X_noisy,X);
%% 3. Reconstruction using GFT
% V contains the eigenvectors, Lambda contains the eigenvalues (vector
% format)
[V,Lambda] = eigs(diag(sum(W,2))-W,n_ROI,'SA',opts);
Lambda = diag(Lambda);
Lambda = Lambda - Lambda(1);
C_GFT = V' * X_noisy;
idx_GFT = 1;
Lambda_cutoff = linspace(0,max(Lambda),300);
for lam = Lambda_cutoff
H = zeros(n_ROI,n_ROI);
for i = 1:n_ROI
if Lambda(i) <= lam
H(i,i) = 1;
else
H(i,i) = 0;
end
end
X_GFT = V * H * C_GFT;
MSE_GFT(:,idx_GFT) = GSP_ComputeMSE(X_GFT,X);
idx_GFT = idx_GFT + 1;
end
TB4 = figure;
set(gcf,'color','w');
plot(Lambda_cutoff,mean(MSE_GFT(S,:),1),'color',[0 0.2 0.8],'LineWidth',2);
xlabel('Cutoff lambda');
ylabel('MSE');
%% 4. Reconstruction using a localised GFT (perfect set of nbodes)
% V contains the eigenvectors, Lambda contains the eigenvalues (vector
% format)
[V_loc,Lambda_loc] = eigs(diag(sum(W(S,S),2))-W(S,S),floor(n_S),'SA',opts);
Lambda_loc = diag(Lambda_loc);
Lambda_loc = Lambda_loc - min(Lambda_loc);
C_LGFT = V_loc' * X_noisy(S,:);
idx_LGFT = 1;
Lambda_cutoff = linspace(0,max(Lambda_loc),300);
MSE_LGFT = nan(n_ROI,length(Lambda_cutoff));
for lam = Lambda_cutoff
H = zeros(n_S,n_S);
for i = 1:n_S
if Lambda_loc(i) <= lam
H(i,i) = 1;
else
H(i,i) = 0;
end
end
X_LGFT = V_loc * H * C_LGFT;
tmp_MSE = nan(n_ROI,1);
tmp_MSE(S) = GSP_ComputeMSE(X_LGFT,X(S,:));
% Error made within the selected area of interest
MSE_LGFT_S(idx_LGFT) = nanmean(tmp_MSE(S));
idx_LGFT = idx_LGFT + 1;
end
TB5 = figure;
set(gcf,'color','w');
plot(Lambda_cutoff,MSE_LGFT_S,'LineWidth',2,'Color','k');
%% 5. Reconstruction using a localised GFT (imperfect set of nbodes)
% Number of iterations
n_iter = 50;
% Fraction of nodes to keep
f_nodes = 0.4;
for n = 1:n_iter
disp(['Iteration ',num2str(n),'...']);
% Set of nodes selected (always imperfect)
tmp = S;
tmp = tmp(randperm(n_S));
LGFT_set = tmp(1:floor(n_S*f_nodes));
% V contains the eigenvectors, Lambda contains the eigenvalues (vector
% format)
[V_loc,Lambda_loc] = eigs(diag(sum(W(LGFT_set,LGFT_set),2))-W(LGFT_set,LGFT_set),floor(n_S*f_nodes),'SA',opts);
Lambda_loc = diag(Lambda_loc);
Lambda_loc = Lambda_loc - min(Lambda_loc);
C_LGFT = V_loc' * X_noisy(LGFT_set,:);
idx_LGFT = 1;
Lambda_cutoff = linspace(0,max(Lambda_loc),300);
MSE_LGFT = nan(n_ROI,length(Lambda_cutoff));
for lam = Lambda_cutoff
H = zeros(length(LGFT_set),length(LGFT_set));
for i = 1:length(LGFT_set)
if Lambda_loc(i) <= lam
H(i,i) = 1;
else
H(i,i) = 0;
end
end
X_LGFT = V_loc * H * C_LGFT;
tmp_MSE = nan(n_ROI,1);
tmp_MSE(LGFT_set) = GSP_ComputeMSE(X_LGFT,X(LGFT_set,:));
% Error made within the selected area of interest
MSE_LGFT_NS(idx_LGFT,n) = nanmean(tmp_MSE(LGFT_set));
idx_LGFT = idx_LGFT + 1;
end
end
TB6 = figure;
set(gcf,'color','w');
plot_ci(1:300,[mean(MSE_LGFT_NS,2),mean(MSE_LGFT_NS,2)-std(MSE_LGFT_NS,[],2),mean(MSE_LGFT_NS,2)+std(MSE_LGFT_NS,[],2)],...
'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ...
'MainLineColor', [1 0.2 0], 'LineWidth', 1, 'LineStyle','-', ...
'LineColor', 'none');
%% 5. Reconstruction using the Slepians (perfect set of nodes, only local frequency filtering)
% Limit to determine that a vector is 'concentrated' (a Slepian is
% concentrated if Xi > epsilon)
epsilon = 1e-5;
is_epsilon = 1;
% Index for the bandwidth
idx_BW = 1;
% Bandwidth values to try
BandWidth = [5 10 15 20 30 40 100 200 n_ROI];
% I want to compare bandwidths
for BW = BandWidth
% For each bandwidth, I compute my Slepians
[Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,S,{X});
% I get my eigenvalues
Xi = diag(Xi);
Xi = Xi - min(Xi);
% Coefficients in the Slepian domain
C_SLEP = Slep' * X_noisy;
% index for the cutoff frequency
idx_SLEP = 1;
% Frequency range to probe
Lambda_cutoff = linspace(0,max(Xi),300);
% Will contain the error
MSE_SLEP = nan(n_ROI,length(Lambda_cutoff),length(BandWidth));
% For all possible cutoff values
for lam = Lambda_cutoff
% Will contain only my coefficients of interest
H = zeros(BW,BW);
% Fills if local AND low frequency
for i = 1:BW
if ~is_epsilon
if Xi(i) <= lam %&& Xi(i) > epsilon
H(i,i) = 1;
end
else
if Xi(i) <= lam && Mu(i) > epsilon
H(i,i) = 1;
end
end
end
% Reconstructed signal
X_SLEP = Slep * H * C_SLEP;
% Error for the present iteration (all nodes)
tmp_MSE = GSP_ComputeMSE(X_SLEP,X);
% Error to use in the end
MSE_SLEP_S1(idx_SLEP,idx_BW) = nanmean(tmp_MSE(S));
idx_SLEP = idx_SLEP + 1;
end
idx_BW = idx_BW + 1;
end
TB7 = figure;
imagesc(MSE_SLEP_S1);
%% 6. Reconstruction using the Slepians (imperfect set of nodes)
% Limit to determine that a vector is 'concentrated' (a Slepian is
% concentrated if Xi > epsilon)
epsilon = 1e-5;
is_epsilon = 1;
% Index for the bandwidth
idx_BW = 1;
% Bandwidth values to try
BandWidth = [5 10 15 20 30 40 100 200 n_ROI];
% Number of iterations
n_iter = 10;
% Feaction of nodes to keep
f_nodes = 0.4;
% I want to compare bandwidths
for BW = BandWidth
for n = 1:n_iter
disp(['Iteration ',num2str(n),'...']);
% Set of nodes selected (always imperfect)
tmp = S;
tmp = tmp(randperm(n_S));
SLEP_set = tmp(1:floor(n_S*f_nodes));
% For each bandwidth, I compute my Slepians
[Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,SLEP_set,{X});
% I get my eigenvalues
Xi = diag(Xi);
Xi = Xi - min(Xi);
% Coefficients in the Slepian domain
C_SLEP = Slep' * X_noisy;
% index for the cutoff frequency
idx_SLEP = 1;
% Frequency range to probe
Lambda_cutoff = linspace(0,max(Xi),300);
% Will contain the error
MSE_SLEP = nan(n_ROI,length(Lambda_cutoff),length(10:10:n_ROI));
% For all possible cutoff values
for lam = Lambda_cutoff
% Will contain only my coefficients of interest
H = zeros(BW,BW);
% Fills if local AND low frequency
for i = 1:BW
if ~is_epsilon
if Xi(i) <= lam %&& Xi(i) > epsilon
H(i,i) = 1;
end
else
if Xi(i) <= lam && Mu(i) > epsilon
H(i,i) = 1;
end
end
end
% Reconstructed signal
X_SLEP = Slep * H * C_SLEP;
% Error for the present iteration (all nodes)
tmp_MSE = GSP_ComputeMSE(X_SLEP,X);
% Error to use in the end
Error_Slepians(idx_SLEP,idx_BW,n) = nanmean(tmp_MSE(S));
idx_SLEP = idx_SLEP + 1;
end
end
idx_BW = idx_BW + 1;
end
TB8 = figure;
imagesc(mean(Error_Slepians,3));
c4science · Help