function [source] = dipolefitting(cfg, data)

% DIPOLEFITTING perform grid search and non-linear fit with one or multiple
% dipoles and try to find the location where the dipole model is best able
% to explain the measured EEG or MEG topography.
%
% This function will initially scan the whole brain with a single dipole on
% a regular coarse grid, and subsequently start at the most optimal location
% with a non-linear search. Alternatively you can specify the initial
% location of the dipole(s) and the non-linear search will start from there.
%
% Use as
%   [source] = dipolefitting(cfg, data)
%
% The configuration has the following general fields
%   cfg.numdipoles  = number, default is 1
%   cfg.nonlinear   = 'yes' or 'no'
%   cfg.gridsearch  = 'yes' or 'no'
%   cfg.symmetry    = 'x' or 'y' symmetry for two dipoles (can be empty)
%   cfg.channel     = cell-array (default is 'all', see CHANNELSELECTION)
%   cfg.hdmfile     = filename of headmodel, or
%   cfg.vol         = headmodel
%   cfg.gradfile    = filename of gradiometer definition for MEG, or
%   cfg.grad        = gradiometer definition
%   cfg.elecfile    = filename of electrode locations for EEG, or
%   cfg.elec        = electrode locations and labels
%
% If you start with a grid search, you should specify
%   cfg.grid        = predefined grid for scanning, see PREPARE_LEADFIELD
% or
%   cfg.xgrid       = xmin:xstep:xmax for initial scanning, can be 'auto'
%   cfg.ygrid       = ymin:ystep:ymax for initial scanning, can be 'auto'
%   cfg.zgrid       = zmin:zstep:zmax for initial scanning, can be 'auto'
%   cfg.resolution  = number, for automatic grid generation
%
% If you do not start with a grid search, you have to give a starting location
% for the nonlinear search
%   cfg.dip.pos     = initial dipole position, matrix of Ndipoles x 3
%
% The conventional approach is to fit dipoles to event-related averages, which
% within fieldtrip can be obtained from the TIMELOCKANALYSIS or from
% the TIMELOCKGRANDAVERAGE function. This has the additional options
%   cfg.latency     = [begin end] or 'all'
%   cfg.model       = 'moving', 'regional' or 'fixed'
% A moving dipole model has a different position (and orientation) for each
% timepoint, or for each component. A regional dipole model has the same
% position for each timepoint or component, and a different orientation. A
% fixed dipole model has the same position and direction for all timepoints.
%
% You can also fit dipoles to the spatial topographies of an independent
% component analysis, obtained from the COMPONENTANALYSIS function.
% This has the additional options
%   cfg.component   = array with numbers (can be empty -> all)
%
% You can also fit dipoles to the spatial topographies that are present
% in the data in the frequency domain, which can be obtained using the
% FREQANALYSIS function. This has the additional options
%   cfg.frequency   = should be a single number (in Hz)
%
% See also SOURCEANALYSIS, PREPARE_LEADFIELD

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: dipolefitting.m,v $
% Revision 1.1  2007/01/04 21:35:06  toby
% Initial revision
%
% Revision 1.23  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.22  2005/12/21 09:26:29  roboos
% implemented support for a fixed dipole model, only dipole moment is computed (i.e. no gridsearch or nonlinear fit is done)
% some changes in whitespace by auto indentation
%
% Revision 1.21  2005/11/01 09:37:46  roboos
% removed local subfunction compute_dipole_fit() and replaced by a
% call to the new and general dipole_fit() function which supports
% both EEG and MEG
%
% Revision 1.20  2005/10/14 16:23:11  roboos
% added support for dipole fitting in the frequency domain (based on fourier or csd input)
% restructured the handling of comp and freq input, now both consistent with timelock input
% changed the channel and latency selection of the topographies
%
% Revision 1.19  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.18  2005/06/01 09:15:36  roboos
% replaced fprintf for each grid point with more controlled progress indicator, using cfg.feedback
%
% Revision 1.17  2005/06/01 08:00:59  roboos
% some small changes to the help documentation
%
% Revision 1.16  2005/05/17 17:50:37  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
%
% Revision 1.15  2005/03/14 10:44:54  roboos
% added informative fprintf about optimization toolbox (no functional change)
%
% Revision 1.14  2005/02/21 07:56:54  roboos
% deleted the subfunction compute_leadfield, it is now a separate function
%
% Revision 1.13  2004/12/08 18:00:14  roboos
% implemented consistent method of selecting a subset of channels for
% forward and inverse computations using cfg.channel and updated the
% ducumentation
%
% Revision 1.12  2004/09/06 07:47:03  roboos
% fixed bug for finding optimium after grid scanning fixed dipole model
%
% Revision 1.11  2004/09/01 17:59:28  roboos
% added copyright statements to all filed
% added cfg.version to all functions that give configuration in their output
% added cfg.previous to all functions with input data containing configuration details
%
% Revision 1.10  2004/08/30 13:18:06  roboos
% fixed 2 bogs for moving dipole model
%
% Revision 1.9  2004/06/21 20:13:25  roberto
% fixed two bugs: one with sens.pnt and the other with inside/outside dipole positions after scanning
%
% Revision 1.8  2004/06/03 15:47:21  roberto
% improved detecton of correct configuration (numdipoles)
%
% Revision 1.7  2004/05/19 15:39:35  roberto
% fixed bug in time, which can be cell-array in case of ICA on raw data
% added complete version details to output configuration
%
% Revision 1.6  2004/05/13 23:35:47  roberto
% fixed multiple bugs related to new implementation moving/regional
%
% Revision 1.5  2004/05/13 22:36:30  roberto
% implemented cfg.method = 'regional' (old behaviour) and 'moving' (additional)
%
% Revision 1.4  2004/03/06 13:04:03  roberto
% fixed some small bugs and added time/component number to the output
%
% Revision 1.3  2004/02/05 12:20:06  roberto
% moved the code that was common with dipolefitting into separate prepare_grid and prepare_sens_vol functions
% fixed multiple small bugs
%
% Revision 1.2  2004/01/27 14:39:05  roberto
% multiple small bug fixes
%
% Revision 1.1  2004/01/22 21:42:52  roberto
% initial version, not tested
% support for single dipole scanning, fitting and symmetric dipole pair
%

% test whether fminunc from the Matlab optimization toolbox is present
% the low-level dipole fitting code will use either fminunc or fminsearch,
% which is present in the default Matlab toolboxes.
    
hasoptim = exist('fminunc');
if hasoptim
  fprintf('using optimization toolbox\n');
else
  fprintf('optimization toolbox not found, using fminsearch\n');
end

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

% set the defaults
if ~isfield(cfg, 'channel'),     cfg.channel = 'all';        end
if ~isfield(cfg, 'component'),   cfg.component = [];         end  % for comp input
if ~isfield(cfg, 'frequency'),   cfg.frequency = [];         end  % for freq input
if ~isfield(cfg, 'latency'),     cfg.latency = 'all';        end
if ~isfield(cfg, 'feedback'),    cfg.feedback = 'text';      end
if ~isfield(cfg, 'gridsearch'),  cfg.gridsearch = 'yes';     end
if ~isfield(cfg, 'nonlinear'),   cfg.nonlinear = 'yes';      end
if ~isfield(cfg, 'resolution'),  cfg.resolution = 1;         end
if ~isfield(cfg, 'symmetry'),    cfg.symmetry = [];          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

% the default for this depends on the data type
if ~isfield(cfg, 'model'),
  if ~isempty(cfg.component)
    % each component is fitted independently
    cfg.model = 'moving';
  elseif ~isempty(cfg.latency)
    % fit the data with a dipole at one location
    cfg.model = 'regional';
  elseif ~isempty(cfg.latency)
    % fit the data with a dipole at one location
    cfg.model = 'regional';
  end
end

if ~isfield(cfg, 'numdipoles')
  if isfield(cfg, 'dip')
    cfg.numdipoles = size(cfg.dip(1).pos,1);
  else
    cfg.numdipoles = 1;
  end
end

% set up the symmetry constraints
if ~isempty(cfg.symmetry)
  if cfg.numdipoles~=2
    error('symmetry constraints are only supported for two-dipole models');
  elseif strcmp(cfg.symmetry, 'x')
    cfg.constr.reduce = [1 2 3];
    cfg.constr.expand = [1 2 3 1 2 3];
    cfg.constr.mirror = [1 1 1 -1 1 1];
  elseif strcmp(cfg.symmetry, 'y')
    cfg.constr.reduce = [1 2 3];
    cfg.constr.expand = [1 2 3 1 2 3];
    cfg.constr.mirror = [1 1 1 1 -1 1];
  else
    error('unrecognized symmetry constraint');
  end
else
  % no constraints have been specified
  cfg.constr = [];
end

if isfield(data, 'topolabel')
  % this looks like a component analysis
  iscomp = 1;
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = comp2timelock(cfg, data);
else
  iscomp = 0;
end

if isfield(data, 'freq')
  % this looks like a frequency analysis
  isfreq = 1;
  % transform the data into a representation on which the timelocked dipole fit can perform its trick
  data = freq2timelock(cfg, data);
else
  isfreq = 0;
end

% prepare the volume conduction model and the sensor array
% this updates the configuration with the appropriate fields
[vol, sens, cfg] = prepare_vol_sens(cfg, data);

% select the desired channels
[seldata, selsens] = match_str(data.label, sens.label);
Vdata = data.avg(seldata, :);

if iscomp
  % select the desired component topographies
  Vdata = Vdata(:, cfg.component);
elseif isfreq
  % the desired frequency topographies have already been selected
  Vdata = Vdata(:, :);
else
  % select the desired latencies
  if ischar(cfg.latency) && strcmp(cfg.latency, 'all')
    cfg.latency = data.time([1 end]);
  end
  tbeg = nearest(data.time, cfg.latency(1));
  tend = nearest(data.time, cfg.latency(end));
  cfg.latency = [data.time(tbeg) data.time(tend)];
  Vdata = Vdata(:, tbeg:tend);
end

nchans = size(Vdata,1);
ntime  = size(Vdata,2);
fprintf('selected %d channels\n', nchans);
fprintf('selected %d topographies\n', ntime);

if nchans<3*cfg.numdipoles
  warning('not enough channels to perform a dipole fit');
end

if ntime<1
  error('no spatial topography selected');
end

% check whether EEG is average referenced
if strcmp(sens.type, 'electrodes')
  if rv(Vdata, avgref(Vdata))>0.01
    warning('the EEG data is not average referenced, correcting this');
  end
  Vdata = avgref(Vdata);
end

if strcmp(cfg.model, 'regional')
  % set to zeros if no initial position was specified for the dipoles
  if ~isfield(cfg, 'dip')
    cfg.dip.pos = zeros(cfg.numdipoles, 3);
  end
  % check the input configuration
  if size(cfg.dip.pos,1)~=cfg.numdipoles
    error('inconsistent number of dipoles in configuration')
  end
elseif strcmp(cfg.model, 'moving')
  if ~isfield(cfg, 'dip')
    for t=1:ntime
      % set to zeros if no initial position was specified for the dipoles
      cfg.dip(t).pos = zeros(cfg.numdipoles, 3);
    end
  end
  for t=1:ntime
    % check the input configuration
    if size(cfg.dip(t).pos,1)~=cfg.numdipoles
      error('inconsistent number of dipoles in configuration')
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the dipole scan, this is usefull for generating an initial guess
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.gridsearch, 'yes')
  % test whether we have a valid configuration for dipole scanning
  if cfg.numdipoles==1
    % this is ok
  elseif cfg.numdipoles==2 && ~isempty(cfg.constr)
    % this is also ok
  else
    error('dipole scanning is only possible for a single dipole or a symmetric dipole pair');
  end

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

  if ~isempty(cfg.constr)
    % expand the number of parameters from one (3) to two dipoles (6)
    grid.pos = grid.pos(:,cfg.constr.expand) .* repmat(cfg.constr.mirror, size(grid.pos,1), 1);
  end
%  progress('init', cfg.feedback, 'scanning grid'); %toby: Commented out because it crashes Matlab 7.3 1/4/2006
  for i=1:length(grid.inside)
%    progress(i/length(grid.inside), 'scanning grid location %d/%d\n', i, length(grid.inside)); %toby 1/4/2006
    indx = grid.inside(i);
    if isfield(grid, 'leadfield')
      % reuse the previously computed leadfield
      lf = grid.leadfield{indx};
    else
      pos = reshape(grid.pos(indx,:), [3 cfg.numdipoles])';
      lf = compute_leadfield(pos, sens, vol);
    end

    % the model is V=lf*mom+noise, therefore mom=pinv(lf)*V estimates the
    % dipole moment this makes the model potential U=lf*pinv(lf)*V and the
    % model error is norm(V-U) = norm(V-lf*pinv(lf)*V) = norm((eye-lf*pinv(lf))*V)
    if strcmp(cfg.model, 'regional')
      % sum over all latencies
      grid.error(indx,1) = sum(sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2));
    elseif  strcmp(cfg.model, 'moving')
      % remember error for every latency independently
      grid.error(indx,:) = sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2);
    elseif strcmp(cfg.model, 'fixed')
      error('gridsearch is not supported for a fixed dipole model');
    end
  end
%  progress('close');   %toby 1/4/2006
  % find the grid point(s) with the minimum error
  [err, indx] = min(grid.error(grid.inside,:), [], 1);
  indx = grid.inside(indx);
  for i=1:length(indx)
    dip(i).pos = grid.pos(indx(i),:);
    fprintf('found minimum after scanning on grid point [%f %f %f]\n', dip(i).pos(1), dip(i).pos(2), dip(i).pos(3));
  end
else
  % use the initial guess supplied in the configuration
  dip = cfg.dip;
end

if ~isfield(dip, 'mom')
  % this must be added, otherwise there is an error in the assignment of
  % the dipole fit result later
  for t=1:length(dip)
    dip(t).mom = [];
  end
end

if strcmp(cfg.nonlinear, 'yes')
  if strcmp(cfg.model, 'regional')
    % perform the non-linear dipole fit for all latencies together
    try
      dip = dipole_fit(dip, Vdata, sens, vol, cfg.constr);
    catch
      % catch errors due to non-convergence
      dip.pos = nan*dipin.pos;
      dip.mom = [];
    end
  elseif strcmp(cfg.model, 'moving')
    % perform the non-linear dipole fit for each latency independently
    for t=1:ntime
      try
        dip(t) = dipole_fit(dip(t), Vdata(:,t), sens, vol, cfg.constr);
      catch
        % catch errors due to non-convergence
        dip(t).pos = nan*dipin.pos;
        dip(t).mom = [];
      end
    end
  elseif strcmp(cfg.model, 'fixed')
    error('non-linear fit is not supported for a fixed dipole model');
  end
end

if strcmp(cfg.model, 'regional')
  % re-compute the leadfield in order to compute the model potential and dipole moment
  lf = compute_leadfield(dip.pos, sens, vol);
  % compute all details of the final dipole model
  dip.mom = pinv(lf)*Vdata;
  dip.pot = lf*dip.mom;
  dip.rv  = rv(Vdata, dip.pot);

elseif strcmp(cfg.model, 'moving')
  for t=1:ntime
    % re-compute the leadfield in order to compute the model potential and dipole moment
    lf = compute_leadfield(dip(t).pos, sens, vol);
    % compute all details of the final dipole model
    dip(t).mom = pinv(lf)*Vdata(:,t);
    dip(t).pot = lf*dip(t).mom;
    dip(t).rv  = rv(Vdata(:,t), dip(t).pot);
  end % ntime

elseif strcmp(cfg.model, 'fixed')
  % re-compute the leadfield in order to compute the model potential and dipole moment
  lf = [];
  for i=1:cfg.numdipoles
    tmp = compute_leadfield(dip.pos(i,:), sens, vol);
    tmp = tmp * dip.mom(i,:)';
    lf = [lf tmp];
  end
  % compute all details of the final dipole model
  dip.mom = pinv(lf)*Vdata;
  dip.pot = lf*dip.mom;
  dip.rv  = rv(Vdata, dip.pot);

end % model

if isfreq && strcmp(cfg.model, 'regional')
  % the matrix with the dipole moment is encrypted and cannot be interpreted straight away
  % reconstruct the frequency representation of the data at the source level
  [dip.pow, dip.csd, dip.fourier] = timelock2freq(dip.mom);
elseif isfreq && strcmp(cfg.model, 'moving')
  % although this is technically possible sofar, it does not make any sense
  warning('a moving dipole model in the frequency domain is not supported');
elseif isfreq && strcmp(cfg.model, 'fixed')
  % although this is technically possible sofar, it does not make any sense
  warning('a fixed dipole model in the frequency domain is not supported');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
source.label  = cfg.channel; % these channels were used in fitting
source.dip    = dip;
source.Vdata  = Vdata;

% assign a latency, frequeny or component axis to the output
if iscomp
  source.component = cfg.component;
elseif isfreq
  source.freq = cfg.frequency;
else
  tbeg = nearest(data.time, cfg.latency(1));
  tend = nearest(data.time, cfg.latency(end));
  source.time = data.time(tbeg:tend);
end

% add the version details of this function call 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: dipolefitting.m,v 1.1 2007/01/04 21:35:06 toby Exp $';
% remember the configuration details of the input data
try, cfg.previous = data.cfg; end
% remember the exact configuration details in the output
source.cfg = cfg;
