%% This function performs within-group ISFC computations between a given
% subject, and all the other subjects from the same group
%
% Inputs:
% - X is the data of all the subjects that we want to consider (cell array 
% with each cell n_ROI x n_TP)
% - Xref is the data of the reference group that we want to use in the
% process (same dimensions of the cells)
% - n_boot is the number of folds over which to randomly select a subset of
% reference subjects
% - W is the window size (in samples)
% - Step is the window step (in samples)
%
% Outputs:
% - Y is the output (cell array with each cell of size n_pairs x n_swtp)
% containing dynamic ISFC measurements for the subjects in X (one subject
% output per cell)
% - CM_stat is a cell array, with each cell containing the n_ROI x n_ROI
% static ISFC information for one subject of X when compared to the
% reference group
%
% Written by Thomas Bolton, last checked on July 24th
function [Y] = ISFC_ToReferenceGroup(X,Xref,W,Step)

    % Number of time points of the data
    n_TP = size(X{1},2);
    
    % Number of subjects
    n_subjects = length(X);
    
    % Current end of the window; we want to iterate the computations until
    % the window end is past the end of the data
    WindowEnd = W;
    WindowStart = 1;
    
    time_w = 1;
    
    full_index = 1;
    
    % Contains our output
    Y = cell(1,n_subjects);
    %CM_stat = cell(1,n_subjects);
    
    % We will do the same for each sliding window...
    while WindowEnd <= n_TP

        % For each subject, we want to compute the estimate for a
        % given region pair comparing this particular region to all
        % the reference group subjects' paired area
        for i = 1:n_subjects

            % tmp is our data for subject i
            tmp = X{i}(:,WindowStart:WindowEnd);

            tmp_corr = 0;
            
            % We will basically compute all pairs of correlation, and
            % average the results
            for j = 1:length(Xref)
                
                to_add = corr(tmp',(Xref{j}(:,WindowStart:WindowEnd))');
                tmp_corr = tmp_corr + ((to_add + to_add')/2);
                
                clear to_add
            end

            tmp_corr = tmp_corr/length(Xref);
            
            % Samples back the upper diagonal (contains all the
            % information, since we symmeterized tmp_corr) and puts it in Y
            Y{i}(:,time_w) = jUpperTriMatToVec(tmp_corr); 

            clear tmp_corr
            clear tmp
        end    

        % We shift the window by the step size
        WindowStart = WindowStart + Step;
        WindowEnd = WindowEnd + Step;
        time_w = time_w + 1;
        full_index = full_index + 1;
    end

    % Now we do something similar, but for the stationary part
%     for i = 1:n_subjects
%         
%         tmp = X{i};
% 
%         tmp_corr_stat = 0;
%         
%         for j = 1:length(Xref)
%                 
%             to_add_stat = corr(tmp',Xref{j}');
%             tmp_corr_stat = tmp_corr_stat + ((to_add_stat + to_add_stat')/2);
% 
%             clear to_add_stat
%         end
% 
%         tmp_corr_stat = tmp_corr_stat/length(Xref);
%             
%         CM_stat{i} = tmp_corr_stat; 
%         
%         clear tmp_corr_stat
%         clear tmp
%     end
end