% BOOTSTAT - accumulate surrogate data to assess significance by permutation of some 
%              measure of two input variables. 
%
%              If 'distfit','on', fits the psd with a 4th-order polynomial using the 
%              data kurtosis, as in Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., 
%              Mykkytka, E.F. "A probability distribution and its uses in fitting data." 
%              Technometrics, 21:201-214, 1979.
% Usage:
%            >> [rsignif,rboot] = bootstat( { arg1 arg2 ...}, formula, varargin ...);
% Inputs:
%    arg1    - [array] 1-D, 2-D or 3-D array of values
%    arg2    - [array] 1-D, 2-D or 3-D array of values
%    formula - [string] formula to compute the given measure. Takes arguments
%                   'arg1', 'arg2' as inputs and 'res' (result, by default) as output.
%                   For data arrays of more than 1 dimension, the formula must be iterative 
%                   so that shuffling can occur at each step while scanning the last 
%                   array dimension.  Examples:
%                   'res = arg1 - arg2'              % difference of two 1-D data arrays
%                   'res = mean( arg1 .* arg2)'      % mean projection of two 1-D data arrays
%                   'res = res + arg1 .* conj(arg2)' % iterative, for use with 2|3-D arrays
% Optional inputs:
%   'boottype '   - ['rand'|'shuffle']
%                   'rand'  = do not shuffle data. Only flip polarity randomly (for real 
%                             number) or phase (for complex numbers).
%                   'shuffle' = shuffle values of first argument (see two options below). 
%                   Default.
%   'shuffledim'  - [integer] indices of dimensions to shuffle. For instance, [1 2] will
%                   shuffle the first two dimensions. Default is to shuffle along 
%                   dimension 2.
%   'shufflemode' - ['swap'|'regular'] shuffle mode. Either swap dimensions (for instance
%                   swap rows then columns if dimension [1 2] are selected) or shuffle 
%                   in each dimension independently (slower). If only one dimension is 
%                   selected for shuffling, this option does not change the result.
%   'randmode'    - ['opposite'|'inverse'] randomize sign (or phase for complex number,
%                   or randomly set half the value to reference.
%   'alpha'       - [real] significance level (between 0 and 1) {default 0.05}.
%   'naccu'       - [integer] number of exemplars to accumulate {default 200}.
%   'bootside'    - ['both'|'upper'] side of the surrogate distribution to
%                   consider for significance. This parameter affects the size
%                   of the last dimension of the accumulation array ('accres') 
%                   (size is 2 for 'both' and 1 for 'upper') {default: 'both'}.
%   'basevect'    - [integer vector] time vector indices for baseline in second dimension.
%                   {default: all time points}.
%   'rboot'       - accumulation array (from a previous call). Allows faster 
%                   computation of the 'rsignif' output {default: none}.
%   'formulaout'  - [string] name of the computed variable {default: 'res'}.
%   'dimaccu'     - [integer] use dimension in result to accumulate data.
%                   For instance if the result array is size [60x50] and this value is 2,
%                   the function will consider than 50 times 60 value have been accumulated.
%
% Fitting distribution:
%   'distfit'     - ['on'|'off'] fit distribution with known function to compute more accurate 
%                   limits or exact p-value (see 'vals' option). The MATLAB statistical toolbox 
%                   is required. This option is currently implemented only for 1-D data.
%   'vals'        - [float array] significance values. 'alpha' is ignored and 
%                   rsignif returns the p-values. Requires 'distfit' (see above).
%                   This option currently implemented only for 1-D data.
%   'correctp'    - [phat pci zerofreq] parameters for correcting for a biased probability 
%                   distribution (requires 'distfit' above). See help of CORRECTFIT.
% Outputs: 
%    rsignif      - significance arrays. 2 values (low high) for each point (use
%                   'alpha' to change these limits).
%    rboot        - accumulated surrogate data values.
%
% Authors: Arnaud Delorme, Bhaktivedcanta Institute, Mumbai, India, Nov 2004
%
% See also: TIMEF

% NOTE: There is an undocumented parameter, 'savecoher', [0|1]
% HELP TEXT REMOVED:  (Ex: Using option 'both', coherence during baseline would be 
%                      ignored since times are shuffled during each accumulation.

% Copyright (C) 9/2002  Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD
%
% 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.

% *************************************
% To fit the psd with as 4th order polynomial using the distribution kurtosis,
% Reference: Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., Mykkytka, E.F. 
% "A probability distribution and its uses in fitting data." 
% Technimetrics, 1979, 21: 201-214.
% *************************************

function [accarrayout, Rbootout, Rbootout2] = bootstat(oriargs, formula, varargin)
%	nb_points, timesout, naccu, baselength, baseboot, boottype, alpha, rboot);
	
if nargin < 2
	help bootstat;
	return;
end

if ~ischar(formula)
	error('The second argument must be a string formula');
end

g = finputcheck(varargin, ...
                { 'dims'          'integer'  []                       []; ...
                  'naccu'         'integer'  [0 10000]                200; ...
                  'bootside'      'string'   { 'both','upper' }       'both'; ...
                  'basevect'      'integer'  []                       []; ...
				  'boottype'      'string'  { 'rand','shuffle' }      'shuffle'; ...
				  'shufflemode'   'string'  { 'swap','regular' }      'swap'; ...
				  'randmode'      'string'  { 'opposite','inverse' }  'opposite'; ...
				  'shuffledim'    'integer'  [0 Inf]                  []; ...
				  'label'         'string'   []                       formula; ...
				  'alpha'         'real'     [0 1]                    0.05; ...
				  'vals'          'real'     []                       []; ...
				  'distfit'       'string'   {'on','off' }            'off'; ...
				  'dimaccu'       'integer'  [1 Inf]                  []; ...
				  'correctp'      'real'     []                       []; ...
				  'rboot'         'real'     []                       NaN	});
if ischar(g)
	error(g);
end
if isempty(g.shuffledim) && strcmpi(g.boottype, 'rand')
    g.shuffledim = []; 
elseif isempty(g.shuffledim)
    g.shuffledim = 2;
end; 
unitname = '';
if 2/g.alpha > g.naccu 
    if strcmpi(g.distfit, 'off') || ~((size(oriarg1,1) == 1 || size(oriarg1,2) == 1) && size(oriarg1,3) == 1)
        g.naccu = 2/g.alpha; 
        fprintf('Adjusting naccu to compute alpha value');
    end
end
if isempty(g.rboot)
	g.rboot = NaN;
end

% function for bootstrap computation
% ----------------------------------
if ~iscell(oriargs) || length(oriargs) == 1, 
    oriarg1 = oriargs;
    oriarg2 = []; 
else 
    oriarg1 = oriargs{1};
    oriarg2 = oriargs{2};
end
[nb_points times trials] = size(oriarg1);
if times == 1, disp('Warning 1 value only for shuffling dimension'); end

% only consider baseline
% ----------------------
if ~isempty(g.basevect)
    fprintf('\nPermutation statistics baseline length is %d (out of %d) points\n', length(g.basevect), times);
    arg1 = oriarg1(:,g.basevect,:);
    if ~isempty(oriarg2)
        arg2 = oriarg2(:,g.basevect,:);
    end
else
    arg1 = oriarg1;
    arg2 = oriarg2;
end

% formula for accumulation array
% ------------------------------
% if g.dimaccu is not empty, accumulate over that dimension
% of the resulting array to speed up computation
formula = [ 'res=' formula ];
g.formulapost = [ 'if index == 1, ' ...
                  '   if ~isempty(g.dimaccu), ' ...
                  '      Rbootout= zeros([ ceil(g.naccu/size(res,g.dimaccu)) size( res ) ]);' ...
                  '   else,' ...
                  '      Rbootout= zeros([ g.naccu size( res ) ]);' ...
                  '   end;' ...
                  'end,' ...
                  'Rbootout(count,:,:,:) = res;' ...
                  'count = count+1;' ...
                  'if ~isempty(g.dimaccu), ' ...
                  '   index = index + size(res,g.dimaccu);' ...
                  '   fprintf(''%d '', index-1);' ...
                  'else ' ...
                  '   index=index+1;' ...
                  '   if rem(index,10)  == 0, fprintf(''%d '', index); end;' ...
                  '   if rem(index,100) == 0, fprintf(''\n''); end;' ...
                  'end;' ];

% **************************
% case 1: precomputed values
% **************************
if ~isnan(g.rboot)
    Rbootout = g.rboot;
% ***********************************
% case 2: randomize polarity or phase
% ***********************************
elseif strcmpi(g.boottype, 'rand') && strcmpi(g.randmode, 'inverse')
    fprintf('Bootstat function: randomize inverse values\n');
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
    
    % compute random array
    % --------------------
    multarray = ones(size(arg1));
    totlen    = prod(size(arg1));
    if isreal(arg1), 
        multarray(1:round(totlen/2)) = 0;
    end
    for shuff = 1:ndims(multarray)
         multarray = supershuffle(multarray,shuff); % initial shuffling
    end
    if isempty(g.shuffledim), g.shuffledim = 1:ndims(multarray); end
    invarg1 = 1./arg1;
    
    % accumulate
    % ----------
    index = 1;
    count = 1;
    while index <= g.naccu
        for shuff = g.shuffledim
            multarray = supershuffle(multarray,shuff);
        end
        tmpinds = find(reshape(multarray, 1, prod(size(multarray))));
        arg1 = oriarg1;
        arg1(tmpinds) = invarg1(tmpinds);
        eval([ formula ';' ]);
        eval( g.formulapost ); % also contains index = index+1
    end
elseif strcmpi(g.boottype, 'rand') % opposite
    fprintf('Bootstat function: randomize polarity or phase\n');
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
    
    % compute random array
    % --------------------
    multarray = ones(size(arg1));
    totlen    = prod(size(arg1));
    if isreal(arg1), 
        multarray(1:round(totlen/2)) = -1;
    else
        tmparray            = exp(j*linspace(0,2*pi,totlen+1));
        multarray(1:totlen) = tmparray(1:end-1);
    end
    for shuff = 1:ndims(multarray)
         multarray = supershuffle(multarray,shuff); % initial shuffling
    end
    if isempty(g.shuffledim), g.shuffledim = 1:ndims(multarray); end
    
    % accumulate
    % ----------
    index = 1;
    count = 1;
    while index <= g.naccu
        for shuff = g.shuffledim
            multarray = supershuffle(multarray,shuff);
        end
        arg1 = arg1.*multarray;
        eval([ formula ';' ]);
        eval( g.formulapost ); % also contains index = index+1
    end
% ********************************************
% case 3: shuffle vector of only one dimension
% ********************************************
elseif length(g.shuffledim) == 1
    fprintf('Bootstat function: shuffling along dimension %d only\n', g.shuffledim);
    fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);

    index = 1;
    count = 1;
    while index <= g.naccu
        arg1 = shuffleonedim(arg1,g.shuffledim);
        eval([ formula ';' ]);
        eval( g.formulapost );
    end
% ***********************************************
% case 5: shuffle vector along several dimensions
% ***********************************************
else 
    if strcmpi(g.shufflemode, 'swap') % swap mode
        fprintf('Bootstat function: shuffling along dimension %s (swap mode)\n', int2str(g.shuffledim));
        fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
        index = 1;
        count = 1;
        while index <= g.naccu
            for shuff = g.shuffledim
                arg1 = supershuffle(arg1,shuff);
            end
            eval([ formula ';' ]);
            eval( g.formulapost );
        end
    else  % regular shuffling
        fprintf('Bootstat function: shuffling along dimension %s (regular mode)\n', int2str(g.shuffledim));
        fprintf('Processing permutation statistics for %s (naccu=%d):', g.label, g.naccu);
        index = 1;
        count = 1;
        while index <= g.naccu
            for shuff = g.shuffledim
                arg1 = shuffleonedim(arg1,shuff);
            end
            eval([ formula ';' ]);
            eval( g.formulapost );
        end
    end
end
Rbootout(count:end,:,:,:) = [];

% **********************
% assessing significance
% **********************

% get accumulation array
% ----------------------
accarray = Rbootout;
if ~isreal(accarray)
    accarray = sqrt(accarray .* conj(accarray)); % faster than abs()
end
% reshape the output if necessary
% -------------------------------
if ~isempty(g.dimaccu)
    if g.dimaccu+1 == 3
        accarray = permute( accarray, [1 3 2]);
    end
    accarray = reshape( accarray, size(accarray,1)*size(accarray,2), size(accarray,3) );
end
if size(accarray,1) == 1, accarray = accarray'; end; % first dim contains g.naccu

% ******************************************************
% compute thresholds on array not fitting a distribution
% ******************************************************
if strcmpi(g.distfit, 'off')
  
    % compute bootstrap significance level
    % ------------------------------------
    accarray  = sort(accarray,1); % always sort on naccu
    Rbootout2 = accarray;
    i         = round(size(accarray,1)*g.alpha);
    accarray1 = squeeze(mean(accarray(size(accarray,1)-i+1:end,:,:),1));
    accarray2 = squeeze(mean(accarray(1:i                     ,:,:),1));
    if abs(accarray(1,1,1) - accarray(end,1,1)) < abs(accarray(1,1,1))*1e-15
        accarray1(:) = NaN;
        accarray2(:) = NaN;
    end

else
    % *******************
    % fit to distribution 
    % *******************
    sizerboot   = size (accarray);
    accarray1   = zeros(sizerboot(2:end));
    accarray2   = zeros(sizerboot(2:end));
    
    if ~isempty(g.vals{index})
        if ~all(size(g.vals{index}) == sizerboot(2:end) )
            error('For fitting, vals must have the same dimension as the output array (try transposing)');
        end
    end
    
    % fitting with Ramberg-Schmeiser distribution
    % -------------------------------------------
    if ~isempty(g.vals{index}) % compute significance for value
        for index1 = 1:size(accarrayout,1)
            for index2 = 1:size(accarrayout,2)
                accarray1(index1,index2) = 1 - rsfit(squeeze(accarray(:,index1,index2)), g.vals{index}(index1, index2));
                if length(g.correctp) == 2
                    accarray1(index1,index2) = correctfit(accarray1, 'gamparams', [g.correctp 0]); % no correction for p=0
                else
                    accarray1(index1,index2) = correctfit(accarray1, 'gamparams', g.correctp);
                end
            end
        end
    else % compute value for significance
        for index1 = 1:size(accarrayout,1)
            for index2 = 1:size(accarrayout,2)
                [p c l chi2] = rsfit(Rbootout(:),0);
                pval = g.alpha;   accarray1(index1,index2) = l(1) + (pval.^l(3) - (1-pval).^l(4))/l(2);
                pval = 1-g.alpha; accarray2(index1,index2) = l(1) + (pval.^l(3) - (1-pval).^l(4))/l(2);        
            end
        end
    end

    % plot results 
    % -------------------------------------
    % figure;
    % hist(abs(Rbootout)); tmpax = axis;
    % hold on; 
    % valcomp = linspace(min(abs(Rbootout(:))), max(abs(Rbootout(:))), 100);
    % normy = normpdf(valcomp, mu, sigma);
    % plot(valcomp, normy/max(normy)*tmpax(4), 'r');
    % return;
end

% set output array: backward compatible
% -------------------------------------
if strcmpi(g.bootside, 'upper'); % only upper significance
    accarrayout = accarray1;
else 
    if size(accarray1,1) ~= 1 && size(accarray1,2) ~= 1
        accarrayout        = accarray2;
        accarrayout(:,:,2) = accarray1;
    else
        accarrayout = [ accarray2(:) accarray1(:) ];
    end
end
accarrayout = squeeze(accarrayout);
if size(accarrayout,1) == 1 && size(accarrayout,3) == 1, accarrayout = accarrayout'; end

% better but not backward compatible
% ----------------------------------
% accarrayout = { accarray1 accarray2 }; 
    
return;

    % fitting with normal distribution (deprecated)
    % --------------------------------
    [mu sigma] = normfit(abs(Rbootout(:)));
    accarrayout = 1 - normcdf(g.vals, mu, sigma); % cumulative density distribution
                                        % formula of normal distribution
                                        % y = 1/sqrt(2) * exp( -(x-mu).^2/(sigma*sigma*2) ) / (sqrt(pi)*sigma);
% % Gamma and Beta fits:
% elseif strcmpi(g.distfit, 'gamma')
%  [phatgam pcigam] = gamfit(abs(Rbootout(:)));
%  gamy = gampdf(valcomp, phatgam(1), pcigam(2))
%  p = 1 - gamcdf(g.vals, phatgam(1), pcigam(2)); % cumulative density distribution
% elseif strcmpi(g.distfit, 'beta')
%  [phatbeta pcibeta] = betafit(abs(Rbootout(:)));
%  betay = betapdf(valcomp, phatbeta(1), pcibeta(1));
%  p = 1 - betacdf(g.vals, phatbeta(1), pcibeta(1)); % cumulative density distribution
% end


    if strcmpi(g.distfit, 'off')
        tmpsort = sort(Rbootout);
        i = round(g.alpha*g.naccu);
        sigval = [mean(tmpsort(1:i)) mean(tmpsort(g.naccu-i+1:g.naccu))];
        if strcmpi(g.bootside, 'upper'), sigval = sigval(2); end
        accarrayout = sigval;
    end

% this shuffling preserve the number of -1 and 1
% for columns and rows (assuming matrix size is multiple of 2
% -----------------------------------------------------------
function array = supershuffle(array, dim)
    if size(array, 1) == 1 || size(array,2) == 1
        array = shuffle(array);
        return;
    end
    if size(array, dim) == 1, return; end

    if dim == 1
        indrows = shuffle(1:size(array,1));
        for index = 1:2:length(indrows)-rem(length(indrows),2) % shuffle rows
            tmparray                    = array(indrows(index),:,:);
            array(indrows(index),:,:)   = array(indrows(index+1),:,:);
            array(indrows(index+1),:,:) = tmparray;
        end
    elseif dim == 2
        indcols = shuffle(1:size(array,2));
        for index = 1:2:length(indcols)-rem(length(indcols),2) % shuffle columns
            tmparray                    = array(:,indcols(index),:);
            array(:,indcols(index),:)   = array(:,indcols(index+1),:);
            array(:,indcols(index+1),:) = tmparray;
        end
    else
        ind3d = shuffle(1:size(array,3));
        for index = 1:2:length(ind3d)-rem(length(ind3d),2) % shuffle columns
            tmparray                  = array(:,:,ind3d(index));
            array(:,:,ind3d(index))   = array(:,:,ind3d(index+1));
            array(:,:,ind3d(index+1)) = tmparray;
        end
    end

% shuffle one dimension, one row/columns at a time
% -----------------------------------------------
function array = shuffleonedim(array, dim)
    if size(array, 1) == 1 || size(array,2) == 1
        array = shuffle(array, dim);
    else
        if dim == 1
            for index1 = 1:size(array,3)
                for index2 = 1:size(array,2)
                    array(:,index2,index1) = shuffle(array(:,index2,index1));
                end
            end
        elseif dim == 2
            for index1 = 1:size(array,3)
                for index2 = 1:size(array,1)
                    array(index2,:,index1) = shuffle(array(index2,:,index1));
                end
            end
        else
            for index1 = 1:size(array,1)
                for index2 = 1:size(array,2)
                    array(index1,index2,:) = shuffle(array(index1,index2,:));
                end
            end
        end
    end
