function Results = lcmvbf(OPTIONS);
%LCMVBF - Spatial filtering solution as called from command line or BST_SOURCEIMAGING.M
% function Results = lcmvbf(OPTIONS);
% -------------------------------------------------------------------------------------------------
% NOTICE:
% This function is not optimized for stand-alone command calls.
% Please use the generic BST_SOURCEIMAGING function to access Min.Norm. source map computations.
% -------------------------------------------------------------------------------------------------
%
%
% INPUTS
%       OPTIONS - a structure of parameters
% OUTPUTS
%       Results - a regular BrainStorm Results structure (see documentation for details)
%
% Details of the OPTIONS structure
%
%   Data Information:
%          .DataFile   : A string actual data file requested for localization. If isempty, look for data array(s) in .Data
%          .Data       : Array of nChannel x nTime array(s) of data.
%          .SegmentIndices: Vector of all time samples (in sample values)
%                           of the latencies specifying the time window for analysis.
%          .Trial      : A scalar specifying the trial from which to extract the data segment
%                        (not empty only when data in .DataFile is stored in native format)
%          .ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
%          .DataTypes  : A scalar or a vector of codes to specidy which is the type of data to use:
%                        1 is MEG, 2 is EEG, 3 is Fusion
%                        (Default is all data type available)
%          .FusionType : Code for the type of fusion to apply between MEG and EEG data
%                        1 is linear, 2 is non-linear
%                        (Defaut is 1)
%          .Data       : Specify some data to be passed independently of any data file
%                            .F : Cell array of nChannel x nTime array(s) of data.
%                                 If not empty: uses this array as input for inverse process.
%                            .ChannelFile : Channel file corresponding to data in .F
%                            .Time : Vector specifying the time instants of data samples in .Data.F
%                        (Default for .Data is empty)
%
%   Source Support Information:
%          .GridLoc    : a string or a 3xN matrix
%                        - if a string : specifies the name of the BrainStorm tessellation file where cortical support is stored
%                        - if a matrix : contains the coordinates of the cortical locations where to compute the cortical current density
%          .iGrid      : when .GridLoc is a file name, specify which imaging grid to use (Default is 1).
%
%   Forward Field Information :
%          .HeadModelFile: Cell array of strings for filename of the BrainStorm headmodel file(s) to use for imaging. If is empty, look for.GainFile field.
%                          There are only 2 options at this point:
%                                - .HeadModelFile is of length 1, a single headmodel file is used for all data sets passed to the bst_sourceimaging.
%                                - .HeadModelFile is of length the number of data sets passed to bst_sourceimaging (i.e. the length of .DataFile).
%          .GainFile   : Cell array of filenames of image gain matrix (.bin) for the selected cortical support. If is empty, look for .Gain field.
%                        Remarks about the length of .HeadModelFile hold for .GainFile too.
%          .Gain       : Cell array containing the gain matrix of the cortical support (either a single cell or a cell array the length the number of data sets passed to function)
%                        Remarks about the length of .HeadModelFile hold for .GainFile too.
%                        (Default is empty).
%
%   Processing Parameters :
%          .Method     : A string that specifies the imaging method
%                        'Minimum-Norm Imaging': MN estimate of current density using MN priors on source amplitudes (Default)
%                        'Recursive Min. Norm.': Recursive MN estimate based on adaptation of FOCUSS (NOT AVAILABLE)
%                        'SPTF' : Scanning technique constrained on cortical surface (see Denis Schwartz) (UNDER DEVELOPMENT)
%                        'ST-MAP' : Non-linear imaging technique based on sparse-focal image models (NOT AVAILABLE)
%                        'RAP-MUSIC' : link to Mosher's GUI (UNDER DEVELOPMENT)
%                        (Default is 'Minimum-Norm Imaging')
%          .Tikhonov   : Hyperparameter used in Tykhonov regularization
%                        (Default is 10)
%          .Rank       : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
%                        (Default is full rank)
%          .ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
%                         (Default is 0: compute full ImageGridAmp array)
%
%   Output Information:
%          .ResultFile : optional string for filename where to write results
%                        if 'default', filename as used in GUI environment is used
%                        (Default is '' (no Result file created))
%          .Verbose    : Turn ON/OFF verbose mode
%                        (Default is 1 (ON))
%          .GUIHandles : A vector of GUI handles to the SOURCEIMAGING GUI
%                        when operating in graphic mode: passed automatically to bst_sourceimaging by the source imaging switchyard
%                        in command line mode: leave this field empty (default).
%          .NNormalization : Apply noise normalization
%                        (Default is 0)

%<autobegin> ---------------------- 12-Oct-2004 12:01:16 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Unknown Category
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\colnorm.m
%   toolbox\good_channel.m
%   toolbox\inorcol.m
%   toolbox\meg4read.m
%   toolbox\norcol.m
%   toolbox\read_gain.m
%   toolbox\regcheck.m
%
% Subfunctions in this file, in order of occurrence in file:
%   varargout = gainmat_covar(varargin)
%
% At Check-in: $Author: Mosher $  $Revision: 3 $  $Date: 10/12/04 10:23a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 12-Oct-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 12-Oct-2004 12:01:16 -----------------------



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% /---Script Author--------------------------------------\
% | *** Esen Kucukaltun-Yildirim                           |
% |  USC-SIPI                                              |
% \------------------------------------------------------/


% Date of creation: October 2004 (adapted from minnorm.m)
%
% ----------------------------- Script History ---------------------------------
% EK  7-Oct-2004  New outputput formats 'Source Power' and Noise Normalized
%                 Filter Outputs are added
% JCM 12-Oct-2004 Header comments
% ----------------------------- Script History ---------------------------------


% General hard-wired parameter values:--------------------------------------------------------------
DataType = {'MEG','EEG'}; % String labels corresponding to OPTIONS.DataTypes flags
OPTIONS.BLOCK = 10000; % process by blocks
OPTIONS.ExpoiG = .4; % Exponential weighting for column normalization
%---------------------------------------------------------------------------------------------------

% Load Headmodel information -----------------------------------------------------------------------

% From file, if passed to function
if ~isempty(OPTIONS.HeadModelFile)
    HeadModel = load(OPTIONS.HeadModelFile,'Function','Param',...
        'GridLoc','Gain'); % Load Basics

    if ~isempty(OPTIONS.GridLoc)
        HeadModel.GridLoc = {OPTIONS.GridLoc};
    end
    if ~isempty(OPTIONS.Gain)
        HeadModel.Gain = OPTIONS.Gain;
    end

    HeadModelFileWhos = whos('-file',OPTIONS.HeadModelFile);

    clear tmp

    % Gain matrix full filename
    gainfile = fullfile(fileparts(OPTIONS.HeadModelFile),HeadModel.Gain{1});

end

% Load Data Information -----------------------------------------------------------------------

if(~isempty(findstr('_results',deblank(lower(OPTIONS.DataFile))))),
    % CHEAT
    % The DataName has the string '_results' in it
    % We are trying to treat independent topographies as min norm data. Eventually this
    %  capability will be moved to a new routine.
    % name has the string '_results' in it
    Results = load(OPTIONS.DataFile);
    % now synthesize a Data structure from this
    [ignorm,Anrm] = colnorm(Results.IndepTopo);

    Data = struct('F',Anrm,'ChannelFlag',Results.ChannelFlag,...
        'Time',[1:size(Results.IndepTopo,2)],... % time is now simply an indexer
        'NoiseCov',Results.NoiseCov,'SourceCov',Results.SourceCov,'Projector',Results.Projector,...
        'Comment',sprintf('Synthesized Data Set from Results: %s',Results.Comment));
else
    if isempty(OPTIONS.Data)
        % assume it is a standard data file
        DataArray = 0; % Correponsing flag
        Data = load(OPTIONS.DataFile,'F','ChannelFlag','Time','Projector','Comment');
    else
        DataArray = 1;
        Data = OPTIONS.Data;
        Data.Projector = [];
        Data.Comment = '';
        Data.ChannelFlag = ones(size(Data.F,1),1);
    end

end

% Data & Channel information ------------------------------------------------------------------
if isfield(OPTIONS,'ChannelFile')
    Channel = load(OPTIONS.ChannelFile);
    OPTIONS.Channel = Channel.Channel; clear Channel
end

% now alter the data according to the bad channels
OPTIONS.GoodChannel = good_channel(OPTIONS.Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality
OPTIONS.EEGChannel = good_channel(OPTIONS.Channel,ones(size(OPTIONS.Channel)),'EEG');
OPTIONS.MEGChannel = good_channel(OPTIONS.Channel,ones(size(OPTIONS.Channel)),'MEG');

if ischar(Data.F) % Data is in native file format
    VERBOSE = 1;
    Data.F = meg4read(Data.F,OPTIONS.Trial,OPTIONS.SegmentIndices,OPTIONS.Verbose);
    Data.FBaseline = meg4read(Data.F,OPTIONS.Trial,OPTIONS.NoiseSegmentIndices,OPTIONS.Verbose);
else
    % Trim the data to the requested segments
    if OPTIONS.NNormalization
        Data.FBaseline = Data.F(:,OPTIONS.NoiseSegmentIndices);
    end
    Data.F = Data.F(:,OPTIONS.SegmentIndices);
end


try
    Data.F = Data.F(OPTIONS.GoodChannel,:); % good channels only
    if OPTIONS.NNormalization
        Data.FBaseline = Data.FBaseline(OPTIONS.GoodChannel,:); % good channels for noise segment too
    end

catch
    errordlg(sprintf('Data array has %d rows while there are %d %s channels. Possible problem is wrong type assignement of reference vs. measurement channels in channel file %s',...
        size(Data.F,1),length(OPTIONS.GoodChannel), DataType{OPTIONS.DataType}, OPTIONS.ChannelFile ),'Inconsistent channel and data array sizes')
    return
end


Time = Data.Time(OPTIONS.SegmentIndices);


% Is there any projector for the data ?
if(~isempty(Data.Projector)),
    Data.Projector = Data.Projector(OPTIONS.GoodChannel,:);
end

% Get HeadModel parmaters for GoodChannels
HeadModel.Param = HeadModel.Param(OPTIONS.GoodChannel);



% Update MEG and EEG channel indices by removing bad channels.
OPTIONS.EEGChannel = intersect(OPTIONS.GoodChannel,OPTIONS.EEGChannel);
OPTIONS.MEGChannel = intersect(OPTIONS.GoodChannel,OPTIONS.MEGChannel);

if OPTIONS.DataType == 2 & ...
        length(OPTIONS.GoodChannel) < length(OPTIONS.EEGChannel) % EEG data with bad channels

    % Adjust average reference contribution from covariance matrix
    if isempty(find(mean(Data.F >eps)))
        % Average Reference (don't use the channel.comment field that may not have been altered following manipulations in the data viewer tool)
        % Load gain matrix
        G = read_gain(gainfile,[]);
        if strcmp(OPTIONS.Col_norm,'y');
            G = G*inorcol(G);
        end

        Gm = mean(G,1)';

        I = ones(length(OPTIONS.GoodChannel),1);
        % Reopen the Gain file that was closed during call to read_gain
        fid = fopen(...
            fullfile(fileparts(OPTIONS.HeadModelFile),HeadModel.Gain{1}),...
            'r','ieee-be'); % HeadModel.Gain{1} - CHEAT Dipoles only
    else
        % Do nothing
    end
end

OPTIONS = regcheck(OPTIONS); % returns OPTIONS.REG set to the appropriate field

% Did the user provide an existing control?
if(isfield(Data,'Projector')),
    %Decompose it:
    if(isempty(Data.Projector)),
        A = [];
        Ua = [];
    else
        A = Data.Projector; % initialize
        Ua = orth(A);
    end
else % user did not give a projector
    A = [];
    Ua = [];
end

if isempty(OPTIONS.Rank)
    OPTIONS.Rank = min(size(Data.F));
end

if(OPTIONS.Rank < min(size(Data.F))), % reduced rank

    % User has asked for an rank RANK subspace search.
    % Form the signal subspace

    if(size(Data.F,1) >= size(Data.F,2)), % skinny matrix
        [Uf,Sf,Vf] = svd(Data.F,0); % economy
    else % wide matrix
        [Vf,Sf,Uf] = svd(Data.F',0);
    end

    % form the reduced rank data
    Data.F = Uf(:,1:OPTIONS.Rank)*(Sf(1:OPTIONS.Rank,1:OPTIONS.Rank)*Vf(:,1:OPTIONS.Rank)');

    % project away if needed
    if(~isempty(Ua)),
        Data.F = Data.F - Ua*(Ua'*Data.F);
    end

    clear Ua Sa Va Uf Sf Vf
end


%                             --------------------------------------------------------------------------------------

% Begin LCMV BF estimate --------------------------------------------------------------------------------------

%                             --------------------------------------------------------------------------------------

Kernel = read_gain(gainfile,OPTIONS.GoodChannel,[]);
nSources = size(Kernel,2);
nChannels = size(Kernel,1);

Cm = cov(Data.F');
%invert data covariance
[U, S, V] = svd(Cm);
% if user select not to regularize .Tikhonov is 0 and it is a direct
% inverse
S = diag(S);
Lambda = S(1)*OPTIONS.Tikhonov/100 ;
Si = diag(1../(S+Lambda)); % filtered inverse
Cm_inv = V*Si*U';   clear U S Si V
%design beamformer
source_power_inv = sum( Kernel.*(Cm_inv*Kernel) );
if get(OPTIONS.GUIHandles.sourcePower, 'Value') ==1 %Calculate source power
    bst_message_window('Calculating the Source Power. . .');
    
    ImageGridAmp = (1 ./ source_power_inv)' ;
    ImagingKernel = [];
    Fsynth = [];
elseif get(OPTIONS.GUIHandles.normalizedOut, 'Value') ==1
    bst_message_window('Calculating the spatial filter. . .');
    spatialFilter = Kernel'*Cm_inv;
    spatialFilter= (spatialFilter) ./ repmat(source_power_inv', 1,size(spatialFilter,2));
    ImageGridAmp=[];
    Fsynth = [];

    tmp = (spatialFilter*Data.FBaseline);
    mean_Baseline= mean(tmp, 2);
    std_Baseline = std(tmp, 0, 2);
    tmp = (spatialFilter);
    %     ImagingKernel = (tmp - repmat(mean_Baseline, 1, size(tmp, 2))) ./  repmat(std_Baseline, 1, size(tmp, 2));
    ImagingKernel = tmp  ./  repmat(std_Baseline, 1, size(tmp, 2));
    ImagingKernel = single(ImagingKernel);   
    
    
elseif OPTIONS.NNormalization  %Calculate the neural index
    bst_message_window('Calculating the Neural activity index. . .');
    Cn = cov(Data.FBaseline');
    [U, S, V] = svd(Cn);S = diag(S);
    Lambda = S(1)*OPTIONS.Tikhonov/100 ;
    Si = diag(1../(S+Lambda)); % filtered inverse
    Cn_inv = V*Si*U';   clear U S Si V
    
    noise_power_inv = sum( Kernel.*(Cn_inv*Kernel) );
    ImageGridAmp = (noise_power_inv ./ source_power_inv)' ;
    ImagingKernel = [];
    Fsynth = [];
elseif OPTIONS.NNormalization==0 %Calculate the current density kernel multiply by the data Data.F to get the reconstruction map
    bst_message_window('Calculating the spatial filter. . .');
    spatialFilter = Kernel'*Cm_inv;
    ImagingKernel= (spatialFilter) ./ repmat(source_power_inv', 1,size(spatialFilter,2));
    ImageGridAmp=[]; clear spatialFilter
    Fsynth = [];
end


if OPTIONS.Verbose
    bst_message_window({...
        'Inverse Transformation -> DONE.'...
        })
end

%
% beamformer estimate - DONE --------------------------------------------------------------------------------------------------------------------------


% Save Results in Result file -----------------------------------------------------------------------------------------------------------------
%
% Fields from Results structure
Channel = OPTIONS.Channel;
ChannelFlag = Data.ChannelFlag;
Comment = sprintf('MIN NORM, rank %.0f',OPTIONS.Rank);
DataFile = OPTIONS.DataFile;
Date = datestr(datenum(now),0);
Time = [OPTIONS.SegmentIndices];
ImageGridTime = Data.Time(Time);
ImageGridAmp = ImageGridAmp;  % monster matrix when ImagingKernel is not computed
Function = mfilename; % name of this calling routine
HeadModelFile = OPTIONS.HeadModelFile;
ModelGain = HeadModel.Gain{1};
Projector = Data.Projector;
SourceLoc = HeadModel.GridLoc{1};
SourceOrder = -1;

NoiseCov = [];
SourceCov = [];
PatchNdx = [];
PatchAmp = [];


% Reduced Results structure for call from bst_sourceimaging.
Results = struct('Comment',Comment,'Function',mfilename,...
    'OPTIONS',OPTIONS,'Date','','Time',Time,'HeadModelFile',HeadModelFile,...
    'Channel',Channel,'ChannelFlag',Data.ChannelFlag,'DataFile',OPTIONS.DataFile,...
    'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
    'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
    'SourceOrientation',[],'TimeSeries',[],'ImagingKernel',ImagingKernel,'ModelGain',ModelGain,...
    'PatchNdx',PatchNdx,'PatchAmp',PatchAmp, ...
    'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);

return

% -----------------------------------------------------------------------------------------------------
%
%
%                                             SubFunctions
%
%
% -----------------------------------------------------------------------------------------------------

function varargout = gainmat_covar(varargin)

% GAINMAT_COVAR Computation of gain matrix covariance
% function [GainCovar] = gainmat_covar(gainfile,OPTIONS)
% Computation of gain matrix covariance
%
% INPUTS
%
% gainfile : a string - filename of the binary gain matrix file (imaging grid)
% OPTIONS  : a structure specifying the options for the computation
%            .Verbose : turn on/off verbose mode (1/0)
%            .DataType: scalar specifying the type of data
%                       1 is MEG
%                       2 is EEG
%                       3 is MEG/EEG fusion
%            .FFNormalization : flag specifying whether forward field normalization is requested
%                       0 off
%                       1 on
%            .GoodChannel : Indices of good channels for current DataType
%            .MEGChannel  : Indices of all MEG channels
%            .EEGChannel  : Indcies of all EEG channels
%
% OUTPUTS
% GainCovar : a structure with the proper gain covar matrix outputs
%
% /---Script Author--------------------------------------\
% | *** John C. Mosher, Ph.D.                            |
% |  Biophysics Group                                    |
% |                                                      |
% | *** Sylvain Baillet, Ph.D.                           |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   |
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% \------------------------------------------------------/
%
% Date of creation: October 2002
% Script History -----------------------------------------------------------------------------------
%---------------------------------------------------------------------------------------------------

% Map entries
gainfile = varargin{1};
OPTIONS = varargin{2};

% Read number of sources and total number of channels
tmp = read_gain(gainfile,OPTIONS.GoodChannel(1),[]);
nsrc = length(tmp); clear tmp % Number of sources
tmp = read_gain(gainfile,[],1);
nchan = length(tmp); clear tmp % Number of channels

% if OPTIONS.Verbose
%     bst_message_window(sprintf('Image support: %.0f grid points . . .',nsrc));
% end

AAt = zeros(nchan); %  Initialize GainCovar matrix
Gain_RowNorm = [];  %  Initialize matrix containing norms of gain matrix rows (useless and therefore empty if DataType ~= 3)

for i = 1:OPTIONS.BLOCK:nsrc,

    if OPTIONS.Verbose & ~isempty(OPTIONS.GUIHandles)
        waitbar(i/nsrc,OPTIONS.hwaitbar);
    end

    ndx = [0:(OPTIONS.BLOCK-1)]+i;
    if(ndx(end) > nsrc),  % last block too long
        ndx = [ndx(1):nsrc];
    end

    % Normalize EEG and MEG forward fields separately
    % Fusion Case
    if OPTIONS.DataType == 3 % Compute gain matrix row norm on first block only
        if i == 1
            Gain_RowNorm = zeros(length(OPTIONS.Channel));
            if ~isempty(OPTIONS.GUIHandles)
                hkk = waitbar(0,'Processing Lead-Field Normalization...');
            end
            for kk = OPTIONS.GoodChannel
                Gain_RowNorm(kk,kk) = 1/norm(read_gain(gainfile,kk,1:nsrc));
                if ~isempty(OPTIONS.GUIHandles)
                    waitbar((kk-OPTIONS.GoodChannel(1)+1)/length(OPTIONS.GoodChannel))
                end

            end

            if ~isempty(OPTIONS.GUIHandles)
                close(hkk)
            end
        end

        temp = Gain_RowNorm([OPTIONS.GoodChannel],[OPTIONS.GoodChannel])*read_gain(gainfile,[OPTIONS.GoodChannel],ndx);
        OPTIONS.iG = sqrt(norcol(temp)'); % NORCOL computes the column norms to the SQUARE
        OPTIONS.iG(OPTIONS.iG == 0)= min(OPTIONS.iG(OPTIONS.iG ~= 0));
        OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
        AAt(OPTIONS.GoodChannel,OPTIONS.GoodChannel) = AAt(OPTIONS.GoodChannel,OPTIONS.GoodChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations

    else % MEG | EEG
        temp = read_gain(gainfile,[OPTIONS.MEGChannel],ndx);
        if ~isempty(OPTIONS.MEGChannel) & ~isempty(temp)
            temp(find(isnan(temp))) = 0;
            if OPTIONS.FFNormalization % Forward field normalization is requested
                OPTIONS.iG = sqrt(norcol(temp)');  % NORCOL computes the column norms to the SQUARE
                OPTIONS.iG(OPTIONS.iG == 0)= min(OPTIONS.iG(OPTIONS.iG ~= 0));
                OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
                AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) = AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations
            else % no FF normalization
                AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) = AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) + temp*temp';  % next chunk of correlations
            end
        end

        clear temp OPTIONS.iG

        temp = read_gain(gainfile,OPTIONS.EEGChannel,ndx);
        if ~isempty(OPTIONS.EEGChannel) & ~isempty(temp)
            temp(find(isnan(temp))) = 0;
            if OPTIONS.FFNormalization
                OPTIONS.iG = sqrt(norcol(temp)');  % NORCOL computes the column norms to the SQUARE
                if min(OPTIONS.iG(:)) == max(OPTIONS.iG(:)) % as it happens when EEG channels stand for ECoG or EMG
                    OPTIONS.iG = eps;
                end

                OPTIONS.iG(OPTIONS.iG == 0) = min(OPTIONS.iG(OPTIONS.iG ~= 0));

                OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
                AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) = AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations
            else
                AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) = AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) + temp*temp';  % next chunk of correlations
            end
        end

        clear temp OPTIONS.iG

    end % datatype

end % for each block of sources

varargout{1} = AAt;
varargout{2} = Gain_RowNorm;



% -----------------------------------------------------------------------------------------------------------------------
