diff --git a/Code/GSP_Script_ISBI_2017_Simulated.m b/Code/GSP_Script_ISBI_2017_Simulated.m index c08b08b..88ef157 100644 --- a/Code/GSP_Script_ISBI_2017_Simulated.m +++ b/Code/GSP_Script_ISBI_2017_Simulated.m @@ -1,298 +1,313 @@ %% 1. Generation of the graph % Number of regions (to be close to HCP) n_ROI = 200; % 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 = 10; 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 = 1.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 + 1e-3; C_GFT = V' * X_noisy; idx_GFT = 1; for Lambda_cutoff = 0:0.1:50 H = zeros(n_ROI,n_ROI); for i = 1:n_ROI if Lambda(i) < Lambda_cutoff 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); % Error made within the selected area of interest MSE_GFT_1(idx_GFT) = mean(MSE_GFT(S1,idx_GFT)); MSE_GFT_2(idx_GFT) = mean(MSE_GFT(S2,idx_GFT)); MSE_GFT_rest(idx_GFT) = mean(MSE_GFT(S_rest,idx_GFT)); idx_GFT = idx_GFT + 1; end TB4 = figure; set(gcf,'color','w'); plot(0:0.1:50,MSE_GFT_1,'color',[0 0.2 0.8],'LineWidth',2); hold on; plot(0:0.1:50,MSE_GFT_2,'color',[1 0.2 0.2],'LineWidth',2); -plot(0:0.1:50,MSE_GFT_rest,'color',[0.6 0.6 0.6],'LineWidth',2); +% plot(0:0.1:50,MSE_GFT_rest,'color',[0.6 0.6 0.6],'LineWidth',2); xlabel('Cutoff lambda'); ylabel('MSE'); legend('First type','Second type','Rest'); %% 4. Reconstruction using a localised GFT for n = 1:100 disp(['Iteration ',num2str(n),'...']); % Set of nodes selected (always imperfect) tmp = S; tmp = tmp(randperm(n_S)); LGFT_set = tmp(1:floor(n_S/2)); % 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/2),'SA',opts); Lambda_loc = diag(Lambda_loc); Lambda_loc = Lambda_loc + 1e-3; C_LGFT = V_loc' * X_noisy(LGFT_set,:); idx_LGFT = 1; Lambda_cutoff = 0:0.1:50; 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_cutoff(i) < lam + 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(LGFT_set,idx_LGFT) = GSP_ComputeMSE(X_LGFT,X(LGFT_set,:)); % Error made within the selected area of interest MSE_LGFT_1(idx_LGFT,n) = nanmean(MSE_LGFT(S1,idx_LGFT)); MSE_LGFT_2(idx_LGFT,n) = nanmean(MSE_LGFT(S2,idx_LGFT)); - MSE_LGFT_rest(idx_LGFT,n) = nanmean(MSE_LGFT(S_rest,idx_LGFT)); + % MSE_LGFT_rest(idx_LGFT,n) = nanmean(MSE_LGFT(S_rest,idx_LGFT)); idx_LGFT = idx_LGFT + 1; end end TB4 = figure; set(gcf,'color','w'); plot_ci(Lambda_cutoff,[mean(MSE_LGFT_1,2),mean(MSE_LGFT_1,2)-std(MSE_LGFT_1,[],2),mean(MSE_LGFT_1,2)+std(MSE_LGFT_1,[],2)],... 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... 'MainLineColor', [1 0.2 0], 'LineWidth', 1, 'LineStyle','-', ... 'LineColor', 'none'); hold on; plot_ci(Lambda_cutoff,[mean(MSE_LGFT_2,2),mean(MSE_LGFT_2,2)-std(MSE_LGFT_2,[],2),mean(MSE_LGFT_2,2)+std(MSE_LGFT_2,[],2)],... 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... 'MainLineColor', [0 0.2 0.8], 'LineWidth', 1, 'LineStyle','-', ... 'LineColor', 'none'); xlabel('Cutoff lambda'); ylabel('MSE'); legend('First type','Second type'); -%% Reconstruction using the Slepians +%% 5. Reconstruction using the Slepians + +% Limit to determine that a vector is 'concentrated' +epsilon = 1e-3; % I want to compare bandwidths -for BW = 10:10:n_ROI +for BW = 10:10:40 for n = 1:100 disp(['Iteration ',num2str(n),'...']); % Set of nodes selected (always imperfect) tmp = S; tmp = tmp(randperm(n_S)); SLEP_set = tmp(1:floor(n_S/2)); % V contains the eigenvectors, Lambda contains the eigenvalues (vector % format) % For each bandwidth, I compute my Slepians [Slep,Slep_coefs,Xi,Mu] = GSP_ComputeSlepians(W,BW,SLEP_set,{X}); Xi = diag(Xi); - C_SLEP = Slep' * X_noisy(SLEP_set,:); + C_SLEP = Slep' * X_noisy; - idx_LGFT = 1; + idx_SLEP = 1; Lambda_cutoff = 0:0.1:50; - MSE_LGFT = nan(n_ROI,length(Lambda_cutoff)); + MSE_SLEP = nan(n_ROI,length(Lambda_cutoff),length(10:10:n_ROI)); for lam = Lambda_cutoff - H = zeros(length(LGFT_set),length(LGFT_set)); + H = zeros(BW,BW); + H2 = zeros(BW,BW); - for i = 1:length(LGFT_set) - if Lambda_cutoff(i) < lam + for i = 1:BW + if Xi(i) < lam H(i,i) = 1; - else - H(i,i) = 0; + + if Xi(i) > epsilon + H2(i,i) = 1; + end end end - X_LGFT = V_loc * H * C_LGFT; + X_SLEP = Slep * H * C_SLEP; + X2_SLEP = Slep * H2 * C_SLEP; - MSE_LGFT(LGFT_set,idx_LGFT) = GSP_ComputeMSE(X_LGFT,X(LGFT_set,:)); + MSE_SLEP(:,idx_SLEP,BW/10) = GSP_ComputeMSE(X_SLEP,X); + MSE2_SLEP(:,idx_SLEP,BW/10) = GSP_ComputeMSE(X2_SLEP,X); % Error made within the selected area of interest - MSE_LGFT_1(idx_LGFT,n) = nanmean(MSE_LGFT(S1,idx_LGFT)); - MSE_LGFT_2(idx_LGFT,n) = nanmean(MSE_LGFT(S2,idx_LGFT)); - MSE_LGFT_rest(idx_LGFT,n) = nanmean(MSE_LGFT(S_rest,idx_LGFT)); + MSE_SLEP_1(idx_SLEP,BW/10,n) = nanmean(MSE_SLEP(S1,idx_SLEP,BW/10)); + MSE_SLEP_2(idx_SLEP,BW/10,n) = nanmean(MSE_SLEP(S2,idx_SLEP,BW/10)); + + MSE2_SLEP_1(idx_SLEP,BW/10,n) = nanmean(MSE2_SLEP(S1,idx_SLEP,BW/10)); + MSE2_SLEP_2(idx_SLEP,BW/10,n) = nanmean(MSE2_SLEP(S2,idx_SLEP,BW/10)); - idx_LGFT = idx_LGFT + 1; + idx_SLEP = idx_SLEP + 1; end end end TB4 = figure; set(gcf,'color','w'); +imagesc(Lambda_cutoff,10:10:n_ROI,mean(MSE_SLEP_1,3)); -plot_ci(Lambda_cutoff,[mean(MSE_LGFT_1,2),mean(MSE_LGFT_1,2)-std(MSE_LGFT_1,[],2),mean(MSE_LGFT_1,2)+std(MSE_LGFT_1,[],2)],... - 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... - 'MainLineColor', [1 0.2 0], 'LineWidth', 1, 'LineStyle','-', ... - 'LineColor', 'none'); - -hold on; - -plot_ci(Lambda_cutoff,[mean(MSE_LGFT_2,2),mean(MSE_LGFT_2,2)-std(MSE_LGFT_2,[],2),mean(MSE_LGFT_2,2)+std(MSE_LGFT_2,[],2)],... - 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... - 'MainLineColor', [0 0.2 0.8], 'LineWidth', 1, 'LineStyle','-', ... - 'LineColor', 'none'); -xlabel('Cutoff lambda'); -ylabel('MSE'); -legend('First type','Second type'); +TB4 = figure; +set(gcf,'color','w'); +imagesc(Lambda_cutoff,10:10:n_ROI,mean(MSE_SLEP_2,3)); + +% plot_ci(Lambda_cutoff,[mean(MSE_LGFT_1,2),mean(MSE_LGFT_1,2)-std(MSE_LGFT_1,[],2),mean(MSE_LGFT_1,2)+std(MSE_LGFT_1,[],2)],... +% 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... +% 'MainLineColor', [1 0.2 0], 'LineWidth', 1, 'LineStyle','-', ... +% 'LineColor', 'none'); +% +% hold on; +% +% plot_ci(Lambda_cutoff,[mean(MSE_LGFT_2,2),mean(MSE_LGFT_2,2)-std(MSE_LGFT_2,[],2),mean(MSE_LGFT_2,2)+std(MSE_LGFT_2,[],2)],... +% 'PatchColor', [0.6 0.6 0.6], 'PatchAlpha', 0.25, 'MainLineWidth', 2, ... +% 'MainLineColor', [0 0.2 0.8], 'LineWidth', 1, 'LineStyle','-', ... +% 'LineColor', 'none'); +% xlabel('Cutoff lambda'); +% ylabel('MSE'); +% legend('First type','Second type');