vb_dipole_search

 Dipole position serach by minimizing error
   [Vout, err, iter] = ...
        vb_dipole_search(Vinit, bexp, pick, Qpick, Wsensor, errfunc, Option, mode)
 --- Output
 Vout : Estimated position
 err  : residual normalized error
 iter : iteration number

 --- INPUT
 errfunc : error function to calculate dipole estimation error
 Vinit      : initial dipole position (Ndipole x 3)
 Vinit(n,:) : dipole position (3-D coordinate) at n-th vertex

 bexp   : MEG-data   (Nsensor ,Time) 

  pick(k, 1:3) : Sensor coil position  : コイル位置, 
 Qpick(k, 1:3)    : Sensor coil direction : コイル方向

 Wsensor(m,n) = n-th coil weight for m-th sensor channel
    basis(channel,dipole) = Wsensor * basis(coil,dipole)

 Option  : Option in optimization
 mode  : Minimization mode 
    =0 : fminsearch
    =1 : fminunc

 2006-12-16 made by M.Sato
 2007-4-3 made by M.Sato
 ---

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Vout, err, iter] = ...
0002     vb_dipole_search(Vinit, bexp, pick, Qpick, Wsensor, errfunc, Option, mode)
0003 % Dipole position serach by minimizing error
0004 %   [Vout, err, iter] = ...
0005 %        vb_dipole_search(Vinit, bexp, pick, Qpick, Wsensor, errfunc, Option, mode)
0006 % --- Output
0007 % Vout : Estimated position
0008 % err  : residual normalized error
0009 % iter : iteration number
0010 %
0011 % --- INPUT
0012 % errfunc : error function to calculate dipole estimation error
0013 % Vinit      : initial dipole position (Ndipole x 3)
0014 % Vinit(n,:) : dipole position (3-D coordinate) at n-th vertex
0015 %
0016 % bexp   : MEG-data   (Nsensor ,Time)
0017 %
0018 %  pick(k, 1:3) : Sensor coil position  : コイル位置,
0019 % Qpick(k, 1:3)    : Sensor coil direction : コイル方向
0020 %
0021 % Wsensor(m,n) = n-th coil weight for m-th sensor channel
0022 %    basis(channel,dipole) = Wsensor * basis(coil,dipole)
0023 %
0024 % Option  : Option in optimization
0025 % mode  : Minimization mode
0026 %    =0 : fminsearch
0027 %    =1 : fminunc
0028 %
0029 % 2006-12-16 made by M.Sato
0030 % 2007-4-3 made by M.Sato
0031 % ---
0032 %
0033 % Copyright (C) 2011, ATR All Rights Reserved.
0034 % License : New BSD License(see VBMEG_LICENSE.txt)
0035 
0036 Vinit = Vinit(:);
0037 
0038 % error function name for dipole search
0039 if ~exist('errfunc','var')
0040     errfunc = 'vb_dipole_error';
0041 end
0042 %
0043 % Set minimization option
0044 %
0045 % Optimization setting
0046 if isfield(Option,'MaxIter') & ~isempty(Option.MaxIter)
0047     MaxIter = Option.MaxIter; 
0048 else
0049     MaxIter = 5000; 
0050 end;
0051 
0052 if isfield(Option,'TolFun') & ~isempty(Option.TolFun)
0053     TolFun = Option.TolFun; 
0054 else
0055     TolFun = 1.e-10; 
0056 end;
0057 
0058 if isfield(Option,'TolX') & ~isempty(Option.TolX)
0059     TolX = Option.TolX; 
0060 else
0061     TolX = 1.e-8; 
0062 end;
0063 
0064 if isfield(Option,'Display') & ~isempty(Option.Display)
0065     Display = Option.Display; 
0066 else
0067     Display = 'off';
0068 end;
0069 
0070 OPTIONS = optimset('MaxIter', MaxIter, ...
0071                    'MaxFunEvals', 2*MaxIter,...
0072                    'TolX',    TolX, ...
0073                    'TolFun',  TolFun, ...
0074                    'Display', Display, ...
0075                    'LargeScale', 'off', ...
0076                    'GradObj','off'); 
0077 
0078 % Minimize the error function
0079 
0080 BB = bexp*bexp';    % BB(Nsensor,Nsensor)
0081 
0082 % Search method
0083 if ~exist('mode','var') | isempty(mode)
0084     mode = vb_tool_check('Optimization');
0085     %     = 0;    % Use 'fminsearch'
0086     %     = 1;    % Use 'fminunc' in Optimization Toolbox
0087 end
0088 
0089 switch    mode
0090 case    0
0091     % Nelder-Mead Simplex Method
0092     disp('Nelder-Mead Simplex Method');
0093     [Vout, err, flag, outputs]= ...
0094         fminsearch(errfunc,Vinit, OPTIONS, BB, pick, Qpick, Wsensor); 
0095 case    1
0096     % Unconstrained optimization
0097     % Optimization toolbox
0098     disp('Unconstrained optimization')
0099     [Vout, err, flag, outputs]= ...
0100         fminunc(errfunc,Vinit, OPTIONS, BB, pick, Qpick, Wsensor); 
0101 case    2
0102     % Unconstrained optimization
0103     disp('Unconstrained optimization')
0104     [Vout, err, flag, outputs]= ...
0105         fminunc1(errfunc,Vinit, OPTIONS, BB, pick, Qpick, Wsensor); 
0106 end
0107 
0108 iter = outputs.iterations;
0109 
0110 NV      = length(Vout(:))/3;        % Number of dipoles
0111 Vout  = reshape(Vout,[NV 3]);