Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/GSP_QuantifyError.m
0e8896a72098master

Last Change
Raw File

GSP_QuantifyError.m
View Options

% This function quantifies how well, compared to a ground truth, different
% approaches fare
%
% Inputs:
% - W is the structural scaffold
% - X is the noisy signal
% - X_GT is the ground truth signal
% - S contains the indices of the nodes that are part of the community of
% interest (on which we focus)
% - n_S is the number of elements in S
% - f_nodes is the fraction of nodes that are sampled, from S, at each
% subset iteration (when assessing local GFT vs Slepians)
% - n_iter is the number of subset trials that are run
% - is_epsilon defines, for Slepians, whether we take all the vectors (0),
% or only the ones that are concentrated within the subset of interest (1)
% - epsilon is the tolerance level past which we say that a vector is not
% concentrated within the chosen subset
% - n_BW is the number of bandwidths to test for Slepians
function [MSE_GFT,MSE_LGFT,MSE_LGFT_SS,MSE_SLEP,MSE_SLEP_SS] = GSP_QuantifyError(W,X_noisy,X_GT,S,n_S,f_nodes,n_iter,epsilon,is_epsilon,n_BW,n_ROI)
n_lambda = 300;
%% 1. Reconstruction using GFT
disp('Running GFT computations...');
% V contains the eigenvectors, Lambda contains the eigenvalues (vector
% format)
opts.issym=1;
opts.isreal=1;
opts.maxit=2500;
opts.disp=0;
[V,Lambda] = eigs(diag(sum(W,2))-W,n_ROI,'SA',opts);
Lambda = diag(Lambda);
Lambda = Lambda - Lambda(1);
% Computing the coefficients (the expression of the signal on the
% graph)
C_GFT = V' * X_noisy;
Lambda_cutoff = linspace(0,max(Lambda),n_lambda);
idx_GFT = 1;
% I gradually increment the number of building blocks that I use
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
% Reconstructed signal
X_GFT = V * H * C_GFT;
% Computation of the error at each node
MSE_GFT(:,idx_GFT) = GSP_ComputeMSE(X_GFT,X_GT);
clear X_GFT
idx_GFT = idx_GFT + 1;
end
%% 2. Reconstruction using a localised GFT (perfect set of nodes)
disp('Running localised GFT computations...');
% 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),n_lambda);
% Initialises at NaN because some nodes will never have any error
% measure
MSE_LGFT = nan(n_ROI,n_lambda);
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;
MSE_LGFT(S,idx_LGFT) = GSP_ComputeMSE(X_LGFT,X_GT(S,:));
idx_LGFT = idx_LGFT + 1;
end
%% 3. Reconstruction using a localised GFT (imperfect set of nodes)
disp('Running localised GFT computations (subset of nodes)...');
% Going to contain the estimation errors for the trials investigated
% here
MSE_LGFT_SS = nan(n_ROI,n_lambda,n_iter);
% In each iteration, I pick a different set of nodes
for n = 1:n_iter
% 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),n_lambda);
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;
MSE_LGFT_SS(LGFT_set,idx_LGFT,n) = GSP_ComputeMSE(X_LGFT,X_GT(LGFT_set,:));
idx_LGFT = idx_LGFT + 1;
end
clear V_loc
clear Lambda_loc
clear Lambda_cutoff
clear LGFT_set
end
%% 4. Reconstruction using the Slepians
disp('Running Slepians computations...');
% Index for the bandwidth
idx_BW = 1;
% Bandwidth values to try (from 5% of the whole set to maximal value)
BandWidth = round(linspace(0.05*n_ROI,n_ROI,n_BW));
% Will contain the error
MSE_SLEP = nan(n_ROI,n_lambda,length(BandWidth));
% I want to compare bandwidths
for BW = BandWidth
% For each bandwidth, I compute my Slepians
[Slep,~,Xi,Mu] = GSP_ComputeSlepians(W,BW,S,{X_noisy});
% I get my eigenvalues
Xi = diag(Xi);
Xi = Xi - min(Xi);
% Coefficients in the Slepian domain
C_SLEP = Slep' * X_noisy;
% Frequency range to probe
Lambda_cutoff = linspace(0,max(Xi),n_lambda);
idx_SLEP = 1;
% 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
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)
MSE_SLEP(:,idx_SLEP,idx_BW) = GSP_ComputeMSE(X_SLEP,X_GT);
idx_SLEP = idx_SLEP + 1;
end
idx_BW = idx_BW + 1;
end
%% 7. Reconstruction using the Slepians (imperfect set of nodes)
% Error...
MSE_SLEP_SS = nan(n_ROI,n_lambda,length(BandWidth),n_iter);
% Bandwidth values to try (from 5% of the whole set to maximal value)
BandWidth = round(linspace(0.05*n_ROI,n_ROI,n_BW));
for n = 1:n_iter
% Set of nodes selected (always imperfect)
tmp = S;
tmp = tmp(randperm(n_S));
SLEP_set = tmp(1:floor(n_S*f_nodes));
idx_BW = 1;
for BW = BandWidth
% For each bandwidth, I compute my Slepians
[Slep,~,Xi,Mu] = GSP_ComputeSlepians(W,BW,SLEP_set,{X_noisy});
% 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),n_lambda);
% 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
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)
MSE_SLEP_SS(:,idx_SLEP,idx_BW,n) = GSP_ComputeMSE(X_SLEP,X_GT);
idx_SLEP = idx_SLEP + 1;
end
idx_BW = idx_BW + 1;
end
end
end
c4science · Help