function [D] = spm_eeg_invert_prepro(D,val)
% Preprocessing for inversion stage
% Includes spatial and temporal dimension reduction.
%
% This version only handles single subject single modality data.
% The removal of many scaling factors makes it easier to compare between
% forward models.
% ReML inversion of multiple forward models for EEG-MEG
% FORMAT [D] = spm_eeg_invert_classic(D)
% ReML estimation of regularisation hyperparameters using the
% spatiotemporal hierarchy implicit in EEG/MEG data
%
% Requires:
% D{i}.inv{val}.inverse:
%
%     inverse.modality - modality to use in case of multimodal datasets
%
%     inverse.trials - D.events.types to invert
%     inverse.type   - 'GS' Greedy search on MSPs
%                      'ARD' ARD search on MSPs
%                      'MSP' GS and ARD multiple sparse priors
%                      'LOR' LORETA-like model
%                      'IID' minimum norm
%                       'EBB' for empirical bayes beamformer
%
%    inverse.priors{} - a cell array of anatomical and functional priors to be
%    considered
%  
%
%     inverse.woi    - time window of interest ([start stop] in ms)
%     inverse.lpf    - band-pass filter - low frequency cut-off (Hz)
%     inverse.hpf    - band-pass filter - high frequency cut-off (Hz)
%     inverse.Han    - switch for Hanning window
%     inverse.xyz    - (n x 3) locations of spherical VOIs
%     inverse.rad    - radius (mm) of VOIs
%
%     inverse.Nm     - maximum number of channel modes
%     inverse.Nr     - maximum number of temporal modes
%     inverse.Np     - number of sparse priors per hemisphere
%     inverse.smooth - smoothness of source priors (0 to 1)
%     inverse.Na     - number of most energetic dipoles
%     inverse.sdv    - standard deviations of Gaussian temporal correlation
%     inverse.Qe     - any sensor error components (e.g. empty-room data)
%     inverse.Qe0    - minimum amount of sensor noise power relative to
%                      signal eg 0.1 would correspond to power SNR of 10.0
%     inverse.A      - predefined spatial modes (Nchans*Nmodes) to project
%                      sensor data through
%
% Evaluates:
%
%     inverse.M      - MAP projector (reduced)
%     inverse.J{i}   - Conditional expectation (i conditions) J = M*U*Y
%     inverse.L      - Lead field (reduced UL := U*L)
%     inverse.qC     - spatial covariance
%     inverse.qV     - temporal correlations
%     inverse.T      - temporal projector
%     inverse.U(j)   - spatial projector (j modalities) - derived from data
%     inverse.A      - pre-specified spatial projector
%     inverse.Y{i}   - reduced data (i conditions) UY = UL*J + UE
%     inverse.Is     - Indices of active dipoles
%     inverse.It     - Indices of time bins
%     inverse.Ic{j}  - Indices of good channels (j modalities)
%     inverse.Nd     - number of dipoles
%     inverse.pst    - peristimulus time
%     inverse.dct    - frequency range
%     inverse.F      - log-evidence
%     inverse.VE     - variance explained in spatial/temporal subspaces (%)
%     inverse.R2     - variance in subspaces accounted for by model (%)
%     inverse.scale  - scaling of data for each of j modalities
%__________________________________________________________________________
%
% This version is for single subject single modality analysis and therefore
% contains none of the associated scaling factors.
% No symmetric priors are used in this implementation (just single patches)
% There is an option for a Beamforming prior : inversion type 'EBB'
% also added new beamforming method- using GS rather than ARD, from Juan 
% David Martinez Vargas ('EBBgs').
%
% The code was used in
% Lopez, J. D., Penny, W. D., Espinosa, J. J., Barnes, G. R. (2012).
% A general Bayesian treatment for MEG source reconstruction incorporating
% lead field uncertainty.
% Neuroimage 60(2), 1194-1204 doi:10.1016/j.neuroimage.2012.01.077.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Jose David Lopez, Gareth Barnes, Vladimir Litvak
% $Id: spm_eeg_invert_prepro.m 6444 2015-05-21 11:15:48Z guillaume $


Nl = length(D);

if Nl > 1
    error('function only defined for a single subject');
else
    D = D{1};
end

% D - SPM data structure
%==========================================================================
if nargin > 1
    D.val = val;
elseif ~isfield(D, 'val')
    D.val = 1;
end


val       = D.val;

inverse   = D.inv{val}.inverse;

% forward model
%--------------------------------------------------------------------------


% defaults
%--------------------------------------------------------------------------

try, type = inverse.type;   catch, type = 'GS';     end
try, s    = inverse.smooth; catch, s    = 0.6;      end
try, Np   = inverse.Np;     catch, Np   = 256;      end
try, Nr   = inverse.Nr;     catch, Nr   = 16;       end % requested number of temporal modes, could be changed depending on svd
try, xyz  = inverse.xyz;    catch, xyz  = [0 0 0];  end
try, rad  = inverse.rad;    catch, rad  = 128;      end
try, hpf  = inverse.hpf;    catch, hpf  = 48;       end % need to one day put these the correct way round
try, lpf  = inverse.lpf;    catch, lpf  = 0;        end
try, sdv  = inverse.sdv;    catch, sdv  = 4;        end
try, Han  = inverse.Han;    catch, Han  = 1;        end
try, woi  = inverse.woi;    catch, woi  = [];       end
try, Nm   = inverse.Nm;     catch, Nm   = [];       end
try, Nt   = inverse.Nt;     catch, Nt   = [];       end % fixed number of temporal modes
try, Ip   = inverse.Ip;     catch, Ip   = [];       end
try, QE   = inverse.QE;     catch, QE   = 1;        end %  empty room noise measurement
try, Qe0  = inverse.Qe0;    catch, Qe0  = exp(-5);  end % set noise floor at 1/100th signal power i.e. assume amplitude SNR of 10
try, inverse.A;             catch, inverse.A = [];  end % orthogonal channel modes



% defaults
%--------------------------------------------------------------------------
%type = inverse.type;    % Type of inversion scheme


% get specified modalities to invert (default to all)

%--------------------------------------------------------------------------
modalities = D.inv{val}.forward.modality;       % MEG in this case

Nmax  = 16;         % max number of temporal modes

% check lead fields and get number of dipoles (Nd) and channels (Nc)
%==========================================================================

fprintf('Checking leadfields')
[L,D] = spm_eeg_lgainmat(D);    % Generate/load lead field
Nd=size(L,2);

% Check gain or lead-field matrices
%--------------------------------------------------------------------------

if size(modalities,1)>1
    error('not defined for multiple modalities');
end
Ic  = setdiff(D.indchantype(modalities), badchannels(D));
Nd   = size(L,2);      % Number of dipoles




%==========================================================================
% Spatial projectors (adjusting for different Lead-fields)
%==========================================================================

fprintf('Optimising and aligning spatial modes ...\n')

% eliminate low SNR spatial modes
%--------------------------------------------------------------------------
disp('FIXING FULL RANK A');
inverse.A=eye(length(Ic)); Nm=size(inverse.A,1);

if isempty(inverse.A) % no spatial modes prespecified
    if isempty(Nm), %% number of modes not specifiedd
        [U,ss,vv]    = spm_svd((L*L'),exp(-16));
        A     = U';                 % spatial projector A
        UL    = A*L;
        
    else % number of modes pre-specified
        [U,ss,vv]    = spm_svd((L*L'),0);
        if length(ss)<Nm,
            disp('number available');
            length(ss)
            error('Not this many spatial modes in lead fields');
            
        end;
        
        ss = ss(1:Nm);
        disp('using preselected number spatial modes !');
        A     = U(:,1:Nm)';                 % spatial projector A
        UL    = A*L;
    end
else %% U was specified in input
    disp('Using pre-specified spatial modes');
    if isempty(Nm),
        error('Need to specify number of spatial modes if U is prespecified');
    end;
    %
    A=inverse.A;
    UL=A*L;
end

Nm    = size(UL,1);         % Number of spatial projectors


clear ss vv

% Report
%--------------------------------------------------------------------------
fprintf('Using %d spatial modes',Nm)


%==========================================================================
% Temporal projector
%==========================================================================
AY    = {};                                      % pooled response for MVB
AYYA  = 0;                                       % pooled response for ReML

% Time-window of interest
%--------------------------------------------------------------------------

if isempty(woi)
    w      = 1000*[min(D.time) max(D.time)];
else
    w=woi; %% in milliseconds
end

It     = (w/1000 - D.timeonset)*D.fsample + 1;
It     = max(1,It(1)):min(It(end), length(D.time));
It     = fix(It);
disp(sprintf('Number of samples %d',length(It)));

    

% Peristimulus time
%--------------------------------------------------------------------------
pst    = 1000*D.time;                   % peristimulus time (ms)
pst    = pst(It);                       % windowed time (ms)
dur    = (pst(end) - pst(1))/1000;      % duration (s)
dct    = (It - It(1))/2/dur;            % DCT frequencies (Hz)
Nb     = length(It);                    % number of time bins

% Serial correlations
%--------------------------------------------------------------------------
K      = exp(-(pst - pst(1)).^2/(2*sdv^2)); %% sdv set to 4 by default
K      = toeplitz(K);
qV     = sparse(K*K'); %% Samples* samples covariance matrix- assumes smooth iid

% Confounds and temporal subspace
%--------------------------------------------------------------------------

T      = spm_dctmtx(Nb,Nb);         % use plot(T) here!

j      = find( (dct >= lpf) & (dct <= hpf) ); %% THis is the wrong way round but leave for nowfor compatibility with spm_eeg_invert
T      = T(:,j);                    % Apply the filter to discrete cosines
dct    = dct(j);                    % Frequencies accepted

% Hanning window
%--------------------------------------------------------------------------

if Han
    W  = sparse(1:Nb,1:Nb,spm_hanning(Nb)); %% use hanning unless specified
else
    W=1;
end



% get trials or conditions
%--------------------------------------------------------------------------
try
    trial = D.inv{D.val}.inverse.trials;
catch
    trial = D.condlist;
end
Ntrialtypes=length(trial);
% get temporal covariance (Y'*Y) to find temporal modes
%==========================================================================

YY    = 0;

N=0;
badtrialind=D.badtrials;
for j = 1:Ntrialtypes,                          % pool over conditions
    c     = D.indtrial(trial{j});     % and trials
    c=setxor(c,badtrialind); %% ignore bad trials
    Nk    = length(c);
    for k = 1:Nk
        Y     = A*D(Ic,It,c(k));
        
        YY    = YY + Y'*Y;
        N     = N + 1;
    end
end
YY=YY./N;


% Apply any Hanning and filtering
%--------------------------------------------------------------------------
YY         = W'*YY*W;     % Hanning
YTY        = T'*YY*T;     % Filter


%==========================================================================

if isempty(Nt), %% automatically assign appropriate number of temporal modes    
    [U,E]  = spm_svd(YTY,exp(-8));          % get temporal modes
    if isempty(U), %% fallback
        warning('nothing found using spm svd, using svd');
        [U,E]  = svd(YTY);          % get temporal modes
    end
    E      = diag(E)/trace(YTY);            % normalise variance
    Nr     = min(length(E),Nmax);           % number of temporal modes
    Nr=max(Nr,1); %% use at least one mode
else %% use predefined number of modes
    [U,E]  = svd(YTY);          % get temporal modes
    E      = diag(E)/trace(YTY);            % normalise variance
    disp('Fixed number of temporal modes');
    Nr=Nt;
end

V      = U(:,1:Nr);                     % temporal modes
VE     = sum(E(1:Nr));                  % variance explained

fprintf('Using %i temporal modes, ',Nr)
fprintf('accounting for %0.2f percent average variance\n',full(100*VE))

% projection and whitening
%--------------------------------------------------------------------------
S      = T*V;                           % temporal projector
Vq     = S*pinv(S'*qV*S)*S';            % temporal precision

%disp('FIXIN TEMP PROJECTOR TO FULL RANK');
%S=eye(size(YY,1));

% get spatial covariance (Y*Y') for Gaussian process model
%==========================================================================
% loop over Ntrialtypes trial types 
%--------------------------------------------------------------------------
UYYU    = 0;
AYYA    = 0;
Nn      = 0;                                 % number of samples
AY      = {};
Ntrials = 0;

for j = 1:Ntrialtypes
    
    UY{j} = sparse(0);
    c     = D.indtrial(trial{j});
    c     = setxor(c,badtrialind);           % ignore bad trials
    Nk    = length(c);
       
    % loop over epochs
    %----------------------------------------------------------------------
    
    for k = 1:Nk
        
        % stack (scaled aligned data) over modalities
        %------------------------------------------------------------------
        
        Y       = D(Ic,It,c(k))*S;           % in temporal subspace
        Y       = A*Y;                       % in spatial subspace
  
        
        % accumulate first & second-order responses
        %------------------------------------------------------------------
        Nn      = Nn + Nr;                   % number of samples
        
        YY      = Y*Y';                      % and covariance
        Ntrials = Ntrials+1;
        
        % accumulate statistics (subject-specific)
        %------------------------------------------------------------------
        UY{j}   = UY{j} + Y;                 % condition-specific ERP
        UYYU    = UYYU + YY;                 % subject-specific covariance
        
        % and pool for optimisation of spatial priors over subjects
        %------------------------------------------------------------------
        AY{end + 1} = Y;                     % pooled response for MVB
        AYYA        = AYYA    + YY;          % pooled response for ReML
        
    end
end

AY = spm_cat(AY);     % goes to MVB/GS algorithm

ID = spm_data_id(AY); % get a unique ID for these filtered data



%======================================================================

inverse.AY     = AY; %% concatenated data in reduced spatial modes 
inverse.AYYA   = AYYA;%% sensor level covariance of all trials, conditions and samples (containing Nn temporal modes in total)
%inverse.type  = type;                  % inverse model
%inverse.Y     = Y;                     % ERP data (reduced)
inverse.UY     = UY;                    % ERP data (reduced)
inverse.L      = UL;                    % Lead-field (reduced)
inverse.Nn=Nn; %% number of independent samples (temporal modes* trials* conditions)

inverse.qV     = Vq;                    % temporal correlations
inverse.T      = S;                     % temporal projector
inverse.U      = {A};                   % spatial projector
%inverse.Is    = Is;                    % Indices of active dipoles

inverse.Nd     = Nd;                    % number of dipoles
inverse.pst    = pst;                   % peristimulus time
inverse.dct    = dct;                   % frequency range

inverse.ID     = ID;                    % data ID
inverse.woi    = w;                     % time-window inverted

inverse.modality = modalities;          % modalities inverted


% save in struct
%----------------------------------------------------------------------
D.inv{val}.inverse = inverse;
D.inv{val}.method  = 'Imaging';

disp('Done invert prepro');
