Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/GSP_GenerateData.m
8a186cd82567master

Last Change
Raw File

GSP_GenerateData.m
View Options

%% This function creates a graph, and a signal on this graph, according to
% a set of parameters
%
% Inputs:
% - n_ROI is the number of regions in the graph
% - Nc is the number of communities to include in the graph
% - min_comm is the minimal number of nodes belonging to a community
% - WD is the world density (typically 1/n_ROI)
% - n_TP is the number of time points to have
% - is_plot defines whether to show plots (1) or not (0)
% - Std is the level of noise (larger values mean that the noise can take
% larger values; i.e., arises from a distribution with larger Std)
% - n_rn is the amount of nodes that co-activate as a noise pattern
% - c_fac defines the density of the smaller scale patterns: a larger
% number means more condensed patterns (use c_fac = 4)
% - Type is the type of signal to create: only two subsets of nodes with a
% paradigm (1), or along with other communities activity (2), or this and
% also smaller scale interactions (3)
function [W,X,X_noisy,MSE_start,S,n_S,S_rest,G] = GSP_GenerateData(n_ROI,Nc,min_comm,WD,n_TP,is_plot,Std,n_rn,c_fac,Type)
%% 1. Generation of the graph itself
% We pick a small enough world density and reasonably high community
% numbers
param.Nc = Nc;
param.min_comm = min_comm;
param.world_density = WD;
% 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:4),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);
idx_others = 1;
for c = 1:(param.Nc)
if c ~= tmp
S_others{idx_others} = find(idx == c);
idx_others = idx_others + 1;
end
end
S_rest = find(idx ~= tmp);
n_S = length(S);
%% 2. Generation of the sets of co-activating graph nodes
% A. Halves of the community of interest
BB1 = zeros(n_ROI,1);
BB2 = zeros(n_ROI,1);
type_vec = randperm(n_S);
S1 = S(type_vec(1:floor(n_S/2)));
S2 = S(type_vec(floor(n_S/2)+1:end));
BB1(S1) = 1;
BB2(S2) = 1;
BB_others = [];
for c = 1:(param.Nc-1)
BB_others(:,c) = zeros(n_ROI,1);
BB_others(S_others{c},c) = 1;
end
% B. Finds the nodes linked to that community (that are outside)
BB = [];
idx_pairs = 1;
% Checking all nodes from the graph
for i = 1:size(W,1)
% If we consider a node part of the community of interest...
if ismember(i,S)
for j = 1:size(W,2)
% If the other node is part of another module and there is a
% connection, we consider
if ismember(j,S_rest) && W(i,j) == 1
BB = [BB,zeros(n_ROI,1)];
BB(i,end) = 1;
BB(j,end) = 1;
Mu_reg(idx_pairs) = 500 + 200*randn;
Sigma_reg(idx_pairs) = 50 + 10*randn;
idx_pairs = idx_pairs + 1;
end
end
end
end
BBbis = [];
tmp_vec = zeros(n_ROI,3);
idx_bis = 1;
for i = 1:size(BB,2)
if mod(i,c_fac) ~= 0
tmp_vec(:,mod(i,c_fac)) = BB(:,i);
else
tmp_vec(:,c_fac) = BB(:,i);
BBbis(:,idx_bis) = sum(tmp_vec,2);
idx_bis = idx_bis + 1;
end
end
BBbis(BBbis > 0) = 1;
if is_plot
% Summary plot
figure;
set(gcf,'color','w');
imagesc([BB1,BB2,BB_others,BB,BBbis]);
end
%% 3. Generation of the timings of activation
% Duration of individual paradigms
Delta = floor(n_TP/(param.Nc+1));
% Creation of the paradigms
time = 1:1:n_TP;
% For both of my region of interest communities
Par1 = zeros(1,n_TP);
Par1(time <= Delta) = 1;
Par2 = zeros(1,n_TP);
Par2(time > Delta & time <= 2*Delta) = 1;
% Other paradigms
for c = 1:(param.Nc-1)
Par_others(c,:) = zeros(1,n_TP);
Par_others(c,time > (c+1)*Delta & time <= (c+2)*Delta) = 1;
end
for i = 1:size(BB,2)
for t = 1:n_TP
Par3(i,t) = normpdf(t,Mu_reg(i),Sigma_reg(i));
end
Par3(i,:) = Par3(i,:)/max(Par3(i,:));
end
if is_plot
figure;
set(gcf,'color','w');
plot(time,Par1,'r');
hold on
plot(time,Par2,'b');
plot(time,Par_others);
%plot(time,Par3);
end
%% 4. Generation of the signal itself
% A. Building blocks
% Initialised at zero
X = zeros(n_ROI,n_TP);
% Adding the subset paradigms
X = X + (BB1 * Par1 + BB2 * Par2);
if Type > 1
% Adding the other modules
for i = 1:size(BB_others,2)
X = X + (BB_others(:,i) * Par_others(i,:));
end
end
if Type > 2
% Adding the smaller scale interactions
for i = 1:size(BBbis,2)
X = X + (BBbis(:,i) * Par3(i,:));
end
end
if is_plot
TB2 = figure;
set(gcf,'color','w');
imagesc(X);
end
% B. Noise
% Construction of the noise and adding it
BB_noise = zeros(n_ROI,n_TP);
for t = 1:n_TP
rvec = randperm(n_ROI);
for i = 1:n_rn
BB_noise(rvec(i),t) = 1;
end
end
X_noisy = X + Std * BB_noise .* randn(n_ROI,n_TP);
if is_plot
TB3 = figure;
set(gcf,'color','w');
imagesc(X_noisy);
end
% Error between X and X_noisy
MSE_start = GSP_ComputeMSE(X_noisy,X);
end
c4science · Help