clustGRP.m

% clustGRP() - Tests the null hypothesis that the grand average voltage
%             of a between-subject difference wave or ERP averaged across two
%             groups is mu (usually 0) using a cluster-based 
%             permutation test and the "cluster mass" statistic 
%             (Bullmore et al., 1999; Maris & Oostenveld, 2007).  Note, mu 
%             is assumed to be 0 by default.  This function requires 
%             individual subject ERPs to be stored in a "GRP" structure 
%             and outputs the test results in a number of graphical 
%             and text formats. For analogous within-subject comparisons use
%             the function clustGND.m.
%               Note, when applied to a bin that is the mean ERP across two
%             groups (i.e., NOT a difference wave), a one-sample test is
%             executed.
%             
% 
% 
% Usage:
%  >> [GRP, prm_pval, data_t]=clustGRP(GRP_or_fname,bin,varargin);
%
% Required Inputs:
%   GRP_or_fname - A GRP structure variable or the filename of a GRP 
%                  structure that has been saved to disk.  A GRP variable 
%                  is based on GND variables. To create a GRP variable from 
%                  GND variables use GNDs2GRP.m.  See Mass Univariate ERP 
%                  Toolbox documentation for detailed information about the 
%                  format of a GRP variable. If you specifiy a filename be 
%                  sure to include the file's path, unless the file is in 
%                  the current working directory.
%   bin          - [integer] The bin to contrast against the mean of the
%                  null hypothesis. Use the function headinfo.m to see what 
%                  bins are stored in a GRP variable.  Use the function 
%                  bin_opGRP.m to create a difference wave between two bins 
%                  whose significance you can test with this function.
%
% Optional Inputs:
%   tail          - [1 | 0 | -1] An integer specifying the tail of the
%                   hypothesis test.  "1" means upper-tailed (i.e., alternative 
%                   hypothesis is that the ERP/difference wave is greater 
%                   than mu).  "0" means two-tailed (i.e., alternative hypothesis
%                   is that the ERP/difference wave is not equal to mu).  
%                   "-1" means lower-tailed (i.e., alternative hypothesis
%                   is that the ERP/difference wave is less than mu).
%                   {default: 0}
%   alpha         - A number between 0 and 1 specifying the family-wise
%                   alpha level of the test. {default: 0.05}
%   thresh_p      - The test-wise p-value threshold for cluster inclusion. If
%                   a channel/time-point has a t-score that corresponds to an
%                   uncorrected p-value greater than thresh_p, it is assigned
%                   a p-value of 1 and not considered for clustering. Note
%                   that thresh_p automatically takes into account the tail of
%                   the test (e.g., you will get positive and negative t-score
%                   thresholds for a two-tailed test).
%   chan_hood     - A scalar or a 2D symmetric binary matrix that indicates
%                   which channels are considered neighbors of other 
%                   channels. E.g., if chan_hood(2,10)=1, then Channel 2 
%                   and Channel 10 are nieghbors. You can produce a 
%                   chan_hood matrix using the function spatial_neighbors.m. 
%                   If a scalar is provided, then all electrodes within that 
%                   distance of a particular electrode are considered 
%                   neighbors. Note, EEGLAB's electrode coordinates assume
%                   the head has a radius of 1. See the help documentation 
%                   of the function spatial_neighbors to see how you could
%                   convert this distance threshold to centimeters. 
%                   {default: 0.61}
%   head_radius   - The radius of the head in whatever units the Cartesian
%                   coordinates in GRP.chanlocs are in. This is used to
%                   convert scalar values of chan_hood into centimeters.
%                   {default: []}
%   n_perm        - The number of permutations to use in the test.  As this
%                   value increases, the test takes longer to compute and 
%                   the results become more reliable.  Manly (1997) suggests 
%                   using at least 1000 permutations for an alpha level of 
%                   0.05 and at least 5000 permutations for an alpha level 
%                   of 0.01. {default: 2500}
%   time_wind     - Pair of time values specifying the beginning and
%                   end of a time window in ms (e.g., [160 180]). Every
%                   single time point in the time window will be individually
%                   tested (i.e., maximal temporal resolution) if mean_wind
%                   option is NOT used. Note, boundaries of time window
%                   may not exactly correspond to desired time window                  
%                   boundaries because of temporal digitization (i.e., you
%                   only have samples every so many ms). {default: 0 ms to
%                   the end of the epoch}
%   mean_wind     - ['yes' or 'no'] If 'yes', the permutation test will be
%                   performed on the mean amplitude within the time window 
%                   specified by time_wind.  This sacrifices temporal 
%                   resolution to increase test power by reducing the number
%                   of comparisons.  If 'no', every single time point within
%                   time_wind's time windows will be tested individually.
%                   {default: 'no'}
%   null_mean     - [number] The mean of the null hypothesis (in units of 
%                   microvolts). {default: 0}
%   exclude_chans - A cell array of channel labels to exclude from the
%                   permutation test (e.g., {'A2','lle','rhe'}).  This option 
%                   sacrifices spatial resolution to increase test power by 
%                   reducing the number of comparisons. Use headinfo.m to see
%                   the channel labels stored in the GRP variable. You cannot
%                   use both this option and 'include_chans' (below).{default: 
%                   not used, all channels included in test}
%   include_chans - A cell array of channel labels to use in the permutation
%                   test (e.g., {'A2','lle','rhe'}).  All other channels will
%                   be ignored. This option sacrifices spatial resolution to 
%                   increase test power by reducing the number of comparisons.
%                   Use headinfo.m to see the channel labels stored in the GRP
%                   variable. You cannot use both this option and 
%                   'exclude_chans' (above). {default: not used, all channels 
%                   included in test}
%   verblevel     - An integer specifiying the amount of information you want
%                   this function to provide about what it is doing during runtime.
%                    Options are:
%                      0 - quiet, only show errors, warnings, and EEGLAB reports
%                      1 - stuff anyone should probably know
%                      2 - stuff you should know the first time you start working
%                          with a data set {default value}
%                      3 - stuff that might help you debug (show all
%                          reports)
%   plot_gui      - ['yes' or 'no'] If 'yes', a GUI is created for
%                   visualizing the results of the permutation test using the 
%                   function gui_erp.m. The GUI vizualizes the grand average 
%                   ERPs in each bin via various stats (uV, t-scores), shows 
%                   topographies at individual time points, and illustrates 
%                   which electrodes significantly differ from the null 
%                   hypothesis.  This option does not work if mean_wind 
%                   option is set to 'yes.' This GUI can be reproduced using
%                   the function gui_erp.m. {default: 'yes'}
%   plot_raster   - ['yes' or 'no'] If 'yes', a two-dimensional (time x channel)
%                   binary "raster" diagram is created to illustrate the
%                   results of the permutation test.  Significant negative and
%                   positive deviations from the null hypothesis are shown
%                   as black and white rectangles respectively. Non-
%                   significant comparisons are shown as gray rectangles. 
%                   Clicking on the rectangles will show you the 
%                   corresponding time and channel label for that
%                   rectangle. This figure can be reproduced with the 
%                   function sig_raster.m. {default: 'yes'}
%   plot_mn_topo  - ['yes' or 'no'] If 'yes', the topography of the mean
%                   voltages/effects in the time window is produced.  More
%                   specifically, two figures are produced: one showing the
%                   topographies in uV the other in t-scores. Significant/
%                   nonsignificant comparisons are shown as white/black 
%                   electrodes. Clicking on electrodes will show the
%                   electrode's name.  This figure can be reproduced with
%                   the function sig_topo.m.  This option has NO effect if
%                   mean_wind option is set to 'no'. {default: 'yes'}
%   output_file   - A string indicating the name of a space delimited text
%                   file to produce containing the p-values of all comparisons 
%                   and the details of the test (e.g., number of permutations, 
%                   family-wise alpha level, etc...). If mean_wind option is
%                   set to 'yes,' t-scores of each comparison are also 
%                   included since you cannot derive them from the t-scores
%                   at each time point/electrode in a simple way. When 
%                   importing this file into a spreadsheet be sure NOT to count
%                   consecutive spaces as multiple delimiters. {default: none}
%   save_GRP      - ['yes' or 'no'] If 'yes', the GRP variable will be
%                   saved to disk after the permutation test is completed 
%                   and added to it. User will first be prompted to verify 
%                   file name and path. {default: 'yes'}
%   reproduce_test- [integer] The number of the permutation test stored in
%                   the GRP variable to reproduce.  For example, if 
%                   'reproduce_test' equals 2, the second t-test 
%                   stored in the GRP variable (i.e., GRP.t_tests(2)) will 
%                   be reproduced.  Reproduction is accomplished by setting
%                   the random number generator used in the permutation test 
%                   to the same initial state it was in when the permutation 
%                   test was first applied.
%
% Outputs:
%   GRP           - GRP structure variable.  This is the same as
%                   the input GRP variable with one addition: the 
%                   field GRP.t_tests will contain the results of the 
%                   permutation test and the test parameters. 
%   prm_pval      - A two-dimensional matrix (channel x time) of the
%                   p-values of each comparison.
%   data_t        - A two-dimensional matrix (channel x time) of the
%                   t-scores of each comparison.
%
%   Note also that a great deal of information about the test is displayed 
%   in the MATLAB command window.  You can easiy record of all this
%   information to a text file using the MATLAB command "diary."
%
% Global Variables:
%   VERBLEVEL = Mass Univariate ERP Toolbox level of verbosity (i.e., tells 
%               functions how much to report about what they're doing during
%               runtime) set by the optional function argument 'verblevel'
%
% Notes:
% -To add a difference wave to a GRP variable, use the function "bin_opGRP".
%
% -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
% are possible (at most the number of possible permutations).  Since the
% number of possible permutations grows rapdily with the number of
% participants, this is only issue for small sample sizes (e.g., 3
% participants in each group).  When you have such a small sample size, the
% limited number of p-values may make the test overly conservative (e.g., 
% you might be forced to use an alpha level of .0286 since it is the biggest
% possible alpha level less than .05).
%
%
% Author:
% David Groppe
% May, 2011
% Kutaslab, San Diego
%
% References:
% Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor, 
% E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory 
% and permutation, for a difference between two groups of structural MR 
% images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42. 
% doi:10.1109/42.750253
%
% Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of 
% EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 
% doi:10.1016/j.jneumeth.2007.03.024
%
% Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
% biology. 2nd ed. Chapmn and Hall, London.

%%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
%
% 12/11/2011-Now uses Amy Guthormsen's recursiveless find_clusters.m.
%


function [GRP, prm_pval, data_t]=clustGRP(GRP_or_fname,bin,varargin)

global VERBLEVEL;

p=inputParser;
p.addRequired('GRP_or_fname',@(x) ischar(x) || isstruct(x));
p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0));
p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('alpha',0.05,@(x) isnumeric(x) && (x>0) && (x<1));
p.addParamValue('thresh_p',0.05,@(x) (length(x)==1) && isnumeric(x) && (x>0) && (x<1));
p.addParamValue('chan_hood',0.61,@isnumeric);
p.addParamValue('head_radius',[],@isnumeric);
p.addParamValue('time_wind',[],@(x) isnumeric(x) && (size(x,2)==2));
p.addParamValue('mean_wind','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
p.addParamValue('n_perm',2500,@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x));
p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x));
p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
p.addParamValue('output_file',[],@ischar);
p.addParamValue('reproduce_test',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('save_GRP','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));

p.parse(GRP_or_fname,bin,varargin{:});

if isempty(p.Results.verblevel),
    if isempty(VERBLEVEL),
        VERBLEVEL=2;
    end
else
   VERBLEVEL=p.Results.verblevel; 
end
    
mean_wind=str2bool(p.Results.mean_wind);

%Load GRP struct
if ischar(GRP_or_fname),
    fprintf('Loading GRP struct from file %s.\n',GRP_or_fname);
    load(GRP_or_fname,'-MAT');
else
    GRP=GRP_or_fname;
    clear GRP_or_fname;
end
[n_chan, n_pt, n_bin]=size(GRP.grands);
n_group=length(GRP.GND_fnames);
VerbReport(sprintf('Experiment: %s',GRP.exp_desc),2,VERBLEVEL);

if (bin>n_bin),
    error('There is no Bin %d in this GRP variable.',bin); 
end
grpA=GRP.bin_info(bin).groupA;
grpB=GRP.bin_info(bin).groupB;


%% Figure out which channels to ignore if any
%But first make sure exclude & include options were not both used.
if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans)
    error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.');
end
if ischar(p.Results.exclude_chans),
    exclude_chans{1}=p.Results.exclude_chans;
elseif isempty(p.Results.exclude_chans)
    exclude_chans=[];
else
    exclude_chans=p.Results.exclude_chans;
end
if ischar(p.Results.include_chans),
    include_chans{1}=p.Results.include_chans;
elseif isempty(p.Results.include_chans)
    include_chans=[];
else
    include_chans=p.Results.include_chans;
end

% exclude and include chan options
if ~isempty(exclude_chans),
    ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem
    ct=0;
    for x=1:length(exclude_chans),
        found=0;
        for c=1:n_chan,
            if strcmpi(exclude_chans{x},GRP.chanlocs(c).labels),
                found=1;
                ct=ct+1;
                ignore_chans(ct)=c;
            end
        end
        if ~found,
            watchit(sprintf('I attempted to exclude %s.  However no such electrode was found in GRP variable.', ...
                exclude_chans{x}));
        end
    end
    ignore_chans=ignore_chans(1:ct);
    use_chans=setdiff(1:n_chan,ignore_chans);
elseif ~isempty(include_chans),
    use_chans=zeros(1,length(include_chans)); %preallocate mem
    ct=0;
    for x=1:length(include_chans),
        found=0;
        for c=1:n_chan,
            if strcmpi(include_chans{x},GRP.chanlocs(c).labels),
                found=1;
                ct=ct+1;
                use_chans(ct)=c;
            end
        end
        if ~found,
            watchit(sprintf('I attempted to include %s.  However no such electrode was found in GRP variable.', ...
                include_chans{x}));
        end
    end
    use_chans=use_chans(1:ct);
else
    use_chans=1:n_chan;
end

% Establish spatial neighborhood matrix for making clusters
if isscalar(p.Results.chan_hood),
    chan_hood=spatial_neighbors(GRP.chanlocs(use_chans),p.Results.chan_hood);
else
    chan_hood=p.Results.chan_hood;
end


%% Find time points
if isempty(p.Results.time_wind),
    time_wind=[0 GRP.time_pts(end)]; %default time window
else
    time_wind=p.Results.time_wind;
end
time_wind=sort(time_wind,2); %first make sure earlier of each pair of time points is first
time_wind=sort(time_wind,1); %next sort time windows from earliest to latest onset
n_wind=size(time_wind,1);
if n_wind>1,
    error('clustGRP.m can only handle a single mean time window at the moment.  Try tmaxGRP.m or tfdrGRP.m if you want to simultaneously test hypotheses in multiple mean time windows.');
end

if mean_wind,
    use_tpts=cell(1,n_wind);
else
    use_tpts=[];
end
for a=1:n_wind,
    VerbReport(sprintf('Time Window #%d:',a),1,VERBLEVEL);
    VerbReport(sprintf('Attempting to use time boundaries of %d to %d ms for hypothesis test.',time_wind(a,1),time_wind(a,2)), ...
        1,VERBLEVEL);
    start_tpt=find_tpt(time_wind(a,1),GRP.time_pts);
    end_tpt=find_tpt(time_wind(a,2),GRP.time_pts);
    if mean_wind,
        use_tpts{a}=[start_tpt:end_tpt];
    else
        use_tpts=[use_tpts [start_tpt:end_tpt]];
    end
    %replace desired time points with closest matches
    time_wind(a,1)=GRP.time_pts(start_tpt);
    time_wind(a,2)=GRP.time_pts(end_tpt);
    VerbReport(sprintf('Exact window boundaries are %d to %d ms (that''s from time point %d to %d).', ...
        time_wind(a,1),time_wind(a,2),start_tpt,end_tpt),1,VERBLEVEL);
end
if ~mean_wind,
    use_tpts=unique(use_tpts); %sorts time points and gets rid of any redundant time points
end


%% Compile data from two groups of participants
%pre-allocate memory
grp_ct=0;
for grp=[grpA grpB],
    grp_ct=grp_ct+1;
    if grp_ct==1,
        %Group A's bin
        ur_bin=GRP.bin_info(bin).source_binA;
    else
        %Group B's bin
        ur_bin=GRP.bin_info(bin).source_binB;
    end
    
    %Load GND variable for this group
    load(GRP.GND_fnames{grp},'-MAT');
    VerbReport(sprintf('Loading individual participant ERPs from Bin %d (%s) of GND variable from %s.', ...
        ur_bin,GND.bin_info(ur_bin).bindesc,GRP.GND_fnames{grp}),2,VERBLEVEL);
    
    %Check to make sure time points are still compatible across GND and GRP
    %variables
    if ~isequal(GND.time_pts,GRP.time_pts)
        error('The time points in the GND variable from file %s are NOT the same as those in your GRP variable.', ...
            GRP.GND_fnames{grp});
    end
    
    %Derive channel indexes for this GND variable in case GND.chanlocs differs from
    %GRP.chanlocs (in order or in identity of channels)
    n_chanGND=length(GND.chanlocs);
    use_chansGND=[];
    for a=use_chans,
        found=0;
        for b=1:n_chanGND,
            if strcmpi(GRP.chanlocs(a).labels,GND.chanlocs(b).labels),
                found=1;
                use_chansGND=[use_chansGND b];
                break;
            end
        end
        if ~found,
            error('GND variable in file %s is missing channel %s.',GRP.GND_fnames{grp},GRP.chanlocs(a).labels);
        end
    end
    
    %Use only subs with data in relevant bin(s)
    use_subs=find(GND.indiv_bin_ct(:,ur_bin));
    n_sub=length(use_subs);
    
    if mean_wind,
        %Take mean amplitude in time blocks and then test
        if grp_ct==1,
            %Group A's ERPs
            n_subA=n_sub;
            use_subsA=use_subs;
            erpsA=zeros(length(use_chans),n_wind,n_sub);
            for a=1:n_wind,
                for sub=1:n_sub,
                    erpsA(:,a,sub)=mean(GND.indiv_erps(use_chansGND,use_tpts{a},ur_bin,use_subs(sub)),2);
                end
            end
        else
            %Group B's ERPs
            n_subB=n_sub;
            use_subsB=use_subs;
            erpsB=zeros(length(use_chans),n_wind,n_sub);
            for a=1:n_wind,
                for sub=1:n_sub,
                    erpsB(:,a,sub)=mean(GND.indiv_erps(use_chansGND,use_tpts{a},ur_bin,use_subs(sub)),2);
                end
            end
        end
    else
        %Use every single time point in time window(s)
        if grp_ct==1,
            %Group A's ERPs
            n_subA=n_sub;
            use_subsA=use_subs;
            n_use_tpts=length(use_tpts);
            erpsA=zeros(length(use_chans),n_use_tpts,n_sub);
            for sub=1:n_sub,
                erpsA(:,:,sub)=GND.indiv_erps(use_chansGND,use_tpts,ur_bin,use_subs(sub));
            end
        else
            %Group B's ERPs
            n_subB=n_sub;
            use_subsB=use_subs;
            n_use_tpts=length(use_tpts);
            erpsB=zeros(length(use_chans),n_use_tpts,n_sub);
            for sub=1:n_sub,
                erpsB(:,:,sub)=GND.indiv_erps(use_chansGND,use_tpts,ur_bin,use_subs(sub));
            end
        end
    end
end
df=n_subA+n_subB-2;

%% Report tail of test & alpha levels
VerbReport(sprintf('Testing null hypothesis that the grand average ERPs in GRP variable''s Bin %d (%s) have a mean of %f microvolts.',bin, ...
    GRP.bin_info(bin).bindesc,p.Results.null_mean),1,VERBLEVEL);
if p.Results.tail==0
    VerbReport(sprintf('Alternative hypothesis is that the ERPs differ from %f (i.e., two-tailed test).',p.Results.null_mean), ...
        1,VERBLEVEL);
elseif p.Results.tail<0,
    VerbReport(sprintf('Alternative hypothesis is that the ERPs are less than %f (i.e., lower-tailed test).',p.Results.null_mean), ...
        1,VERBLEVEL);
else
    VerbReport(sprintf('Alternative hypothesis is that the ERPs are greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ...
        1,VERBLEVEL);
end

%% Optionally reset random number stream to reproduce a previous test
if isempty(p.Results.reproduce_test),
    seed_state=[];
else
    if p.Results.reproduce_test>length(GRP.t_tests),
        error('Value of argument ''reproduce_test'' is too high.  You only have %d permutation tests stored with this GRP variable.',length(GRP.t_tests));
    else
        if isnan(GRP.t_tests(p.Results.reproduce_test).n_perm)
            error('t-test set %d is NOT a permutation test. You don''t need to seed the random number generator to reproduce it.', ...
                p.Results.reproduce_test);
        else
            seed_state=GRP.t_tests(p.Results.reproduce_test).seed_state;
        end
    end
end

%% Compute the permutation test
if strcmpi(GRP.bin_info(bin).op,'(A+B)/n'),
    %one sample t-test
    VerbReport('Performing one sample/repeated measures t-tests.',1,VERBLEVEL);
    erpsAB=zeros(length(use_chans),size(erpsA,2),n_subA+n_subB);
    erpsAB(:,:,1:n_subA)=erpsA-p.Results.null_mean;
    erpsAB(:,:,(n_subA+1):(n_subA+n_subB))=erpsB-p.Results.null_mean;
    [prm_pval, data_t, clust_info, seed_state, est_alpha]=clust_perm1(erpsAB, ...
        chan_hood,p.Results.n_perm,p.Results.alpha,p.Results.tail,p.Results.thresh_p, ...
        VERBLEVEL,seed_state,0);
else
    %independent samples t-test
    VerbReport('Performing independent samples t-tests.',1,VERBLEVEL);
    [prm_pval, data_t, clust_info, seed_state, est_alpha]=clust_perm2(erpsA-p.Results.null_mean, ...
        erpsB,chan_hood,p.Results.n_perm,p.Results.alpha,p.Results.tail,p.Results.thresh_p,VERBLEVEL,seed_state,0);
end

%% Command line summary of results
if p.Results.tail>=0,
    %upper or two-tailed test
    n_pos=length(clust_info.pos_clust_pval);
    fprintf('# of positive clusters found: %d\n',n_pos);
    fprintf('# of significant positive clusters found: %d\n',sum(clust_info.pos_clust_pval<est_alpha));
    fprintf('Positive cluster p-values range from %g to %g.\n',min(clust_info.pos_clust_pval), ...
        max(clust_info.pos_clust_pval));
end
if p.Results.tail<=0,
    %lower or two-tailed test
    n_neg=length(clust_info.neg_clust_pval);
    fprintf('# of negative clusters found: %d\n',n_neg);
    fprintf('# of significant negative clusters found: %d\n',sum(clust_info.neg_clust_pval<est_alpha));
    fprintf('Negative cluster p-values range from %g to %g.\n',min(clust_info.neg_clust_pval), ...
        max(clust_info.neg_clust_pval));
end

sig_ids=find(prm_pval<p.Results.alpha);
if isempty(sig_ids),
    fprintf('ERPs are NOT significantly different from zero (alpha=%f) at any time point/window analyzed.\n', ...
        p.Results.alpha);
    fprintf('All p-values>=%f\n',min(min(prm_pval)));
else
    [dummy min_t_id]=min(abs(data_t(sig_ids)));
    min_t=data_t(sig_ids(min_t_id));
    VerbReport(['Smallest significant t-score(s):' num2str(min_t)],1,VERBLEVEL);
    if p.Results.tail
        %one-tailed test
        tw_alpha=1-cdf('t',max(abs(min_t)),n_sub-1);
    else
        %two-tailed test
        tw_alpha=(1-cdf('t',max(abs(min_t)),n_sub-1))*2;
    end
    
    VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL);
    VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.alpha/(size(prm_pval,1)* ...
        size(prm_pval,2))),1,VERBLEVEL);
    sig_tpts=find(sum(prm_pval<p.Results.alpha));
    
    fprintf('Significant differences from zero (in order of earliest to latest):\n');
    max_sig_p=0;
    min_sig_p=2;
    for t=sig_tpts,
        if mean_wind
            %time windows instead of time points
            fprintf('%d to %d ms, electrode(s): ',GRP.time_pts(use_tpts{t}(1)), ...
                GRP.time_pts(use_tpts{t}(end)));
        else
            fprintf('%d ms, electrode(s): ',GRP.time_pts(use_tpts(t)));
        end
        sig_elec=find(prm_pval(:,t)<p.Results.alpha);
        ct=0;
        for c=sig_elec',
            ct=ct+1;
            if prm_pval(c,t)>max_sig_p,
                max_sig_p=prm_pval(c,t);
            end
            if prm_pval(c,t)<min_sig_p,
                min_sig_p=prm_pval(c,t);
            end
            
            if ct==length(sig_elec),
                fprintf('%s.\n',GRP.chanlocs(use_chans(c)).labels);
            else
                fprintf('%s, ',GRP.chanlocs(use_chans(c)).labels);
            end
        end
    end
    fprintf('All significant corrected p-values are between %f and %f\n',max_sig_p,min_sig_p);
end


%Add permutation results to GRP struct
n_t_tests=length(GRP.t_tests);
neo_test=n_t_tests+1;
GRP.t_tests(neo_test).bin=bin;
GRP.t_tests(neo_test).time_wind=time_wind;
GRP.t_tests(neo_test).used_tpt_ids=use_tpts;
n_use_chans=length(use_chans);
include_chans=cell(1,n_use_chans);
for a=1:n_use_chans,
   include_chans{a}=GRP.chanlocs(use_chans(a)).labels; 
end
GRP.t_tests(neo_test).include_chans=include_chans;
GRP.t_tests(neo_test).used_chan_ids=use_chans;
GRP.t_tests(neo_test).mult_comp_method='cluster mass perm test';
GRP.t_tests(neo_test).n_perm=p.Results.n_perm;
GRP.t_tests(neo_test).desired_alphaORq=p.Results.alpha;
GRP.t_tests(neo_test).estimated_alpha=est_alpha;
GRP.t_tests(neo_test).null_mean=p.Results.null_mean;
if mean_wind,
    GRP.t_tests(neo_test).data_t=data_t;
    GRP.t_tests(neo_test).mean_wind='yes';
else
    GRP.t_tests(neo_test).data_t='See GRP.grands_t';
    GRP.t_tests(neo_test).mean_wind='no';
end
GRP.t_tests(neo_test).crit_t=NaN;
GRP.t_tests(neo_test).df=df;
GRP.t_tests(neo_test).adj_pval=prm_pval;
GRP.t_tests(neo_test).fdr_rej=NaN;
GRP.t_tests(neo_test).seed_state=seed_state;
GRP.t_tests(neo_test).clust_info=clust_info;
GRP.t_tests(neo_test).chan_hood=chan_hood;

if strcmpi(p.Results.plot_raster,'yes'),
    sig_raster(GRP,neo_test,'verblevel',0,'use_color','rgb');
end

if mean_wind,
    if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo),
        sig_topo(GRP,neo_test,'units','t','verblevel',0); %t-score topographies
        sig_topo(GRP,neo_test,'units','uV','verblevel',0); %microvolt topographies
    end
else
    if strcmpi(p.Results.plot_gui,'yes'),
        gui_erp('initialize','GNDorGRP',GRP,'t_test',neo_test,'stat','t', ...
            'verblevel',1);
    end
end

if ~isempty(p.Results.output_file)
    [fid msg]=fopen(p.Results.output_file,'w');
    if fid==-1,
        error('Cannot create file %s for writing.  According to fopen.m: %s.', ...
            p.Results.file,msg);
    else
        %Write header column of times
        % Leave first column blank for channel labels   
        if mean_wind,
            for t=1:n_wind
                fprintf(fid,' %d-%d',GRP.time_pts(use_tpts{t}(1)), ...
                    GRP.time_pts(use_tpts{t}(end)));
            end
            
            %write a couple spaces and then write header for t-scores
            fprintf(fid,'  ');
            for t=1:n_wind
                fprintf(fid,' %d-%d',GRP.time_pts(use_tpts{t}(1)), ...
                    GRP.time_pts(use_tpts{t}(end)));
            end
        else
            for t=use_tpts,
                fprintf(fid,' %d',GRP.time_pts(t));
            end
        end
        fprintf(fid,' Milliseconds\n');
        
        % Write channel labels and p-values
        chan_ct=0;
        for c=use_chans,
            chan_ct=chan_ct+1;
            fprintf(fid,'%s',GRP.chanlocs(c).labels);
            for t=1:length(use_tpts),
                fprintf(fid,' %f',prm_pval(chan_ct,t));
            end
            fprintf(fid,' p-value');
            
            if mean_wind,
                %write a couple spaces and then write t-scores if mean amp
                %in time windows used
                fprintf(fid,' ');
                for t=1:n_wind
                    fprintf(fid,' %f',data_t(chan_ct,t));
                end
                fprintf(fid,' t-score \n');
            else
                fprintf(fid,'\n');
            end
        end
        
        % Write permutation test details
        fprintf(fid,'Experiment: %s\n',GRP.exp_desc);
        fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals: %f\n',bin,p.Results.null_mean);
        fprintf(fid,'#_of_time_windows: %d\n',n_wind);
        fprintf(fid,'#_of_permutations: %d\n',p.Results.n_perm);
        fprintf(fid,'Tail_of_test: ');
        if ~p.Results.tail,
            fprintf(fid,'Two_tailed\n');
        elseif p.Results.tail>0
            fprintf(fid,'Upper_tailed\n');
        else
            fprintf(fid,'Lower_tailed\n');
        end
        fprintf(fid,'Degrees_of_freedom: %d\n',df);
        fprintf(fid,'Alpha_level: %f\n',p.Results.alpha);
        fprintf(fid,'Cluster_inclusion_p_value_threshold: %f\n',p.Results.thresh_p);
        if isscalar(p.Results.chan_hood)
            fprintf(fid,'Max_spatial_neighbor_distance: %f\n',p.Results.chan_hood);
        end
        
        for grp=[grpA grpB],
            % # of participants and filenames
            if grp==grpA,
                use_subs=use_subsA;
                fprintf(fid,'GND_fname_GroupA: %s\n',GRP.GND_fnames{grp});
                fprintf(fid,'#_of_participants_GroupA: %d\n',n_subA);
                fprintf(fid,'Participant_names_GroupA: \n');
            else
                use_subs=use_subsB;
                fprintf(fid,'GND_fname_GroupB: %s\n',GRP.GND_fnames{grp});
                fprintf(fid,'#_of_participants_GroupB: %d\n',n_subB);
                fprintf(fid,'Participant_names_GroupB: \n');
            end
            for s=1:length(use_subs),
                fprintf(fid,'%s\n',GRP.indiv_subnames{grp}{use_subs(s)});
            end
        end
    end
    fclose(fid);
end

if ~strcmpi(p.Results.save_GRP,'no'),
    GRP=save_matmk(GRP);
end

%
%% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function bool=str2bool(str)
%function bool=str2bool(str)

if ischar(str),
    if strcmpi(str,'yes') || strcmpi(str,'y')
        bool=1;
    else
        bool=0;
    end
else
   bool=str; 
end