%% Cluster-based permutation test - Example
function [Sig, p_stat] = CBPtest(data1, data2, time, c_alpha, p_alpha, randNum)
% This function requires FieldTrip toolbox
%

NoS  = size(data1, 1); % Number of NoSect

if isempty(data2)
	data2 = zeros(size(data1, 1), size(data1, 2));
end

%% Parameters
Cond1(:,1,:) = data1; % Subject x Channel x Time
Cond2(:,1,:) = data2;

if exist('c_alpha')
	if ~isempty(c_alpha)
		clustAlpha = c_alpha;
	else
		clustAlpha = 0.05;
	end
else
	clustAlpha = 0.05;
end

if exist('p_alpha')
	if ~isempty(p_alpha)
		pAlpha = p_alpha;
	else
		pAlpha = 0.025;
	end
else
	pAlpha = 0.025;
end

if exist('randNum')
	if ~isempty(randNum)
		bootNum = randNum;
	else
		bootNum = 5000;
	end
else
	bootNum = 5000;
end

clear Cond1Mat
Cond1Mat.label      = cell(1,1);
Cond1Mat.label{1}   = 'MLC11'; % Any channel name
Cond1Mat.fsample    = 100;
Cond1Mat.individual = Cond1;
Cond1Mat.avg        = mean(Cond1,1);
Cond1Mat.time       = time;
Cond1Mat.dimord     = 'subj_chan_time';

clear Cond2Mat
Cond2Mat.label      = cell(1,1);
Cond2Mat.label{1}   = 'MLC11'; % Any channel name
Cond2Mat.fsample    = 100;
Cond2Mat.individual = Cond2;
Cond2Mat.avg        = mean(Cond2,1);
Cond2Mat.time       = time;
Cond2Mat.dimord     = 'subj_chan_time';

cfg                  = [];
cfg.channel          = 1;
cfg.latency          = 'all'; % Time of interest
cfg.parameter        = 'individual';
cfg.method           = 'montecarlo';
cfg.statistic        = 'depsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = clustAlpha;
cfg.clusterstatistic = 'maxsum';
cfg.clusterthreshold = 'nonparametric_individual';
cfg.neighbours       = [];
cfg.clustertail      = 0;
cfg.alpha            = pAlpha;
cfg.correcttail      = 'alpha';
cfg.tail             = 0;
cfg.numrandomization = bootNum;

% Design matrix
design = zeros(2,2*NoS);
for i = 1:NoS
    design(1,i) = i;
end
for i = 1:NoS
    design(1,NoS+i) = i;
end
design(2,1:NoS)       = 1;
design(2,NoS+1:2*NoS) = 2;

cfg.design = design;
cfg.uvar   = 1;
cfg.ivar   = 2;

[stat] = ft_timelockstatistics(cfg, Cond1Mat, Cond2Mat);

%% Visualization

switch cfg.tail
    case 0
        if isempty(stat.posclusters) && isempty(stat.negclusters)
            return
        end
        
        if ~isempty(stat.posclusters) && isempty(stat.negclusters)
            
            sigMat_pos = [];
            for i = 1:size(stat.posclusters, 2)
                if stat.posclusters(i).prob < 0.05
                    sigMat_pos = [sigMat_pos i];
                end
            end
            
            sigWind = [];
            
            for j = 1:length(sigMat_pos)
                sigWind = [sigWind stat.time(find(stat.posclusterslabelmat==sigMat_pos(j)))];
            end
            
        end
        
        if  size(stat.negclusters,2) ~=0 && size(stat.posclusters,2) == 0
            sigMat_neg = [];
            for i = 1:size(stat.negclusters, 2)
                if stat.negclusters(i).prob < 0.05
                    sigMat_neg = [sigMat_neg i];
                end
            end
            
            sigWind = [];
            
            for j = 1:length(sigMat_neg)
                sigWind = [sigWind stat.time(find(stat.negclusterslabelmat==sigMat_neg(j)))];
            end
        end
        
        if size(stat.posclusters,2) ~= 0 && size(stat.negclusters,2) ~=0
            sigMat_pos = [];
            for i = 1:size(stat.posclusters, 2)
                if stat.posclusters(i).prob < 0.05
                    sigMat_pos = [sigMat_pos i];
                end
            end
            
            sigMat_neg = [];
            for i = 1:size(stat.negclusters, 2)
                if stat.negclusters(i).prob < 0.05
                    sigMat_neg = [sigMat_neg i];
                end
            end
            
            sigWind = [];
            
            for j = 1:length(sigMat_pos)
                sigWind = [sigWind stat.time(find(stat.posclusterslabelmat==sigMat_pos(j)))];
            end
            
            for j = 1:length(sigMat_neg)
                sigWind = [sigWind stat.time(find(stat.negclusterslabelmat==sigMat_neg(j)))];
            end
        end
        
    case 1
        
        if ~isempty(stat.posclusters)
            sigMat_pos = [];
            for i = 1:size(stat.posclusters, 2)
                if stat.posclusters(i).prob < 0.05
                    sigMat_pos = [sigMat_pos i];
                end
            end
            
            sigWind = [];
            
            for j = 1:length(sigMat_pos)
                sigWind = [sigWind stat.time(find(stat.posclusterslabelmat==sigMat_pos(j)))];
            end
            
        end
        
    case -1
        if  size(stat.negclusters,2) ~=0 
            sigMat_neg = [];
            for i = 1:size(stat.negclusters, 2)
                if stat.negclusters(i).prob < 0.05
                    sigMat_neg = [sigMat_neg i];
                end
            end
            
            sigWind = [];
            
            for j = 1:length(sigMat_neg)
                sigWind = [sigWind stat.time(find(stat.negclusterslabelmat==sigMat_neg(j)))];
            end
        end
end
        
% 
Sig = sigWind;
p_stat = stat;
