function [source] = sourceanalysis(cfg, data, baseline);

% SOURCEANALYSIS performs beamformer dipole analysis on EEG or MEG data
% after preprocessing and a timelocked or frequency analysis
%
% Use as either
%   [source] = sourceanalysis(cfg, freq)
%   [source] = sourceanalysis(cfg, timelock)
%
% where the data in freq or timelock should be organised in a structure
% as obtained from the FREQANALYSIS or TIMELOCKANALYSIS function. The
% configuration "cfg" is a structure containing information about
% source positions and other options.
% 
% The different source reconstruction algorithms that are implemented 
% for timelocked data are
%   cfg.method     = 'lcmv'    linear constrained minimum variance beamformer
%   cfg.method     = 'mne'     minimum norm estimation
%   cfg.method     = 'loreta'  minimum norm estimation with smoothness constraint
%   cfg.method     = 'rv'      scan residual variance with single dipole
%   cfg.method     = 'music'   multiple signal classification
%   cfg.method     = 'pcc'     partial cannonical correlation/coherence
%
% The different source reconstruction algorithms that are implemented 
% for frequency domain data are
%   cfg.method     = 'dics'    dynamic imaging of coherent sources
%   cfg.method     = 'pcc'     partial cannonical correlation/coherence
% 
% The positions of the sources are either specified as regular orthogonal grid  
%   cfg.xgrid      = something like -20:1:20 or 'auto'
%   cfg.ygrid      = something like -20:1:20 or 'auto'
%   cfg.zgrid      = something like   0:1:20 or 'auto'
%   cfg.resolution = number (e.g. 1 cm) required for automatic grid generation
% or as a selected number of dipoles at locations of interest that can be
% obtained from an anatomical or functional MRI
%   cfg.grid.pos    = Ndipoles x 3 matrix with [x y z] values
% or as a grid with precomputed leadfields from a previous run of
% sourecanalysis or from the PRECOMPUTE_LEADFIELDS function.
%   cfg.grid       = precomputed source structure
%
% The following strategies are supported to obtain statistics for the source parameters using 
% multiple trials in the data, either directly or through a resampling-based approach
%   cfg.singletrial   = 'no' or 'yes'   construct filter from average, apply to single trials
%   cfg.rawtrial      = 'no' or 'yes'   construct filter from single trials, apply to single trials
%   cfg.jacknife      = 'no' or 'yes'   jacknife resampling of trials
%   cfg.pseudovalue   = 'no' or 'yes'   pseudovalue resampling of trials
%   cfg.bootstrap     = 'no' or 'yes'   bootstrap resampling of trials
%   cfg.numbootstrap  = number of bootstrap replications (e.g. number of original trials)
% If none of these options is specified, the average over the trials will
% be computed prior to computing the source reconstruction.
%
% To obtain statistics over the source parameters between two conditions, you
% can also use a resampling procedure that reshuffles the trials over both
% conditions. In that case, you should call the function with two datasets
% containing single trial data like
%   [source] = sourceanalysis(cfg, freqA, freqB)
%   [source] = sourceanalysis(cfg, timelockA, timelockB)
% and you should specify
%   cfg.randomization      = 'no' or 'yes'
%   cfg.permutation        = 'no' or 'yes'
%   cfg.numrandomization   = number, e.g. 500
%   cfg.numpermutation     = number, e.g. 500 or 'all'
%
% You should specify the volume conductor model and, if the sensor information
% is not contained in the data itself, you should also specify the sensor
% information using
%   cfg.hdmfile       = file containing the volume conduction model (required)
%   cfg.gradfile      = string (optional, can already be present in data)
%   cfg.elecfile      = string (optional, can already be present in data)
% or alternatively
%   cfg.vol           = structure with volume conduction model
%   cfg.grad          = structure with gradiometer definition
%   cfg.elec          = structure with electrode definition
% 
% Other configuration options are 
%   cfg.channel       = cell-array (default is 'all', see CHANNELSELECTION)
%   cfg.frequency     = single number in Hz
%   cfg.latency       = single number in seconds, for time-frequency analysis
%   cfg.lambda        = number or empty for automatic default
%   cfg.refchan       = reference channel label (for coherence)
%   cfg.refdip        = reference dipole location (for coherence)
%   cfg.supchan       = suppressed channel label(s)
%   cfg.supdip        = suppressed dipole location(s)
%   cfg.keeptrials    = 'no' or 'yes'
%   cfg.keepleadfield = 'no' or 'yes'
%   cfg.projectnoise  = 'no' or 'yes'
%   cfg.keepfilter    = 'no' or 'yes'
%   cfg.keepcsd       = 'no' or 'yes'
%   cfg.keepmom       = 'no' or 'yes'
%   cfg.feedback      = 'gui', 'text' or 'no'
%
% See also SOURCESTATISTICS, PRECOMPUTE_LEADFIELD

% Undocumented general configuration options are
%   cfg.trialweight        = 'equal' or 'proportional'
%   cfg.normalize          = for leadfields in beamformer methods
%   cfg.reducerank         = for leadfields in beamformer methods
%   cfg.powmethod          = 'lambda1' or 'trace'

% Undocumented configuration options for parallellization on linux cluster are
%   cfg.parallel           = 'no', 'beowulf' or 'mentat' (using different backends)
%   cfg.wakewulf           = 'yes' or 'no' (automatically start cluster, default is yes)
%   cfg.killwulf           = 'yes' or 'no' (automatically stop cluster, default is yes)

% Undocumented options for cortex segmentation
%   cfg.mri           = can be filename, MRI structure or segmented volume structure
%   cfg.sourceunits   = 'mm' or 'cm' (default is 'cm')
%   cfg.mriunits      = 'mm' or 'cm' (default is 'mm')
%   cfg.threshold     = 0.1, relative to the maximum value in the segmentation
%   cfg.smooth        = 5, smoothing in voxels

% Copyright (c) 2003-2004, Robert Oostenveld, F.C. Donders Centre
%
% $Log: sourceanalysis.m,v $
% Revision 1.94  2006/02/24 16:41:39  roboos
% changed foi and toi into freq and time for frequency data
%
% Revision 1.93  2006/02/23 10:28:16  roboos
% changed dimord strings for consistency, changed toi and foi into time and freq, added fixdimord where neccessary
%
% Revision 1.92  2006/02/07 22:21:01  roboos
% changed the xgrid/ygrid/zgrid and dim in the output source structure (from cfg to grid), they can be different from the cfg since prepare_dipole_grid will make the box tight
%
% Revision 1.91  2006/02/07 20:08:23  roboos
% changed all occurences of a dimord with chancmb (was previous sgncmb) into chan
%
% Revision 1.90  2006/02/01 12:26:00  roboos
% made all uses of dimord consistent with the common definition of data dimensions, see the fixdimord() function
%
% Revision 1.89  2006/01/24 20:51:57  roboos
% fixed obvious typo
%
% Revision 1.88  2006/01/24 20:02:50  roboos
% renamed dics_cohrefdip and cohrefchan into cfg.method=dics (backwards compatible)
% added extra local variable "submethod" that indicates the details for dics
%
% Revision 1.87  2006/01/24 14:28:33  roboos
% small change in documentation
%
% Revision 1.86  2005/11/24 16:25:16  roboos
% added average or single trial ERF to the input of PCC beamformer in timedomain
%
% Revision 1.85  2005/11/16 09:07:03  roboos
% added option cfg.powmethod to determine whether lambda1 or trace is used in beamformer
%
% Revision 1.84  2005/11/08 11:07:13  roboos
% added support for normalize and reducerank, by passing these options on to the beamformer subfunction
%
% Revision 1.83  2005/09/29 00:45:00  roboos
% added multiple signal classification as reconstruction technique for timelocked data
%
% Revision 1.82  2005/10/25 08:42:56  roboos
% changed the detection of either frequency or timelock data using flags (isfreq, iscomp, istimelock) instead of by renaming the input data object to freq/timelock
% renamed the freq/timelock variable throughout the code into "data"
% added support for freq/comp data to the non-beamforming timelocked source reconstruction methods
%
% Revision 1.81  2005/10/14 15:47:03  roboos
% add an identity covariance matrix to timelocked data if not present
% some cosmetic changes
%
% Revision 1.80  2005/10/05 11:14:11  roboos
% updated documentation
% added support for LORETA (still experimental)
% removed all code that referred to unimplemented parallelization for some of the reconstruction methods
%
% Revision 1.79  2005/09/05 06:36:56  jansch
% keep cumtapcnt of freq-input in output
%
% Revision 1.78  2005/08/16 13:15:55  jansch
% *** empty log message ***
%
% Revision 1.77  2005/08/16 12:43:03  jansch
% included possibility to have 'rpttap_sgncmb_frq' as input
%
% Revision 1.76  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.75  2005/06/17 09:01:13  roboos
% moved implementation of keepmom and keepcsd over to beamformer subfunction
%
% Revision 1.74  2005/05/17 17:56:38  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
% implemented support for projecting the single trial fft matrix through the filters of beamformer_pcc
%
% Revision 1.73  2005/04/25 16:05:11  roboos
% includedupdated options for cortex segmentation based grid in the documentation (not in main documentation)
%
% Revision 1.72  2005/03/03 13:57:29  roboos
% implemented pseudovalue resampling, i.e. a combination of jacknife and the complete average
%
% Revision 1.71  2005/03/03 10:57:43  roboos
% Renamed the source analysis methods, the submethods for dics are now
% recognized by their additional required fields (refchan/refdip). Added
% support for experimental pcc method (not yet implemented in a clean way),
% pcc is now a re3cognised submethod of the standard beamformer function.
% Included mne method in the documentation.
%
% Revision 1.70  2005/02/14 21:53:14  roboos
% changed assignment of empty struct to dip into 'clear dip' (for mentat-style parallelization)
%
% Revision 1.69  2005/01/26 08:38:50  roboos
% extended and fixed the experimental partcanoncorr code
%
% Revision 1.68  2005/01/11 12:43:03  roboos
% added temporary hack to support partcanoncorr as a beamformer method
%
% Revision 1.67  2005/01/10 17:26:46  roboos
% added support for 2nd method of parallelizing the code, separated
% two methods by the name of the backend toolbox (beowulf or mentat)
%
% Revision 1.66  2004/12/08 18:00:13  roboos
% implemented consistent method of selecting a subset of channels for
% forward and inverse computations using cfg.channel and updated the
% ducumentation
%
% Revision 1.65  2004/10/21 09:01:34  roboos
% added numcondition=2 for averaging two conditions while resampling
%
% Revision 1.64  2004/09/28 08:01:03  roboos
% added residual variance scanning with single dipole (cfg.method=rv)
% added 'yes' as allowed value for cfg.parallel
%
% Revision 1.63  2004/09/21 14:30:07  roboos
% fixed bug in DICS/power for two conditions (Nbaseline~=Ntrials)
%
% Revision 1.62  2004/09/21 13:22:15  roboos
% in the previous revision (1.61) initial support for MNE was added but accidentally not logged in CVS
% this revision only adds a single comment line with the undocumented method=mne option
%
% Revision 1.61  2004/09/21 13:15:32  roboos
% *** empty log message ***
%
% Revision 1.60  2004/09/21 12:49:00  roboos
% fixed a bug in averaging of 2 input conditions using prepare_resampled_data (using cfg2)
%
% Revision 1.59  2004/09/14 16:28:49  roboos
% incorporated the changes on the "exp-resampling-parallel" branch
% new implementation of resampling/randomization
% new implementation of parallellization
%
% Revision 1.58.4.5  2004/09/08 16:36:00  roboos
% fixed many bugs that stopped execution through rigorous testing
% correctness of the output still has to be evaluated
% non-proportional trial weighing is not implemented yet
%
% Revision 1.58.4.4  2004/09/08 12:25:13  roboos
% auto-indented code using matlab editor, fixed one incorrect "end"
%
% Revision 1.58.4.3  2004/09/08 09:35:49  roboos
% implemented prepare_resampled_data and new evalwulf for DICS, some related changes to LCMV
%
% Revision 1.58.4.2  2004/09/06 16:32:28  roboos
% implemented new prepare_resampled_data for a few resampling strategies and lcmv only, untested
%
% Revision 1.58.4.1  2004/09/06 16:10:26  roboos
% disabled parallelization over grid (still in code)
% moved documentation for paralellization to hidden location
% fixed bug in jacknife+lcmv
%
% Revision 1.58  2004/08/16 12:14:33  roboos
% fixed small bug that was introduced in the previous update
%
% Revision 1.57  2004/08/16 10:53:09  roboos
% fixed bug for lcmv together with randomization (trial assignment was oubviously wrong)
%
% Revision 1.56  2004/08/16 10:01:23  roboos
% updated help, fixed bug in output assignment with permutation, changed handling of keepleadfield and precomputed leadfields
%
% Revision 1.55  2004/08/06 12:23:04  roboos
% fixed small bug in default setting for cfg.permutation
%
% Revision 1.54  2004/08/05 07:14:25  roboos
% implemented permutation analogous to randomization, updated and restructured help
%
% Revision 1.53  2004/08/05 06:37:29  roboos
% extra removal of leadfield in output source and cfg, removed old revision comments
%
% Revision 1.52  2004/08/04 11:53:03  roboos
% updated the help
%
% Revision 1.51  2004/06/09 10:07:19  roberto
% added version information and input data configuration to output cfg
%
% Revision 1.50  2004/03/01 11:01:55  roberto
% added precomputed filters to grid and using these, implemented singletrial for dics
%
% Revision 1.49  2004/02/05 12:19:36  roberto
% moved the code that was common with dipolefitting into separate prepare_grid and prepare_sens_vol functions
% also renamed internal subfunction freq_to_matrices to prepare_freq_matrices and moved into separate function
%
% Revision 1.48  2004/02/02 10:48:11  roberto
% fixed two bugs in randomization for lcmv & timelocked data
%
% Revision 1.47  2004/01/26 09:04:54  roberto
% fixed serious bug in selection of spheres for multi-sphere volume model
% selection was based on wrong indices, causing reference magnetometer spheres
% and gradiometer spheres to be mixed up
%
% Revision 1.46  2004/01/22 10:50:55  roberto
% fixed bug in cfg.keepmom
% improved handling of pre-computation of leadfield without scanning
%
% Revision 1.45  2004/01/21 13:07:52  roberto
% added a fprintf-line to indicate the progress in randomization (which takes quite some time)
%
% Revision 1.44  2004/01/20 09:09:45  roberto
% added computation of leadfield in case not provided and multiple trials have to
% be scanned
%
% Revision 1.43  2004/01/14 14:29:43  roberto
% changed handling of sgncmbindx etc inside freq_to_matrices subfunction
% mainly cosmetical changes, but also fixed bug in refindx for pow_refchan
%
% Revision 1.42  2004/01/14 08:50:05  roberto
% changed the name of output variables for randomization
%
% Revision 1.41  2004/01/08 09:25:35  roberto
% canged cfg.label into freq.sgn inside freq_to_matrices subfunction to fix bug
%
% Revision 1.40  2004/01/08 09:14:49  roberto
% added option keepmom to save memory, default is 'yes'

% set a timer to determine how long the sourceanalysis takes in total
tic;

% for backward compatibility with old data structures
data = fixdimord(data);

% determine the type of input data
if isfield(data, 'freq')
  isfreq     = 1;
  iscomp     = 0;
  istimelock = 0;
elseif isfield(data, 'topo')
  isfreq     = 0;
  iscomp     = 1;
  istimelock = 0;
elseif isfield(data, 'time')
  iscomp     = 0;
  isfreq     = 0;
  istimelock = 1;
else
  error('input data is not recognized');
end

% set the defaults
if ~isfield(cfg, 'method') && istimelock, cfg.method = 'lcmv';      end
if ~isfield(cfg, 'method') && isfreq,     cfg.method = 'dics';      end
if ~isfield(cfg, 'keeptrials')        cfg.keeptrials  = 'no';     end
if ~isfield(cfg, 'keepfilter')        cfg.keepfilter = 'no';      end
if ~isfield(cfg, 'keepleadfield')     cfg.keepleadfield = 'no';   end
if ~isfield(cfg, 'keepcsd')           cfg.keepcsd     = 'no';     end
if ~isfield(cfg, 'keepmom')           cfg.keepmom     = 'yes';    end
if ~isfield(cfg, 'projectnoise')      cfg.projectnoise = 'yes';   end
if ~isfield(cfg, 'trialweight')       cfg.trialweight = 'equal';  end
if ~isfield(cfg, 'jacknife'),         cfg.jacknife    = 'no';     end
if ~isfield(cfg, 'pseudovalue'),      cfg.pseudovalue = 'no';     end
if ~isfield(cfg, 'bootstrap'),        cfg.bootstrap   = 'no';     end
if ~isfield(cfg, 'singletrial'),      cfg.singletrial = 'no';     end
if ~isfield(cfg, 'rawtrial'),         cfg.rawtrial    = 'no';     end
if ~isfield(cfg, 'randomization'),    cfg.randomization = 'no';   end
if ~isfield(cfg, 'numrandomization'), cfg.numrandomization = 100; end
if ~isfield(cfg, 'permutation'),      cfg.permutation = 'no';     end
if ~isfield(cfg, 'numpermutation'),   cfg.numpermutation = 100;   end
if ~isfield(cfg, 'parallel'),         cfg.parallel    = 'no';     end
if ~isfield(cfg, 'wakewulf'),         cfg.wakewulf    = 'yes';    end
if ~isfield(cfg, 'killwulf'),         cfg.killwulf    = 'yes';    end
if ~isfield(cfg, 'feedback'),         cfg.feedback    = 'text';   end
if ~isfield(cfg, 'xgrid'),            cfg.xgrid = 'auto';         end
if ~isfield(cfg, 'ygrid'),            cfg.ygrid = 'auto';         end
if ~isfield(cfg, 'zgrid'),            cfg.zgrid = 'auto';         end
if ~isfield(cfg, 'supdip'),           cfg.supdip = [];            end
if ~isfield(cfg, 'lambda'),           cfg.lambda = [];            end
if ~isfield(cfg, 'powmethod'),        cfg.powmethod = [];         end
if ~isfield(cfg, 'channel'),          cfg.channel = 'all';        end
if ~isfield(cfg, 'normalize'),        cfg.normalize = 'no';       end
if ~isfield(cfg, 'reducerank'),       cfg.reducerank = 'no';      end

convertfreq = 0;
convertcomp = 0;
if ~istimelock && (strcmp(cfg.method, 'mne') || strcmp(cfg.method, 'loreta') || strcmp(cfg.method, 'rv') || strcmp(cfg.method, 'music'))
  % these timelock methods are also supported for frequency or component data
  if isfreq
    [data, cfg] = freq2timelock(cfg, data);
    convertfreq = 1;  % flag indicating that the data was converted
  elseif iscomp
    [data, cfg] = comp2timelock(cfg, data);
    convertcomp = 1;  % flag indicating that the data was converted
  end
  istimelock = 1;     % from now on the data can be treated as timelocked
  isfreq     = 0;
  iscomp     = 0;
end

% for backwards compatibility
if isfield(cfg, 'refchannel')
  cfg.refchan = cfg.refchannel;
  cfg = rmfield(cfg, 'refchannel');
end

% these should be cell-arrays
if isfield(cfg, 'refchan') && ischar(cfg.refchan)
  cfg.refchan = {cfg.refchan};
end
if isfield(cfg, 'supchan') && ischar(cfg.supchan)
  cfg.supchan = {cfg.supchan};
end

% for backwards compatibility, these are all submethods for dics
% the appropriate submethod will be selected, depending on the other cfg settings
if strcmp(cfg.method, 'power')
  cfg.method = 'dics';
elseif strcmp(cfg.method, 'coh_refchan')
  cfg.method = 'dics';
elseif strcmp(cfg.method, 'coh_refdip')
  cfg.method = 'dics';
elseif strcmp(cfg.method, 'dics_cohrefchan')
  cfg.method = 'dics';
elseif strcmp(cfg.method, 'dics_cohrefdip')
  cfg.method = 'dics';
end

if nargin>2 && (strcmp(cfg.randomization, 'no') && strcmp(cfg.permutation, 'no'))
  error('input of two conditions only makes sense if you want to randomize or permute');
elseif nargin<3 && (strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
  error('randomization or permutation requires that you give two conditions as input');
end

if isfield(cfg, 'latency') && istimelock
  error('specification of cfg.latency is only required for time-frequency data');
end

if sum([strcmp(cfg.jacknife, 'yes'), strcmp(cfg.bootstrap, 'yes'), strcmp(cfg.pseudovalue, 'yes'), strcmp(cfg.singletrial, 'yes'), strcmp(cfg.rawtrial, 'yes'), strcmp(cfg.randomization, 'yes'), strcmp(cfg.permutation, 'yes')])>1
  error('jacknife, bootstrap, pseudovalue, singletrial, rawtrial, randomization and permutation are mutually exclusive');
end

if isfreq
  if ~strcmp(data.dimord, 'chan_freq') && ...
     ~strcmp(data.dimord, 'chan_freq_time') && ...
     ~strcmp(data.dimord, 'rpt_chan_freq') && ...
     ~strcmp(data.dimord, 'rpt_chan_freq_time') && ...
     ~strcmp(data.dimord, 'rpttap_chan_freq')
    error('dimord of input frequency data is not recognized');
  end
end
if isempty(cfg.method)
  % parallellization is not required since only leadfield computation is requested
  cfg.parallel = 'no';
end

% check cfg.parallel for backward compatibility with old scripts
if strcmp(cfg.parallel, 'grid')
  error('parallellization over grid is not supported any more (cfg.parallel)')
elseif strcmp(cfg.parallel, 'jacknife') || strcmp(cfg.parallel, 'bootstrap') || strcmp(cfg.parallel, 'yes')
  % these options were supported in an older version, offer backward compatibility
  % since the parallellization is the same as for single trial analysis
  warning('the options for parallelization have changed, please contact Robert Oostenveld');
  warning('assuming default type of parallelization (''beowulf'') over trials');
  cfg.parallel = 'beowulf';
end

% check for the presence of these two optional paralellization toolboxes
hasbeowulf = (exist('evalwulf')==2 & exist('evalwulf')==2 & exist('evalwulf')==2);
hasmentat  = (exist('pcompile')==2 & exist('pfor')==2 & exist('peval')==2);

if strcmp(cfg.parallel, 'beowulf') && ~hasbeowulf
  error('beowulf parallelization toolbox was not found');
end

if strcmp(cfg.parallel, 'mentat') && ~hasmentat
  error('mentat parallelization toolbox was not found');
end

% start the beowulf cluster if parallelization is required
if strcmp(cfg.parallel, 'beowulf') && strcmp(cfg.wakewulf, 'yes')
  fprintf('atempting to start beowulf cluster from within sourceanalysis\n');
  wakewulf;
end

% collect and preprocess the electrodes/gradiometer and head model
[vol, sens, cfg] = prepare_vol_sens(cfg, data);

% construct the grid on which the scanning will be done
[grid, cfg] = prepare_dipole_grid(cfg, vol, sens);

if ~isfield(grid, 'leadfield') &&    ( ...
    strcmp(cfg.permutation, 'yes')   || ...
    strcmp(cfg.randomization, 'yes') || ...
    strcmp(cfg.bootstrap, 'yes')     || ...
    strcmp(cfg.jacknife, 'yes')      || ...
    strcmp(cfg.pseudovalue, 'yes')   || ...
    strcmp(cfg.singletrial, 'yes')   || ...
    strcmp(cfg.rawtrial, 'yes')      )
  % also precompute the leadfields if multiple trials have to be processed
  fprintf('precomputing leadfields for efficient handling of multiple trials\n');
  grid = prepare_leadfield(cfg, data);
  fprintf('precomputing leadfield ready\n');
elseif strcmp(cfg.keepleadfield, 'yes')
  % precompute the leadfields upon the users request
  grid = prepare_leadfield(cfg, data);
end

% select only those channels that have a sensor
cfg.channel = channelselection(cfg.channel, sens.label);
Nchans      = length(cfg.channel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do frequency domain source reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfreq && (strcmp(cfg.method, 'dics') || strcmp(cfg.method, 'pcc'))

  if strcmp(cfg.method, 'pcc')
    % HACK: requires some extra defaults
    if ~isfield(cfg, 'refdip'), cfg.refdip = []; end
    if ~isfield(cfg, 'supdip'), cfg.supdip = []; end

    % HACK: use some experimental code
    if nargin>2
      error('not supported')
    end
    tmpcfg         = cfg;
    tmpcfg.refchan = [];	% prepare_freq_matrices should not know explitely about the refchan
    tmpcfg.channel = cfg.channel(:)';
    if isfield(cfg, 'refchan')
      % add the refchan implicitely
      tmpcfg.channel = [tmpcfg.channel cfg.refchan(:)'];
    end
    if isfield(cfg, 'supchan')
      % add the supchan implicitely
      tmpcfg.channel = [tmpcfg.channel cfg.supchan(:)'];
    end

    % select the data in the channels and the frequency of interest
    [Cf, Cr, Pr, Ntrials, tmpcfg] = prepare_freq_matrices(tmpcfg, data);

    if isfield(cfg, 'refchan') && ~isempty(cfg.refchan)
      refchanindx = match_str(tmpcfg.channel,cfg.refchan);
    else
      refchanindx = [];
    end
    if isfield(cfg, 'supchan') && ~isempty(cfg.supchan)
      supchanindx = match_str(tmpcfg.channel,cfg.supchan);
    else
      supchanindx = [];
    end
    Nchans = length(tmpcfg.channel); % update the number of channels 

    % if the input data has a complete fourier spectrum, project it through the filters
    if isfield(data, 'fourierspctrm')
      chansel = match_str(data.label, tmpcfg.channel);
      fbin = nearest(data.freq, cfg.frequency);
      if strcmp(data.dimord, 'chan_freq')
        avg = data.fourierspctrm(chansel, fbin);
      elseif strcmp(data.dimord, 'rpt_chan_freq') || strcmp(data.dimord, 'rpttap_chan_freq'),
        avg = data.fourierspctrm(:, chansel, fbin)';
      elseif strcmp(data.dimord, 'chan_freq_time')
        tbin = nearest(data.time, cfg.latency);
        avg = data.fourierspctrm(chansel, fbin, tbin);
      elseif strcmp(data.dimord, 'rpt_chan_freq_time') || strcmp(data.dimord, 'rpttap_chan_freq_time'),
        tbin = nearest(data.time, cfg.latency);
        avg = data.fourierspctrm(:, chansel, fbin, tbin)';
      end
    else
      avg = [];
    end

  else
    % HACK: use the default code
    % convert the input data, so that Cf, Cr and Pr contain either the average over all trials (if Ntrials==1)
    % or the individual cross-spectral-densities/powers of each individual trial (if Ntrials>1)
    [Cf, Cr, Pr, Ntrials, cfg] = prepare_freq_matrices(cfg, data);
  end

  % assign a descriptive name to each of the dics sub-methods, the default is power only
  if strcmp(cfg.method, 'dics') && isfield(cfg, 'refdip') && ~isempty(cfg.refdip);
    submethod = 'dics_refdip';
  elseif strcmp(cfg.method, 'dics') && isfield(cfg, 'refchan') && ~isempty(cfg.refchan);
    submethod = 'dics_refchan';
  else
    submethod = 'dics_power';
  end  

  % fill these with NaNs, so that I dont have to treat them separately
  if isempty(Cr), Cr = nan*zeros(Ntrials, Nchans, 1); end
  if isempty(Pr), Pr = nan*zeros(Ntrials, 1, 1); end

  if nargin>2
    % repeat the conversion for the baseline condition
    [bCf, bCr, bPr, Nbaseline, cfg] = prepare_freq_matrices(cfg, baseline);
    % fill these with NaNs, so that I dont have to treat them separately
    if isempty(bCr), bCr = nan*zeros(Nbaseline, Nchans, 1); end
    if isempty(bPr), bPr = nan*zeros(Nbaseline, 1, 1); end
    % rename the active condition for convenience
    aCf = Cf;
    aCr = Cr;
    aPr = Pr;
    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
    % this is required for randomizing/permuting 2 conditions using prepare_resampled_data
    cfg.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is neccessary
  if (Ntrials<=1) && (strcmp(cfg.jacknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
    error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random permutation
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.randomization, 'yes') 
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random resampling
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.jacknife, 'yes')
    % compute the cross-spectral density matrix with jacknife resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the cross-spectral density matrix with bootstrap resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes')
    % compute the cross-spectral density matrix with pseudovalue resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    error('this option contains a bug, and is therefore not supported at the moment');
    % The filter to use on singletrial Cfs has either been provided or
    % it has been computed above based on the average Cf so
    % keep all the individual trials, do not average them
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % keep all the individual trials, do not average them
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    % compute the average from the individual trials
    Cf = reshape(sum(Cf, 1) / Ntrials, [Nchans Nchans]);
    Cr = reshape(sum(Cr, 1) / Ntrials, [Nchans 1]);
    Pr = reshape(sum(Pr, 1) / Ntrials, [1 1]);
    Nrepetitions = 1;

  elseif Ntrials==1
    % no rearrangement of trials is neccesary, the data already represents the average
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cf  = reshape(Cf , [1 Nchans Nchans]);
    Cr  = reshape(Cr , [1 Nchans 1]);
    Pr  = reshape(Pr , [1 1 1]);
  end

  % start the beamformer scanning for all repetitions
  if strcmp(cfg.parallel, 'beowulf')
    fprintf('scanning %d repetitions in parallel\n', Nrepetitions);
    if strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_power')
      dip = evalwulf('i=1:Nrepetitions', 'dip(i)', 'beamformer', 'grid, sens, vol, [], squeeze(Cf(i,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd');
    elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refdip')
      dip = evalwulf('i=1:Nrepetitions', 'dip(i)', 'beamformer', 'grid, sens, vol, [], squeeze(Cf(i,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd, ''refdip'', cfg.refdip');
    elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refchan')                                                                                                                                                                                                                                                                                             
      dip = evalwulf('i=1:Nrepetitions', 'dip(i)', 'beamformer', 'grid, sens, vol, [], squeeze(Cf(i,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd, ''Cr'', Cr(i,:), ''Pr'', Pr(i)');
    else
      error(sprintf('method ''%s'' is unsupported for source reconstruction in the frequency domain using parallelization', cfg.method));
    end

  elseif strcmp(cfg.parallel, 'mentat')
    fprintf('scanning %d repetitions in parallel using mentat binaries\n', Nrepetitions);
    clear dip
    if strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_power')
      pfor(1:Nrepetitions, 'dip(%d) = beamformer(grid, sens, vol, [], squeeze(Cf(%d,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd)');
    elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refdip')
      pfor(1:Nrepetitions, 'dip(%d) = beamformer(grid, sens, vol, [], squeeze(Cf(%d,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd, ''refdip'', cfg.refdip)');
    elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refchan')
      pfor(1:Nrepetitions, 'dip(%d) = beamformer(grid, sens, vol, [], squeeze(Cf(%d,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepfilter'', cfg.keepfilter, ''supdip'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod, ''keepmom'', cfg.keepmom, ''keepcsd'', cfg.keepcsd, ''Cr'', Cr(%d,:), ''Pr'', Pr(%d))');
    else
      error(sprintf('method ''%s'' is unsupported for source reconstruction in the frequency domain using parallelization', cfg.method));
    end
  
  else
    for i=1:Nrepetitions
      fprintf('scanning repetition %d\n', i);
      if strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_power')
        dip(i) = beamformer(grid, sens, vol, [], squeeze(Cf(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepfilter', cfg.keepfilter, 'supdip', cfg.supdip, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod, 'keepmom', cfg.keepmom, 'keepcsd', cfg.keepcsd);
      elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refdip')
        dip(i) = beamformer(grid, sens, vol, [], squeeze(Cf(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepfilter', cfg.keepfilter, 'supdip', cfg.supdip, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod, 'keepmom', cfg.keepmom, 'keepcsd', cfg.keepcsd, 'refdip', cfg.refdip);
      elseif strcmp(cfg.method, 'dics') && strcmp(submethod, 'dics_refchan')
        dip(i) = beamformer(grid, sens, vol, [], squeeze(Cf(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepfilter', cfg.keepfilter, 'supdip', cfg.supdip, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod, 'keepmom', cfg.keepmom, 'keepcsd', cfg.keepcsd, 'Cr', Cr(i,:), 'Pr', Pr(i));
      elseif strcmp(cfg.method, 'pcc')
        % experimental implementation of partial canonical correlation/coherence beamformer method
        dip(i) = beamformer(grid, sens, vol, avg, squeeze(Cf(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepfilter', cfg.keepfilter, 'supdip', cfg.supdip, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod, 'keepmom', cfg.keepmom, 'keepcsd', cfg.keepcsd, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supchan', supchanindx);
      else
        error(sprintf('method ''%s'' is unsupported for source reconstruction in the frequency domain', cfg.method));
      end
    end
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do time domain source reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif istimelock && (strcmp(cfg.method, 'lcmv') || strcmp(cfg.method, 'mne') || strcmp(cfg.method, 'loreta') || strcmp(cfg.method, 'rv') || strcmp(cfg.method, 'music') || strcmp(cfg.method, 'pcc')) 

  % determine the size of the data
  Nsamples = size(data.avg,2);
  if isfield(data, 'cov') && length(size(data.cov))==3
    Ntrials = size(data.cov,1);
  elseif isfield(data, 'trial') && length(size(data.trial))==3
    Ntrials = size(data.trial,1);
  else
    Ntrials = 1;
  end

  if ~isfield(data, 'cov')
    % add a identity covariance matrix, this simplifies the handling of the different source reconstruction methods
    % since the covariance is only used by some reconstruction methods and might not allways be present in the data
    if Ntrials==1
      data.cov = eye(length(data.label));
    else
      for i=1:Ntrials
        data.cov(i,:,:) = eye(length(data.label));
      end
    end
  end


  if strcmp(cfg.method, 'pcc')
    % HACK: requires some extra defaults
    if ~isfield(cfg, 'refdip'), cfg.refdip = []; end
    if ~isfield(cfg, 'supdip'), cfg.supdip = []; end

    % HACK: experimental code
    if nargin>2
      error('not supported')
    end
    tmpcfg = [];
    tmpcfg.channel = cfg.channel(:)';
    if isfield(cfg, 'refchan')
      tmpcfg.channel = [tmpcfg.channel cfg.refchan(:)'];
    end
    if isfield(cfg, 'supchan')
      tmpcfg.channel = [tmpcfg.channel cfg.supchan(:)'];
    end

    % select the data in the channels of interest
    seldata = match_str(data.label, tmpcfg.channel);
    if Ntrials==1
      data.avg = data.avg(seldata,:);
      data.cov = data.cov(seldata,seldata);
    else
      data.avg   = data.avg(seldata,:);
      data.cov   = data.cov(:,seldata,seldata);
      data.trial = data.trial(:,seldata,:);
    end
    data.label = data.label(seldata);

    if isfield(cfg, 'refchan') && ~isempty(cfg.refchan)
      refchanindx = match_str(data.label,cfg.refchan);
    else
      refchanindx = [];
    end
    if isfield(cfg, 'supchan') && ~isempty(cfg.supchan)
      supchanindx = match_str(data.label,cfg.supchan);
    else
      supchanindx = [];
    end
    Nchans = length(tmpcfg.channel); % update the number of channels 

  else
    % HACK: use the default code
    % select the channels of interest
    seldata = match_str(data.label, cfg.channel);
    if Ntrials==1
      data.avg = data.avg(seldata,:);
      data.cov = data.cov(seldata,seldata);
    else
      data.avg   = data.avg(seldata,:);
      data.cov   = data.cov(:,seldata,seldata);
      data.trial = data.trial(:,seldata,:);
    end
    data.label = data.label(seldata);

  end

  if nargin>2
    % baseline and active are only available together for resampling purposes,
    % hence I assume here that there are multiple trials in both
    baseline.avg   = baseline.avg(seldata,:);
    baseline.cov   = baseline.cov(:,seldata,seldata);
    baseline.trial = baseline.trial(:,seldata,:);
    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is neccessary
  if (strcmp(cfg.jacknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes')) && ~(Ntrials>1)
    error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random permutation
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.randomization, 'yes') 
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random resampling
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.jacknife, 'yes') 
    % compute the jacknife repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the bootstrap repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes') 
    % compute the pseudovalue repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    error('this option contains a bug, and is therefore not supported at the moment');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FIXME, the following implements proportional weighing, which is not consistently implemented sofar
    % Cy = zeros(Nchans,Nchans);
    % for trial=1:Ntrials
    %   Cy = Cy + squeeze(data.cov(trial,:,:)) * data.numcovsamples(trial);
    % end
    % Cy = Cy / sum(data.numcovsamples);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % average the individual trial covariance matrices
    Cy = squeeze(mean(data.cov,1));
    % make a copy of the average covariance matrix for every individual trial
    tmp = [];
    for i=1:Ntrials
      tmp(i,:,:) = Cy(:,:);
    end
    Cy = tmp;
    % keep the single trials, rename them to avg for convenience
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % do not do any resampling, keep the single trial covariances
    Cy(:,:,:) = data.cov;
    % do not do any resampling, keep the single trials (rename them to avg for convenience)
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FIXME, the following implements proportional weighing, which is not consistently implemented sofar
    % Cy = zeros(Nchans,Nchans);
    % for trial=1:Ntrials
    %   Cy = Cy + squeeze(data.cov(trial,:,:)) * data.numcovsamples(trial);
    % end
    % Cy = Cy / sum(data.numcovsamples);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % average the covariance for each trial and select the already averaged ERF/ERP
    Cy  = reshape(mean(data.cov,1), [Nchans Nchans]);
    avg = data.avg;
    % this reduces the number of trials to one
    Nrepetitions = 1;

  elseif Ntrials==1
    % no processing neccessary, just take the average covariance and ERF/ERP
    Cy  = data.cov;
    avg = data.avg;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cy  = reshape(Cy , [1 Nchans Nchans]);
    avg = reshape(avg, [1 Nchans Nsamples]);
  end

  % collect the extra arguments for the beamformer scan in one convenient cell-array
  % optargs = {'feedback' cfg.feedback 'keepfilter' cfg.keepfilter 'supdip' cfg.supdip 'lambda' cfg.lambda};

  if strcmp(cfg.method, 'lcmv')
    % start the beamformer scanning for all repetitions
    if strcmp(cfg.parallel, 'beowulf')
      fprintf('scanning %d repetitions in parallel\n', size(Cy,1));
      dip = evalwulf('i=1:Nrepetitions', 'dip(i)', 'beamformer', 'grid, sens, vol, squeeze(avg(i,:,:)), squeeze(Cy(i,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepmom'', cfg.keepmom, ''keepfilter'', cfg.keepfilter, ''sup'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod');
    elseif strcmp(cfg.parallel, 'mentat')
      fprintf('scanning %d repetitions in parallel using mentat binaries\n', size(Cy,1));
      clear dip
      pfor(1:Nrepetitions, 'dip(%d) = beamformer(               grid, sens, vol, squeeze(avg(%d,:,:)), squeeze(Cy(%d,:,:)), ''method'', cfg.method, ''feedback'', cfg.feedback, ''projectnoise'', cfg.projectnoise, ''keepmom'', cfg.keepmom, ''keepfilter'', cfg.keepfilter, ''sup'', cfg.supdip, ''lambda'', cfg.lambda, ''normalize'', cfg.normalize, ''reducerank'', cfg.reducerank, ''powmethod'', cfg.powmethod)');
    else
      for i=1:Nrepetitions
        fprintf('scanning repetition %d\n', i);
        dip(i) = beamformer(grid, sens, vol, squeeze(avg(i,:,:)), squeeze(Cy(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepmom', cfg.keepmom, 'keepfilter', cfg.keepfilter, 'sup', cfg.supdip, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod);
      end
    end
  elseif strcmp(cfg.method, 'pcc')
    for i=1:Nrepetitions
      fprintf('estimating residual variance for repetition %d\n', i);
      % experimental implementation of partial canonical correlation/coherence beamformer method
      dip(i) = beamformer(grid, sens, vol, squeeze(avg(i,:,:)), squeeze(Cy(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback, 'projectnoise', cfg.projectnoise, 'keepmom', cfg.keepmom, 'keepfilter', cfg.keepfilter, 'supdip', cfg.supdip, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supchan', supchanindx, 'lambda', cfg.lambda, 'normalize', cfg.normalize, 'reducerank', cfg.reducerank, 'powmethod', cfg.powmethod);
    end
  elseif strcmp(cfg.method, 'mne')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      dip(i) = minimumnormestimate(grid, sens, vol, squeeze(avg(i,:,:)), 'method', cfg.method, 'feedback', cfg.feedback);
    end
  elseif strcmp(cfg.method, 'loreta')
    for i=1:Nrepetitions
      fprintf('estimating LORETA reconstruction for repetition %d\n', i);
      dip(i) = loreta(grid, sens, vol, squeeze(avg(i,:,:)), 'feedback', cfg.feedback);
    end
  elseif strcmp(cfg.method, 'rv')
    for i=1:Nrepetitions
      fprintf('estimating residual variance for repetition %d\n', i);
      dip(i) = residualvariance(grid, sens, vol, squeeze(avg(i,:,:)), 'feedback', cfg.feedback);
    end
  elseif strcmp(cfg.method, 'music')
    for i=1:Nrepetitions
      fprintf('computing multiple signal classification for repetition %d\n', i);
      dip(i) = music(grid, sens, vol, squeeze(avg(i,:,:)), 'feedback', cfg.feedback, 'numcomponents', cfg.numcomponents);
    end
  else
    error(sprintf('method ''%s'' is unsupported for source reconstruction in the time domain', cfg.method));
  end

end % if freq or timelock data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up and collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% kill the beowulf cluster if parallelization was required
if strcmp(cfg.parallel, 'beowulf') && strcmp(cfg.killwulf, 'yes')
  killwulf;
end

if isfield(grid, 'xgrid')
  % the dipoles are placed on a regular grid that is aligned with the carthesian axes
  % copy the description of the grid axes
  source.xgrid = grid.xgrid;
  source.ygrid = grid.ygrid;
  source.zgrid = grid.zgrid;
  source.dim    = [length(source.xgrid) length(source.ygrid) length(source.zgrid)];
end

source.vol = vol;
if exist('grad', 'var')
  source.grad = grad;
elseif exist('elec', 'var')
  source.elec = elec;
end

if istimelock
  % add the time axis to the output
  source.time = data.time;
elseif iscomp
  % FIXME, add the component numbers to the output
elseif isfreq
  % add the frequency axis to the output
  cfg.frequency    = data.freq(nearest(data.freq, cfg.frequency));
  source.frequency = cfg.frequency;
  if isfield(data, 'time') && isfield(cfg, 'latency')
    cfg.latency    = data.time(nearest(data.time, cfg.latency));
    source.latency = cfg.latency;
  end
  if isfield(data, 'cumtapcnt'),
    source.cumtapcnt = data.cumtapcnt;
  end
end

if exist('dip', 'var')
  % do some cleaning up, keep the dipole positions etc. in the global structure and not in each trial
  source.pos     = dip(1).pos;
  source.inside  = dip(1).inside;
  source.outside = dip(1).outside;
  dip = rmfield(dip, 'pos');
  dip = rmfield(dip, 'inside');
  dip = rmfield(dip, 'outside');
  if isfield(dip(1), 'leadfield')
    source.leadfield = dip(1).leadfield;
    dip = rmfield(dip, 'leadfield');
  end
elseif exist('grid', 'var')
  % no scanning has been done, probably only the leadfield has been computed
  try, source.pos       = grid.pos;       end
  try, source.inside    = grid.inside;    end
  try, source.outside   = grid.outside;   end
  try, source.leadfield = grid.leadfield; end
end

if strcmp(cfg.keepleadfield, 'yes') && ~isfield(source, 'leadfield')
  % add the precomputed leadfields to the output source
  source.leadfield = grid.leadfield;
elseif strcmp(cfg.keepleadfield, 'no') && isfield(source, 'leadfield')
  % remove the precomputed leadfields from the output source
  source = rmfield(source, 'leadfield');
end

% remove the precomputed leadfields from the cfg regardless of what keepleadfield is saying
if isfield(cfg, 'grid') && isfield(cfg.grid, 'leadfield')
  cfg.grid = rmfield(cfg.grid, 'leadfield');
end

if convertfreq
  % FIXME, convert the source reconstruction back to a frequency representation
elseif convertcomp
  % FIXME, convert the source reconstruction back to a component representation
end

if strcmp(cfg.jacknife, 'yes')
  source.method = 'jacknife';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.bootstrap, 'yes')
  source.method = 'bootstrap';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.pseudovalue, 'yes')
  source.method = 'pseudovalue';
  source.trial = dip;
elseif strcmp(cfg.singletrial, 'yes')
  source.method = 'singletrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.rawtrial, 'yes')
  source.method = 'rawtrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.randomization, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'randomization';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numrandomization));
  source.trialB = dip(2+2*(1:cfg.numrandomization));
elseif strcmp(cfg.permutation, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'permutation';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numpermutation));
  source.trialB = dip(2+2*(1:cfg.numpermutation));
elseif exist('dip', 'var')
  % it looks like beamformer analysis was done on an average input, keep the average source reconstruction
  source.method = 'average';
  source.avg = dip;
else
  % apparently no computations were performed
end

if (strcmp(cfg.jacknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes')) && strcmp(cfg.keeptrials, 'yes')
  % keep the source reconstruction for each repeated or resampled trial
  source.trial	= dip;
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: sourceanalysis.m,v 1.94 2006/02/24 16:41:39 roboos Exp $';
% remember the configuration details of the input data
if nargin==2
  try, cfg.previous    = data.cfg;     end
elseif nargin==3
  cfg.previous = [];
  try, cfg.previous{1} = data.cfg;     end
  try, cfg.previous{2} = baseline.cfg; end
end
% remember the exact configuration details in the output 
source.cfg = cfg;

fprintf('total time in sourceanalysis %.1f seconds\n', toc);

