%% DEMO file to compute statistics between a voxel inter-subject correlation matrix 
% computed with the ISC toolbox with option "PublicParams.calcCorMatrices = 1;" (i.e. save
% full correlation matrices) and a "model" inter-subject matrix computed with corr/pdist
% from behavioural individual data.
%
% This analysis is also called Mantel test: see bramila_mantel.m for a reusable function
% with code similar to the below.
%
% If you use this analysis please cite: 
%
% Nummenmaa, L., Glerean, E., Viinikainen, M., Jaaskelainen, I. P., Hari, R., & Sams, M.
% (2012). Emotions promote social interaction by synchronizing brain activity across 
% individuals. Proceedings of the National Academy of Sciences, 109(24), 9599–604. 
% doi:10.1073/pnas.1206095109
%
% (c) Enrico Glerean 2013 - Brain and Mind Laboratory Aalto University http://becs.aalto.fi/bml/

%% STEP 1: DEFINE THE MODEL MATRIX

% example A: you have continuous ratings for each subject stored in the matrix behav_data 
% with T rows (num of time points) and N columns (num of subjects)

model = corr(behav_data);   % model is now an intersubject correlation matrix

% example B: you have behavioural scores for each subject (e.g. age) stored in the matrix
% behav_data with N columns (num of subjects)

temp = squareform(pdist(behav_data','euclidean')); % note the transpose operator
temp = 2*temp/max(abs(temp(:))); % normalize the DISSIMILARITY matrix
model = 1-temp; % not we have a SIMILARITY matrix with values between -1 and 1

%% STEP 2: LOAD DATA COMPUTED WITH ISC TOOLBOX
load ISCresults/memMaps.mat
NC=size(memMaps.cormatMap.whole.band0.Session1.cor.Data.xyzc,4); % number of comparisons
data=memMaps.cormatMap.whole.band0.Session1.cor.Data.xyzc;
data=reshape(data,[],NC);
data=data'; % every column is a voxel, every row is a pair of subjects
voxels=intersect(find(var(data)~=0),find(~isnan(var(data)))); % which voxels contain data
iscIDs=find(triu(ones(N),1));   % how the data stored by the ISC toolbox 
                                % remaps to an ISC matrix

%% STEP 3: COMPUTE STATISTICS
% Statistics are stored as spearman correlation values between the top triangle of the 
% model ISC matrix with the single voxel

tempdata=data(:,voxels); % let's do the math only for where we have data
NITER = 5000; % Number of permutations

corr_volume=zeros(91,109,91); % where we stored the correlation values
corr_volume(voxels)=corr(model(iscIDs),tempdata,'type','spearman');

pvals_volume=zeros(91,109,91); % where we stored the estimated p values

% create NITER permutations of the model
surro=zeros(length(iscIDs),NITER);
for iter=1:NITER;
    pe=randperm(size(model,1));
    tempmodel=model(pe,pe);
    surro(:,iter)=tempmodel(iscIDs);
end

% estimate the p value for each voxel
% the code below can be highly optimized if you have large amount of RAM or multiple CPUs
for vID=1:length(voxels)
    voxel=voxels(vID); % which voxel we are calculating
    surroout=corr(surro,tempdata(:,vID),'type','spearman');
    if(isnan(sum(surroout))) continue; end
    if(mod(vID,1000)==1) disp(num2str(vID/length(voxels))); end
    [fi xi]=ksdensity(surroout,'function','cdf','npoints',200);
    pvals_volume(voxel)=1-interp1([-1 xi 1],[0 fi 1],corr_volume(voxel),'spline');
end

%% STEP 4: store results and apply multiple comparisons correction
% store the corr_volume and store the pvals_volume
save_nii(make_nii(corr_volume),'filename_corr.nii');
save_nii(make_nii(pvals_volume),'filename_pvals.nii');

% apply standard FDR with
q=mafdr(pvals_volume(voxels),'BHFDR', TF);

% or correct for multiple comparisons using cluster correction.
% Another option is to store the max statistic from the surrogate data and take the 95th
% percentile as corrected p<0.05.





