function stat = clusterTest(cfg0, X, smooth_factor)
% Modify by alexis to include a smoothing effect in each iteration as in
% the main data.

%% Initialization
% Capture current state of rng before randomizing it
oldRng = rng;
rng('shuffle');

% Default cfg0 settings
if ~isfield(cfg0, 'tFeedback')
    tFeedback = 2;
else
    tFeedback = cfg0.tFeedback;
end

% Initialize variables
tNullPos = nan(cfg0.nPerm, 1);
tNullNeg = nan(cfg0.nPerm, 1);

nSubj = length(X);

stat = [];

%% Do permutations
tStart = tic;
% compute matrix size
matrix_size =   size(X{1});
 
for iPerm = 1:cfg0.nPerm

    
    % Flip half of the accuracies
   iFlip = (rand(matrix_size) < .5);
   
   Xperm = {};
    for a = 1:nSubj
    % randomly flip the values of half of the accuracies/Discriminant
    % channel valuies
    Xmat =  X{a};
    Xmat(iFlip) =  -Xmat(iFlip);
    Xmat =  imgaussfilt( Xmat,smooth_factor); % apply smoothing factor to data to equate the phenomenon in the real data
    Xperm{a} =    Xmat;

    end

    % Convert to matrix
    Xperm = conv2mat(Xperm);

    % Perform univariate tests
    [~, pPerm, ~, tstats] = ttest(Xperm, 0, 'dim', length(size(Xperm)));
    tPerm = tstats.tstat;
   
    % Find clusters, separately for negative and positive ones
    posClusters = feval(cfg0.clusterFun, cfg0.clusterFunCfg, (pPerm < cfg0.pThreshold) & (tPerm >= 0));
    negClusters = feval(cfg0.clusterFun, cfg0.clusterFunCfg, (pPerm < cfg0.pThreshold) & (tPerm < 0));
    
    % Find maximum test statistic, seperately for negative and positive clusters
    for iCluster = 1:length(posClusters)
        % Calculate test-statistic
        tCluster = sum(tPerm(posClusters{iCluster}));
        
        % Only keep largest
        tNullPos(iPerm) = max(tNullPos(iPerm), tCluster);
    end

    for iCluster = 1:length(negClusters)
        % Calculate test-statistic
        tCluster = sum(tPerm(negClusters{iCluster}));
        
        % Only keep lowest value
        tNullNeg(iPerm) = min(tNullNeg(iPerm), tCluster);
    end
    
    % Feedback
    if (toc(tStart) > tFeedback)
        fprintf('Finished permutation %g/%g\n', iPerm, cfg0.nPerm);
        tStart = tic;
    end
end

%% Calculate clusters and test-statistics for real data
% Convert to matrix
X = conv2mat(X);

% Perform univariate tests
[~, pReal, ~, tstats] = ttest(X, 0, 'dim', length(size(X)));
tReal = tstats.tstat;

% Find clusters, separately for negative and positive ones
stat.posClusters = feval(cfg0.clusterFun, cfg0.clusterFunCfg, (pReal < cfg0.pThreshold) & (tReal >= 0));
stat.negClusters = feval(cfg0.clusterFun, cfg0.clusterFunCfg, (pReal < cfg0.pThreshold) & (tReal < 0));

nPosCluster = length(stat.posClusters);
nNegCluster = length(stat.negClusters);

% Calculate p-value per cluster
stat.pPos = nan(nPosCluster, 1);
stat.pNeg = nan(nNegCluster, 1);

for iCluster = 1:nPosCluster
    % Calculate test-statistic
    tCluster = sum(tReal(stat.posClusters{iCluster}));

    % Calculate p-value
    stat.pPos(iCluster) = mean(tNullPos > tCluster);    
end

for iCluster = 1:nNegCluster
    % Calculate test-statistic
    tCluster = sum(tReal(stat.negClusters{iCluster}));

    % Calculate p-value
    stat.pNeg(iCluster) = mean(tNullNeg < tCluster);    
end

%% Return and clean-up
% Sort according to p-value
[~, iPos] = sort(stat.pPos);
[~, iNeg] = sort(stat.pNeg);

stat.pPos = stat.pPos(iPos);
stat.pNeg = stat.pNeg(iNeg);

stat.posClusters = stat.posClusters(iPos);
stat.negClusters = stat.negClusters(iNeg);

% Multiply p-values by two, cause we're testing two-sided
stat.pPos = stat.pPos*2;
stat.pPos(stat.pPos > 1) = 1;

stat.pNeg = stat.pNeg*2;
stat.pNeg(stat.pNeg > 1) = 1;

% Restore old state of rng
rng(oldRng);

%% Save to disk
if isfield(cfg0, 'outputFile')
    save(cfg0.outputFile, 'stat', 'cfg0');
end       

end