% STATCOND  - compare two or more data conditions statistically using 
%               standard parametric or nonparametric permutation-based ANOVA 
%               (1-way or 2-way) or t-test methods. Parametric testing uses 
%               FCDF from the Matlab Statistical Toolbox.
% Usage:
%          >> [stats, df, pvals, surrog] = statcond( data, 'key','val'... );
%
% Inputs:
%   data       = one-or two-dimensional cell array of data matrices. 
%                   For nonparametric, permutation-based testing, the 
%                last dimension of the data arrays (which may be of up to 
%                4 dimensions) is permuted across conditions, either in
%                a 'paired' fashion (not changing the, e.g., subject or
%                trial order in the last dimension) or in an umpaired 
%                fashion (not respecting this order). If the number of 
%                elements in the last dimension is not the same across 
%                conditions, the 'paired' option is turned 'off'.  Note: 
%                All other dimensions MUST be constant across conditions. 
%                   For example, consider a (1,3) cell array of matrices 
%                of size (100,20,x) each holding a (100,20) time/frequency 
%                transform from each of x subjects. Only the last dimension 
%                (here x, the number of subjects) may differ across the 
%                three conditions. 
%                   The test used depends on the size of the data array input.
%                When the data cell array has 2 columns and the data are 
%                paired, a paired t-test is performed; when the data are 
%                unpaired, an unpaired t-test is performed. If 'data' 
%                has only one row (paired or unpaired) and more than 2 
%                columns, a one-way ANOVA is performed. If the data cell 
%                array contains several rows and columns, and the data is
%                paired, a two-way repeated measure ANOVA is performed. 
%                NOTE THAT IF THE DATA is unpaired, EEGLAB will use a 
%                balanced 1 or 2 way ANOVA and parametric results might not 
%                be meaningful (bootstrap and permstatcondutation should be fine).
%
% Optional inputs:
%   'paired'   = ['on'|'off'] pair the data array {default: 'on' unless 
%                the last dimension of data array is of different lengths}.
%   'method'   = ['perm'|'bootstrap'|'param'] method for computing the p-values:
%                 'param' or 'parametric' = parametric testing (standard ANOVA
%                                           or t-test); 
%                 'perm' or 'permutation' = non-parametric testing using 
%                                           surrogate data
%                 'bootstrap' = non-parametric bootstrap 
%                  made by permuting the input data {default: 'param'}
%   'naccu'    = [integer] Number of surrogate data copies to use in 'perm' 
%                 or 'bootstrap' method estimation (see above) {default: 200}.
%   'verbose'  = ['on'|'off'] print info on the command line {default: 'on'}.
%   'variance' = ['homegenous'|'inhomogenous'] this option is exclusively
%                for parametric statistics using unpaired t-test. It allows
%                to compute a more accurate value for the degree of freedom
%                using the formula for inhomogenous variance (see
%                ttest2_cell function). Default is 'inhomegenous'.
%   'surrog'   = surrogate data array (see output).
%   'stats'    = F- or T-value array (see output).
%   'tail'     = ['one'|'two'] run one-tailed (F-test) or two tailed
%                (T-test). This option is only relevant when using the
%                'surrog' input. Otherwise it is ignored.
%   'forceanova' = ['on'|'off'] force the use of ANOVA calculation even
%                for 2x1 designs. Default is 'off'.
%   'alpha'    = [float] p-value threshold value. Allow returning
%                confidence intervals and mask (requires structoutput below).
%   'structoutput' = ['on'|'off'] return an output structure instead of 
%                the regular output. Allow to output mask and confidence
%                intervals.
%   'cluster'  = ['on'|'off'] cluster correction for multiple comparison.
%                Only functional when alpha is set (if alpha is NaN, it
%                sets it to 0.05) and for 1-way Anova or t-test. Note that
%                this type of cluster correction only work for matrices,
%                not for channels. When correcting in channel space, use
%                the statcondfieldtrip function instead.
%
% Legacy parameters:
%   'threshold' - now 'alpha'
%   'mode'      - now 'method'
%
% Outputs:
%   stats      = F- or T-value array of the same size as input data without 
%                the last dimension. A T value is returned only when the data 
%                includes exactly two conditions.
%   df         = degrees of freedom, a (2,1) vector, when F-values are returned
%   pvals      = array of p-values. Same size as input data without the last
%                data dimension. All returned p-values are two-tailed.
%   surrog     = surrogate statistic values (same size as stats output with the last 
%                dimension filled with a number ('naccu') of surrogate data sets.
%
% Important note: When a two-way ANOVA is performed, outputs are cell arrays
%                 with three elements: output(1) = row effects; 
%                 output(2) = column effects; output(3) = interactions
%                 between rows and columns.
%
% Examples:
%      >> a = { rand(1,10) rand(1,10)+0.5 }; % pseudo 'paired' data vectors
%         [t df pvals] = statcond(a);        % perform paired t-test
%           pvals =                  
%              5.2807e-04 % standard t-test probability value
%         % Note: for different RAND outputs, results will differ.
%
%         [t df pvals surog] = statcond(a, 'method', 'perm', 'naccu', 2000); 
%           pvals =
%              0.0065 % nonparametric t-test using 2000 permuted data sets
%
%         a = { rand(2,11) rand(2,10) rand(2,12)+0.5 }; % pseudo 'unpaired' 
%         [F df pvals] = statcond(a); % perform an unpaired ANOVA 
%           pvals =
%              0.00025 % p-values for difference between columns 
%              0.00002 % for each data row
%
%         a = { rand(3,4,10) rand(3,4,10) rand(3,4,10); ...
%               rand(3,4,10) rand(3,4,10) rand(3,4,10)+0.5 }; 
%         % pseudo (2,3)-condition data array, each entry containing 
%         %                                    ten (3,4) data matrices
%         [F df pvals] = statcond(a);  % perform a paired 2-way ANOVA 
%         % Output:
%           pvals{1} % a (3,4) matrix of p-values; effects across rows
%           pvals{2} % a (3,4) matrix of p-values; effects across columns 
%           pvals{3} % a (3,4) matrix of p-values; interaction effects
%                                      % across rows and columns
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005-
%         With thanks to Robert Oostenveld for fruitful discussions 
%         and advices on this function.
%
% See also: ANOVA1_CELL, ANOVA2_CELL, ANOVA2RM_CELL, FCDF

% perform a paired t-test
% -----------------------
% a = { rand(2,10) rand(2,10) };
% [t df pval] = statcond(a); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p
% [h p t stat] = ttest( a{1}(2,:), a{2}(2,:)); p
%
% compare significance levels
% --------------------------
% a = { rand(1,10) rand(1,10) }; 
% [F df pval] = statcond(a, 'method', 'perm', 'naccu', 200); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ ori_vals, df, pvals, surrogval ] = statcond( data, varargin );
    
    if nargin < 1
        help statcond;
        return;
    end
    try, warning('off', 'MATLAB:divideByZero'); catch, end;    
    
    if exist('finputcheck')
        g = finputcheck( varargin, { 'naccu'      'integer'   [1 Inf]             200;
                                     'method'     'string'    { 'param','parametric','perm','permutation','bootstrap' }  'param';
                                     'mode'       'string'    { }                 '';
                                     'paired'     'string'    { 'on','off','auto' }      'auto'; 
                                     'surrog'     { 'real','cell' }      []       []; 
                                     'stats'      { 'real','cell' }      []       []; 
                                     'structoutput' 'string'  { 'on','off' }      'off'; 
                                     'forceanova'   'string'  { 'on','off' }      'off'; 
                                     'arraycomp'  'string'    { 'on','off' }      'on'; 
                                     'cluster'    'string'    { 'on','off' }      'off'; 
                                     'alpha'      'real'      []                  NaN;
                                     'tail'       'string'    { 'one','both','upper','lower'}    'both'; 
                                     'variance'   'string'    { 'homogenous','inhomogenous' }    'inhomogenous'; 
                                     'returnresamplingarray' 'string'    { 'on','off' }      'off'; 
                                     'verbose'    'string'    { 'on','off' }      'on' }, 'statcond');
        if ischar(g), error(g); end
    else
        g = struct(varargin{:});
        if ~isfield(g, 'naccu'),     g.naccu = 200; end
        if ~isfield(g, 'method'),    g.method  = 'param'; end
        if ~isfield(g, 'paired'),    g.paired = 'auto'; end
        if ~isfield(g, 'surrog'),    g.surrog = []; end
        if ~isfield(g, 'orivals'),   g.orivals = []; end
        if ~isfield(g, 'arraycomp'), g.arraycomp = 'on'; end
        if ~isfield(g, 'verbose'),   g.verbose = 'on'; end
        if ~isfield(g, 'tail'),      g.tail = 'both'; end
        if ~isfield(g, 'variance'),  g.variance = 'homogenous'; end
        if ~isfield(g, 'structoutput'), g.structoutput = 'on'; end
        if ~isfield(g, 'returnresamplingarray'),   g.returnresamplingarray = 'off'; end
    end
    if ~isempty(g.mode), g.method = g.mode; end
    if size(data,2) == 1, data  = transpose(data); end % cell array transpose
    
    % other settings
    % --------------
    if strcmpi(g.method, 'parametric'), g.method = 'param'; end
    if strcmpi(g.method, 'permutation'), g.method = 'perm'; end
    if strcmpi(g.verbose, 'on'), verb = 1; else verb = 0; end
    if strcmp(g.method, 'param' ) && exist('fcdf') ~= 2
      myfprintf('on',['statcond(): parametric testing requires fcdf() \n' ...
               '            from the Matlab StatsticaL Toolbox.\n' ...
               '            Running nonparametric permutation tests\n.']);
      g.method = 'perm';
    end
    g.naccu = round(g.naccu);
    
    % pairing
    % -------
    paired{1} = 'on';
    if size(data,1) == 1 && length(unique(cellfun('size', data, ndims(data{1}) ))) > 1
        paired{1} = 'off';
    else
        paired{1} = 'on';
        for iCol = 1:size(data,2)
            if length(unique(cellfun('size', data(:,iCol), ndims(data{1}) ))) > 1
                paired{1} = 'off';
            end
        end
        paired{2} = 'on';
        for iRow = 1:size(data,1)
            if length(unique(cellfun('size', data(iRow,:), ndims(data{1}) ))) > 1
                paired{2} = 'off';
            end
        end
    end
    if length(paired) > 1
        if (strcmpi(paired{1}, 'off') && strcmpi(paired{2}, 'on')) || ...
                (strcmpi(paired{1}, 'on') && strcmpi(paired{2}, 'off'))
            myfprintf(g.verbose, 'Possible mixed paired and unpaired independent variables, using balanced Anova (which assumes all unpaired)\n');
            paired{1} = 'off';
        end
    end
    if strcmpi(g.paired, 'auto')
        g.paired = paired{1};
    else
        if strcmpi(g.paired, 'on') && strcmpi(paired{1}, 'off')
            myfprintf(g.verbose, 'You set to use paired statistics but the number of cases differs\n');
            g.paired = 'off';
        end
    end
            
    % reshape matrices
    % ----------------
    nd = size(data{1});
    nd = nd(1:end-1);
    for index = 1:prod(size(data))
        data{index} = reshape(data{index}, [prod(nd) size(data{index},myndims(data{index}))]);
    end
    
    if ~strcmpi(g.method, 'param') && isempty(g.surrog)
         tmpsize   = size(data{1});
         surrogval = zeros([ tmpsize(1:end-1) g.naccu ], 'single');
    else surrogval = [];
    end
    
    % check for NaNs or Inf
    % ---------------------
    for iDat = 1:length(data(:))
        if any(isnan(reshape(data{iDat}, prod(size(data{iDat})),1))) || ...
                any(isinf(reshape(data{iDat}, prod(size(data{iDat})),1)))
            error('Statcond: One of the input array contains NaNs or Infinite values');
        end
    end
        
    % bootstrap flag
    % --------------
    if strcmpi(g.method, 'bootstrap'), bootflag = 1;
    else                               bootflag = 0;
    end
    
    if isempty(g.surrog)
        
        % return resampling array
        % -----------------------
        if strcmpi(g.returnresamplingarray, 'on')
            if strcmpi(g.arraycomp, 'on')
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
            else
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
            end
            return;
        end
        
        % text output
        % -----------
        myfprintf(verb,'%d x %d, ', size(data,1), size(data,2));
        if strcmpi(g.paired, 'on')
             myfprintf(verb,'paired data, ');
        else myfprintf(verb,'unpaired data, ');
        end
        if size(data,1) == 1 && size(data,2) == 2
             myfprintf(verb,'computing T values\n');
        else myfprintf(verb,'computing F values\n');
        end
        if size(data,1) > 1 
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 2-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 2-way ANOVA (not suitable for parametric testing, only bootstrap)\n');
            end
        elseif size(data,2) > 2
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 1-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 1-way ANOVA (equivalent to Matlab anova1)\n');
            end
        else
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using paired t-test\n');
            else myfprintf(verb,'Using unpaired t-test\n');
            end
        end
        if ~strcmpi(g.method, 'param')
            if bootflag, myfprintf(verb,'Bootstraps (of %d):', g.naccu);
            else         myfprintf(verb,'Permutations (of %d):', g.naccu);
            end
        end
    end
    
    tail = g.tail;
    if isempty(g.surrog)
        if size(data,1) == 1 % only one row

            if size(data,2) == 2 && strcmpi(g.forceanova, 'off')

                % paired t-test (very fast)
                % -------------
                [ori_vals, df] = ttest_cell_select(data, g.paired, g.variance);

                if strcmpi(g.method, 'param')
                    
                    % Check if exist tcd.m file from the Statistics Toolbox (Bug 1352 )
                    if exist('tcdf','file') == 2  && license('test', 'Statistics_Toolbox')
                        pvals = 2*tcdf(-abs(ori_vals), df);
                    else
                        pvals = 2*mytcdf(-abs(ori_vals), df);
                    end
                    
                    pvals = reshape(pvals, size(pvals));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = ttest_cell_select( res, g.paired, g.variance);
                        catch
                           lasterr
                           myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computation');
                           g.arraycomp = 'off';
                        end
                    end
                    if strcmpi(g.arraycomp, 'off')
                        [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            res = surrogdistrib( {}, 'precomp', precomp);
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end
                            if myndims(res{1}) == 1
                                 surrogval(index)     = ttest_cell_select(res, g.paired, g.variance);
                            else surrogval(:,index)   = ttest_cell_select(res, g.paired, g.variance);
                            end
                        end
                    end
                end
            else
                % one-way ANOVA (paired) this is equivalent to unpaired t-test
                % -------------
                tail = 'one';
                [ori_vals, df] = anova1_cell_select( data, g.paired );
                if strcmpi(g.method, 'param')
                    pvals = 1-fcdf(ori_vals, df(1), df(2));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');                        
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = anova1_cell_select( res, g.paired );
                        catch,
                            myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computing');
                            g.arraycomp = 'off';
                        end
                    end
                    if strcmpi(g.arraycomp, 'off')
                        [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end

                            res = surrogdistrib( {}, 'precomp', precomp);
                            if myndims(data{1}) == 1
                            	 surrogval(index)     = anova1_cell_select( res, g.paired );
                            else surrogval(:,index)   = anova1_cell_select( res, g.paired );
                            end
                        end
                    end
                end
            end
        else
            % two-way ANOVA (paired or unpaired)
            % ----------------------------------
            tail = 'one';
            [ ori_vals{1}, ori_vals{2}, ori_vals{3}, df{1}, df{2}, df{3} ] = anova2_cell_select( data, g.paired );
            if strcmpi(g.method, 'param')
                pvals{1} = 1-fcdf(ori_vals{1}, df{1}(1), df{1}(2));
                pvals{2} = 1-fcdf(ori_vals{2}, df{2}(1), df{2}(2));
                pvals{3} = 1-fcdf(ori_vals{3}, df{3}(1), df{3}(2));
            else
                surrogval = { surrogval surrogval surrogval };
                dataori   = data;
                if strcmpi(g.arraycomp, 'on')
                    try
                        myfprintf(verb,'...');
                        res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                        [ surrogval{1}, surrogval{2}, surrogval{3} ] = anova2_cell_select( res, g.paired );
                    catch
                        myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computing');
                        g.arraycomp = 'off';
                    end
                end
                if strcmpi(g.arraycomp, 'off')
                    [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                    for index = 1:g.naccu
                        if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                        if mod(index, 100) == 0, myfprintf(verb,'\n'); end

                        res = surrogdistrib( {}, 'precomp', precomp);
                        if myndims(data{1}) == 1
                         	 [ surrogval{1}(index),     surrogval{2}(index),     surrogval{3}(index)     ] = anova2_cell_select( res, g.paired );
                        else [ surrogval{1}(:,index),   surrogval{2}(:,index),   surrogval{3}(:,index)   ] = anova2_cell_select( res, g.paired );
                        end
                    end
                end
            end
        end
        myfprintf(verb,'\n');
    else
        surrogval = g.surrog;
        ori_vals  = g.stats;
        df        = [];
    end
    
    % compute p-values
    % ----------------
    if ~strcmpi(g.method, 'param')
        if iscell( surrogval )
            pvals{1} = stat_surrogate_pvals(surrogval{1}, ori_vals{1}, tail);
            pvals{2} = stat_surrogate_pvals(surrogval{2}, ori_vals{2}, tail);
            pvals{3} = stat_surrogate_pvals(surrogval{3}, ori_vals{3}, tail);
        else
            if strcmpi(g.cluster, 'on')
                if isnan(g.alpha)
                    g.alpha = 0.05;
                    disp('Alpha value set to 0.05 automatically');
                end
                if size(surrogval,3) > 1 || size(surrogval,4) > 1 || isnan(g.alpha)
                    error('Cluster method not implemented for alpha NaN or for more than 2 dims');
                else                    
                    tmpPvals = 2*tcdf(-abs(surrogval), size(surrogval,1)) < g.alpha;
                    largestCluster = zeros(1, size(surrogval,2));
                    signifPos = zeros(1, size(surrogval,1));
                    for iSurog = 1:size(surrogval,2)
                        signifPos = tmpPvals(:,iSurog);
                        while any(signifPos)
                            tmpInd = find(signifPos);
                            currentInd = tmpInd(1);
                            currentCount = 0;
                            while currentInd+currentCount <= length(signifPos) && signifPos(currentInd+currentCount) == 1
                                signifPos(currentInd+currentCount) = 0;
                                currentCount = currentCount+1;
                            end
                            if currentCount > largestCluster(iSurog), largestCluster(iSurog) = currentCount; end
                        end
                    end
                    largestCluster = sort(largestCluster);
                    thresholdVal = largestCluster(round(length(largestCluster)*(1-g.alpha)));
                    
                    % assign p-val status
                    signifPos = 2*tcdf(-abs(ori_vals), length(ori_vals)) < g.alpha;
                    pvals     = ones(size(signifPos));
                    while any(signifPos)
                        tmpInd = find(signifPos);
                        currentInd = tmpInd(1);
                        currentCount = 0;
                        while currentInd+currentCount <= length(signifPos) && signifPos(currentInd+currentCount) == 1
                            signifPos(currentInd+currentCount) = 0;
                            currentCount = currentCount+1;
                        end
                        if currentCount >= thresholdVal, pvals(currentInd:currentInd+currentCount-1) = 0; end
                    end
                    
                end
            else
                pvals = stat_surrogate_pvals(surrogval, ori_vals, tail);
            end
        end
        try, warning('on', 'MATLAB:divideByZero'); catch, end
    end

    [ ori_vals, pvals ] = reshape_results( nd, ori_vals, pvals);
    [ surrogval ]       = reshape_results( [nd g.naccu], surrogval);
    
    % confidence intervals
    % --------------------
    if ~isnan(g.alpha)
        outputstruct.ci = stat_surrogate_ci(surrogval, g.alpha, tail);
        if strcmpi(g.structoutput, 'off')
            disp('Warning: returning confidence interval requires an output structure');
        end
        if iscell(pvals)
            for ind = 1:length(pvals)
                outputstruct.mask{ind} = pvals{ind} < g.alpha;
            end
        else
            outputstruct.mask = pvals < g.alpha;
        end
    end
    
    % create a structure for outputting values
    % ---------------------------------------
    if strcmpi(g.structoutput, 'on')
        outputstruct.method = g.method;
        outputstruct.pval   = pvals;
        outputstruct.df     = df;
        outputstruct.surrog = surrogval;
        if length(data(:)) == 2
             outputstruct.t = ori_vals;
        else outputstruct.f = ori_vals;
        end
        outputstruct.stat   = ori_vals;
        ori_vals = outputstruct;
    end
       
% compute ANOVA 2-way
% -------------------
function [f1, f2, f3, df1, df2, df3] = anova2_cell_select( res, paired)
    if strcmpi(paired,'on')
        [f1, f2, f3, df1, df2, df3] = anova2rm_cell( res );
    else
        [f1, f2, f3, df1, df2, df3] = anova2_cell( res );
    end
    
% compute ANOVA 1-way
% -------------------
function [f, df] = anova1_cell_select( res, paired)
    if strcmpi(paired,'on')
        [f, df] = anova1rm_cell( res );
    else
        [f, df] = anova1_cell( res );
    end

% compute t-test
% -------------------
function [t, df] = ttest_cell_select( res, paired, homogenous)
    if strcmpi(paired,'on')
        [t, df] = ttest_cell( res{1}, res{2});
    else
        [t, df] = ttest2_cell( res{1}, res{2}, homogenous);
    end

% function to compute the number of dimensions
% --------------------------------------------
function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1
            val = 2;
        elseif size(a,2) == 1
            val = 1;
        else
            val = 2;
        end
    end

% function for verbose messages
% -----------------------------
function myfprintf(verb, varargin)
    if verb
        fprintf(varargin{:});
    end

% function to replace tcdf
% ------------------------
function p = mytcdf(x,v)

if length(v) == 1
    v = repmat(v, size(x));
end

x2 = x.^2;
inds1 = (v < x2);
inds2 = (v >= x2);
if any(inds1(:)), p(inds1) = betainc(v(inds1) ./ (v(inds1) + x2(inds1)), v(inds1)/2, 0.5, 'lower') / 2; end
if any(inds2(:)), p(inds2) = betainc(x2(inds2) ./ (v(inds2) + x2(inds2)), 0.5, v(inds2)/2, 'upper') / 2; end
inds = (x > 0); 
if any(inds)
    p(inds) = 1 - p(inds);
end

inds = (v > 1e7);
if any(inds(:)), p(inds) = normcum(x(inds)); end

p(x == 0) = 0.5;
if isempty(p)
    p = ones(size(x));
else
    p = reshape(p, size(x));
end
function [p] = normcum(z)
p = 0.5 * erfc(-z ./ sqrt(2));

% reshape results
% ---------------
function varargout = reshape_results(nd, varargin)
    if length(varargin) > 1
        for index = 1:length(varargin)
            varargout{index} = reshape_results(nd, varargin{index});
        end
    elseif iscell(varargin{1})
        for index = 1:length(varargin{1})
            varargout{1}{index} = reshape_results(nd, varargin{1}{index});
        end
    else
        if ~isempty(varargin{1})
            if length(nd) == 1, nd = [ nd 1 ]; end
            varargout{1} = reshape(varargin{1}, nd);
        else varargout{1} = [];
        end
    end
