eegmp_preproc/code/a6_automatic_artifact_correction/H_STSW_automatic_artifact_correction_180123_noTardis.m

function H_STSW_automatic_artifact_correction_180123_noTardis()

%% H_SS_automatic_artifact_correction_20170922

% 180326 | instead of equalizing to smallest amount of samples, interpolate
%           for the same timing across subjects. Note that this does not
%           counteract any natural jitter in the paradigm and should only
%           equalize the size of matrices (fill with zeros! at ends).
%           Artefact detection does NOT currently work with potential NaN
%           values that would naturally arise from the interpolation.

%% initialize

%restoredefaultpath;
%clear all; close all; pack; clc;

%% pathdef

if ismac
    pn.study        = '/Volumes/LNDG/Projects/StateSwitch/';
    pn.eeg_root     = [pn.study, 'dynamic/data/eeg/task/A_preproc/SA_preproc_study/'];
    pn.dynamic_In   = [pn.eeg_root, 'B_data/B_EEG_ET_ByRun/'];
    pn.triggerTiming= [pn.eeg_root, 'C_figures/D_TriggerTiming/'];
    pn.EEG          = [pn.eeg_root, 'B_data/C_EEG_FT/']; mkdir(pn.EEG);
    pn.History      = [pn.eeg_root, 'B_data/D_History/']; mkdir(pn.History);
    pn.logs         = [pn.eeg_root, 'Y_logs/H_ArtCorr/'];
    % add ConMemEEG tools
    pn.MWBtools     = [pn.eeg_root, 'T_tools/fnct_MWB/'];           addpath(genpath(pn.MWBtools));
    pn.THGtools     = [pn.eeg_root, 'T_tools/fnct_THG/'];           addpath(genpath(pn.THGtools));
    pn.commontools  = [pn.eeg_root, 'T_tools/fnct_common/'];        addpath(genpath(pn.commontools));
    pn.fnct_JQK     = [pn.eeg_root, 'T_tools/fnct_JQK/'];           addpath(genpath(pn.fnct_JQK));
    pn.FT           = [pn.eeg_root, 'T_tools/fieldtrip-20170904/']; addpath(pn.FT); ft_defaults;
    pn.helper       = [pn.eeg_root, 'A_scripts/helper/'];           addpath(pn.helper);
else
    pn.root         = '/home/mpib/LNDG/StateSwitch/WIP_eeg/SA_preproc_study/';
    pn.EEG          = [pn.root, 'B_data/C_EEG_FT/'];
    pn.History      = [pn.root, 'B_data/D_History/'];
    pn.triggerTiming= [pn.root, 'C_figures/D_TriggerTiming/'];
    pn.logs         = [pn.root, 'Y_logs/H_ArtCorr/'];
end

%% define Condition & IDs for preprocessing

condEEG = 'dynamic';

%% define IDs for segmentation

% N = 47 YAs + 53 OAs;
IDs = {'1117';'1118';'1120';'1124';'1126';'1131';'1132';'1135';'1136';'1138';...
    '1144';'1151';'1158';'1160';'1163';'1164';'1167';'1169';'1172';'1173';...
    '1178';'1182';'1215';'1216';'1219';'1221';'1223';'1227';'1228';'1233';...
    '1234';'1237';'1239';'1240';'1243';'1245';'1247';'1250';'1252';'1257';...
    '1261';'1265';'1266';'1268';'1270';'1276';'1281';...
    '2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';...
    '2130';'2131';'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';...
    '2149';'2157';'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';...
    '2211';'2213';'2214';'2215';'2216';'2217';'2219';'2222';'2224';'2226';...
    '2227';'2236';'2237';'2238';'2241';'2244';'2246';'2248';'2250';'2251';...
    '2252';'2258';'2261'};

%%  loop IDs

ID_unavailable = cell(length(IDs),1);
for id = 1:length(IDs)
    try
        display(['processing ID ' num2str(IDs{id})]);
        for iRun = 1:4
        
            % load config
            load([pn.History, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_config'],'config')
            
            if ~isfield(config, 'ArtDect')
            
                % load data
                load([pn.EEG, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_EEG_Rlm_Fhl_rdSeg'],'data_EEG');

                dataICA = load([pn.EEG, IDs{id}, '_', condEEG, '_EEG_Rlm_Fhl_Ica'],'data');

                data = data_EEG; clear data_EEG;
                data.elec = dataICA.data.elec;
                data.chanlocs = dataICA.data.chanlocs; clear dataICA;

               %% ------------------ ARTIFACT DETECTION - PREPARATION ----------------- %%

                config.elec = data.elec;

                %% ICA (from weights)

                % ica config
                cfg.method           = 'runica';
                cfg.channel          = {'all','-ECG','-A2'};
                cfg.trials           = 'all';
                cfg.numcomponent     = 'all';
                cfg.demean           = 'yes';
                cfg.runica.extended  = 1;

                % ICA solution
                cfg.unmixing     = config.ica1.unmixing;
                cfg.topolabel    = config.ica1.topolabel;

                % components
                comp = ft_componentanalysis(cfg,data);

                % clear cfg
                clear cfg data

                %% remove components

                % get IC labels
                iclabels = config.ica1.iclabels.manual;

                % cfg for rejecting components (reject: blinks, eye movements,
                % ecg, ref); JQK: added artifacts and emg
                cfg.component = sortrows([iclabels.bli(:); iclabels.mov(:); iclabels.hrt(:); iclabels.ref(:); iclabels.art(:); iclabels.emg(:)])';
                cfg.demean    = 'yes';

                % reject components
                data = ft_rejectcomponent(cfg,comp);

                % clear cfg
                clear cfg comp

                %% remove eye & reference channels

                cfg.channel     = {'all','-IOR','-LHEOG','-RHEOG','-A1'};
                cfg.demean      = 'yes';

                % remove channels
                tmpdat = ft_preprocessing(cfg,data);

                % clear cfg & data variable
                clear cfg data


                %% ------------------------- ARTIFACT DETECTION ------------------------ %%

                % open log file
                fid = fopen([pn.logs, 'log_' IDs{id}, '_r',num2str(iRun), '_' condEEG '_ArtCorr.txt'],'a');

                % write log
              % fprintf(fid,['*********************  ' BLOCK '  *********************\n']);                         % Undefined function or variable 'BLOCK'. but see: %% H_prep_data_for_analysis_20141218

                fprintf(fid,['*********************   BLOCK   *********************\n']);
                fprintf(fid,['function: H_SS_automatic_artifact_correction_20170922.m \n\n']);                
                fprintf(fid,['eeg file = ' config.data_file '\n\n']);

                n_trl = length(tmpdat.trial);

                %%  get artifact contaminated channels by kurtosis, low & high frequency artifacts

                cfg.criterion = 3;
                cfg.recursive = 'no';

                [index0 parm0 zval0] = cm_MWB_channel_x_epoch_artifacts_20170922(cfg,tmpdat);

                % write log
                tmp_log = '';
                for j = 1:length(index0.c)
                    tmp_log = [tmp_log num2str(index0.c(j)) ' '];
                end; clear j
                tmp_log = [tmp_log(1:end-1) '\n'];
                fprintf(fid,'(1) automatic bad channel detection:\n');
                fprintf(fid,['MWB:          channel(s) ' tmp_log]);

                % clear cfg
                clear cfg tmp_log

                %%  get artifact contaminated channels by FASTER

                cfg.criterion = 3;
                cfg.recursive = 'no';

                [index1 parm1 zval1] = THG_FASTER_1_channel_artifacts_20140302(cfg,tmpdat);

                % write log
                tmp_log = '';
                for j = 1:length(index1)
                    tmp_log = [tmp_log num2str(index1(j)) ' '];
                end; clear j
                tmp_log = [tmp_log(1:end-1) '\n'];
                fprintf(fid,['FASTER:      channel(s) ' tmp_log]);

                % clear cfg
                clear cfg tmp_log

                %%  interpolate artifact contaminated channels

                % collect bad channels
                badchan = unique([index0.c; index1]);

                fprintf(fid,['--> ' num2str(length(badchan)) ' channels interpolated\n\n']);

                cfg.method     = 'spline';
                cfg.badchannel = tmpdat.label(badchan);
                cfg.trials     = 'all';
                cfg.lambda     = 1e-5; 
                cfg.order      = 4; 
                cfg.elec       = config.elec;

                % interpolate
                tmpdat = ft_channelrepair(cfg,tmpdat);

                % clear cfg
                clear cfg

                %% interpolate data to requested time (necessary due to acquisition jitter)
                
                % 2160 has only 62 trials, i.e. the first two trials were
                % not recorded.
                
                timeRequested = -1.5:1/500:9.566;
                
                for indTrial = 1:size(tmpdat.trial,2)
                    recTime = tmpdat.time{indTrial};
                    recData = tmpdat.trial{indTrial};
                    for indChannel = 1:size(recData,1)
                        interpData(indChannel,:) = interp1(recTime,recData(indChannel,:),timeRequested,'linear');
                    end
                    interpData(isnan(interpData)) = 0;
                    tmpdat.trial{indTrial} = interpData;
                    tmpdat.time{indTrial} = timeRequested;
                end
                
                %figure; plot(timeRequested,interpData, 'k'); hold on; plot(recTime,recData(1,:), 'r')
                
                %%  get artifact contaminated epochs & exclude epochs
                % includes: - THG_MWB_channel_x_epoch_artifacts_20140311
                %           - THG_FASTER_2_epoch_artifacts_20140302
                % recursive epoch exclusion!

                [tmpdat index] = THG_automatic_artifact_correction_trials_20170922(tmpdat);

                % write log
                fprintf(fid,'(2) automatic recursive bad epoch detection:\n');
                fprintf(fid,['MWB & FASTER: ' num2str(n_trl-length(index)) ' bad epoch(s) detected\n\n']);

                %%  get channel x epoch artifacts

                % cfg
                cfg.criterion = 3;
                cfg.recursive = 'yes';

                [index3 parm3 zval3] = THG_FASTER_4_channel_x_epoch_artifacts_20140302(cfg,tmpdat);

                % write log
                fprintf(fid,'(3) automatic single epoch/channel detection:\n');
                fprintf(fid,['FASTER:       ' num2str(sum(sum(index3))) ' channel(s) x trial(s) detected\n\n']);

                % clear cfg
                clear cfg

                %%  collect and save detected artifacts & artifact correction infos

                % include ArtDect.parameters

                % bad channels
                ArtDect.channels = badchan;

                % bad trials
                tmp  = zeros(n_trl,1); tmp(index,1) = 1;
                badtrl = find(tmp==0);
                ArtDect.trials  = badtrl; 
                clear tmp

                % bad single epoch(s)/channel(s) - after exclusion of bad epochs
                ArtDect.channels_x_trials = index3;

                % overview
                ind = [1:n_trl];
                ind(badtrl) = [];
                tmp = ones(length(tmpdat.label),n_trl);
                tmp(:,ind) = index3;
                tmp(badchan,:) = 1;
                ArtDect.channels_x_trials_all = tmp;
                clear tmp

                % save version
                ArtDect.version = 'JQK 20180123';

                % add to config
                config.ArtDect = ArtDect;

                % save config
                save([pn.History, IDs{id}, '_r',num2str(iRun), '_', condEEG, '_config'],'config')

                %%  finalize log

                % log
                fprintf(fid,'Artifact detection completed.\n');
                fprintf(fid,['Information saved to: ' IDs{id} '_config.mat\n\n']);
                fclose(fid);
                
                
            else
                disp([IDs{id}, ' Run ',num2str(iRun), ' has already run through automatic artifact detection']);
            end % if data was not processed yet
        end % run
    catch ME
        warning('Error occured. Please check.');
        rethrow(ME)
        fprintf('\n ID not availble! Skip! \n')
        ID_unavailable{id,1} = (IDs{id});
    end
    
end; clear id

end