%%% PREPROCESSING  %%%
close all
clear
addpath('/home/common/matlab/fieldtrip-20170630/');
ft_defaults


subject_options;

% Nsb =27;%:2;%1:4;%input('Subjects number(s): '); %% which subject(s) to process
% Nsb =[2 3 6 12 14];
Nsb =1:numel(subjname);
sb = Nsb;

cwd = pwd;


for sb = Nsb
    do_contrej     = 0;
    do_trialdef     = 1;
    do_musclejump   = 0;
    do_visremove    = 0;
    do_calcica      = 0;
    do_compremove   = 0;
    do_eogremove    = 0;
    do_localizer    = 0;
    do_restore_old  = 0;
    
    
    trialdef.interval = [-0.5, 1.5]; % cut out a trial between 500 ms pre-stim and 1500 ms
    trialdef.prestim=0.5;
    trialdef.poststim=1.5;
    cwd = pwd;
    
    
    fprintf('Preprocessing subject %s... \n', subjname{sb});

    if do_localizer        
        trialdef.trig.eventnum = [11:13];
        infofile = strcat(infopath,subjname{sb},'_info_localizer2.mat'); 
    else        
        trialdef.trig.eventnum = [21:26,71:76];
        infofile = strcat(infopath,subjname{sb},'_info_maintask2.mat');  
    end
    
    clear info;
    cfg = [];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Trial definition
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if do_contrej
        cfg = [];
        cfg.dataset = fullfile(['../Raw/',dataset{sb}]);
        cfg.channel = 'MEG';
        data = ft_preprocessing(cfg);
        % open the browser and page through the trials
        cfg=[];
        cfg.channel = 'MEG';
        artf=ft_databrowser(cfg,data);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Trial definition
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if do_trialdef
        if ~exist(infofile)
            answer = 'y'; % infofile can be created
        elseif exist(infofile)
            answer = input('Warning: infofile already exists. Overwrite (y/n)?','s');
        end
        
        if isequal(answer,'y')
            cfg.dataset = fullfile(['../Raw/',dataset{sb}]);
            cfg.trialdef=trialdef;
            cfg.trialdef.eventtype = 'gui';
            cfg.trialfun='trialfun_biao';
            cfg.trialdef.offset=0;
            cfg = ft_definetrial(cfg);
            info.preproc.cfg = cfg;
            save(infofile,'info');
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Muscle and jump artifact rejection
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if do_musclejump
        %     load(infofile);
        
        cfg = info.preproc.cfg;
        cfg.continuous = 'yes';
        cfg.channel = 'MEG';
        
        cfg.artfctdef.reject = 'complete';
        cfg.artfctdef.type = {'jump','muscle'};
        cfg.artfctdef.jump.cutoff   = 20;
        cfg.artfctdef.muscle.cutoff = 20;
        cfg.artfctdef.jump.feedback = 'yes';
        cfg.artfctdef.muscle.feedback = 'yes';
        
        cfg = ft_rejectartifact(cfg);
        
        info.preproc.cfg = cfg;
        info.preproc.cfg.version.do_musclejump = 'yes';
        save(infofile,'info','-APPEND');
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Visual inspection of cleaned up data to manually remove suspicious trials
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if do_visremove
        load(infofile);
        cfg = info.preproc.cfg;
        cfg.channel = 'MEG';
        cfg.dataset = fullfile(['../Raw/',dataset{sb}]);    
        data = ft_preprocessing(cfg);
%         matlabmail('raintwoto@gmail.com','ok','ok')
        
        data2 = ft_rejectvisual([],data);
        
        info.preproc.cfg.trlold = info.preproc.cfg.trl;
        
        tmp = info.preproc.cfg.trlold(1,3); %
        info.preproc.cfg.trl = [data2.sampleinfo tmp*ones(length(data2.sampleinfo),1) data2.trialinfo];
        info.preproc.cfg.version.do_visremove = 'yes';
        
        save(infofile,'info','-APPEND');
        
    end
    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate ICA components
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if do_calcica
        load(infofile);
        cfg = info.preproc.cfg;
        
        cfg.blc = 'yes';  % apply baseline correction; removes the mean of the signal
        cfg.detrend ='no';
        cfg.continuous = 'yes';
        cfg.channel = 'MEG';
        cfg.checkmaxfilter=1;
        data = ft_preprocessing(cfg);
        info.label = data.label;
        
        % resample to speed up the ICA decomposition, especially useful for 1200Hz MEG data
        cfg            = [];
        cfg.detrend    = 'no';
        cfg.resamplefs = 150;
        data           = ft_resampledata(cfg, data);
        
        % perform the independent component analysis, i.e. decompose the downsampled original data
        % this results in a mixing matrix and in the component timecourses
        warning('OFF','MATLAB:dispatcher:nameConflict'); eeglab_dir = '/home/action/flodlan/groupshare/matlab/eeglab'; addpath(genpath(eeglab_dir),'-END');
        cfg            = [];
        cfg.channel    = 'MEG';
        
        comp           = ft_componentanalysis(cfg, data);
        info.comp.topo = comp.topo;
        info.comp.topolabel = comp.topolabel;
        info.preproc.cfg.version.do_calcica = 'yes';
        clear comp;
        
        save(infofile,'info','-APPEND'); % save the indices of strongest ICA EOG and ECG components
        
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Confirm the bad ICA components
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if do_compremove
        load(infofile);
        load([matpath,subjname{sb},'_data_clean_noica.mat']);
        cfg = [];
        cfg.topo      = info.comp.topo;
        cfg.topolabel = info.comp.topolabel;
        comp_data     = ft_componentanalysis(cfg, data);
        
        cfg           = [];
        cfg.component = [1:20];       % specify the component(s) that should be plotted
        cfg.layout    = 'CTF275.lay'; % specify the layout file that should be used for plotting
        cfg.comment   = 'no';
        figure; ft_topoplotIC(cfg, comp_data)
        
        
        comp_artifact = input('type the array of all artifact components: ');
        info.comp.artifact = comp_artifact;
        
        save(infofile,'info','-APPEND');
        
    end
    %%
    if do_eogremove
        
        
        load(infofile);
        cfg = info.preproc.cfg;
        
             
        newcfg=cfg;
        
        
        % the ft_rejectartifact has bug
        % from 0 to 0.5 s
        for i=1:size(newcfg.trl,1)
            newcfg.trl(i,1)=newcfg.trl(i,1)+600;
            newcfg.trl(i,2)=newcfg.trl(i,2)-700;
        end
        
        [~, artifact_EOG] = ft_artifact_zvalue(newcfg);
        
        
        
        cfg = info.preproc.cfg;
        
        %    cfg=[];
        if cfg.artfctdef.type{4}
            cfg.artfctdef.type(4)=[];
        end
        
        
        cfg.artfctdef.eog.artifact    = artifact_EOG;
        
        cfg.artfctdef.reject = 'complete';
        cfg = ft_rejectartifact(cfg);
        
        info.preproc.cfg = cfg;
        info.preproc.cfg.version.do_eog = 'yes';
        save(infofile,'info','-APPEND');
    end
    %%
    %%
    if do_restore_old
        
        
        load(infofile);
        info.preproc.cfg.version.do_eog = 'no';
        info.preproc.cfg.trl = info.preproc.cfg.trlold;        
        save(infofile,'info','-APPEND');
    end
    
end

cd(cwd);

