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
%__________________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_firstlevel.m 6146 2014-09-02 11:27:43Z vladimir $


SVNrev = '$Rev: 6146 $';

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

D        = spm_eeg_load(char(S.sess.D));

statdir = D.path;

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

job{1}.spm.stats.fmri_design = [];
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);

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;
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;

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')
            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);
            
            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 isequal(S.timing.units, 'secs')
                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
                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
        
        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

eventchanind = D.indchannel('event');
if ~isempty(eventchanind) && nc
    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');
    
    [dummy, job] = spm_jobman('harvest', job);
    spm_run_fmri_spec(job{1}.spm.stats.fmri_design);
    
    load(fullfile(statdir, 'SPM.mat'));
    
    X = SPM.xX.X*SPM.xBF.dt;
    X = X(:, 1:(end-1));
    
    aX = any(X', 1);
    caX = cumsum(aX);
    
    if isTF
        Y = squeeze(mean(D(eventchanind, :, :), 2));
    else
        Y = D(eventchanind, :);
    end
    
    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')
        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

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


U = struct([]);

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);

if ~isequal(S.sess.multi_conv_reg, {''})
    temp = load(S.sess.multi_conv_reg{1});
    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

job{1}.spm.stats.fmri_design.timing.fmri_t        = S.timing.utime;
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));

[dummy, job] = spm_jobman('harvest', job);
spm_run_fmri_spec(job{1}.spm.stats.fmri_design);

clear job

load(fullfile(statdir, 'SPM.mat'));

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

nb = size(SPM.xBF.bf, 1);
nr = size(SPM.xBF.bf, 2);

T = round(1/(SPM.xBF.dt*D.fsample));

if ~isempty(U)
    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
        u = U(i).u;
        
        u = u(:);
        
        if length(u)~= D.nsamples
            error('Convolution regressor should be the same length as the data.');
        end
           
        u = [u(:); zeros(pad, 1)];
        
        if T~=1
            u = interp1(1:length(u), u, 1:1/T:length(u));
        end
        
        U(i).u = u(:);
        
        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, :);
    
    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')                           %-#
save(fullfile(statdir, 'SPM.mat'), 'SPM', spm_get_defaults('mat.format'));
fprintf('%30s\n','...SPM.mat saved')                                    %-#

if S.sess.hpf
    K = [];
    K.RT     = dt;
    K.row    = SPM.Sess.row;
    K.HParam = S.sess.hpf;
else
    K = 1;
end

xX    = spm_sp('Set', spm_filter(K, X)); % filter design
if ~isfield(xX,'pX')
    xX.pX = spm_sp('x-',xX);
end

if ~isempty(SPM.Sess.U)            
    label = {};
    for i = 1:numel(SPM.Sess.U)
        label = [label SPM.Sess.U(i).name];
    end
    ne = numel(label);
    
    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');
    
    
    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
        
        W = ones(nt, 1);
        W(D.badsamples(channels(c), ':', 1)) = exp(-256);
        W = spdiags(W, 0, nt, nt);
        
        axX    = spm_sp('Set', spm_filter(K, W*X));
        if ~isfield(axX,'pX')
            axX.pX = spm_sp('x-',axX);
        end
        
        Y = reshape(D(channels(c), :, :, :), nf, nt);
        
        Y = spm_filter(K, W*Y');
        
        B = axX.pX*Y;
        
        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

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
    
    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');
    
    
    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
        
        W = ones(nt, 1);
        W(D.badsamples(channels(c), ':', 1)) = exp(-256);
        W = spdiags(W, 0, nt, nt);
        
        axX    = spm_sp('Set', spm_filter(K, W*X));
        if ~isfield(axX,'pX')
            axX.pX = spm_sp('x-',axX);
        end
        
        Y = reshape(D(channels(c), :, :, :), nf, nt);
        
        Y = spm_filter(K, W*Y');
        
        B = axX.pX*Y;
        
        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