% EEG_EEGREJ - reject porition of continuous data in an EEGLAB 
%                dataset
%
% Usage:
%   >> EEGOUT = eeg_eegrej( EEGIN, regions );
%
% Inputs:
%   INEEG      - input dataset
%   regions    - array of regions to suppress. number x [beg end]  of 
%                regions. 'beg' and 'end' are expressed in term of points
%                in the input dataset. Size of the array is
%                number x 2 of regions.
%
% Outputs:
%   INEEG      - output dataset with updated data, events latencies and 
%                additional boundary events.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 8 August 2002
%
% See also: EEGLAB, EEGPLOT, POP_REJEPOCH

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, com] = eeg_eegrej( EEG, regions)

com = '';
if nargin < 2
	help eeg_eegrej;
	return;
end
if nargin<3
    probadded = [];
end
if isempty(regions)
	return;
end

% regions = sortrows(regions,3); % Arno and Ramon on 5/13/2014 for bug 1605

% Ramon on 5/29/2014 for bug 1619
if size(regions,2) > 2
    regions = sortrows(regions,3);
else
    regions = sortrows(regions,1);
end

try
    % For AMICA probabilities...Temporarily add model probabilities as channels
    %-----------------------------------------------------
    if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
        if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
            if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

                EEG = eeg_formatamica(EEG);
                %-------------------------------------------

                [EEG com] = eeg_eegrej(EEG,regions);

                %-------------------------------------------

                EEG = eeg_reformatamica(EEG);
                EEG = eeg_checkamica(EEG);
                return;
            else
                disp('AMICA probabilities not compatible with size of data, probabilities cannot be epoched')
                disp('Load AMICA components before extracting epochs')
                disp('Resuming rejection...')
            end
        end

    end
    % ------------------------------------------------------
catch
    warnmsg = strcat('your dataset contains amica information, but the amica plugin is not installed.  Continuing and ignoring amica information.');
    warning(warnmsg)
end

% handle regions from eegplot
% ---------------------------
if size(regions,2) > 2, regions = regions(:, 3:4); end
regions = combineregions(regions);

% remove events within regions
% ----------------------------
if ~isempty(EEG.event) && isfield(EEG.event, 'latency')
    allEventLatencies = [ EEG.event.latency];
    allEventFlag      = zeros(1,length(allEventLatencies));
    for iRegion = 1:size(regions,1)
        allEventFlag = allEventFlag | ( allEventLatencies >= regions(iRegion,1) & allEventLatencies <= regions(iRegion,2));
    end
    boundaryIndices = eeg_findboundaries(EEG);
    allEventFlag(boundaryIndices) = false; % do not remove boundary events
    EEG.event(allEventFlag) = [];
end

% reject data
% -----------
[EEG.data, EEG.xmax, event2, boundevents] = eegrej( EEG.data, regions, EEG.xmax-EEG.xmin, EEG.event);
oldEEGpnts = EEG.pnts;
oldEEGevents = EEG.event;
EEG.pnts   = size(EEG.data,2);
EEG.xmax   = EEG.xmax+EEG.xmin;
if length(event2) > 1 && event2(1).latency == 0, event2(1) = []; end
if length(event2) > 1 && event2(end).latency == EEG.pnts, event2(end) = []; end
if length(event2) > 2 && event2(end).latency == event2(end-1).latency
    if isfield(event2, 'type') && isequal(event2(end).type, event2(end-1).type)
        event2(end) = []; 
    end
end

% add boundary events
% -------------------
[ EEG.event ] = eeg_insertbound(EEG.event, oldEEGpnts, regions);
EEG = eeg_checkset(EEG, 'eventconsistency');
if ~isempty(EEG.event) && EEG.trials == 1 && EEG.event(end).latency-0.5 > EEG.pnts
    EEG.event(end) = []; % remove last event if necessary
end

% double check event latencies
% the function that insert boundary events and recompute latency is
% delicate so we do it twice using different methods and check
% the results. It takes longer, but accuracy is paramount.
eeglab_options;
warnflag = false;
if isfield(EEG.event, 'latency') && length(EEG.event) < 3000
    alllats = [ EEG.event.latency ];
    if ~isempty(event2)
        otherlatencies = [event2.latency];
        if ~isequal(alllats, otherlatencies)
            warning([ 'Discrepancy when checking event latencies using legacy method.' 10 'Often the discrepancy is minor and the new method (used here) is correct' 10 'still, try to reproduce the problem and send us your dataset' ]);
            warnflag = true;
        end
    end
end

% double check boundary event latencies
if ~isempty(EEG.event) && length(EEG.event) < 3000 && ischar(EEG.event(1).type) && isfield(EEG.event, 'duration') && isfield(event2, 'duration')
    try
        if ~isempty(EEG.event) && ischar(EEG.event(1).type)
            indBound1 = find(cellfun(@(x)strcmpi(num2str(x), 'boundary'), { EEG.event(:).type }));
            indBound2 = find(cellfun(@(x)strcmpi(num2str(x), 'boundary'), { event2(:).type }));
        else
            indBound1 = find([ EEG.event(:).type ] == -99);
            indBound2 = find([ event2(:).type ] == -99);
        end
        duration1 = [EEG.event(indBound1).duration]; duration1(isnan(duration1)) = [];
        duration2 = [event2(indBound2).duration]; duration2(isnan(duration2)) = [];
        if ~isequal(duration1, duration2)
            duration1(duration1 == 0) = [];
            if ~isequal(duration1, duration2) && ~warnflag
                warning([ 'Inconsistency in boundary event duration using legacy method.' 10 'Often the discrepancy is minor and the new method (used here) is correct' 10 'still, try to reproduce the problem and send us your dataset' ]);
            end
        end
    catch, warning('Unknown error when checking event latency - please send us your dataset');
    end
end

% debugging code below
% regions, n1 = 1525; n2 = 1545; n = n2-n1+1;
% a = zeros(1,n); a(:) = 1; a(strmatch('boundary', { event2(n1:n2).type })') = 8; 
% [[n1:n2]' alllats(n1:n2)' [event2(n1:n2).latency]' alllats(n1:n2)'-[event2(n1:n2).latency]' otherorilatencies(n1:n2)' a']
% figure; ev = 17; range = [-1000:1000]; plot(EEG.data(1,EEG.event(ev).latency+range)); hold on; plot(tmpdata(1,tmpevent(EEG.event(ev).urevent).latency+range+696), 'r'); grid on;

com = sprintf('EEG = eeg_eegrej( EEG, %s);', vararg2str({ regions })); 

% combine regions if necessary
% it should not be necessary but a 
% bug in eegplot makes that it sometimes is
% ----------------------------
% function newregions = combineregions(regions)
% newregions = regions;
% for index = size(regions,1):-1:2
%     if regions(index-1,2) >= regions(index,1)
%         disp('Warning: overlapping regions detected and fixed in eeg_eegrej');
%         newregions(index-1,:) = [regions(index-1,1) regions(index,2) ];
%         newregions(index,:)   = [];
%     end
% end

function res = issameevent(evt1, evt2)

res = true;
if isequal(evt1,evt2)
    return;
else
    if isfield(evt1, 'type') && isnumeric(evt2.type) && ~isnumeric(evt1.type) 
        evt2.type = num2str(evt2.type);
        if isequal(evt1,evt2)
            return;
        end
    end
    if isfield(evt1, 'duration') && isfield(evt2, 'duration')
        if isnan(evt1.duration) && isnan(evt2.duration)
            evt1.duration = 1;
            evt2.duration = 1;
        end
        if abs( evt1.duration - evt2.duration) < 1e-10
            evt1.duration = 1;
            evt2.duration = 1;
        end
        if isequal(evt1,evt2)
            return;
        end
    end
    if isfield(evt1, 'latency') && isfield(evt2, 'latency')
        if abs( evt1.latency - evt2.latency) < 1e-10
            evt1.latency = 1;
            evt2.latency = 1;
        end
        if isequal(evt1,evt2)
            return;
        end
    end
end
res = false;
return;

function newregions = combineregions(regions)
% 9/1/2014 RMC
regions = sortrows(sort(regions,2)); % Sorting regions
allreg = [ regions(:,1)' regions(:,2)'; ones(1,numel(regions(:,1))) -ones(1,numel(regions(:,2)')) ].';
allreg = sortrows(allreg,1); % Sort all start and stop points (column 1),

mboundary = cumsum(allreg(:,2)); % Rationale: regions will start always with 1 and close with 0, since starts=1 end=-1
indx = 0; count = 1;

while indx ~= length(allreg) 
    newregions(count,1) = allreg(indx+1,1);
    [tmp,I]= min(abs(mboundary(indx+1:end)));
    newregions(count,2) = allreg(I + indx,1);
    indx = indx + I ;
    count = count+1;
end

% Verbose
if size(regions,1) ~= size(newregions,1)
    disp('Warning: overlapping regions detected and fixed in eeg_eegrej');
end
