function [tfce_score,thresholded_maps] = limo_tfce_handling(varargin)

% routine to create tfce files commensurate to boostrapped files
%
% FORMAT  [tfce_score,thresholded_maps] = limo_tfce_handling(filename,'checkfile','yes')
%
% INPUTS filename is the stat file that need to be tfced (if H0 exist it is done too)
%        'checkfile' is 'yes' by default - if 'no' and tfce files already exist,
%                     it overwrites without asking otherwise user is prompted
%
% OUTPUTS tfce_* files are saved on the drive in a tfce folder
%         H0_tfce_* files are saved on the drive in the H0 folder
%         tfce_score is the tfce_score for the data of the file in
% ------------------------------------------------------------------
% Cyril R. Pernet - Copyright (C) LIMO Team 2020


% FOR TFCE VALIDATION - get all the maps before integrating cluster mass
% ----------------------------------------------------------------------
% FORMAT  [tfce_score,thresholded_maps] = limo_tfce_handling(filename,'checkfile','yes')
% OUTPUT  thresholded_maps is all the thresholded maps i.e. for each dh
%                          value in sum(extent(h)^E*height^H*dh)
%         if a boostrap file exist, thresholded_maps is a 1 * 2 cell array
%         with {1} the maps for the observed data and {2} a cell of cells
%         with the maps of each boostrap

%% Check inputs

file = varargin{1};
[filepath,filename,ext] = fileparts(file);
if isempty(filepath); filepath = pwd; end
if isempty(ext)
    file = dir(fullfile(filepath,[filename '*']));
    filename = file.name;
end

if exist(fullfile(filepath,[filename ext]),'file')
    if ~exist(fullfile(filepath,'LIMO.mat'),'file')
        error('no LIMO.mat found next to %s',filename)
    else
        filename = [filename ext];
        LIMO = load(fullfile(filepath,'LIMO.mat'));
        LIMO = LIMO.(cell2mat(fieldnames(LIMO)));
    end
else
    error('can''t find %s',varargin{1})
end

checkfile = 'yes';
for i=2:2:nargin
    if contains(varargin{i},'checkfile')
        checkfile = varargin{i+1};
    end
end

%% quick user check
% ----------------
[fpath,fname,ext] = fileparts(filename);
if isempty(fpath)
    filename = fullfile(LIMO.dir,fname);
end

% files to create
if LIMO.Level == 1
    subname      = limo_get_subname([filename ext]);
    stat_type    = extractAfter(fname,'desc-');
    tfce_file    = fullfile(LIMO.dir,['tfce' filesep subname '_desc-' stat_type 'tfce']);
    H0_tfce_file = fullfile(LIMO.dir,['H0'   filesep subname '_desc-' stat_type 'tfceH0']);
    H0filename   = fullfile(LIMO.dir,['H0'   filesep subname '_desc-' stat_type 'H0.mat']); % .mat as exists 'on disk'
else
    tfce_file    = fullfile(LIMO.dir,['tfce' filesep fname '_desc-tfce']);
    H0_tfce_file = fullfile(LIMO.dir,['H0'   filesep fname '_desc-tfceH0']);
    H0filename   = fullfile(LIMO.dir,['H0'   filesep fname '_desc-H0.mat']);
end

if strcmpi(checkfile,'yes')
    if exist(tfce_file,'file')
        answer = questdlg('tfce file already exist - overwrite?','data check','No','Yes','Yes');
        if strcmp(answer,'Yes')
            LIMO.design.tfce = 1;
            save(fullfile(LIMO.dir,'LIMO.mat'),'LIMO')
        else
            LIMO.design.tfce = 0;
            return
        end
    else
        LIMO.design.tfce = 1;
        if exist(LIMO.dir,'dir')
            save(fullfile(LIMO.dir,'LIMO.mat'),'LIMO','-v7.3')
            if ~exist(fullfile(LIMO.dir,'tfce'),'dir')
                mkdir(fullfile(LIMO.dir,'tfce'));
            end
        else
            save(fullfile(pwd,'LIMO.mat'),'LIMO','-v7.3')
            if ~exist(fullfile(pwd,'tfce'),'dir')
                mkdir(fullfile(pwd,'tfce'));
            end
        end
    end
end

% check if there is a neighbouring matrix
% (since TFCE integrates over clusters)
if ~isfield(LIMO.data,'neighbouring_matrix')
    warning('no neighbouring matrix found, this is required for TFCE')
    [~, LIMO.data.neighbouring_matrix] = limo_expected_chanlocs;
    if isempty(LIMO.data.neighbouring_matrix)
        return
    else
        save(fullfile(LIMO.dir,'LIMO.mat'),'LIMO')
    end
end

% create tfce folder
if ~exist(fullfile(LIMO.dir,'tfce'),'dir')
    mkdir(fullfile(LIMO.dir,'tfce'))
end
fprintf('Thresholding %s using TFCE \n',filename);

% -------------------------------------------------------------------------------
if contains(filename,'R2') || ...
        contains(filename,'semi_partial') % these files last dimension is R2, F, p
    % ----------------------------------------------------------------------------
    if ~isfile(tfce_file)
        R2 = load(filename);
        R2 = R2.(cell2mat(fieldnames(R2)));
        if size(R2,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score(1,:,:),thresholded_maps] = limo_tfce(2, squeeze(R2(:,:,:,2)),[]); % no neighbouring, time-freq cluster
            else
                [tfce_score(1,:),thresholded_maps]   = limo_tfce(1, squeeze(R2(:,:,2)),LIMO.data.neighbouring_matrix);
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score,thresholded_maps] = limo_tfce(3, squeeze(R2(:,:,:,2)),LIMO.data.neighbouring_matrix);
            else
                [tfce_score,thresholded_maps] = limo_tfce(2, squeeze(R2(:,:,2)),LIMO.data.neighbouring_matrix);
            end
        end
        save(tfce_file,'tfce_score','-v7.3'); clear R2 ;
    end

    if exist(H0filename,'file')
        fprintf('Applying TFCE to null data (this may take a while)... \n')
        H0_R2 = load(H0filename);
        H0_R2 = H0_R2.(cell2mat(fieldnames(H0_R2)));
        nboot = size(H0_R2,ndims(H0_R2));
        tfce_H0_thmaps = cell(1,nboot);
        if size(H0_R2,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                tfce_H0_score  = NaN(1,size(H0_R2,2),size(H0_R2,3),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(1,:,:,b),tfce_H0_thmaps] = limo_tfce(2,squeeze(H0_R2(:,:,:,2,b)),[],0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score  = NaN(1,size(H0_R2,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(1,:,b),tfce_H0_thmaps] = limo_tfce(1,squeeze(H0_R2(:,:,2,b)),neighbouring_matrix,0); %#ok<*PFTUSW,*PFOUS>
                end
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score  = NaN(size(H0_R2,1),size(H0_R2,2),size(H0_R2,3),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,:,b),tfce_H0_thmaps] = limo_tfce(3,squeeze(H0_R2(:,:,:,2,b)),neighbouring_matrix,0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score  = NaN(size(H0_R2,1),size(H0_R2,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,b),tfce_H0_thmaps] = limo_tfce(2,squeeze(H0_R2(:,:,2,b)),neighbouring_matrix,0);
                end
            end
        end
        save(H0_tfce_file,'tfce_H0_score','-v7.3'); clear H0_R2 tfce_H0_score;
        tmp                 = thresholded_maps;     clear thresholded_maps;
        thresholded_maps{1} = tmp;                  clear tmp
        thresholded_maps{2} = tfce_H0_thmaps;       clear tfce_H0_thmaps;
    end
    
    % -------------------------------------------------------------------------------------------------------------
elseif contains(filename,'con') || ...
        contains(LIMO.design.name,'One sample','IgnoreCase',true) || ...
        contains(LIMO.design.name,'Two samples','IgnoreCase',true) || ...
        contains(LIMO.design.name,'Paired','IgnoreCase',true)  % these file last dimension is mean, se, df, t and p
    % --------------------------------------------------------------------------------------------------------------
    if ~isfile(tfce_file)
        tval = load(filename);
        tval = tval.(cell2mat(fieldnames(tval)));
        if size(tval,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score(1,:,:),thresholded_maps] = limo_tfce(2, squeeze(tval(:,:,:,end-1)),[]); % no neighbouring, time-freq cluster
            else
                [tfce_score(1,:),thresholded_maps]   = limo_tfce(1, squeeze(tval(:,:,end-1)),LIMO.data.neighbouring_matrix);
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score,thresholded_maps] = limo_tfce(3, squeeze(tval(:,:,:,end-1)),LIMO.data.neighbouring_matrix);
            else
                [tfce_score,thresholded_maps] = limo_tfce(2, squeeze(tval(:,:,4)),LIMO.data.neighbouring_matrix);
            end
        end
        save(tfce_file,'tfce_score','-v7.3'); clear tval ;
    end

    if exist(H0filename,'file')
        fprintf('Applying TFCE to null data (this may take a while)... \n')
        H0_tval = load(H0filename);
        H0_tval = H0_tval.(cell2mat(fieldnames(H0_tval)));
        nboot   = size(H0_tval,ndims(H0_tval));
        tfce_H0_thmaps = cell(1,nboot);
        if size(H0_tval,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                tfce_H0_score  = NaN(1,size(H0_tval,2),size(H0_tval,3),nboot);
                H0_tval        = squeeze(H0_tval(:,:,:,end-1,:));
                parfor b=1:nboot
                    [tfce_H0_score(1,:,:,b),tfce_H0_thmaps{b}] = limo_tfce(2,H0_tval(:,:,:,b),[],0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                H0_tval             = squeeze(H0_tval(:,:,end-1,:));
                tfce_H0_score       = NaN(1,size(H0_tval,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(1,:,b),tfce_H0_thmaps{b}] = limo_tfce(1,H0_tval(:,:,b),neighbouring_matrix,0);
                end
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                H0_tval             = squeeze(H0_tval(:,:,:,end-1,:));
                tfce_H0_score       = NaN(size(H0_tval,1),size(H0_tval,2),size(H0_tval,3),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,:,b),tfce_H0_thmaps{b}] = limo_tfce(3,H0_tval(:,:,:,b),neighbouring_matrix,0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                H0_tval             = squeeze(H0_tval(:,:,end-1,:));
                tfce_H0_score       = NaN(size(H0_tval,1),size(H0_tval,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,b),tfce_H0_thmaps{b}] = limo_tfce(2,H0_tval(:,:,b),neighbouring_matrix,0);
                end
            end
        end
        save(H0_tfce_file,'tfce_H0_score','-v7.3'); clear H0_tval tfce_H0_score;
        tmp                 = thresholded_maps;     clear thresholded_maps;
        thresholded_maps{1} = tmp;                  clear tmp
        thresholded_maps{2} = tfce_H0_thmaps;       clear tfce_H0_thmaps;
    end
    
    % ------------------------------------------
else % anything else last dimension is F and p
    % ------------------------------------------
    if ~isfile(tfce_file)
        Fval = load(filename);
        Fval = Fval.(cell2mat(fieldnames(Fval)));
        if contains(filename,'ess')
            Fval = Fval(:,:,end-1:end);
        end

        if size(Fval,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score(1,:,:),thresholded_maps] = limo_tfce(2,squeeze(Fval(:,:,:,1)),[]);
            else
                [tfce_score(1,:),thresholded_maps] = limo_tfce(1,squeeze(Fval(:,:,1)),LIMO.data.neighbouring_matrix);
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                [tfce_score,thresholded_maps] =  limo_tfce(3,squeeze(Fval(:,:,:,1)),LIMO.data.neighbouring_matrix);
            else
                [tfce_score,thresholded_maps] =  limo_tfce(2,squeeze(Fval(:,:,1)),LIMO.data.neighbouring_matrix);
            end
        end
        save(tfce_file,'tfce_score','-v7.3'); clear Fval;
    end

    if exist(H0filename,'file')
        fprintf('Applying TFCE to null data (this may take a while)... \n')
        H0_Fval = load(H0filename);
        H0_Fval = H0_Fval.(cell2mat(fieldnames(H0_Fval)));
        nboot = size(H0_Fval,ndims(H0_Fval));
        tfce_H0_thmaps = cell(1,nboot);
        if size(H0_Fval,1) == 1
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score = NaN(1,size(H0_Fval,2),size(H0_Fval,3),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(1,:,:,b),tfce_H0_thmaps{b}] = limo_tfce(1,squeeze(H0_Fval(:,:,:,1,b)),neighbouring_matrix,0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score = NaN(1,size(H0_Fval,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(1,:,b),tfce_H0_thmaps{b}] = limo_tfce(1,squeeze(H0_Fval(:,:,1,b)),neighbouring_matrix,0);
                end
            end
        else
            if strcmpi(LIMO.Analysis,'Time-Frequency')
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score = NaN(size(H0_Fval,1),size(H0_Fval,2),size(H0_Fval,3),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,:,b),tfce_H0_thmaps{b}] = limo_tfce(2,squeeze(H0_Fval(:,:,:,1,b)),neighbouring_matrix,0);
                end
            else
                neighbouring_matrix = LIMO.data.neighbouring_matrix;
                tfce_H0_score = NaN(size(H0_Fval,1),size(H0_Fval,2),nboot);
                parfor b=1:nboot
                    [tfce_H0_score(:,:,b),tfce_H0_thmaps{b}] = limo_tfce(2,squeeze(H0_Fval(:,:,1,b)),neighbouring_matrix,0);
                end
            end
        end
        save(H0_tfce_file,'tfce_H0_score','-v7.3'); clear H0_Fval tfce_H0_score
    end
    
    tmp                 = thresholded_maps;     clear thresholded_maps;
    thresholded_maps{1} = tmp;                  clear tmp
    if exist('tfce_H0_thmaps','var') % no bootstrap
        thresholded_maps{2} = tfce_H0_thmaps;       clear tfce_H0_thmaps;
    end
end