function varargout = process_evt_detect_seeg_es( varargin )
% PROCESS_EVT_DETECT_SEEG_ES: Epileptic spikes detection for a group of recordings file
%
% USAGE:  OutputFiles = process_evt_detect_seeg_es('Run', sProcess, sInputs)

% @=============================================================================
% This function aims to reduce potential ictal artifacts in SEEG data according to the following study:
% Hirshorn, E. a., Li, Y., Ward, M. J., Richardson, R. M., Fiez, J. a., & Ghuman, A. S. (2016). Decoding and disrupting
% left midfusiform gyrus activity during word reading. PNAS, 113(29), 201604126. https://doi.org/10.1073/pnas.1604126113
% 
% Copyright (c)2000-2020 University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Shuai Wang, 2020

eval(macro_method);
end


%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Detect potential epileptic artifacts (Hirshorn 2016)';
    sProcess.FileTag     = 'ES';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 46;  % the same Index as PROCESS_EVT_DETECT_BADSEGMENT
    sProcess.Description = 'See SI in Hirshorn et al., 2016, PNAS.';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    
%     % Time window
%     sProcess.options.timewindow.Comment   = 'Time window: ';
%     sProcess.options.timewindow.Type      = 'timewindow';
%     sProcess.options.timewindow.Value     = [];
    % Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'SEEG';    
    % Threshold1 - SDs from the average amplitude
    sProcess.options.threshold1.Comment = {'x SDs from the Mean: '};
    sProcess.options.threshold1.Type    = 'value';
    sProcess.options.threshold1.Value   = {5, 'SDs', 0};
    % Threshold2 - absolute upper limit
    sProcess.options.threshold2.Comment = {'The upper limit of peak amplitude (abs.): '};
    sProcess.options.threshold2.Type    = 'value';
    sProcess.options.threshold2.Value   = {350, 'uV', []};    
    % Threshold3 - variation between consecutive time points
    sProcess.options.threshold3.Comment = {'The threshold of variation between consecutive points (abs.): '};
    sProcess.options.threshold3.Type    = 'value';
    sProcess.options.threshold3.Value   = {25, 'uV', []};     

end


%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end


%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % ===== GET OPTIONS ===== 
    ThreshSDs   = sProcess.options.threshold1.Value{1};
    ThreshPeak  = sProcess.options.threshold2.Value{1} / 10e5;  % convert microvolts to volts
    ThreshVar   = sProcess.options.threshold3.Value{1} / 10e5;  % convert microvolts to volts  
    SensorTypes = sProcess.options.sensortypes.Value;    
%     % Time window to process
%     if isfield(sProcess.options, 'timewindow') && isfield(sProcess.options.timewindow, 'Value') && iscell(sProcess.options.timewindow.Value) && ~isempty(sProcess.options.timewindow.Value)
%         TimeWindow = sProcess.options.timewindow.Value{1};
%     else
%         TimeWindow = [];
%     end 
   
    % Option structure for function in_fread()
    ImportOptions = db_template('ImportOptions');
    ImportOptions.ImportMode = 'Time';
    ImportOptions.RemoveBaseline = 'no';
    % Get current progressbar position
    progressPos = bst_progress('get');

    % For each file
    iOk = false(1, length(sInputs));
    for iFile = 1:length(sInputs)
        % ===== GET DATA =====
        % Progress bar
        bst_progress('text', 'Reading trial to process...');
        % select good SEEG channels
        ChannelMat = in_bst_channel(sInputs(iFile).ChannelFile);
        sChannel = channel_find(ChannelMat.Channel, SensorTypes);
        ChannelName = {ChannelMat.Channel(:).Name}';
        % Load epochs        
        DataES = in_bst_data(sInputs(iFile).FileName, 'F', 'Time', 'ChannelFlag', 'Comment');
        goodChannels = logical(DataES.ChannelFlag == 1);
        DataES.goodChannels = ChannelName(sChannel(goodChannels(sChannel)));
        DataES.goodF = DataES.F(sChannel(goodChannels(sChannel)), :);
        
        % ===== automatic detection Hirshorn 2016 =====
        DataES.nTime = length(DataES.Time);
        [DataES.peakF, DataES.peakIdx] = max(abs(DataES.goodF), [], 2);
        % peak amplitude x SDs away from the Mean across the rest of the trial
        DataES.restF = DataES.goodF(logical(repmat(1:DataES.nTime, length(DataES.goodChannels), 1) ~= DataES.peakIdx));
        DataES.restF = reshape(DataES.restF, [], DataES.nTime - 1);
        DataES.restMean = mean(DataES.restF, 2);
        DataES.stdF = std(DataES.restF, 0, 2);
        DataES.checkThreshSDs = logical(DataES.peakF > ThreshSDs * abs(DataES.stdF) + abs(DataES.restMean));
        % peak amplitude greater than the upper limit (350uV by default)
        DataES.checkThreshPeak = logical(DataES.peakF > ThreshPeak);
        % variation threshold between consecutive points (25uV by default)
        DataES.diffF = diff(DataES.goodF, 1, 2);
        DataES.checkThreshVar = any(DataES.diffF > ThreshVar, 2);
        % count number of channels with potential ES
        DataES.nES = sum((DataES.checkThreshSDs + DataES.checkThreshPeak + DataES.checkThreshVar) ~= 0);
        
        % ===== SAVE RESULT =====
        % Progress bar
        bst_progress('text', 'Saving results...');
        bst_progress('set', progressPos + round(3 * iFile / length(sInputs) / 3 * 100));
        % Only save changes if something was detected
        if any(DataES.checkThreshSDs) || any(DataES.checkThreshPeak) || any(DataES.checkThreshVar)
            % Save file definition
            save(file_fullpath(sInputs(iFile).FileName), 'DataES', '-append');
            % Report number of detected events
            bst_report('Info', sProcess, sInputs(iFile), sprintf('For trial %s: potential ES detected in %d channels', DataES.Comment, DataES.nES));
            fprintf('For trial %s: potential ES were detected in %d channels.\n', DataES.Comment, DataES.nES)
        else
            bst_report('Info', sProcess, sInputs(iFile), 'No ES detected');
            fprintf('For trial %s: No ES detected.\n', DataES.Comment)
            iOk(iFile) = true;
        end
        
    end
    % Return all the input files
    OutputFiles = {sInputs(iOk).FileName};
%     % Reference OutputFile in the database:
%     db_add_data(sInput.iStudy, OutputFile, DataES);
end







