% - This function is part of ERPLAB Toolbox -
%
% >> A = areaerp(data, fs,latsam, op, coi)
%
% or
%
% >> [A L] = areaerp(data, fs,latsam, op, coi)
%
% where
%
% data     - input array of data
% fs       - sampling frequency
% latsam   - integration limits in samples. Two values for fixed integration limits.
%            One or two values for automatic interation limit seeding.
% op       - 'auto' means look for automatic interation limits. Otherwise is fixed.
% coi      - Only for op='auto'. When an ERP waveform is not delimited by 2 zero-crossing latencies,
%            but the first one (lower integration limit), the second latency (upper integration limit)
%            will be the offset (if any) of a second or higher ERP waveform.
%            Then coi=1 will set the first zero crossing as the lower integration limit, and the
%            latency of the minimum values between the first and the second
%            ERP as the upper integration limit.
%            Then coi=2 will set the latency of the minimum values between the first and the second
%            ERP as the lower integration limit, and the second zero crossing as the upper integration limit.
%
% Author: Javier Lopez-Calderon & Steven Luck
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

%b8d3721ed219e65100184c6b95db209bb8d3721ed219e65100184c6b95db209b
%
% ERPLAB Toolbox
% Copyright © 2007 The Regents of the University of California
% Created by Javier Lopez-Calderon and Steven Luck
% Center for Mind and Brain, University of California, Davis,
% javlopez@ucdavis.edu, sjluck@ucdavis.edu
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = areaerp(data, fs,latsam, op, coi, frac, fracarearep)

if nargin<7
      fracarearep = 0; %NaN;  in case frac area lat is not found.
end
if nargin<6
      frac = 0.5; % for 50% fractional area latency, by default
end
if nargin<5
      coi = 1; % component of interest
end
if nargin<4
      op = 'integral';
end
if length(latsam)==1 && ~ismember({op},{'auto','autot','autop','auton'})
      error('ERPLAB says: error at areaerp(). You must enter 2 latencies for non-automatic limits')
end
if frac<0 || frac>1
      error('ERPLAB says: error at areaerp(). Fractional area value must be between 0 and 1')
end

fsn = 2000;     % new sample rate 10Khz
Tsn = 1/fsn;     % new sample period
Fk  = fsn/fs;    % oversampling factor

if Fk>1 % interpolate data
      nsam = length(data);
      aa = 1; bb = nsam;
      newsam   = linspace(aa,bb,Fk*(bb-aa+1));
      data     = spline(aa:bb, data(aa:bb),newsam); % interpoled data    
      latsam   = round(latsam*Fk)-1;
else
      Fk=1;
end
if ismember({op},{'auto','autot','autop','auton'})
      if length(latsam)==2
            
            %
            % estimates a seed latency to search integration limits, according with
            % the weighted average of the two given latencies.
            %
            vl1 = abs(data(latsam(1)));
            vl2 = abs(data(latsam(2)));
            k1  = vl1/(vl1 + vl2);
            k2  = vl2/(vl1 + vl2);
            latsam = round(k1*latsam(1) +  k2*latsam(2));
      end
      
      %
      % Limit of integrations (automatic)
      %
      vinit = data(latsam);   % voltage at latency
      vsign = sign(vinit);
      b     = find(sign(data(latsam:end))~= vsign,1,'first') + latsam -1; % in sample
      a     = find(sign(data(1:latsam))  ~= vsign,1,'last'); % in sample
      
      if isempty(b)
            b = length(data);
      end
      if isempty(a)
            a = 1;
      end      
      if coi>0
            %
            % Test for detection of overlaped components
            %
            data(1:a)   = 0;
            data(b:end) = 0;
            data = abs(data); % rectification
            ndata   = length(data);
            [datamax imax] = max(data) ;
            ion     = find(data(1:imax)>0.2*datamax, 1,'first'); % onset, at 10% max
            ioff    = find(data((imax+1):end)>0.2*datamax, 1,'last') + imax; % onset, at 10% max
            onset   = data(ion);
            offset  = data(ioff);
            x  = [a ion imax ioff b];
            y  = [0 onset datamax offset 0];
            xx = linspace(a,b,b-a+1);
            simerp = interp1(x, y,xx,'cubic');
            AT1 = single(trapz(data(a:b)));
            AT2 = single(trapz(simerp));
            ATrate = AT1/AT2;
            
            if ATrate<0.9
                  fprintf('***************************************\n')
                  fprintf('overlapping component detected!\n')
                  fprintf('***************************************\n\n')
                  
                  if coi==1
                        fprintf('Correcting upper integration limit...\n');
                  elseif coi==2
                        fprintf('Correcting lower integration limit...\n');
                  end
                  
                  step = 5;
                  
                  for i = (1+step):ndata-step                        
                        leftp  = data(i-step:i-1);
                        rigthp = data(i+1:i+step);
                        targp  = data(i);
                        conL   = length(find(leftp>targp));
                        conR   = length(find(rigthp>=targp));
                        
                        if conL==5 && conR==5                              
                              if coi==1
                                    fprintf('old upper limit: %g\n', b)
                                    fprintf('corrected upper limit: %g\n', i)
                                    b = i;
                              elseif coi==2
                                    fprintf('old lower limit: %g\n', a)
                                    fprintf('corrected lower limit: %g\n', i)
                                    a = i;
                              end
                              break
                        end
                  end
            end
      end
else
      try
            a = latsam(1); % in sample
            b = latsam(2); % in sample
      catch
            error('ERPLAB says: 2 latency values are required.')
      end
end
if isempty(a) || isempty(b)
      A = 0;   % either of the integration limits failed
else      
      if ismember({op},{'autot','total'})
            data = abs(data);
      elseif ismember({op},{'autop','positive'})
            data(data<0) = 0;
      elseif ismember({op},{'auton','negative'})
              data = -data;
              data(data<0) = 0; % negative values get positive, and the previously positive ones get zeroed.
      end
      
      %
      %  Calculates area
      %
      A = single(Tsn*trapz(data(a:b))); % April 21, 2011
end
if nargout==1
      varargout{1} = A;
elseif nargout==2 || nargout==3
      
      varargout{1} = A;
      
      if A~=0
            
            %
            % Finds fractional (50% by default) area latency  (in samples)
            %
            Ao = single(A*frac); % Fractional area
            
            if Ao~=0
                  %
                  % (convergent) binary search algorithm
                  %
                  plow  = a;
                  phigh = b;
                  
                  while plow < phigh                        
                        pmid = round((plow + phigh) / 2);
                        Ax   = single(Tsn*trapz(data(a:pmid)));
                        
                        if Ax > Ao
                              phigh = pmid - 1;
                        elseif Ax < Ao
                              plow = pmid + 1;
                        else
                              break     % It was found!
                        end
                  end                  
                  L = pmid;  % Frac latency in (non-integer) samples                  
            else
                  L = a; % if Ao=0 then L=lower area limit
            end
      else
            if fracarearep==0 % Fractional area latency replacement
                  L = NaN; % if Ao=0 then fractional latency = NaN
            else
                  error('ERPLAB:zeroarea',['Sorry, fractional area latency was not found because area value is zero.\n'...
                  'You may use fractional area latency replacement option.'])
            end
      end
      
      %
      % Frac latency
      %
      varargout{2} = L/Fk; % rescale. samples
      
      if nargout==3
            
            %
            % Integration limits
            %
            varargout{3} = [a b]/Fk; %in samples
      end      
elseif nargout>3
      error('ERPLAB says: error at areaerp(). Too many output arguments!')
end