diff --git a/LogisticRegressionFramework/Plotting/cbrewer/cbrewer.m b/LogisticRegressionFramework/Plotting/cbrewer/cbrewer.m new file mode 100644 index 0000000..26be891 --- /dev/null +++ b/LogisticRegressionFramework/Plotting/cbrewer/cbrewer.m @@ -0,0 +1,128 @@ +function [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% CBREWER - This function produces a colorbrewer table (rgb data) for a +% given type, name and number of colors of the colorbrewer tables. +% For more information on 'colorbrewer', please visit +% http://colorbrewer2.org/ +% +% The tables were generated from an MS-Excel file provided on the website +% http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_updates.html +% +% +% [colormap]=cbrewer(ctype, cname, ncol, interp_method) +% +% INPUT: +% - ctype: type of color table 'seq' (sequential), 'div' (diverging), 'qual' (qualitative) +% - cname: name of colortable. It changes depending on ctype. +% - ncol: number of color in the table. It changes according to ctype and +% cname +% - interp_method: interpolation method (see interp1.m). Default is "cubic" ) +% +% A note on the number of colors: Based on the original data, there is +% only a certain number of colors available for each type and name of +% colortable. When 'ncol' is larger then the maximum number of colors +% originally given, an interpolation routine is called (interp1) to produce +% the "extended" colormaps. +% +% Example: To produce a colortable CT of ncol X 3 entries (RGB) of +% sequential type and named 'Blues' with 8 colors: +% CT=cbrewer('seq', 'Blues', 8); +% To use this colortable as colormap, simply call: +% colormap(CT) +% +% To see the various colormaps available according to their types and +% names, simply call: cbrewer() +% +% This product includes color specifications and designs developed by +% Cynthia Brewer (http://colorbrewer.org/). +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 06.12.2011 +% ------------------------------ +% 18.09.2015 Minor fixes, fixed a bug where the 'spectral' color table did not appear in the preview + + +% load colorbrewer data +load('colorbrewer.mat') +% initialise the colormap is there are any problems +colormap=[]; +if (~exist('interp_method', 'var')) + interp_method='cubic'; +end + +% If no arguments +if (~exist('ctype', 'var') | ~exist('cname', 'var') | ~exist('ncol', 'var')) + disp(' ') + disp('[colormap] = cbrewer(ctype, cname, ncol [, interp_method])') + disp(' ') + disp('INPUT:') + disp(' - ctype: type of color table *seq* (sequential), *div* (divergent), *qual* (qualitative)') + disp(' - cname: name of colortable. It changes depending on ctype.') + disp(' - ncol: number of color in the table. It changes according to ctype and cname') + disp(' - interp_method: interpolation method (see interp1.m). Default is "cubic" )') + + disp(' ') + disp('Sequential tables:') + z={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'Spectral'}; + disp(z') + + disp('Divergent tables:') + z={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn'}; + disp(z') + + disp(' ') + disp('Qualitative tables:') + %getfield(colorbrewer, 'qual') + z={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + disp(z') + + plot_brewer_cmap + return +end + +% Verify that the input is appropriate +ctype_names={'div', 'seq', 'qual'}; +if (~ismember(ctype,ctype_names)) + disp('ctype must be either: *div*, *seq* or *qual*') + colormap=[]; + return +end + +if (~isfield(colorbrewer.(ctype),cname)) + disp(['The name of the colortable of type *' ctype '* must be one of the following:']) + getfield(colorbrewer, ctype) + colormap=[]; + return +end + +if (ncol>length(colorbrewer.(ctype).(cname))) +% disp(' ') +% disp('----------------------------------------------------------------------') +% disp(['The maximum number of colors for table *' cname '* is ' num2str(length(colorbrewer.(ctype).(cname)))]) +% disp(['The new colormap will be extrapolated from these ' num2str(length(colorbrewer.(ctype).(cname))) ' values']) +% disp('----------------------------------------------------------------------') +% disp(' ') + cbrew_init=colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))}; + colormap=interpolate_cbrewer(cbrew_init, interp_method, ncol); + colormap=colormap./255; + return +end + +if (isempty(colorbrewer.(ctype).(cname){ncol})) + + while(isempty(colorbrewer.(ctype).(cname){ncol})) + ncol=ncol+1; + end + disp(' ') + disp('----------------------------------------------------------------------') + disp(['The minimum number of colors for table *' cname '* is ' num2str(ncol)]) + disp('This minimum value shall be defined as ncol instead') + disp('----------------------------------------------------------------------') + disp(' ') +end + +colormap=(colorbrewer.(ctype).(cname){ncol})./255; + +end \ No newline at end of file diff --git a/LogisticRegressionFramework/Plotting/cbrewer/cbrewer_preview.jpg b/LogisticRegressionFramework/Plotting/cbrewer/cbrewer_preview.jpg new file mode 100644 index 0000000..bd2830a Binary files /dev/null and b/LogisticRegressionFramework/Plotting/cbrewer/cbrewer_preview.jpg differ diff --git a/LogisticRegressionFramework/Plotting/cbrewer/change_jet.m b/LogisticRegressionFramework/Plotting/cbrewer/change_jet.m new file mode 100644 index 0000000..b8d4ecb --- /dev/null +++ b/LogisticRegressionFramework/Plotting/cbrewer/change_jet.m @@ -0,0 +1,64 @@ +% This script help produce a new 'jet'-like colormap based on other RGB reference colors + +% ------- I WAS ASKED --------------- +% "is there a chance that you could add a diverging map going from blue to green to red as in jet, +% but using the red and blue from your RdBu map and the third darkest green from your RdYlGn map?"" +% +% ANSWER: +% You should construct the new colormap based on the existing RGB values of 'jet' +% but projecting these RGB values on your new RGB basis. +% ----------------------------------- + +% load colormaps +jet=colormap('jet'); +RdBu=cbrewer('div', 'RdBu', 11); +RdYlGn=cbrewer('div', 'RdYlGn', 11); + +% Define the new R, G, B references (p stands for prime) +Rp=RdBu(1,:); +Bp=RdBu(end, :); +Gp=RdYlGn(end-2, :); +RGBp=[Rp;Gp;Bp]; + +% construct the new colormap based on the existing RGB values of jet +% Project the RGB values on your new basis +newjet = jet*RGBp; + +% store data in a strcuture, easier to handle +cmap.jet=jet; +cmap.newjet=newjet; +cnames={'jet', 'newjet'}; + +% plot the RGB values +fh=figure(); +colors={'r', 'g', 'b'}; +for iname=1:length(cnames) + subplot(length(cnames),1,iname) + dat=cmap.(cnames{end-iname+1}); + for icol=1:size(dat,2) + plot(dat(:,icol), 'color', colors{icol}, 'linewidth', 2);hold on; + end % icol + title([' "' cnames{end-iname+1} '" in RGB plot']) +end + +% plot the colormaps +fh=figure(); +for iname=1:length(cnames) + F=cmap.(cnames{iname}); + ncol=length(F); + fg=1./ncol; % geometrical factor + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + hold all + end % icol + text(-0.1, mean(Y), cnames{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + axis off + set(gcf, 'color', [1 1 1]) + ylim([0.1 1.05.*max(Y)]); + end % iname + diff --git a/LogisticRegressionFramework/Plotting/cbrewer/colorbrewer.mat b/LogisticRegressionFramework/Plotting/cbrewer/colorbrewer.mat new file mode 100644 index 0000000..ec59ef4 Binary files /dev/null and b/LogisticRegressionFramework/Plotting/cbrewer/colorbrewer.mat differ diff --git a/LogisticRegressionFramework/Plotting/cbrewer/interpolate_cbrewer.m b/LogisticRegressionFramework/Plotting/cbrewer/interpolate_cbrewer.m new file mode 100644 index 0000000..e8b5e21 --- /dev/null +++ b/LogisticRegressionFramework/Plotting/cbrewer/interpolate_cbrewer.m @@ -0,0 +1,36 @@ +function [interp_cmap]=interpolate_cbrewer(cbrew_init, interp_method, ncolors) +% +% INTERPOLATE_CBREWER - interpolate a colorbrewer map to ncolors levels +% +% INPUT: +% - cbrew_init: the initial colormap with format N*3 +% - interp_method: interpolation method, which can be the following: +% 'nearest' - nearest neighbor interpolation +% 'linear' - bilinear interpolation +% 'spline' - spline interpolation +% 'cubic' - bicubic interpolation as long as the data is +% uniformly spaced, otherwise the same as 'spline' +% - ncolors=desired number of colors +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + +% just to make sure, in case someone puts in a decimal +ncolors=round(ncolors); + +% How many data points of the colormap available +nmax=size(cbrew_init,1); + +% create the associated X axis (using round to get rid of decimals) +a=(ncolors-1)./(nmax-1); +X=round([0 a:a:(ncolors-1)]); +X2=0:ncolors-1; + +z=interp1(X,cbrew_init(:,1),X2,interp_method); +z2=interp1(X,cbrew_init(:,2),X2,interp_method); +z3=interp1(X,cbrew_init(:,3),X2, interp_method); +interp_cmap=round([z' z2' z3']); + +end \ No newline at end of file diff --git a/LogisticRegressionFramework/Plotting/cbrewer/plot_brewer_cmap.m b/LogisticRegressionFramework/Plotting/cbrewer/plot_brewer_cmap.m new file mode 100644 index 0000000..a5cab9e --- /dev/null +++ b/LogisticRegressionFramework/Plotting/cbrewer/plot_brewer_cmap.m @@ -0,0 +1,50 @@ +% Plots and identifies the various colorbrewer tables available. +% Is called by cbrewer.m when no arguments are given. +% +% Author: Charles Robert +% email: tannoudji@hotmail.com +% Date: 14.10.2011 + + + +load('colorbrewer.mat') + +ctypes={'div', 'seq', 'qual'}; +ctypes_title={'Diverging', 'Sequential', 'Qualitative'}; +cnames{1,:}={'BrBG', 'PiYG', 'PRGn', 'PuOr', 'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral'}; +cnames{2,:}={'Blues','BuGn','BuPu','GnBu','Greens','Greys','Oranges','OrRd','PuBu','PuBuGn','PuRd',... + 'Purples','RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd'}; +cnames{3,:}={'Accent', 'Dark2', 'Paired', 'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3'}; + +figure('position', [314 327 807 420]) +for itype=1:3 + + %fh(itype)=figure(); + subplot(1,3,itype) + + for iname=1:length(cnames{itype,:}) + + ncol=length(colorbrewer.(ctypes{itype}).(cnames{itype}{iname})); + fg=1./ncol; % geometrical factor + + X=fg.*[0 0 1 1]; + Y=0.1.*[1 0 0 1]+(2*iname-1)*0.1; + F=cbrewer(ctypes{itype}, cnames{itype}{iname}, ncol); + + for icol=1:ncol + X2=X+fg.*(icol-1); + fill(X2,Y,F(icol, :), 'linestyle', 'none') + text(-0.1, mean(Y), cnames{itype}{iname}, 'HorizontalAlignment', 'right', 'FontWeight', 'bold', 'FontSize', 10, 'FontName' , 'AvantGarde') + xlim([-0.4, 1]) + hold all + end % icol + %set(gca, 'box', 'off') + title(ctypes_title{itype}, 'FontWeight', 'bold', 'FontSize', 16, 'FontName' , 'AvantGarde') + axis off + set(gcf, 'color', [1 1 1]) + end % iname + ylim([0.1 1.05.*max(Y)]); +end %itype + +set(gcf, 'MenuBar', 'none') +set(gcf, 'Name', 'ColorBrewer Color maps') \ No newline at end of file diff --git a/LogisticRegressionFramework/PostProcessing/compute_clustering_accuracy.m b/LogisticRegressionFramework/PostProcessing/compute_clustering_accuracy.m new file mode 100644 index 0000000..a7faf75 --- /dev/null +++ b/LogisticRegressionFramework/PostProcessing/compute_clustering_accuracy.m @@ -0,0 +1,48 @@ + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Author: +% Xavier Bresson, xavier.bresson@epfl.ch +% +% Course: +% A Network Tour of Data Science +% a.k.a. Graph Algorithms in Data Science +% EPFL, Fall 2015 +% +% Version: +% September 23, 2015 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Clustering accuracy can be defined with the purity measure, defined here: +% Yang-Hao-Dikmen-Chen-Oja, Clustering by nonnegative matrix factorization +% using graph random walk, 2012. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Usages: +% accuracy = compute_clustering_accuracy(C_computed,C_grndtruth,R) +% +% Notations: +% n = nb_data +% +% Input variables: +% C_computed = Computed clusters. Size = n x 1. Values in [1,2,...,R]. +% C_grndtruth = Ground truth clusters. Size = n x 1. Values in [1,2,...,R]. +% R = Number of clusters. +% +% Output variables: +% accuracy = Clustering accuracy of computed clusters. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function accuracy = compute_clustering_accuracy(C_computed,C_grndtruth,R) + +n = length(C_grndtruth); +nb_of_dominant_points_in_class = zeros(R,1); +for k=1:R + v = C_grndtruth(C_computed==k); + [~,nb_of_dominant_points_in_class(k)] = mode(v); +end +accuracy = sum(nb_of_dominant_points_in_class)/ n; + +end \ No newline at end of file diff --git a/RSCHMM_Computations.m b/RSCHMM_Computations.m index 419a6a8..d9f00db 100644 --- a/RSCHMM_Computations.m +++ b/RSCHMM_Computations.m @@ -1,228 +1,263 @@ %% This script shows how to compute the framework's coefficients, and the % ensuing results, for an example simulation. The paper's content is -% essentially re-using similar code across simulation sub-types, with some -% plotting +% essentially re-using similar code across simulation sub-types, across +% methods that retrieve a coefficient set or another, and also for +% experimental data investigations -%% 0. Loading of the example data +%% 1. Loading of the example data and addition of the needed paths to the +% framework functions and plotting utilities +load('ExampleData.mat'); +addpath('./LogisticRegressionFramework/Framework'); +addpath('./LogisticRegressionFramework/PostProcessing'); +addpath(genpath('./LogisticRegressionFramework/Plotting')); +% The loaded data contains: +% - The ground truth parameters for co-activations and causal couplings, +% respectively in Gamma_GT_regions (R x R matrix with network +% assignments)/NetworkAssignment_GT (R x 1 vector with same information), +% and in Beta_GT_networks (N x N causal modulations across +% networks)/Beta_GT_regions (same in R x R format, across regions) +% - The example simulated data, for the training set of S = 50 subjects (in +% TC) and the cross-validation set of S = 30 subjects (in CV) +% - The number of regions R = 35 and number of networks N = 7 -%% 1. If needed, separate TC into TC (training set) and CV -% (cross-validation set) -% Number of training subjects -n_training=50; - -CV = TC(n_training+1:end); -TC = TC(1:n_training); - - - -%% 2. Binarize the time courses +%% 2. Binarizes the time courses TC_bin = RSCHMM_Binarize(TC); CV_bin = RSCHMM_Binarize(CV); clear TC clear CV %% 3. Selects the time points % This specifies the size of each individual structure containing data % points; it does not impact the outputs from the framework, but it plays a % role in computational speed. 5000 was empirically selected as a % satisfying value limit = 5000; % Selection of the time points: X_coact and X_causal are the activity % levels of the "non r" regions at times t+1 and t, respectively, while % Y denotes whether region r has transited or not at time t [X_coact,X_causal,Y] = RSCHMM_Select_WorkingSet_Optimized(TC_bin,limit); [X_coact_CV,X_causal_CV,Y_CV] = RSCHMM_Select_WorkingSet_Optimized(CV_bin,limit); clear TC_bin clear CV_bin %% 4. Computes the coefficients % Tolerance threshold used for convergence epsilon = 1e-2; % Number of iterations n_iter = 5; % Currently unused, but EN enables to perform elastic net regularisation % instead of purely L1 regularisation-based solutions as now. Keep EN = 1 % not to use elastic net EN = 1; % Range of probed lambda values for simulations lambda_range = [10000,9000,8000,7000,6000,5000,4600,4200,3800,... 3200,2800,2600,2400,2200,2000,1800,1600,1400,1200,1000,... 900,800,700,650,600,550,500,450,400,350,300,250,200,180,160,140,120,100,90,80,70,60,... 50,40,30,25,20,18,16,14,12,10,9.6,9.2,8.8,8.4,8,7.6,7.2,6.8,6.4,6,... 5.6,5.2,4.8,4.4,4,3.6,3.2,2.8,2.4,2,1.8,1.6,1.2,0.8,0.4,... 0.2,0.1,0]; % The commented range below is that used in experimental data applications % in the paper. It started at larger values because otherwise, % max(lambda_range) did not yield only null coefficients % lambda_range = [50000,48000,46000,44000,42000,40000,38000,36000,34000,... % 32000,30000,28000,26000,24000,22000,20000,18000,16000,14000,13000,... % 12000,11000,10000,9000,8600,8200,7800,7400,7000,6800,6600,6400,6200,... % 6000,5800,5600,5400,5200,5000,4800,4600,4400,4200,4000,3800,3600,... % 3400,3200,3000,2800,2600,2400,2200,2000,1800,1600,1400,1200,1000,... % 900,800,700,600,500,400,300,200,180,160,140,120,100,90,80,70,60,... % 50,40,30,20,18,16,14,12,10,9.6,9.2,8.8,8.4,8,7.6,7.2,6.8,6.4,6,... % 5.6,5.2,4.8,4.4,4,3.6,3.2,2.8,2.4,2,1.8,1.6,1.4,1.2,1,0.8,0.6,0.4,... % 0.2,0.1,0.06,0.02,0]; % lambda_range = downsample(lambda_range,2); % lambda_range = [lambda_range,0]; n_lambda = length(lambda_range); % Declaration of the coefficient matrices and times for each framework step Alpha = NaN(n_lambda,2,R,5); Beta = NaN(n_lambda,R,2,R,5); Gamma = NaN(n_lambda,R,2,R,5); t1 = NaN(n_lambda,R,5); t2 = NaN(n_lambda,R,5); t3 = NaN(n_lambda,R,5); % Solving per se to retrieve the coefficients. Note that parallel % computations are run across regions with the "parfor" syntax (remove if % you do not have the toolbox, and replace with a simple for loop) idx_rho = 1; % Five values for rho are probed only in the paper; it could be interesting % to examine a larger range, but it will make computations longer for rho = [0,0.25,0.5,0.75,1] parfor r = 1:R [Alpha(:,:,r,idx_rho),Beta(:,:,:,r,idx_rho),Gamma(:,:,:,r,idx_rho),... ~,t1(:,r,idx_rho),t2(:,r,idx_rho),t3(:,r,idx_rho)] = ... RSCHMM_Optimized_ROI_Deluxe(lambda_range,... epsilon,rho,EN,X_coact{r},X_causal{r},Y{r},r,n_iter); end idx_rho = idx_rho + 1; end clear idx_rho clear X_coact clear X_causal clear Y %% 5. Evaluates the log-likelihood on the cross-validation dataset % Initialises the counter for rho values iter_rho = 1; % Declares the log-likelihood matrix LL_RSCHMM = NaN(5,n_lambda,2,R); % Computes the log-likelihood across rho and lambda values on the % cross-validated data for rho = [0,0.25,0.5,0.75,1] iter_regul = 1; for lambda = lambda_range for r = 1:R for i = 1:2 LL_RSCHMM(iter_rho,iter_regul,i,r) = RSCHMM_ComputeLq_Optimized_Deluxe(Alpha(iter_regul,i,r,iter_rho),... Beta(iter_regul,(1:R)~=r,i,r,iter_rho)',Gamma(iter_regul,(1:R)~=r,i,r,iter_rho)',... X_causal_CV{r}(:,i),X_coact_CV{r}(:,i),Y_CV{r}(:,i)); end end iter_regul = iter_regul + 1; end iter_rho = iter_rho + 1; end clear iter_rho clear iter_regul clear X_causal_CV clear X_coact_CV clear Y_CV % Determines the optimal indices (where the log-likelihood is maximal) for % rho and for lambda idx_RSCHMM = NaN(2,R); idx1_RSCHMM = NaN(2,R); idx2_RSCHMM = NaN(2,R); for r = 1:R for i = 1:2 tmp_LL = squeeze(LL_RSCHMM(:,:,i,r)); [~,idx_RSCHMM(i,r)] = max(tmp_LL(:)); [idx1_RSCHMM(i,r),idx2_RSCHMM(i,r)] = ind2sub(size(squeeze(LL_RSCHMM(:,:,i,r))),idx_RSCHMM(i,r)); end end clear tmp_LL %% 6. Puts the coefficients into matrices as a function of coefficient % type (co-activation vs causal) and start state of activity for probed % region (baseline vs active) % Declares the coefficient matrices Beta_final_baseline = NaN(R,R); Beta_final_active = NaN(R,R); Gamma_final_baseline = NaN(R,R); Gamma_final_active = NaN(R,R); % Fills them for r = 1:R Gamma_final_baseline(:,r) = Gamma(idx2_RSCHMM(1,r),:,1,r,idx1_RSCHMM(1,r)); Beta_final_baseline(:,r) = Beta(idx2_RSCHMM(1,r),:,1,r,idx1_RSCHMM(1,r)); Gamma_final_active(:,r) = Gamma(idx2_RSCHMM(2,r),:,2,r,idx1_RSCHMM(2,r)); Beta_final_active(:,r) = Beta(idx2_RSCHMM(2,r),:,2,r,idx1_RSCHMM(2,r)); end % Also creates the same matrix for alpha: Alpha_final has baseline % coefficients as first column, and active ones as second column Alpha_final = RSCHMM_PickAlpha(Alpha,idx1_RSCHMM,idx2_RSCHMM); %% 7. Computes the resulting probabilistic couplings % Separate computations for the baseline and the active state; DP_gamma and % DP_beta outputs are R x R matrices [DP_gammaB,DP_betaB] = RSCHMM_ComputeDeltaP_Optimized(Alpha_final(:,1),Beta_final_baseline,Gamma_final_baseline); [DP_gammaA,DP_betaA] = RSCHMM_ComputeDeltaP_Optimized(Alpha_final(:,2),Beta_final_active,Gamma_final_active); % Subtraction of both state coefficients (used everywhere except in Figure % 5C) DP_gamma = DP_gammaB - DP_gammaA; DP_beta = DP_betaB - DP_betaA; %% 8. Computes the metrics of interest % From the co-activation matrix, determines network assignment assuming N % networks (hierarchical clustering through Ward's linkage analysis) -[T_RSCHMM_combo,Z] = RSCHMM_DetermineCommunities(DP_gamma,N); +[NetworkAssignment_SLR,Z_SLR] = RSCHMM_DetermineCommunities(DP_gamma,N); + +% Accuracy of the network assignments, knowing the ground truth number of +% networks N +Purity_SLR = compute_clustering_accuracy(NetworkAssignment_SLR,NetworkAssignment_GT,N); -Acc_RSCHMM_combo = compute_clustering_accuracy(T_RSCHMM_combo,Figure3_GT,N); -SimG_RSCHMM_combo = corr([jUpperTriMatToVec(Gamma_GT);jUpperTriMatToVec(Gamma_GT')],... +% Similarity of the retrieved co-activation coefficients with the ground +% truth: we consider the upper diagonal and the lower diagonal together, +% because the set of retrieved coefficients is actually not enforced to be +% symmetric by the framework +Sim_Gamma_SLR = corr([jUpperTriMatToVec(Gamma_GT_regions);jUpperTriMatToVec(Gamma_GT_regions')],... [jUpperTriMatToVec(DP_gamma);jUpperTriMatToVec((DP_gamma)')]); -SimB_RSCHMM_combo = corr([jUpperTriMatToVec(Beta_GT_regions);jUpperTriMatToVec(Beta_GT_regions')],... +% Similarity of the retrieved causal coefficients with the ground truth, +% computed in a similar manner +Sim_Beta_SLR = corr([jUpperTriMatToVec(Beta_GT_regions);jUpperTriMatToVec(Beta_GT_regions')],... [jUpperTriMatToVec(DP_beta);jUpperTriMatToVec((DP_beta)')]); -[Graph] = RSCHMM_ConstructPopulationGraph(Beta_GT_regions,Figure3_GT,N); +% Construction of the graph representation associated to the retrieved +% causal outputs, knowing the ground truth networks +[Graph_SLR] = RSCHMM_ConstructPopulationGraph(Beta_GT_regions,NetworkAssignment_GT,N); + +% Computation of sensitivity and specificity +[Sensitivity_SLR,Specificity_SLR] = RSCHMM_ComputeSenSpec(Graph_SLR,Beta_GT_networks); + + + +%% 9. Plots example outputs from the framework + + +% Red/Blue colormap +CM = cbrewer('div','RdBu',1001); + +% Plot of the regional time courses at hand for 200 samples (on one +% subject) +figure; +set(gcf,'color','w'); +imagesc(TC{1}(:,1:200)); +colormap(flipud(CM)); +caxis([-6 6]); +set(gca,'Box','off'); -[Beta_Sensitivity_MAR(idx_noise,idx_subj),Beta_Specificity_MAR(idx_noise,idx_subj),Beta_Accuracy_MAR(idx_noise,idx_subj)] = RSCHMM_ComputeSenSpec(Graph_MAR,Graph); \ No newline at end of file