%% This script runs the ISFC computations and ISFC transient extraction on
% a set of typically developing subjects, for a movie-watching paradigm
%
% Code written by Thomas Bolton, last checked on July 24th 2018

%% 1. Data loading, paths generation and parameters

% Adds the path towards all utilities
addpath(genpath(fullfile(pwd,'Utilities')));

%%%%%%%%%%%%%%%%%% To be entered by the user %%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Name of the atlas file
Atlas_name = 'Craddock_300.nii';

% Name of the data file
Data_name = 'JOVE_Full_Data.mat';

% Name of the file containing the output data
Save_name = 'JOVE_Results';

% TR for the considered dataset
TR=2.0;

% Subjects that were included in the atlasing step a priori (i.e., name of
% the related subject folders)
Subjects_OI = {'S1','S9','S10','S14','S16','S17','S19','S20','S21','S23','S25',...
    'S26','S27','S31','S32'};

Ses_prefix = 'SES';

% Threshold to use for motion selection (scrubbing), according to Power's
% criterion; advised value is 0.5 [mm]
Mot_Threshold = 0.5;

% Percentage frame that can maximally be scrubbed before the subject is
% taken out; advised value is 10 [%], can go up to 30 [%]
TScrub = 10;

% Size of the sliding windows (in TRs); advised value is W = 10 [TR]
W = 10;

% Step between two windows (in TRs); advised value is Step = 1 [TR]
Step = 1;

% Indices of the sessions, from the ones to analyze, combining RS and
% movie-watching subparts
combined_sessions = [1,2];

% Timing (in samples from the data matrix) of the movie-watching and RS
% subparts for those sessions
MW_timing = 3:178;
RS_timing = 194:340;

% Number of bootstrapping folds to use for ISFC computations
n_boot = 250;

% Number of reference subjects to use for the bootstraping process
n_refsubj = 6;

% alpha-value to use for ISFC thresholding
Alpha = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Path towards the atlas
PathToAtlas = fullfile(pwd,'Utilities',Atlas_name);

% Path towards the folder will all my subjects
PathToSubjFolder = fullfile(pwd,'Subjects');

% Loads the data (called 'TCS' from the previous script)
load(fullfile(pwd,'Results',Data_name));

% Number of different sessions acquired per subject
n_sessions = size(TCS,1);

% For each session, we assign a label
Subject_Label = repmat(Subjects_OI,size(TCS,1),1);
Subject_Label = Subject_Label(:);

% Please switch to true if desiring to plot some figures of the process
is_plot = 1;

% Please switch to true if desiring to plot some text about the process
is_text = 1;

%% 2. Removal of bad quality sessions

% Number of subjects
n_subj = length(Subjects_OI);

% Label for session (1 or 2)
Session_Label = repmat(1:n_sessions,1,length(Subjects_OI));

% Computation of motion for all subjects
[~,~,PScrub] = ComputeFD(PathToSubjFolder,Subjects_OI,Mot_Threshold,Ses_prefix);

% We keep subjects if they show less than 30% of time points with large
% framewise displacement
ToKeep = (PScrub <= TScrub);

if is_text
    disp([num2str(sum(~ToKeep)),'/',num2str(length(ToKeep)),' sessions to discard...']);
end

if is_plot
    figure;
    set(gcf,'color','w');
    bar(1:length(PScrub),PScrub);
    set(gca,'Box','off');
    xlabel('Sessions');
    ylabel('Percentage scrubbed frames [%]');
    xlim([0.5 length(PScrub)+0.5]);
    ylim([0 50]);
    hold on;
    plot([1 length(PScrub)],[TScrub TScrub],'k--');
end

Session_Label = Session_Label(ToKeep);
Subject_Label = Subject_Label(ToKeep);

% I convert my 3 x 30 data into a 1 x 90 format in order to keep only the
% sessions of interest
TCS_all = reshape(TCS,1,max(Session_Label)*size(TCS,2));
TCS_all = TCS_all(ToKeep);

clear PScrub

% Removes the regions in which at least one of the subjects has NaN values
% (i.e. the region did not exist). We use all the sessions for this step
[TCS_ip,idx_removedROIs] = RemoveNaNRegions(TCS_all);

if is_text
   disp([num2str(length(idx_removedROIs)),' regions removed from the atlas...']); 
end

% Filtering of the time courses at cutoff 1/W [s]
TCS_rear = cell(1,length(TCS_ip));
for i = 1:length(TCS_ip)
    TCS_rear{i} = y_IdealFilter(TCS_ip{i}', TR, [1/(W*TR) 0]);
    TCS_rear{i} = TCS_rear{i}';
end

%% 3. Atlas-related details

% We want a vector indexing the ROIs, going from 1 to the number of regions
% kept (not NaN)
ROI_vec = 1:size(TCS_all{1},1);
ROI_vec(idx_removedROIs) = [];
ROI_vec = 1:length(ROI_vec);

% Number of regions and of connections
n_ROI = length(ROI_vec);
n_conn = n_ROI*(n_ROI-1)/2;

if is_text
    disp(['Total number of regions: ',num2str(n_ROI),'...']);
    disp(['Total number of connections: ',num2str(n_conn),'...']);
end

% Converts the region vector into a connection vector (an example entry is
% '1-13' between ROI 1 and ROI 13)
Con_vec = GetConVec(ROI_vec);

% I gather the labels of the regions involved in each connection (for the
% above example, C1(i) = 1 and C2(i) = 13)
[C1,C2] = ParseConVec(Con_vec);

% Creation of the atlas codebook and of 3D center of gravity information,
% with COG of size 3 x n_ROI
[CodeBook,COG] = CreateCodeBook(PathToAtlas,ROI_vec,idx_removedROIs);

% Euclidean distance between all connections considered
for i = 1:n_conn
    D_euclidean(i,1) = sqrt(sum((squeeze(COG(:,C1(i)))-squeeze(COG(:,C2(i)))).^2));
end

if is_plot
    figure;
    set(gcf,'color','w');
    hist(D_euclidean,1000);
    xlim([0 200]);
    set(gca,'Box','off');
    xlabel('Distance between regions [mm]');
    ylabel('Number of connections [-]');
end

%% 4. Extraction of resting-state and movie-watching portions

idx_RSMW = 1;
idx_RS = 1;

% We iterate through all the preprocessed data
for i = 1:length(TCS_rear)
    
    % tmp contains the time courses (n_ROI x n_TPs) for one subject
    tmp = TCS_rear{i};
    
    % If the current recording is a combined stimulus + resting-state case,
    % then we sample both subparts
    switch (ismember(Session_Label(i),combined_sessions))
        case 0
            TC_RS2{idx_RS} = tmp;
            idx_RS = idx_RS + 1;
            
        case 1
            TC_RS{idx_RSMW} = tmp(:,RS_timing);
            TC_MW{idx_RSMW} = tmp(:,MW_timing);

            idx_RSMW = idx_RSMW + 1;
    end
    
    clear tmp
end

if is_text
    disp(['Number of retained pure resting-state sessions: ',num2str(length(TC_RS2)),'...']);
    disp(['Number of retained stimulus-related sessions: ',num2str(length(TC_RS)),'...']);
end

%% 5. Computing ISFC estimates on the data

% First, we wish to solely perform ISFC compuations on the movie-watching
% subparts; we create the ISFC_MW cell array, each cell about to contain
% the outputs for one subject (n_conn x n_slidingwindow_timepoint); then,
% we do the same for the resting-state subparts from movie sessions, and
% finally for the 'pure resting-state' sessions

if is_text
    disp('Computing ISFC on stimulus-related sessions (stimulus segments)...');
end
[ISFC_MW] = Perform_Bootstrapping(TC_MW,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);

if is_text
    disp('Computing ISFC on stimulus-related sessions (resting-state segments)...');
end
[ISFC_RS] = Perform_Bootstrapping(TC_RS,Session_Label(ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);

if is_text
    disp('Computing ISFC on purely resting-state sessions...');
end
[ISFC_RS2] = Perform_Bootstrapping(TC_RS2,Session_Label(~ismember(Session_Label,combined_sessions)),W,Step,n_boot,n_refsubj);

%% 6. Significance thresholding of connectivity pairs
if is_text
    disp('Thresholding ISFC time courses...');
end

% First, I want to determine when in time, and for what connections, I
% observe a significant excursion above from the mean connectivity
% SignMat has size n_conn x n_TP x n_sessions, and contains:
% -1 if there is a negative excursion
% 0 if there is no excursion
% +1 if there is a positive excursion
SignMat = FindExcursions(ISFC_MW,ISFC_RS,ISFC_RS2,Alpha);

% Saves the data
save(fullfile(pwd,'Results',Save_name),'CodeBook','TC_RS','TC_MW',...
    'TC_RS2','ISFC_MW','ISFC_RS','ISFC_RS2','SignMat','Alpha','W','Step',...
    'n_boot','n_refsubj','TR','Mot_Threshold','TScrub','-v7.3');