% ins_temporalClusterPermuTest_WithinBetween.m
function result = ins_temporalClusterPermuTest_WithinBetween(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
% 
% Nonparametric MEG data analysis with two factors (groups and conditions)
% on two levels each. Permutation test based on temporal clusters on the
% data grouped into regions of interest.
% [Analyse non paramétrique de données MEG avec deux facteurs (groupes et
% conditions) à deux niveaux chacun. Test par permutation basé sur des
% clusters temporels sur les données regroupées en régions d'intérêt.]
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_WithinBetween(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% data              = array <1xN struct>
%                   where N = person x roi x cond x other criterium (side)
%                   example of struct : person = 'sujet01'
%                                       group = 'control'
%                                       condition = 'condition_A'
%                                       roi = 'roi_3'
%                                       side = 'left'
%                                       time = <1x193 double>
%                                       avgabs = <1x193 double>
% def               = structure arrays with correspondance between waited
%                   field names and data field names, and list of levels
%                   for each categorical field. Waited field names are :
%                   subject, group, cond, roi, other (other selection criterium),
%                   time, value.
% cfg               = configuration parameters (cf. example).
% structTemplate	= template to reconstruct the data structure waited by
%                   the FieldTrip function ft_timelockgrandaverage.
%                   (cf. http://www.fieldtriptoolbox.org/reference/ft_timelockgrandaverage)
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Oct/2015 (first version)
% Oct/2015 (this version)
% http://ins.univ-amu.fr
 
%% global parameters
cfg.method = 'montecarlo' ;       % use the Monte Carlo Method to calculate the significance probability
cfg.correctm = 'cluster';         % Apply multiple-comparison correction
 
% Interaction between groups and conditions
if def.stats{1}
    disp('*** 1 - Interaction between groups and conditions ***');
    resultat = ins_temporalClusterPermuTest_interaction(data, def, cfg, structTemplate);
    disp(['=== ' resultat ' ===']);
end
 
% Main effect of conditions (within)
if def.stats{2}
    disp('*** 2 - Main effect of conditions (within) ***');
    resultat = ins_temporalClusterPermuTest_within(data, def, cfg, structTemplate);
    disp(['=== ' resultat ' ===']);
end
 
% Main effect of groups (between)
if def.stats{3}
    disp('*** 3 - Main effect of groups (between) ***');
    resultat = ins_temporalClusterPermuTest_between(data, def, cfg, structTemplate);
    disp(['=== ' resultat ' ===']);
end
 
% Simple effect of conditions by group
if def.stats{4}
    disp('*** 4 - Simple effect of conditions by group ***');
    resultat = ins_temporalClusterPermuTest_within_simple(data, def, cfg, structTemplate);
    disp(['=== ' resultat ' ===']);
end
 
% Simple effect of groups by condition
if def.stats{5}
    disp('*** 5 - Simple effect of groups by condition ***');
    resultat = ins_temporalClusterPermuTest_between_simple(data, def, cfg, structTemplate);
    disp(['=== ' resultat ' ===']);
end
 
result = 'END of script ins_temporalClusterPermuTest_WithinBetween.m';

% RunningExample_TCPTWB.m
% Example of utilization of the function
%    ins_temporalClusterPermuTest_WithinBetween.m
% Bernard Giusiano & Sophie Chen - oct 2015
 
%% directory of scripts and example data :
cd('D:\_Recherche\DEV_temporalClusterPermuTest_WithinBetween\scripts ins_TCPTWB');
% cd('C:\Users\dr00245\Dropbox\_RECHERCHE\INS_MEG_MCI3\clustersTemporels');
load('structurePorteuseDataMEG.mat');
% -> variable 'structurePorteuseDataMEG' (-> paramètre structTemplate)
% data loading
load('mydata.mat');
% makink directory for output
currentFolder = pwd;
outputFolder = sprintf('%s\\output', currentFolder); %%% maj dans les fonctions %%%%
if ~exist(outputFolder, 'dir')
  mkdir(outputFolder);
end
 
%% definition of the experiment
mydef.subject = 'person';
mydef.list.subject = {'01','02','03','04','05','06','07','11',...
    '12','13','14','15','16','17','18','19'};
mydef.group = 'group';
mydef.list.group = {'subject','subject','subject','subject','subject','subject','subject','subject',...
    'control','control','control','control','control','control','control','control'};
mydef.cond = 'condition';
mydef.list.cond =  {'condition_A','condition_B'};
mydef.roi = 'roi';
% mydef.list.roi = {'roi_1','roi_2','roi_3','roi_4','roi_5','roi_6','roi_7'};
mydef.list.roi = {'roi_1','roi_5','roi_6','roi_7'}; % select 4 roi only
mydef.other = 'side';
% mydef.list.other = {'L','R'};
mydef.list.other = {'L'};                           % select left side only
mydef.time = 'time';
mydef.value = 'avgabs';
 
mydef.stats = {...  % select statistical calculations to be made (true) or not (false):
    true,...    % Interaction between groups and conditions
    false,...    % Main effect of conditions (within)
    true,...    % Main effect of groups (between)
    true,...    % Simple effect of conditions by group
    false,...    % Simple effect of groups by condition
    };
 
%% global parameters
% for other parameters values cf. http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
cfg = [];
cfg.latency = [-0.3 0.9];         % time interval over which the experimental 
                                  % conditions must be compared (in seconds) or 'all' (default = 'all')
cfg.numrandomization = 10000;     % Number of draws from the permutation distribution
cfg.clusteralpha = 0.05;          % alpha level of the sample-specific statistic that will be used for thresholding
                                  % the values used to construct clusters
cfg.tail = 0;                     % -1, 1 or 0 (default = 0); one-sided or two-sided test
cfg.clusterstatistic = 'wcm';  % How to combine the single samples that belong to a cluster
                                  % -> test statistic that will be evaluated under the permutation distribution
%   cfg.clusterstatistic = 'maxsum', 'maxsize', 'wcm' (default = 'maxsum')
%                          option 'wcm' refers to 'weighted cluster mass',
%                          a statistic that combines cluster size and
%                          intensity; see Hayasaka & Nichols (2004) NeuroImage
%                          for details
% When cfg.clusterstatistic = 'wcm':
      cfg.wcm_weight = 0.5;      % cf. Hayasaka & Nichols (2004) NeuroImage (page 62)
cfg.clustertail = 0;              % -1, 1 or 0 (default = 0 : two-tail test)
cfg.correcttail = 'alpha';        % correct p-values or alpha-values when doing a two-sided test
                                  % cf. http://www.fieldtriptoolbox.org/faq/why_should_i_use_the_cfg.correcttail_option_when_using_statistics_montecarlo
cfg.alpha = 0.05;  % 0.05/10;     % alpha level of the permutation test / 10 tests
% ATTENTION:
% alpha SHOULD BE CORRECTED based on the number of tests
% example of correction for 7 roi, 2 groups, 2 conditions, 2 sides:
%   - interactions: 7 roi * 2 sides = 14 tests
%   - main effects: 2 factors (groups and conditions) * 7 roi * 2 sides = 28 tests
%   - simple effects: 2 groups * 2 cond * 7 roi * 2 sides = 56 tests
% => 14 + 28 + 56 = 98 tests => Bonferroni's correction: 0.05/98 = 0.0005
% => cfg.alpha = 0.0005;  for a global alpha threshold really equal to 0.05
 
% Other parameters:
% cfg.clusterthreshold = 'nonparametric_common'; % method for single-sample threshold
 
cfg.avgoverchan = 'yes';   % average over channels => neighbours no needed
cfg.neighbours  = {};      % cf. http://mailman.science.ru.nl/pipermail/fieldtrip/2013-July/006853.html
 
%% trace %%%%%%%%%%%%%%%%
% traces are saved during execution of script ins_temporalClusterPermuTest_graph.m
log_filename = ['./output/logProbaMin_' datestr(now, 'yyyymmdd_HHMM') '.txt'];
cfg.fid = fopen(log_filename, 'wt');
%%%%%%%%%%%%%%%%%%%%%%%%%
%%
resultat = ins_temporalClusterPermuTest_WithinBetween(mydata, mydef, cfg, structurePorteuseDataMEG);
disp(['### ' resultat ' ###']);
 
disp('*** END of script RunningExample_temporalClusterPermuTest_WithinBetween.m ***');
 
%% trace %%%%%%%%%%%%%%%%
fclose(cfg.fid);
%%%%%%%%%%%%%%%%%%%%%%%%%

% ins_temporalClusterPermuTest_interaction.m
function result = ins_temporalClusterPermuTest_interaction(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of interaction between two factors (groups and conditions)
% on two levels each. Permutation test based on temporal clusters on the
% data grouped into regions of interest.
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_interaction(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_WithinBetween.m
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% * Eric Maris. How to test an interaction effect using cluster-based permutation tests?
%   http://www.fieldtriptoolbox.org/faq/how_can_i_test_an_interaction_effect_using_cluster-based_permutation_tests
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Oct/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
groupDesc = tabulate(def.list.group);
nGroupLevel1 = groupDesc{1,2};
nGroupLevel2 = groupDesc{2,2};
 
sizeSamples = size(data(1,1).(def.time));
nSamples = sizeSamples(2);
 
cfgGM = [];
cfgGM.keepindividual = 'yes';
cfgGM.latency = cfg.latency;
 
%% parameters specific to this function
cfg.statistic = 'indepsamplesT';
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design
design = zeros(1,nGroupLevel1 + nGroupLevel2);
design(1,1:nGroupLevel1) = 1;
design(1,(nGroupLevel1 + 1):(nGroupLevel1 + nGroupLevel2))= 2;
cfg.design = design;            % design matrix
cfg.ivar = 1;                   % number or list with indices indicating the independent variable(s)
 
%% loop on criteria
for roi_ind = 1:length(def.list.roi)
    roi = char(def.list.roi(roi_ind));
    fig = figure;
    set(fig,'Units','pixels');
    set(fig,'Position',[10 415 1419 399]);
    % AxesTitleFont = ['arial'];
    nplot = 1;
 
    for other_ind = 1:length(def.list.other)
        other = char(def.list.other(other_ind));
 
        component = cell(2,2,1);  % 2 groups and 2 conditions
							      % 3° dim -> length(def.list.subject) but initialized to 1  
        ind_component = zeros(2,2);
 
        for subject_ind = 1:length(def.list.subject)
            subject = char(def.list.subject(subject_ind));
 
            group = char(def.list.group(subject_ind));
            group_ind = find(ismember(groupDesc(:,1),group));
 
            for cond_ind = 1:length(def.list.cond)
                cond = char(def.list.cond(cond_ind));
 
                disp([roi,' ',cond,' ',other,' ',group,' ',subject]);
                object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)...
                    & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject));
 
				objectstat = [group cond subject];
                eval([objectstat ' = structTemplate;']);
                eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']);
                eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']);
                eval([objectstat '.fsample = nSamples;']);
                eval([objectstat '.sampleinfo = [1 nSamples];']);
                eval([objectstat '.trialinfo = 1;']);
                ind_component(group_ind,cond_ind) = ind_component(group_ind,cond_ind) + 1;
                component{group_ind,cond_ind,ind_component(group_ind,cond_ind)} = eval(objectstat);
            end
        end
 
        % grand average by group and condition
        % GA11: group level 1 - cond level 1, GA12: group level 1 - cond level 2
        % GA21: group level 2 - cond level 1, GA22: group level 2 - cond level 2
        GA11 = ft_timelockgrandaverage(cfgGM,component{1,1,1:ind_component(1,1)});
        GA12 = ft_timelockgrandaverage(cfgGM,component{1,2,1:ind_component(1,2)});
        GA21 = ft_timelockgrandaverage(cfgGM,component{2,1,1:ind_component(2,1)});
        GA22 = ft_timelockgrandaverage(cfgGM,component{2,2,1:ind_component(2,2)});
 
        set1.GM = GA11;
        set1.GM.individual = GA11.individual - GA12.individual;
        set2.GM = GA11;
        set2.GM.individual = GA21.individual - GA22.individual;
 
        set1.mean = squeeze(mean(set1.GM.individual,1));
        set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nGroupLevel1));
        set2.mean = squeeze(mean(set2.GM.individual,1));
        set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nGroupLevel2));
 
        % test calculus
        stat_dif1_dif2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
        disp([roi,' ',cond,' ',other,' mask: ',num2str(sum(stat_dif1_dif2.mask))]);
 
        % figures
        subplot(1,2,nplot);
        nplot = nplot + 1;
 
        set1.legend = [def.list.cond{1} '-' def.list.cond{2} ' ' groupDesc{1,1}];
        set2.legend = [def.list.cond{1} '-' def.list.cond{2} ' ' groupDesc{2,1}];
        stat_dif1_dif2.title = [roi,' ',other];
        stat_dif1_dif2.ylabel = ['difference between averages (' def.value ')'];
        res = ins_temporalClusterPermuTest_graph(set1, set2, stat_dif1_dif2, def, cfg);
 
    end
    filename = ['./output/f_interaction_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_interaction.m';

% ins_temporalClusterPermuTest_within.m
function result = ins_temporalClusterPermuTest_within(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects within two conditions
% on two levels (paired samples). Permutation test based on temporal clusters on the
% data grouped into regions of interest.
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_within(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_WithinBetween.m
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Oct/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
sizeSamples = size(data(1,1).(def.time));
nSamples = sizeSamples(2);
 
cfgGM = [];
cfgGM.keepindividual = 'yes';
cfgGM.latency = cfg.latency;
 
%% parameters specific to this function
cfg.statistic = 'depsamplesT';    % paired samples
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design
nSubjects = length(def.list.subject);
design = zeros(2,2*nSubjects);
for i = 1:nSubjects
    design(1,i) = i;
end
for i = 1:nSubjects
    design(1,nSubjects+i) = i;
end
design(2,1:nSubjects)             = 1;
design(2,nSubjects+1:2*nSubjects) = 2;
cfg.design = design;            % design matrix
cfg.uvar  = 1;                  % indice of the unit variable
cfg.ivar  = 2;                  % number or list with indices indicating the independent variable(s)
 
%% loop on criteria
for roi_ind = 1:length(def.list.roi)
    roi = char(def.list.roi(roi_ind));
    fig = figure;
    set(fig,'Units','pixels');
    set(fig,'Position',[10 415 1419 399]);
    nplot = 1;
 
    for other_ind = 1:length(def.list.other)
        other = char(def.list.other(other_ind));
 
        % components of grand average
        component = cell(2,1);  % 2 conditions. 2° dim -> length(def.list.subject) but initialized to 1  
        ind_component = zeros(2);
 
        for subject_ind = 1:length(def.list.subject)
            subject = char(def.list.subject(subject_ind));
 
            for cond_ind = 1:length(def.list.cond)
                cond = char(def.list.cond(cond_ind));
 
                disp([roi,' ',cond,' ',other,' ',subject]);
                object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)...
                    & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject));
 
                objectstat = [cond subject];
                eval([objectstat ' = structTemplate;']);
                eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']);
                eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']);
                eval([objectstat '.fsample = nSamples;']);
                eval([objectstat '.sampleinfo = [1 nSamples];']);
                eval([objectstat '.trialinfo = 1;']);
                ind_component(cond_ind) = ind_component(cond_ind) + 1;
                component{cond_ind,ind_component(cond_ind)} = eval(objectstat);
            end
        end
 
        % grand average by condition
        % set1.GM: cond level 1, set2.GM: cond level 2
        set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1:ind_component(1)});
        set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1:ind_component(2)});
 
        set1.mean = squeeze(mean(set1.GM.individual,1));
        set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nSubjects));
        set2.mean = squeeze(mean(set2.GM.individual,1));
        set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nSubjects));
 
        % test calculus
        stat_cond1_cond2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
        disp([roi,' ',cond,' ',other,' ',subject,' mask: ',num2str(sum(stat_cond1_cond2.mask))]);
 
        % figures
        subplot(1,2,nplot);
        nplot = nplot + 1;
 
        set1.legend = def.list.cond{1};
        set2.legend = def.list.cond{2};
        stat_cond1_cond2.title = [roi,' ',other];
        stat_cond1_cond2.ylabel = def.value;
        res = ins_temporalClusterPermuTest_graph(set1, set2, stat_cond1_cond2, def, cfg);
    end
    filename = ['./output/f_within_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_within.m';

% ins_temporalClusterPermuTest_between.m
function result = ins_temporalClusterPermuTest_between(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects between two groups
% on two levels. Permutation test based on temporal clusters on the
% data grouped into regions of interest.
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_between(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_WithinBetween.m
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Oct/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
groupDesc = tabulate(def.list.group);
nGroupLevel1 = groupDesc{1,2};
nGroupLevel2 = groupDesc{2,2};
groupLevels = unique(def.list.group);
 
sizeSamples = size(data(1,1).(def.time));
nSamples = sizeSamples(2);
 
cfgGM = [];
cfgGM.keepindividual = 'yes';
cfgGM.latency = cfg.latency;
 
%% parameters specific to this function
cfg.statistic = 'indepsamplesT';    % independant samples
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design
design = zeros(1,nGroupLevel1*2 + nGroupLevel2*2);
design(1,1:nGroupLevel1*2) = 1;
design(1,(nGroupLevel1*2 + 1):(nGroupLevel1*2 + nGroupLevel2*2))= 2;
cfg.design = design;             % design matrix
cfg.ivar  = 1;                   % number or list with indices indicating the independent variable(s)
 
%% loop on criteria
for roi_ind = 1:length(def.list.roi)
    roi = char(def.list.roi(roi_ind));
    fig = figure;
    set(fig,'Units','pixels');
    set(fig,'Position',[10 415 1419 399]);
    nplot = 1;
 
    for other_ind = 1:length(def.list.other)
        other = char(def.list.other(other_ind));
 
        component = cell(2,2,1);  % 2 groups and 2 conditions (each subject has 2 cond)
							      % 3° dim -> length(def.list.subject) but initialized to 1  
        ind_component = zeros(2,2);
 
        for subject_ind = 1:length(def.list.subject)
            subject = char(def.list.subject(subject_ind));
 
            group = char(def.list.group(subject_ind));
            group_ind = find(ismember(groupDesc(:,1),group));
 
            for cond_ind = 1:length(def.list.cond)
                cond = char(def.list.cond(cond_ind));
 
                disp([roi,' ',cond,' ',other,' ',group,' ',subject]);
                object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)...
                    & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject));
 
                objectstat = [group cond subject];
                eval([objectstat ' = structTemplate;']);
                eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']);
                eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']);
                eval([objectstat '.fsample = nSamples;']);
                eval([objectstat '.sampleinfo = [1 nSamples];']);
                eval([objectstat '.trialinfo = 1;']);
                ind_component(group_ind,cond_ind) = ind_component(group_ind,cond_ind) + 1;
                component{group_ind,cond_ind,ind_component(group_ind,cond_ind)} = eval(objectstat);
            end
        end
 
        % grand average by group
        % set1.GM: group level 1, set2.GM: group level 2
        set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1,1:ind_component(1,1)},...
            component{1,2,1:ind_component(1,2)});
        set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1,1:ind_component(2,1)},...
            component{2,2,1:ind_component(2,2)});
 
        set1.mean = squeeze(mean(set1.GM.individual,1));
        set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nGroupLevel1));
        set2.mean = squeeze(mean(set2.GM.individual,1));
        set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nGroupLevel2));
 
        % test calculus
        stat_group1_group2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
        disp([roi,' ',other,' ',group,' ',subject,' mask: ',num2str(sum(stat_group1_group2.mask))]);
 
        % figure
        subplot(1,2,nplot);
        nplot = nplot + 1;
 
        set1.legend = groupLevels{1};
        set2.legend = groupLevels{2};
        stat_group1_group2.title = [roi,' ',other];
        stat_group1_group2.ylabel = def.value;
        res = ins_temporalClusterPermuTest_graph(set1, set2, stat_group1_group2, def, cfg);
    end
    filename = ['./output/f_between_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_between.m';

% ins_temporalClusterPermuTest_within_simple.m
function result = ins_temporalClusterPermuTest_within_simple(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects within two conditions
% on two levels (paired samples) for each group. Permutation test based on
% temporal clusters on the data grouped into regions of interest.
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_within_simple(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_WithinBetween.m
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Nov/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
groupLevels = unique(def.list.group);
 
sizeSamples = size(data(1,1).(def.time));
nSamples = sizeSamples(2);
 
cfgGM = [];
cfgGM.keepindividual = 'yes';
cfgGM.latency = cfg.latency;
 
%% parameters specific to this function
cfg.statistic = 'depsamplesT';    % paired samples
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design : in the loop on criteria because it may be specific for each group
 
%% loop on criteria
for roi_ind = 1:length(def.list.roi)
    roi = char(def.list.roi(roi_ind));
    fig = figure;
    set(fig,'Units','pixels');
    set(fig,'Position',[10 74 1419 740]);
    nplot = 1;
 
    for group_ind = 1:length(groupLevels)
        group = char(groupLevels(group_ind));
        groupSubjects = def.list.subject(strcmp(def.list.group,group));
 
        % design
        nSubjects = length(groupSubjects);
        design = zeros(2,2*nSubjects);
        for i = 1:nSubjects
            design(1,i) = i;
        end
        for i = 1:nSubjects
            design(1,nSubjects+i) = i;
        end
        design(2,1:nSubjects)             = 1;
        design(2,nSubjects+1:2*nSubjects) = 2;
        cfg.design = design;            % design matrix
        cfg.uvar  = 1;                  % indice of the unit variable
        cfg.ivar  = 2;                  % number or list with indices indicating the independent variable(s)
        % end design
 
        for other_ind = 1:length(def.list.other)
            other = char(def.list.other(other_ind));
 
            % components of grand average
            component = cell(2,1);  % 2 conditions. 2° dim -> length(def.list.subject) but initialized to 1
            ind_component = zeros(2);
 
            for subject_ind = 1:length(groupSubjects)
                subject = char(groupSubjects(subject_ind));
 
                for cond_ind = 1:length(def.list.cond)
                    cond = char(def.list.cond(cond_ind));
 
                    disp([roi,' ',cond,' ',other,' ',group,' ',subject]);
                    object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)...
                        & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject));
 
                    objectstat = [cond subject];
                    eval([objectstat ' = structTemplate;']);
                    eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']);
                    eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']);
                    eval([objectstat '.fsample = nSamples;']);
                    eval([objectstat '.sampleinfo = [1 nSamples];']);
                    eval([objectstat '.trialinfo = 1;']);
                    ind_component(cond_ind) = ind_component(cond_ind) + 1;
                    component{cond_ind,ind_component(cond_ind)} = eval(objectstat);
                end
            end
 
            % grand average by condition
            % set1.GM: cond level 1, set2.GM: cond level 2
            set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1:ind_component(1)});
            set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1:ind_component(2)});
 
            set1.mean = squeeze(mean(set1.GM.individual,1));
            set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nSubjects));
            set2.mean = squeeze(mean(set2.GM.individual,1));
            set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nSubjects));
 
            % test calculus
            stat_cond1_cond2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
            disp([roi,' ',cond,' ',other,' ',subject,' mask: ',num2str(sum(stat_cond1_cond2.mask))]);
 
            % figure
            subplot(2,2,nplot);
            nplot = nplot + 1;
 
            set1.legend = def.list.cond{1};
            set2.legend = def.list.cond{2};
            stat_cond1_cond2.title = [roi,' ',group,' ',other];
            stat_cond1_cond2.ylabel = def.value;
            res = ins_temporalClusterPermuTest_graph(set1, set2, stat_cond1_cond2, def, cfg);
        end
    end
    filename = ['./output/f_within_simple_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_within_simple.m';

% ins_temporalClusterPermuTest_between_simple.m
function result = ins_temporalClusterPermuTest_between_simple(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects between two groups
% on two levels for each condition. Permutation test based on temporal
% clusters on the data grouped into regions of interest.
%
% DEPENDANCES:
% This function uses the FieldTrip MATLAB software toolbox
% (http://www.fieldtriptoolbox.org/)
%
% USAGE:
% result = ins_temporalClusterPermuTest_between_simple(
%                               data,
%                               def,
%                               cfg,
%                               structTemplate)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_WithinBetween.m
%
% OUTPUTS:
% result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% REFERENCES:
% * Cluster-based permutation tests on event related fields.
%   http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Nov/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
groupDesc = tabulate(def.list.group);
nGroupLevel1 = groupDesc{1,2};
nGroupLevel2 = groupDesc{2,2};
groupLevels = unique(def.list.group);
 
sizeSamples = size(data(1,1).(def.time));
nSamples = sizeSamples(2);
 
cfgGM = [];
cfgGM.keepindividual = 'yes';
cfgGM.latency = cfg.latency;
 
%% parameters specific to this function
cfg.statistic = 'indepsamplesT';  % independant samples
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design
design = zeros(1,nGroupLevel1 + nGroupLevel2);
design(1,1:nGroupLevel1) = 1;
design(1,(nGroupLevel1 + 1):(nGroupLevel1 + nGroupLevel2))= 2;
cfg.design = design;             % design matrix
cfg.ivar  = 1;                   % number or list with indices indicating the independent variable(s)
 
%% loop on criteria
for roi_ind = 1:length(def.list.roi)
    roi = char(def.list.roi(roi_ind));
    fig = figure;
    set(fig,'Units','pixels');
    set(fig,'Position',[10 74 1419 740]);
    nplot = 1;
 
    for cond_ind = 1:length(def.list.cond)
        cond = char(def.list.cond(cond_ind));
 
        for other_ind = 1:length(def.list.other)
            other = char(def.list.other(other_ind));
 
            % components of grand average
            component = cell(2,1);  % 2 conditions. 2° dim -> length(def.list.subject) but initialized to 1
            ind_component = zeros(2);
 
            for subject_ind = 1:length(def.list.subject)
                subject = char(def.list.subject(subject_ind));
 
                group = char(def.list.group(subject_ind));
                group_ind = find(ismember(groupDesc(:,1),group));
 
                disp([roi,' ',cond,' ',other,' ',group,' ',subject]);
                object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)...
                    & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject));
 
                objectstat = [group subject];
                eval([objectstat ' = structTemplate;']);
                eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']);
                eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']);
                eval([objectstat '.fsample = nSamples;']);
                eval([objectstat '.sampleinfo = [1 nSamples];']);
                eval([objectstat '.trialinfo = 1;']);
                ind_component(group_ind) = ind_component(group_ind) + 1;
                component{group_ind,ind_component(group_ind)} = eval(objectstat);
            end
 
            % grand average by group
            % set1.GM: group level 1, set2.GM: group level 2
            set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1:ind_component(1)});
            set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1:ind_component(2)});
 
            set1.mean = squeeze(mean(set1.GM.individual,1));
            set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nGroupLevel1));
            set2.mean = squeeze(mean(set2.GM.individual,1));
            set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nGroupLevel2));
 
            % test calculus
            stat_group1_group2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
            disp([roi,' ',other,' ',group,' ',subject,' mask: ',num2str(sum(stat_group1_group2.mask))]);
 
            % figure
            subplot(2,2,nplot);
            nplot = nplot + 1;
 
            set1.legend = groupLevels{1};
            set2.legend = groupLevels{2};
            stat_group1_group2.title = [roi,' ',cond,' ',other];
            stat_group1_group2.ylabel = def.value;
            res = ins_temporalClusterPermuTest_graph(set1, set2, stat_group1_group2, def, cfg);
        end
    end
    filename = ['./output/f_between_simple_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_between_simple.m';

% ins_temporalClusterPermuTest_graph.m
function result = ins_temporalClusterPermuTest_graph(set1, set2, stat_set1_set2, def, cfg)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Specific function to draw graphics of functions
% ins_temporalClusterPermuTest_XXX
%
% DEPENDANCES:
% This function uses the function 'jbfill' to shade area between two curves
% (http://www.mathworks.com/matlabcentral/fileexchange/13188-shade-area-between-two-curves/content/jbfill.m)
%
% USAGE:
% result = ins_temporalClusterPermuTest_graph(
%                               set1,
%                               set2,
%                               stat_set1_set2,
%                               def,
%                               cfg)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_XXX.m
%
% OUTPUTS:
% result = true or false if error
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Nov/2015 (this version)
% http://ins.univ-amu.fr
%
result = false;
 
%% trace %%%%%%%%%%%%%%%%
p = min(stat_set1_set2.prob);
st=dbstack;
fprintf(cfg.fid, '%s;%s;%6.4f\n', st(2,1).name, stat_set1_set2.title, p);
%%%%%%%%%%%%%%%%%%%%%%%%%
 
plot(set1.GM.time,set1.mean,'b','linewidth',2);
hold on;
plot(set2.GM.time,set2.mean,'r','linewidth',2);
 
plot(set1.GM.time,stat_set1_set2.stat,'g','linewidth',0.5);
line(cfg.latency,[0 0],'Color',[0.1 0.1 0.1],'linestyle','--');  % xlim
 
y1 = stat_set1_set2.mask;
y2(~y1) = nan ;
y2(y1) = 0.0 ;
plot(set1.GM.time,y2,'Color',[1.0,0.4,0.0],'linewidth',10); % clusters significatifs en orange
 
legend(set1.legend,set2.legend,'stat diff','Location','NorthWest');
 
x = jbfill(set1.GM.time,(set1.mean+1.96*set1.stdmean).',(set1.mean-1.96*set1.stdmean).','b','b',0,0.05);
x = jbfill(set2.GM.time,(set2.mean+1.96*set2.stdmean).',(set2.mean-1.96*set2.stdmean).','r','r',0,0.05);
 
title(stat_set1_set2.title,'FontWeight','bold','FontSize',14)
xlabel(def.time);ylabel(stat_set1_set2.ylabel);
hold off;
 
%%
result = true;