%% Funtion to read, segment in trials and remove artifacts
% Triggers codes using BITSI in my experiment
% X0 Block number X (e.g. 10)
% X1 Fixation point in block number X (e.g. 11, 21, 31 for 1st, 2nd and 3d block respectively)
% 1X1 1X2 1sr Pedestal (e.g. 111 or 112 center or periphery respectively)
% X2 X3 1st Cue presentation (e.g. 12 or 13, Cue and non cue respectively; important for oddball task)
% 1X3 1X4 Target (e.g. 113 or 114 signal or noise respectively) In noise events actually nothing happens
% X4 X5 2nd Cue presentation (e.g. 14 or 15,  Cue and non cue respectively; important for oddball task)
% X6 Response phase starts (e.g. 16,)
% Not recorded anymore -> X7 X8 Sound presentation (e.g. 17,18, audivisual and visual trials respectively)

clear all
close all

%%first read in and segment the data


% Define globals
global G;
G = [];


%% General paths
exp_main_path = '/project/3018029.06';

SUBJS = {'034'};

curSubj = SUBJS{1}

addpath(genpath( [exp_main_path , filesep 'MEG_analysis' filesep 'AVE_MEG' filesep 'scripts'])); % general processing scripts in this exp.

G.subjPath = [exp_main_path filesep 'MEG' filesep curSubj filesep]; % main subject folder
G.dataPath = [exp_main_path filesep 'raw', filesep 'sub-' curSubj  filesep]; % subjecf MEG data
G.outputPath = [G.subjPath ,'func', filesep]; % subject saving artifact and preprocessing matrixes
G.cfgbehavPath = [exp_main_path filesep 'Behavioral' filesep curSubj filesep] ; % behavioral data of participant


subjscriptsPath = [G.subjPath ,'analysis', filesep]; %  participant specific scripts

% add to path scripts for pre-processing this participant (mainly trial
% function, where the trials manually rejected are selected manually)
addpath(genpath(subjscriptsPath));

inputPath = [G.dataPath  'ses-meg01' filesep ];

% Find datafiles
dirs = dir([inputPath, '*.ds']);
if (length(dirs) ~= 1)
    warning(['No or more than one possible MEG-file found for: %s' curSubj]);
end


%% Segment the MEG data in trials according to the trigger channel (UPPT001)

%% Trial sequece is as follows: 
% fixation (duration random time 0.75/1) -> pedestal (duration random time 1.5/1) -> second pedestal (fix time 1.25) -> 
% responsetime (fix duration 1.5) -> feedback (fix time 0.25) -> next fixation...

cfg = [];
cfg.trialdef.prestim = 2.0;
cfg.trialdef.poststim = 1.5;
cfg.trialdef.eventtype = 'UPPT001';  % [17 18 27 28 37 38 47 48 57 58 67 68]; % 
cfg.trialdef.eventvalue =  [113:114,123:124,133:134,143:144,153:154,163:164 ]; % [11,21,31,41,51,61,71]; % signal or noise for each block of 6
cfg.trialfun = 'trialfun_AVE';
cfg.dataset = [inputPath, dirs(1).name];
cfg.PMtrialfun.inputPath = inputPath;
cfg.PMtrialfun.curSubj = curSubj;
cfg.pathBehav = G.cfgbehavPath;

% lets generate the timings for each trial
cfg = ft_definetrial(cfg);

% remove preproc path
rmpath( subjscriptsPath);


%% Reading behavioral data to append in trialinfo
dirs = dir([cfg.pathBehav, '*detect.mat']);
if (length(dirs) ~= 1)
    warning(['Not exactly one MEG-file found for Subj', num2str(curSubj), '!']);
end

load([cfg.pathBehav, dirs.name])

% concatenating blocks
behav_matrix = [];
for i = 1 : size(Experimentalsession,2)
    behav_matrix = [behav_matrix ;   Experimentalsession(i).detect.responses_vector];
end

% free the space
Experimentalsession = [];

% lets take the relevant variables for the analyses:


% behav_matrix
%% session, trial cond, modality, visualfield, contrast, keypress,
% triggeresp, correct, RT

% 1 2 8 9 10 13 16 17 18 19  (+4 after appending to cfg.trl)

behav_coded_matrix = zeros(size(behav_matrix,1),10);
% recoding
behav_coded_matrix(:,1) = cell2mat(behav_matrix(:,1));
behav_coded_matrix(:,2) = cell2mat(behav_matrix(:,2));
behav_coded_matrix(find(strcmp(behav_matrix(:,8) , 'signal')),3) = 1; % 0 is noise 1 is signal
behav_coded_matrix(find(strcmp(behav_matrix(:,9) , 'AV')),4) = 1; % 0 is V, 1 is AV
behav_coded_matrix(find(strcmp(behav_matrix(:,10) , 'surround')),5) = 1; % 0 is center, 1 is surround
behav_coded_matrix(:,6) = cell2mat(behav_matrix(:,13));
behav_coded_matrix(find(strcmp(behav_matrix(:,16) , 'Y')),7) = 1; % 0 is NO, 1 is YEs and, 3 is NoResponse
behav_coded_matrix(find(strcmp(behav_matrix(:,16) , 'NoResponse')),7) = 3;
behav_coded_matrix(:,8) = cell2mat(behav_matrix(:,17));
behav_coded_matrix(:,9) = cell2mat(behav_matrix(:,18));
behav_coded_matrix(:,10) = cell2mat(behav_matrix(:,19));
% Validate trialinfo

% now lets compare behavioral labels with triggers

%% Validate trialinfo (Lets see if triggers in trl and corresponding behavioral matrix condition match).

% first I recode my relevant triggers to the behavioral values
recoded_triggers = 'signal';
recoded_triggers = num2str(cfg.trl(:,4));
recoded_triggers = recoded_triggers(:,3);
recoded_triggers = str2num(recoded_triggers);

nTrials = size(recoded_triggers,1);
recoded_vector = repmat({'noise'},nTrials ,1);
recoded_vector(find(recoded_triggers == 3)) = {'signal'};

% now lets compare behavioral labels with triggers

if isequal(behav_matrix(:,8),  recoded_vector)
    warning('Trialinfo succesfully validated ');
else
    error('Trialinfo invalid for AVE data,');
end

%% Append behav to trialinfo
cfg.trl = [cfg.trl, behav_coded_matrix]; % if you open the matrix some values can look as dissapeared. This is because the scale


%% Removing no response trials
noresp = (cfg.trl(:, 12) == 0); % in triggerresponse code, 0 is no buttong press
cfg.trl(noresp, :) = [];

%data_preproc.trialinfo = cfg.trl;


%% Read the MEG data

% always check whether there are 576 trials minus the no response

cfg.continuous    = 'yes'
cfg.demean        = 'yes';
data_preproc = ft_preprocessing(cfg);

%cfg = [];
%cfg.offset = 1.25 * data_preproc.hdr.Fs;
%data_preproc = ft_redefinetrial(cfg, data_preproc)


%% artifact correction
%% Identify artifacts

time = data_preproc.time{1};
numT = length(data_preproc.time{1});
numN = length(data_preproc.trial);

cfg = [];
cfg.channel = 'MEG';
% cfg.channel = [cfg.channel, {'-MLF12', '-MLF14', '-MLF25', '-MLT31'}];
tmp_data = ft_selectdata(cfg, data_preproc);
tmp_data.trialinfo = [tmp_data.trialinfo, (1:numN)'];



%% Overall artifacts
cfg = [];
cfg.renderer = 'painters';
tmp_data_overall = ft_rejectvisual(cfg, tmp_data);


removed_n_overall = setdiff(1:numN, tmp_data_overall.trialinfo(:, end));
channels_removed_n_overall = setdiff(tmp_data.label, tmp_data_overall.label);

clear tmp_data_overall;

%% Muscle artifacts
cfg = [];
cfg.preproc.hpfilter = 'yes';
cfg.preproc.hpfreq = 100;
cfg.layout = 'CTF275.lay';
tmp_data_muscle = ft_rejectvisual(cfg, tmp_data);

removed_n_muscle = setdiff(1:numN, tmp_data_muscle.trialinfo(:, end));

clear tmp_data_muscle;
% 
% 
%%  Sounds presentation (Check the timings of the stimulus)
cfg = [];
cfg.channel = {'UADC003'}; % this is the sound channel
tmp_data = ft_selectdata(cfg, data_preproc);
tmp_data.trialinfo = (1:numN)' ;
X = reshape(cell2mat(tmp_data.trial), [1, length(tmp_data.time{1}), length(tmp_data.trial)]);
X_sound = squeeze(X(1, :, :))';

imagesc(time,time, X_sound(find( data_preproc.trialinfo(:,5) == 1 ),:));
sound = mean(X_sound(find( data_preproc.trialinfo(:,5) == 1),:));
imagesc(time, time, sound); colorbar;
plot(time,(sound));

% clear tmp_data;

%%  Blinks during stimulus
% I only read from the eyetracker (no EOG data)
stimOn = [-0.1, 0.1];

cfg = [];
cfg.channel = {'UADC007'};
tmp_data = ft_selectdata(cfg, data_preproc);
tmp_data.trialinfo = (1:numN)';

% not sure what this reshape is doing
X = reshape(cell2mat(tmp_data.trial), [1, length(tmp_data.time{1}), length(tmp_data.trial)]);
X_pupDil = squeeze(X(1, :, :))';
%  X_EOGv = squeeze(X(1, :, :))';
%  X_EOGv = (X_EOGv - nanmean(X_EOGv(:))) ./ nanstd(X_EOGv(:));

X_pupDil = (X_pupDil - nanmean(X_pupDil(:))) ./ nanstd(X_pupDil(:));

%cutoff_z_EOGv = 5;
cutoff_z_pupDil = 5;
stimOnBuffer = .1;

% blinkMask = (abs(X_EOGv) > cutoff_z_EOGv) | (abs(X_pupDil) > cutoff_z_pupDil);
% stimOnMask = repmat(inrange(tmp_data.time{1}, stimOn + [-1 1]*stimOnBuffer), [length(tmp_data.trial), 1]);

blinkMask =  (abs(X_pupDil) > cutoff_z_pupDil);
stimOnMask = repmat(inrange(tmp_data.time{1}, stimOn + [-1 1]*stimOnBuffer), [length(tmp_data.trial), 1]);

%{ Visualize blinks
figure;
%  subplot(4, 1, 1); imagesc(1:numN, time, X_EOGv'); colorbar;
subplot(3, 1, 1); imagesc(1:numN, time, X_pupDil'); colorbar;
subplot(3, 1, 2); imagesc(1:numN, time, blinkMask'); colorbar;
subplot(3, 1, 3); imagesc(1:numN, time, stimOnMask'); colorbar;
%}

% removing those trials where the participant blinke during stimulus
% presentation
removed_n_blinks = tmp_data.trialinfo(any(stimOnMask & blinkMask, 2), end);

clear tmp_data;


% Inspect potentially contaminated trials
cfg = [];
cfg.channel = {'MEG'}; % , 'EEG057', 'EEG058'
%cfg.channel = [cfg.channel, {'-MLF12', '-MLF14', '-MLF25', '-MLT31'}];
% cfg.megscale = 1e8;
cfg.artfctdef.overall.artifact = data_preproc.sampleinfo(removed_n_overall, :);
cfg.artfctdef.muscle.artifact = data_preproc.sampleinfo(removed_n_muscle, :);
cfg.artfctdef.blinks.artifact = data_preproc.sampleinfo(removed_n_blinks, :);
cfg.preproc.hpfilter = 'no';
cfg.preproc.hpfreq = 100;
cfg.renderer = 'painters';
cfgart = ft_databrowser(cfg, data_preproc);

% Store artifact file
save([G.outputPath, 'cfgartAVE_', curSubj, '.mat'], 'cfgart');

% Store reject channels
rejchan =  channels_removed_n_overall; % these channels had a lof of noise in this participant

save([G.outputPath,'rejchan_', curSubj, '.mat'], 'rejchan');

%% 1.2. Reject artifacts and keep data in new variable

cfgartFile = [G.outputPath, 'cfgartAVE_', curSubj, '.mat']; load(cfgartFile);

% Remove artifacts
data_preproc_clean = ft_rejectartifact(cfgart, data_preproc);

% Remove bad channels, if available % consider apply this before ICA
if exist([G.outputPath,'rejchan_',  curSubj, '.mat'], 'file')
    load([G.outputPath,'rejchan_',  curSubj, '.mat']);
    % Select good channels
    cfg = [];
    cfg.channel = setdiff( data_preproc_clean.label, rejchan);
    data_preproc_clean  = ft_selectdata(cfg,  data_preproc_clean);
end


%% downsample for ICA
cfg                 = [];
cfg.resamplefs = 150;
cfg.detrend    = 'no';
data_preproc_clean_ds = ft_resampledata(cfg, data_preproc_clean);


%% perform the independent component analysis in downsampled data (i.e., decompose the data)
cfg        = [];
cfg.method = 'runica';
cfg.channel = {'MEG'}; % if you read more than MEG channels it reports NaN errors..
%cfg.runica.maxsteps = 50;
comp = ft_componentanalysis(cfg, data_preproc_clean_ds);


save([G.outputPath,'comp_ICA_', curSubj, '.mat'], 'comp');

load([G.outputPath,'comp_ICA_', curSubj, '.mat']);

% decompose the original data as it was prior to downsampling to 150Hz
cfg           = [];
cfg.unmixing  = comp.unmixing;
cfg.topolabel = comp.topolabel;
cfg.channel = {'MEG'};
comp_orig = ft_componentanalysis(cfg, data_preproc_clean);
%comp_orig = comp;

%% Visualize components
cfg          = [];
cfg.channel  = [1:10]; % components to be plotted
%cfg.blocksize = 30; % plot time in sec
cfg.viewmode = 'component';
cfg.compscale = 'local';
cfg.layout   = 'CTF275.lay'; % specify the layout file that should be used for plotting
ft_databrowser(cfg, comp)


%% remove the bad components and backproject the data

cfg = [];
cfg.component = [1 2 9]; % to be removed component(s)
cfg.channel = 'MEG';
data_preproc_clean  = ft_rejectcomponent(cfg,  comp_orig, data_preproc_clean) % 

        
    
    
%% Reinspect the data
cfg                 = [];
cfg.channel = 'MEG';
cfg.method =  'summary';
cfg.layout   = 'CTF275.lay';
ft_rejectvisual(cfg,data_preproc_clean)
%%


% inspect data before and after ICA (are the blinks and channel artifacts
% removed?)
cfg          = [];
cfg.renderer = 'painters';
cfg.compscale = 'local';
cfg.channel = 'MEG';
%cfg.layout   = 'CTF275.lay'; % specify the layout file that should be used for plotting
%ft_databrowser(cfg, data_preproc)
ft_databrowser(cfg, data_preproc_clean)



%% if damaged/missing channels, repair them...
% only works with original data (before or after ICA?)
if exist([G.outputPath,'rejchan_',  curSubj, '.mat'], 'file')
    load([G.outputPath,'rejchan_',  curSubj, '.mat']);
    cfg = [];
    cfg.channel = 'MEG';
    cfg.method = 'template';
    neighbours = ft_prepare_neighbours(cfg,data_preproc); % you need a complete template before channels are removed
    
    cfg = [];
    cfg.missingchannel = rejchan;
    cfg.neighbours = neighbours;
    cfg.channel = 'MEG';
    cfg.method = 'weighted'; % if I use weigthed it fails. I think that the problem is that there are two broken channels close in space (weight dont work properly)
    data_preproc_clean = ft_channelrepair(cfg, data_preproc_clean);
end
    
data_preproc_clean.cfg.previous = [];

%% Good! saving repaired channels
save([G.outputPath 'data_preproc_clean'],'data_preproc_clean','-v7.3')






