function D = spm_eeg_firstlevel(S)
% Convolution modelling for M/EEG data
% FORMAT D = spm_eeg_firstlevel(S)
%
%   S - job structure generated by spm_cfg_eeg_firstlevel
%__________________________________________________________________________
% Reference:
% Litvak V, Jha A, Flandin G, Friston K. Convolution models for induced
% electromagnetic responses. Neuroimage. 2013, 64:388-98
%__________________________________________________________________________

% Vladimir Litvak
% Copyright (C) 2013-2022 Wellcome Centre for Human Neuroimaging


%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename);
spm('FigName','M/EEG first level');

% Load the M/EEG .mat file 
% ------------------------
D        = spm_eeg_load(char(S.sess.D));

statdir = D.path;

try
    delete(fullfile(statdir, 'SPM.mat'));
end

% Configure the job for using with conventional fMRI GLM design
% =============================================================

% Create a new field for fMRI analysis
job{1}.spm.stats.fmri_design = [];

% Set basic settings
job{1}.spm.stats.fmri_design.dir            = {statdir};
job{1}.spm.stats.fmri_design.timing.units   = 'secs';
job{1}.spm.stats.fmri_design.timing.fmri_t0 = 1;
job{1}.spm.stats.fmri_design.volt           = S.volt;
job{1}.spm.stats.fmri_design.cvi            = 'none';

isTF     = strncmpi(D.transformtype,'TF',2); 

% Set recording settings 
time     = D.time;
dt       = mean(diff(time));
nt       = D.nsamples;
channels = D.selectchannels(S.channels);
nchan    = numel(channels);
nf       = D.nfrequencies;

job{1}.spm.stats.fmri_design.timing.RT  = dt; % Time resolution
job{1}.spm.stats.fmri_design.sess       = rmfield(S.sess, {'D', 'cond', 'convregress', 'multi_conv_reg', 'savereg'});
job{1}.spm.stats.fmri_design.sess.hpf   = Inf;
job{1}.spm.stats.fmri_design.sess.nscan = nt; % Number of rows

% Process the events for each condition
% -------------------------------------
% details: 
%    .sess
%       .cond
%          .name
%          .define:  How to define events (field 'manual'/'event', see cfg)
%          .tmod:    Time modulation - order of a polynomial basis in time
%          .pmod:    Parameter modulation (polynomial basis) - see cfg
%          .orth:    Whether to orthogonalise modulations
% NB: fields name, tmod, pmod, and orth are not processed, simply copied 

nc = numel(S.sess.cond);

if nc == 0
    job{1}.spm.stats.fmri_design.sess.cond = rmfield(S.sess.cond, 'define');
else
    for j = 1:nc
        conditionlabels = {};
        trigsample      = {};
        duration        = {};
        durationsample  = {};
 
        if isfield(S.sess.cond(j).define, 'event')
            % If events are defined from dataset

            % Configure structure to get event info from dataset
            S1              = [];
            S1.D            = D;
            S1.reviewtrials = 0;
            S1.save         = 0;
            S1.trialdef     = S.sess.cond(j).define.event;
            S1.timewin      = [0 0];
            [S1.trialdef(:).conditionlabel] = deal(S.sess.cond(j).name);
            
            % Get event info from dataset
            trl                       = spm_eeg_definetrial(S1);
            trigsample{j}             = trl(:, 1);
            trig{j}                   = time(trigsample{j});
            duration{j}               = 0*trig{j};
            durationsample{j}         = 0*trig{j};
        elseif isfield(S.sess.cond(j).define, 'manual')
            % If events are defined manually 
            
            if isequal(S.timing.units, 'secs')
                % Populate sample index info from physical time 
                trig{j}         = S.sess.cond(j).define.manual.onset;
                trigsample{j}   = D.indsample(trig{j});
                duration{j}     = S.sess.cond(j).define.manual.duration;
                if isempty(duration)
                    duration{j} = 0*trig{j};
                end
                durationsample{j} = round(duration{j}*D.fsample);
            else
                % Populate physical time fields from sample index info
                trigsample{j}     = S.sess.cond(j).define.manual.onset;
                trig{j}           = time(trigsample{j});
                durationsample{j} = S.sess.cond(j).define.manual.duration;
                if isempty(durationsample)
                    durationsample{j} = 0*trig{j};
                end
                duration{j}       = durationsample{j}/D.fsample;
            end
            conditionlabels{j}    = repmat({S.sess.cond(j).name}, length(trig{j}), 1);
        else
            error('Unsupported option');
        end
        
        % Setup condition in fMRI design
        job{1}.spm.stats.fmri_design.sess.cond(j).name     = S.sess.cond(j).name;
        job{1}.spm.stats.fmri_design.sess.cond(j).onset    = trig{j} - D.timeonset;
        job{1}.spm.stats.fmri_design.sess.cond(j).duration = duration{j};
        job{1}.spm.stats.fmri_design.sess.cond(j).tmod     = S.sess.cond(j).tmod;
        job{1}.spm.stats.fmri_design.sess.cond(j).pmod     = S.sess.cond(j).pmod;
        job{1}.spm.stats.fmri_design.sess.cond(j).orth     = S.sess.cond(j).orth;
    end
end

% Try to get the index of the event channel in the dataset
eventchanind = D.indchannel('event');
if ~isempty(eventchanind) && nc
    % Use the event channel to compute the timing shift (dataset vs constructed)
    % --------------------------------------------------------------------------
    sess =  job{1}.spm.stats.fmri_design.sess;
    job{1}.spm.stats.fmri_design.sess.cond = ...
        rmfield(sess.cond, {'tmod', 'pmod', 'orth'});
    job{1}.spm.stats.fmri_design.timing.fmri_t    = 1;
    job{1}.spm.stats.fmri_design.bases.fir.length = dt;
    job{1}.spm.stats.fmri_design.bases.fir.order  = 1;
    job{1}.spm.stats.fmri_design.sess.regress     = [];
    job{1}.spm.stats.fmri_design.sess =...
        rmfield(job{1}.spm.stats.fmri_design.sess, 'regress');
    
    % Configure default fields for the job
    [dummy, job] = spm_jobman('harvest', job);

    % Build and load the design matrix (populates SPM.mat) 
    spm_run_fmri_spec(job{1}.spm.stats.fmri_design);
    load(fullfile(statdir, 'SPM.mat'));
    
    % NB: for info about fields in SPM.mat: 'help spm_fmri_spm_ui' and
    % 'help spm_fmri_design'
    % xBF:  Basis function structure
    % xX.X: Design matrix
    X = SPM.xX.X*SPM.xBF.dt;

    % Remove the intercept 
    X = X(:, 1:(end-1));
    
    % check whether rows are populated
    aX = any(X', 1);
    caX = cumsum(aX);
    
    % average over trials if time-frequency data
    if isTF
        Y = squeeze(mean(D(eventchanind, :, :), 2));
    else
        Y = D(eventchanind, :);
    end
    
    % Use lagged cross-correlation to find a shift index/time
    Y = Y-min(Y); 
    Y = Y/max(Y); 
    
    [c, lags] = xcorr(aX, Y, 500, 'coeff'); 
    
    [dummy, ind] = max(c); 
    shiftcorr = lags(ind); 
    
    if ~spm('CmdLine')
        % Plot timing realignment
        % -----------------------

        spm_figure('GetWin','Timing check');
        clf
        subplot(2, 3, 1:3);
        plot(lags, c);
        xlim(100*[-1 1]);
        
        ntoplot = 5;
        i1 = find(caX == 1);
        i2 = find(caX == ntoplot);
        
        i3= find((caX-caX(floor(nt/2))) == 1);
        i4= find((caX-caX(floor(nt/2))) == ntoplot);
        
        i5= find((caX-caX(end)) == -ntoplot);
        i6= find((caX-caX(end)) == 0);
        
        subplot(2, 3, 4);
        plot(time(i1(1):i2(1)), Y(i1(1):i2(1)), 'k');
        hold on
        plot(time(i1(1):i2(1)), X(i1(1):i2(1), :));
        
        subplot(2, 3, 5);
        plot(time(i3(1):i4(1)), Y(i3(1):i4(1)), 'k');
        hold on
        plot(time(i3(1):i4(1)), X(i3(1):i4(1), :));
        
        subplot(2, 3, 6);
        plot(time(i5(1):i6(1)), Y(i5(1):i6(1)), 'k');
        hold on
        plot(time(i5(1):i6(1)), X(i5(1):i6(1), :));
    end
    
    delete(fullfile(statdir, 'SPM.mat'));
    job{1}.spm.stats.fmri_design.bases  = [];
    job{1}.spm.stats.fmri_design.sess   = sess;
else
    shiftcorr = 0;
end

% Shift to align onsets with data
% -------------------------------
for c = 1:nc
    job{1}.spm.stats.fmri_design.sess.cond(c).onset = ...
        job{1}.spm.stats.fmri_design.sess.cond(c).onset ... 
        + 1e-3*S.timing.timewin(1) - shiftcorr*dt;
end


% Build the convolution regressors
% ================================
U = struct([]);

% Handle regressors provided by the user
for i = 1:numel(S.sess.convregress) 
    U(i).name = S.sess.convregress(i).name;
    U(i).u    = S.sess.convregress(i).val;
end

ncr = numel(U);

% Check whether convolution regressors where provided through a .mat file
% (see spm_cfg_eeg_firstlevel.m for details on the file structure) 
if ~isequal(S.sess.multi_conv_reg, {''})
    temp = load(S.sess.multi_conv_reg{1});

    % Add each regressor and name it ('R1', 'R2', ...) if necessary 
    for i = 1:size(temp.R, 2)
        U(ncr+i).u = temp.R(:, i); 
        if isfield(temp, 'names')
            U(ncr+i).name = temp.names{i};
        else
            U(ncr+i).name = ['R' num2str(i)];
        end
    end
end

ncr = numel(U);

% Create dummy conditions to replace with continuous regressors later
for i = 1:ncr
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).name     = U(i).name;
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).onset    = D.time(round(D.nsamples/2));
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).duration = 0;
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).tmod     = 0;
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).pmod     = struct([]);
    job{1}.spm.stats.fmri_design.sess.cond(nc+i).orth     = 0;
end

% Set microtime resolution
job{1}.spm.stats.fmri_design.timing.fmri_t        = S.timing.utime;

% Set bases function (field 'fourier', 'fourier_han', or 'fir')
bname                                             = char(fieldnames(S.bases));
job{1}.spm.stats.fmri_design.bases                = S.bases;
job{1}.spm.stats.fmri_design.bases.(bname).length = 1e-3*diff(sort(S.timing.timewin));

% Construct the fMRI GLM design
% =============================

% Configure default fields for the job
[dummy, job] = spm_jobman('harvest', job);

% Build the design (populates SPM.mat)
spm_run_fmri_spec(job{1}.spm.stats.fmri_design);

% Make some space
clear job 

% Load the design 
load(fullfile(statdir, 'SPM.mat'));

% Scale stick functions with sampling interval 
X = SPM.xX.X*SPM.xBF.dt; 

% Unpack the shape of the basis set
% nb: number of time points      (xBF.length / xBF.dt + 1) 
% nr: number of basis functions  (xBF.order * 2 + 1)
nb = size(SPM.xBF.bf, 1);
nr = size(SPM.xBF.bf, 2);

% Compute the up- or downsampling factor 
% (basis functions' dt vs dataset sampling rate) 
T = round(1/(SPM.xBF.dt*D.fsample)); 

if ~isempty(U)
    % Compute the padding necessary for time windows
    if S.timing.timewin(1)<0
        pad = round(abs(1e-3*S.timing.timewin(1)*D.fsample));
    else
        pad = 0;
    end
    
    for i = 1:ncr
        % Unpack convolution regressor
        u = U(i).u;
        u = u(:);
        
        % Check dimensions
        if length(u)~= D.nsamples
            error('Convolution regressor should be the same length as the data.');
        end
           
        % Pad the data at the end 
        u = [u(:); zeros(pad, 1)];
        
        % Resample if necessary (linear interpolation)
        if T~=1
            u = interp1(1:length(u), u, 1:1/T:length(u));
        end
        
        % Repack
        U(i).u = u(:);
        
        % Specify new dt
        U(i).dt = SPM.xBF.dt;
        
        if ~iscell(U(i).name)
            U(i).name = {U(i).name};
        end
    end
    
    
    %-Convolve stimulus functions with basis functions
    %-------------------------------------------------
    [Xu] = spm_Volterra(U, SPM.xBF.bf);
    
    %-Resample regressors 
    %--------------------
    Xu = Xu((0:(nt + pad - 1))*T+1, :);       
    
    Xu = Xu((1+pad):end, :);
    
    % nc  = number of conditions
    % ncr = number of convolution regressors
    % nr  = size of basis set
    X(:, (nc*nr+1):((nc+ncr)*nr)) = Xu*SPM.xBF.dt; 
end

SPM.xX.X = X./SPM.xBF.dt; 

%-Save SPM.mat
%-------------
fprintf('%-40s: ','Saving SPM configuration')                           %-#
fmt = spm_get_defaults('mat.format');
s = whos('SPM');
if s.bytes > 2147483647, fmt = '-v7.3'; end
save(fullfile(statdir, 'SPM.mat'), 'SPM', fmt);
fprintf('%30s\n','...SPM.mat saved')                                    %-#

% Configure the high-pass filter if requested
if S.sess.hpf
    K = [];
    K.RT     = dt;
    K.row    = SPM.Sess.row;
    K.HParam = S.sess.hpf;
else
    K = 1;
end

% Compute SVD of filtered design matrix
xX    = spm_sp('Set', spm_filter(K, X)); 

% Compute the pseudo-inverse in reduced space
if ~isfield(xX,'pX')
    xX.pX = spm_sp('x-',xX); 
end

% Process and save design inputs, compute betas
% ---------------------------------------------
if ~isempty(SPM.Sess.U)            
    label = {};
    for i = 1:numel(SPM.Sess.U)
        label = [label SPM.Sess.U(i).name];
    end
    ne = numel(label);
    
    % Prepare the structure for the output dataset 
    % --------------------------------------------
    if isTF
        Dout = clone(D, spm_file(D.fname, 'prefix', S.prefix), [nchan nf nb ne]);
    else
        Dout = clone(D, spm_file(D.fname, 'prefix', S.prefix), [nchan nb ne]);
    end
    
    Dout = fsample(Dout, 1/SPM.xBF.dt);
    Dout = timeonset(Dout, 1e-3*S.timing.timewin(1));
    Dout = chanlabels(Dout, ':', D.chanlabels(channels));
    Dout = conditions(Dout, ':', label);
    Dout = type(Dout, 'evoked');

    % Process each channel
    % --------------------
    spm_progress_bar('Init', nchan, 'channels done');
    if nchan > 100, Ibar = floor(linspace(1, nchan, 100));
    else Ibar = 1:nchan; end

    for c = 1:nchan

        % Set sample weight for channel 
        W = ones(nt, 1);
        W(D.badsamples(channels(c), ':', 1)) = exp(-256);
        
        % Compute svd with sample weight for that channel
        axX    = spm_sp('Set', spm_filter(K, W.*X));

        % Compute pinv(X'X)*X'
        if ~isfield(axX,'pX')
            axX.pX = spm_sp('x-',axX);
        end
        
        % Reshape and filter M/EEG for that channel
        % -----------------------------------------
        Y = reshape(D(channels(c), :, :, :), nf, nt);
        Y = spm_filter(K, W.*Y');
        
        % Compute beta coefficients (B =(X'X)^{-1} X'Y)
        % ---------------------------------------------
        B = axX.pX*Y;
        
        % Recompose the beta coefficients for the labels 
        % ----------------------------------------------
        % For now, betas are weighting each component of the basis set
        % individualy. Project on the basis to construct the deconvolved
        % response. 
        % i.e. xY(t, i) = sum_k basis(t, k) * beta(k, i) 
        for i = 1:ne
            xY = SPM.xBF.bf*B((i-1)*nr+(1:nr), :);
            if isTF
                Dout(c, :, :, i) = shiftdim(xY', -1); 
            else
                Dout(c, :, i) = xY';
            end
        end
        
        if ismember(c, Ibar), spm_progress_bar('Set', c); end
    end
    
    spm_progress_bar('Clear');
    
    %-Save
    %----------
    save(Dout);
else
    Dout = [];
    nr   = 0;
    ne   = 0;
end

% Process and save regressors if requested
% ----------------------------------------
if ~isempty(SPM.Sess.C.C) && S.sess.savereg
    rlabel = {};
    for i = 1:numel(SPM.Sess.C)
        rlabel = [rlabel SPM.Sess.C(i).name];
    end
    
    ng = numel(rlabel);
    
    if ~isempty(Dout)
        preprefix = 'R';
    else
        preprefix = '';
    end
    
    % Prepare the structure for the output dataset
    % --------------------------------------------
    if isTF
        Dr = clone(D, spm_file(D.fname, 'prefix', [S.prefix preprefix]), [nchan nf 1 ng]);
    else
        Dr = clone(D, spm_file(D.fname, 'prefix', [S.prefix preprefix]), [nchan 1 ng]);
    end
    
    Dr = fsample(Dr, 0);
    Dr = timeonset(Dr, 0);
    Dr = chanlabels(Dr, ':', D.chanlabels(channels));
    Dr = conditions(Dr, ':', rlabel);
    Dr = type(Dr, 'evoked');
    
    % Process each channel
    % --------------------
    spm_progress_bar('Init', nchan, 'channels done');
    if nchan > 100, Ibar = floor(linspace(1, nchan, 100));
    else Ibar = 1:nchan; end
    for c = 1:nchan

        % Compute sample weights for the channel
        % --------------------------------------
        W = ones(nt, 1);
        W(D.badsamples(channels(c), ':', 1)) = exp(-256);
        W = spdiags(W, 0, nt, nt);
        
        % Compute svd with sample weight for that channel
        axX    = spm_sp('Set', spm_filter(K, W*X));
        
        % Compute pinv(X'X)*X'
        if ~isfield(axX,'pX')
            axX.pX = spm_sp('x-',axX); 
        end
        
        % Reshape and filter M/EEG for that channel
        % -----------------------------------------
        Y = reshape(D(channels(c), :, :, :), nf, nt);
        Y = spm_filter(K, W*Y');
        
        % Compute beta coefficients (B = pinv(X'X)*X'*Y)
        % ---------------------------------------------
        B = axX.pX*Y;
        
        % Set data structure to beta
        % --------------------------
        for i = 1:ng
            xY = B(ne*nr+i, :);
            if isTF
                Dr(c, :, :, i) = xY;
            else
                Dr(c, :, i) = xY;
            end
        end
        
        if ismember(c, Ibar), spm_progress_bar('Set', c); end
    end
    
    spm_progress_bar('Clear');
    
    %-Save
    %-----
    save(Dr);
else
    Dr = [];
end

if ~isempty(Dout)
    D = Dout;
else
    D = Dr;
end