function [stat] = sourcestatistics_randomization(cfg, source)

% SOURCESTATISTICS_RANDOMIZATION performs statistics on the output of SOURCEANALYSIS,
% provided the randomization-method has been applied. It will be called from within the
% SOURCESTATISTICS function, when cfg.method is set to 'randomization'. See also SOURCESTATISTICS
% for additional information.
%
% Use as
%  [stat] = sourcestatistics(cfg,source)
% where cfg is a structure containing the configuration details, and source is a structure containing
% the fields trialA and trialB, obtained with the 'randomization'-method in sourceanalysis.
%
% The function performs a non-parametric randomization test for comparing the 
% means of two randomized source reconstructions A and B. 
%
% The following parameters are mandatory:
%    cfg.method           = 'randomization'
%    cfg.parameter        = can be 'coh', 'pow', 'nai', 'noise'
%
% The following parameters are optional:
% The computation of a one- or two-tailed test can be specified with
%   cfg.tail = 0 (test A~=B, default), 1 (test A>B), -1 (test A<B)
% The algorithm can correct for multiple comparisons by comparing against the
% distribution of maximum differences over all sources generated by randomization
%   cfg.correctm = 'yes' or 'no' (default)
% The algorithm z-transforms the tested parameter based on the
% rondomization variance estimate for each voxel:
%   cfg.ztransform = 'yes' or 'no' (default)
% The algorithm either tests the difference between A and B ('difference') or 
% the difference between A and B normalized by B ('relchange') for each source:
%   cfg.comparestat = 'difference' or 'relchange' (default: 'difference')
%
% See also SOURCESTATISTICS, SOURCEANALYSIS 

% Copyright (C) 2003, Robert Oostenveld
% Copyright (C) 2004, Markus Siegel
%
% $Log: sourcestatistics_randomization.m,v $
% Revision 1.10  2006/01/30 12:51:38  roboos
% fixed declaration of function name
%
% Revision 1.9  2005/06/29 12:46:29  roboos
% the following changes were done on a large number of files at the same time, not all of them may apply to this specific file
% - added try-catch around the inclusion of input data configuration
% - moved cfg.version, cfg.previous and the assignment of the output cfg to the end
% - changed help comments around the configuration handling
% - some changes in whitespace
%
% Revision 1.8  2005/06/03 08:58:07  roboos
% transfer homogenous transformation matrix from input to output (if present)
% added an extra method of counting the number of dipoles/voxels in the input source structure
%
% Revision 1.7  2005/05/12 12:17:23  roboos
% added a reshape around the source parameter to ensure that it is a row vector
% if source.pos is not available, determine the number of dipoles from source.dim
%
% Revision 1.6  2005/05/09 14:14:46  roboos
% made the assignment of the output stat fields that describe the grid (like pos, xgrid, dim ...) optional, so that it won't fail when one of these fields is not present
%
% Revision 1.5  2005/03/18 14:37:51  roboos
% changed order of preprocessing: first z-transform, then relchange
% fixed small bug due to latest update
%
% Revision 1.4  2005/03/18 11:16:21  roboos
% fixed bug in two sided testing by replacing the relevant code by the RANDSTATPROB private subfunction
%
% Revision 1.3  2004/10/19 12:04:10  roboos
% fixed  bug in relchange if average power was not row-vector
%
% Revision 1.2  2004/10/13 15:27:41  roboos
% added version information to output cfg
% added configuration of input object to output cfg
% renamed single input argument from varargin to source
%
% Revision 1.1  2004/08/26 10:45:32  marsie
% initial CVS release
%
 
% set the defaults
if ~isfield(cfg, 'method'), cfg.method = 'randomization'; end

% these are the old defaults, they should be moved to the individual core functions
if ~isfield(cfg, 'tail'), cfg.tail = 0; end
if ~isfield(cfg, 'correctm'), cfg.correctm = 'no'; end
if ~isfield(cfg, 'ztransform'), cfg.ztransform = 'no'; end
if ~isfield(cfg, 'comparestat'), cfg.comparestat = 'difference'; end

% check the input variables
if ~strcmp(cfg.method, 'randomization')
  error('unsupported method for this function');
end
if nargin ~= 2
  error('invalid number of arguments, randomization based statistics require a single source structure');
end
if ~isfield(source,'avgA')
  error('no randomzation data included in source structure');
end  
if ~isfield(source.avgA,cfg.parameter)
  error('parameter to test is not included in source structure');
end

if isfield(source, 'pos')
  % count the number of dipole positions, which can be either on a regular or an irregular grid
  nDipole = size(source.pos,1);
elseif isfield(source, 'dim')
  % individual dipole positions are not available after spatial normalisation to a template anatomical MRI
  nDipole = prod(source.dim);
else
  % count the number of functional values, which should be the same as the number of dipoles
  dum = getfield(source.avgA,cfg.parameter);
  nDipole = prod(size(dum));
end

% get the data to work on in nice arrays
nTrial  = size(source.trialA,2);
strialA = zeros(nTrial,nDipole);
strialB = zeros(nTrial,nDipole);
for j = 1:nTrial
 % reshape to ensure that the functional parameter for all voxels is arranged in a row-vector
 strialA(j,:) = reshape(getfield(source.trialA,{j},cfg.parameter), [1 nDipole]);
 strialB(j,:) = reshape(getfield(source.trialB,{j},cfg.parameter), [1 nDipole]);
end
% reshape to ensure that the functional parameter for all voxels is arranged in a row-vector
savgA = reshape(getfield(source.avgA,cfg.parameter), [1 nDipole]);
savgB = reshape(getfield(source.avgB,cfg.parameter), [1 nDipole]);

% restrict computations to sources that are inside the brain
inside  = source.inside;
outside = source.outside;

% transform to z-scores if requested
if strcmp(cfg.ztransform,'yes')
  meanA = repmat(mean(strialA(:,inside)),size(strialA,1),1);
  meanB = repmat(mean(strialB(:,inside)),size(strialB,1),1);
  stdA = repmat(std(strialA(:,inside)),size(strialA,1),1);
  stdB = repmat(std(strialB(:,inside)),size(strialB,1),1);
  strialA(:,inside) = (strialA(:,inside) - meanA) ./ stdA; 
  strialB(:,inside) = (strialB(:,inside) - meanB) ./ stdB;
  savgA(1,inside) = (savgA(1,inside) - meanA(1,:)) ./ stdA(1,:);
  savgB(1,inside) = (savgB(1,inside) - meanB(1,:)) ./ stdB(1,:);
end

% normalize by condition B if requested
if strcmp(cfg.comparestat,'difference'),
   randobs     = strialA-strialB;
   realobs     = savgA-savgB;
elseif strcmp(cfg.comparestat,'relchange'),
   randobs     = (strialA-strialB)./strialB;
   realobs     = (savgA-savgB)./savgB;
end;

% compute p-values by comparing the real observations with those obtained after randomization
stat.prob(inside,:)  = randstatprob(randobs(:,inside)', realobs(:,inside)', cfg.tail, strcmp(cfg.correctm, 'yes'));
stat.prob(outside,:) = nan;

% remember the details of the source reconstruction
try, stat.dim       = source.dim;       end
try, stat.xgrid     = source.xgrid;     end
try, stat.ygrid     = source.ygrid;     end
try, stat.zgrid     = source.zgrid;     end
try, stat.inside    = source.inside;    end
try, stat.outside   = source.outside;   end
try, stat.pos       = source.pos;       end
try, stat.transform = source.transform; end

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: sourcestatistics_randomization.m,v 1.10 2006/01/30 12:51:38 roboos Exp $';
% remember the configuration details of the input data
try, cfg.previous = source.cfg; end
% remember the exact configuration details in the output 
stat.cfg = cfg;

