% this script does the preprocessing and computes event-related fields

definetrial = true;
readdata    = true;
dotimelock  = true;
doplot      = true;
writeflag   = true;

if ~exist('subj', 'var')
  error('specify a subject by calling subj = datainfo_subject(<number>)');
end

%% create the trial definition for each run

if definetrial

  subj.trl = cell(6,1);
  for run_nr = 1:6

    cfg = [];
    cfg.dataset = subj.megfile{run_nr};
    cfg.readbids = 'no'; % it looks like BIDS, but is not complete

    % this is what we can see in the BIDS events.tsv file
    Famous      = [5 6 7];
    Unfamiliar  = [13 14 15];
    Scrambled   = [17 18 19];

    cfg.trialfun = 'ft_trialfun_general';
    cfg.trialdef.eventtype = 'STI101';
    cfg.trialdef.eventvalue = [Famous Unfamiliar Scrambled];
    cfg.trialdef.prestim = 0.5;
    cfg.trialdef.poststim = 1.2;

    cfg = ft_definetrial(cfg);

    % remember the trial definition for each run
    subj.trl{run_nr} = cfg.trl;

  end % for each run

end

%% read the data

if readdata

  rundata = cell(1,6);
  for run_nr = 1:6

    cfg         = [];
    cfg.dataset = subj.megfile{run_nr};
    cfg.trl     = subj.trl{run_nr};

    % MEG specific settings
    cfg.channel = 'MEG';
    cfg.demean  = 'yes';
    cfg.coilaccuracy = 0;
    cfg.bpfilter = 'yes';
    cfg.bpfilttype = 'firws';
    cfg.bpfreq  = [1 40];
    cfg.padding = 3;
    data_meg    = ft_preprocessing(cfg);

    % EEG specific settings
    cfg.channel    = {'EEG' '-EEG061' '-EEG062' '-EEG063' '-EEG064'}; % exclude EOG/ECG/etc hard coded assumed to be this list
    cfg.demean     = 'yes';
    cfg.reref      = 'yes';
    cfg.refchannel = 'all'; % average reference
    data_eeg       = ft_preprocessing(cfg);

    % this is what we can see in the BIDS channels.tsv file
    % EEG061	HEOG	V	horizontaleog
    % EEG062	VEOG	V	verticaleog
    % EEG063	ECG  	V	cardiac

    % settings for EOG and ECG channels
    cfg.channel = {'EEG061' 'EEG062' 'EEG063'};
    cfg.demean  = 'yes';
    cfg.reref   = 'no';
    cfg.bpfilter = 'no';
    % use a one-to-one montage to rename the channels, see FT_APPLY_MONTAGE
    cfg.montage.labelold = {'EEG061' 'EEG062' 'EEG063'};
    cfg.montage.labelnew = {'HEOG' 'VEOG' 'ECG'};
    cfg.montage.tra      = eye(3);
    data_exg    = ft_preprocessing(cfg);

    % settings for all other channels
    cfg.channel = {'all', '-MEG', '-EEG'};
    cfg.montage = []; % remove the previously applied montage
    cfg.demean  = 'no';
    cfg.reref   = 'no';
    cfg.bpfilter = 'no';
    data_other  = ft_preprocessing(cfg);

    cfg            = [];
    cfg.resamplefs = 300;
    data_meg       = ft_resampledata(cfg, data_meg);
    data_eeg       = ft_resampledata(cfg, data_eeg);
    data_exg       = ft_resampledata(cfg, data_exg);
    data_other     = ft_resampledata(cfg, data_other);

    % append the different channel sets into a single structure
    rundata{run_nr} = ft_appenddata([], data_meg, data_eeg, data_exg, data_other);
    clear data_meg data_eeg data_exg data_other
  end % for each run

  data = ft_appenddata([], rundata{:});
  clear rundata;


  if writeflag
    filename = fullfile(subj.outputpath, 'raw2erp', subj.name, sprintf('%s_data', subj.name));
    save(filename, 'data');
  end

end

%% compute the averaged ERF

if dotimelock

  if ~readdata
    % get the precomputed data
    filename = fullfile(subj.outputpath, 'raw2erp', subj.name, sprintf('%s_data', subj.name));
    load(filename, 'data');
  end

  % this is what we can see in the BIDS events.tsv file
  Famous      = [5 6 7];
  Unfamiliar  = [13 14 15];
  Scrambled   = [17 18 19];

  cfg        = [];
  cfg.preproc.demean = 'yes';
  cfg.preproc.baselinewindow = [-0.1 0];

  cfg.trials = ismember(data.trialinfo(:,1), Famous);
  avg_famous = ft_timelockanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Unfamiliar);
  avg_unfamiliar = ft_timelockanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Scrambled);
  avg_scrambled = ft_timelockanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), [Famous Unfamiliar]);
  avg_faces  = ft_timelockanalysis(cfg, data);

  if writeflag
    filename = fullfile(subj.outputpath, 'raw2erp', subj.name, sprintf('%s_timelock', subj.name));
    save(filename, 'avg_famous', 'avg_unfamiliar', 'avg_scrambled', 'avg_faces');
  end

end

%% make figures

if doplot

  if ~dotimelock
    filename = fullfile(subj.outputpath, 'raw2erp', subj.name, sprintf('%s_timelock', subj.name));
    load(filename, 'avg_famous', 'avg_unfamiliar', 'avg_scrambled', 'avg_faces');
  end

  % visualise the magnetometer data
  cfg        = [];
  cfg.layout = 'neuromag306mag_helmet.mat';
  figure; ft_multiplotER(cfg, avg_famous, avg_unfamiliar, avg_scrambled);

  % combine planar gradients and visualise the gradiometer data
  cfg              = [];
  avg_faces_c      = ft_combineplanar(cfg, avg_faces);
  avg_famous_c     = ft_combineplanar(cfg, avg_famous);
  avg_unfamiliar_c = ft_combineplanar(cfg, avg_unfamiliar);
  avg_scrambled_c  = ft_combineplanar(cfg, avg_scrambled);

  cfg        = [];
  cfg.layout = 'neuromag306cmb_helmet.mat';
  figure; ft_multiplotER(cfg, avg_famous_c, avg_unfamiliar_c, avg_scrambled_c);

  % create an EEG channel layout on-the-fly and visualise the eeg data
  cfg      = [];
  cfg.elec = avg_faces.elec;
  layout_eeg = ft_prepare_layout(cfg);

  cfg        = [];
  cfg.layout = layout_eeg;
  figure; ft_multiplotER(cfg, avg_famous, avg_unfamiliar, avg_scrambled);

  % alternatively, it is possible to show the different channel types in a
  % single figure, this allows for interacting simultaneously with the
  % different representation. For the colormap to work, the data needs to
  % be scaled per type of data.

  cfg        = [];
  cfg.layout = 'neuromag306mag_helmet.mat'; % magnetometers
  layout_mag = ft_prepare_layout(cfg);

  cfg        = [];
  cfg.layout = 'neuromag306cmb_helmet.mat'; % combined planar gradiometers
  layout_cmb = ft_prepare_layout(cfg);

  % in order for this to work, the positions should be in the same order of magnitude
  shiftval = min(layout_eeg.pos(1:70,:),[],1);
  layout_eeg.pos = layout_eeg.pos - repmat(shiftval, numel(layout_eeg.label), 1);
  layout_eeg.mask{1} = layout_eeg.mask{1} - repmat(shiftval, size(layout_eeg.mask{1},1), 1);
  for k = 1:numel(layout_eeg.outline)
    layout_eeg.outline{k} = layout_eeg.outline{k} - repmat(shiftval, size(layout_eeg.outline{k},1), 1);
  end

  scaleval = max(layout_eeg.pos(1:70,:),[],1)./500;
  layout_eeg.pos = layout_eeg.pos ./ repmat(scaleval, numel(layout_eeg.label), 1);
  layout_eeg.mask{1} = layout_eeg.mask{1} ./ repmat(scaleval, size(layout_eeg.mask{1},1), 1);
  for k = 1:numel(layout_eeg.outline)
    layout_eeg.outline{k} = layout_eeg.outline{k} ./ repmat(scaleval, size(layout_eeg.outline{k},1), 1);
  end

  layout_eeg.width(:)  = 64;
  layout_eeg.height(:) = 48;

  cfg = [];
  cfg.distance = 180;
  layout = ft_appendlayout(cfg, ft_appendlayout([], layout_mag, layout_cmb), layout_eeg);

  cfg = [];
  cfg.layout    = layout;
  cfg.gridscale = 150;
  cfg.magscale  = 0.25e14;
  cfg.gradscale = 1e12;
  cfg.eegscale  = 1e6;
  figure; ft_multiplotER(cfg, avg_famous_c, avg_unfamiliar_c, avg_scrambled_c);

end


%%

avg          = avg_famous.avg(367,:);
avg_plus_sd  = avg + sqrt(avg_famous.var(367,:));
avg_minus_sd = avg - sqrt(avg_famous.var(367,:));

figure
plot(avg_famous.time, avg, 'b')
hold on
plot(avg_famous.time, avg_plus_sd, 'b:')
plot(avg_famous.time, avg_minus_sd, 'b:')

% there are 295 trials in this condition
avg_plus_sem  = avg + sqrt(avg_famous.var(367,:)) ./ sqrt(avg_famous.dof(367,:));
avg_minus_sem = avg - sqrt(avg_famous.var(367,:)) ./ sqrt(avg_famous.dof(367,:));

figure
plot(avg_famous.time, avg, 'b')
hold on
plot(avg_famous.time, avg_plus_sem, 'b:')
plot(avg_famous.time, avg_minus_sem, 'b:')

%%

figure
ft_plot_vector(avg_famous.time, [avg_plus_sem; avg_minus_sem], 'highlightstyle', 'difference')
ft_plot_vector(avg_famous.time, avg, 'color', 'k', 'linewidth', 2)
axis tight
grid on

%%

significant = abs(avg) > 2*sqrt(avg_famous.var(367,:))./sqrt(avg_famous.dof(367,:));

figure
ft_plot_vector(avg_famous.time, avg, 'highlight', significant, 'highlightstyle', 'thickness', 'color', 'b');
axis tight
grid on

%%

cfg = [];
cfg.keeptrials = 'yes';
avg_all = ft_timelockanalysis(cfg, data);

figure
subplot(2,1,1); imagesc(squeeze(avg_all.trial(:,367,:)))
subplot(2,1,2); plot(avg_all.time, mean(squeeze(avg_all.trial(:,367,:)),1)); axis tight


