function [n_j,kadd]=local_extrema(x,Deltan,varargin)
%LOCAL_EXTREMA Determines the local extrema of a given temporal scale
% n_j = The positions of the local extrema of a partition of scale Deltan
%       as defined at p. 82 in the book [ATE] C. Vamos and M. Craciun,
%       Automatic Trend Estimation, Springer 2012.
%       The positions of the maxima are positive and those of the minima
%       are negative.
% kadd= n_j(kadd) are the local extrema with time scale smaller than Deltan
%       which are added to the partitionsuch that an alternation of maxima
%       and minima is obtained.
% [n_j,kadd]=local_extrema(x,Deltan,ni,nf) determines the partition of
%       scale Deltan of the part x(ni:nf) of the time series x. The
%       algorithm is presented in Appendix E of [ATE].
% [n_j,kadd]=local_extrema(x,Deltan) determines the partition of
%       scale Deltan of the entire time series x
% [n_j,kadd]=local_extrema(x) determines the partition of scale Deltan=1
%       of the entire time series x which contains all the local extrema

N=length(x);
if nargin==1
    Deltan=1;end
if nargin<=2
    ni=1;nf=N;
else 
    ni=varargin{1}; nf=varargin{2};
end

% recognize the plateaus of the time series x defined in [ATE] p. 85
% [n1(n),n2(n)] is the interval with the constant value equal with x(n)
% if x(n) is not contained in a plateau, then n1(n)=n2(n)=n
n1=[1:N]; n2=[1:N];
dx=[1,diff(x)];
k0=find(dx==0); dimk0=length(k0);
if dimk0>0
    dk0=diff(k0);
    k1=[find(dk0>1),dimk0]; dimp=length(k1);
    nk1=1;
    for p=1:dimp
        n01=k0(nk1)-1;
        nk2=k1(p); n02=k0(nk2);
        n1(n01:n02)=n01; n2(n01:n02)=n02;
        nk1=nk2+1;
    end
end

% the iterative algorithm presented in Appendix E of [ATE]
nc=ni-1; % the time step up to which the time series has been acomplished ([ATE] p. 127)
Jest=0; % the number of the local extrema of the partition of scale DeltaN
flagante=0;
iadd=0; % number of the additional local extrema
kadd=[]; % the order number of the additional local extrema between all the local extrema
while 1
    if nc>=nf break; end
    
    nlim=nc+Deltan;
    if nlim>=nf nlim=nf; end
    % [nc+1,nlim] is the interval where the next extremum is searched
    
    [xmin nmin]=min(x(nc+1:nlim)); nmin=nmin+nc;
    nlim1=n1(nmin)-Deltan; if nlim1<ni nlim1=ni; end
    nlim2=n2(nmin)+Deltan; if nlim2>nf nlim2=nf; end
    [xminn nminn]=min(x(nlim1:nlim2)); nminn=nminn+nlim1-1;
    if nminn==nmin flagmin=1;
    else flagmin=0;
    end
    % if flagmin=1 then the minimum at nmin satisfies condition (6.1)
    
    [xmax nmax]=max(x(nc+1:nlim)); nmax=nmax+nc;
    nlim1=n1(nmax)-Deltan; if nlim1<ni nlim1=ni; end
    nlim2=n2(nmax)+Deltan; if nlim2>nf nlim2=nf; end
    [xminn nmaxx]=max(x(nlim1:nlim2)); nmaxx=nmaxx+nlim1-1;
    if nmaxx==nmax flagmax=1;
    else flagmax=0;
    end
    % if flagmax=1 then the minimum at nmax satisfies condition (6.1)
    
    % the extremum closest to nc is kept for analysis
    if flagmin==1 && flagmax==1
        if nmin<nmax flagmax=0;
        else flagmin=0;
        end
    end
    
    if flagante==0 % there is no anterior extremum introduced in the partition
        if flagmax==1
            nc=n1(nmax); flagante=1; Jest=Jest+1;
            n_j(Jest)=floor((n1(nmax)+n2(nmax))/2);
        elseif flagmin==1
            nc=n1(nmin); flagante=-1; Jest=Jest+1;
            n_j(Jest)=-floor((n1(nmin)+n2(nmin))/2);
        else
            nc=nc+Deltan;
        end
    elseif flagante==-1 % the anterior extremum is an minimum
        if flagmax==1 % the current extremum is a maximum
            tminante=abs(n_j(Jest)); xminante=x(tminante);
            if xminante<xmax
                nc=n1(nmax); flagante=1; Jest=Jest+1;
                n_j(Jest)=floor((n1(nmax)+n2(nmax))/2);
            else % the current maximum is smaller than the anterior minimum
                 % an additional maximum is added ([ATE] p. 82 and 83)
                [xmaxx nmaxx]=max(x(tminante:nmax));
                nmaxx=nmaxx+tminante-1;
                nc=n1(nmaxx); flagante=1; Jest=Jest+1;
                n_j(Jest)=floor((n1(nmaxx)+n2(nmaxx))/2);
                iadd=iadd+1; kadd(iadd)=Jest;
            end
        elseif flagmin==1 % the current extremum is also an minimum
                          %  an additional maximum is added ([ATE] p. 82)
            nc=n1(nmin)-1; flagante=1;
            tminante=abs(n_j(Jest));
            [xmax nmax]=max(x(tminante:nc));
            nmax=nmax+tminante-1;
            Jest=Jest+1; n_j(Jest)=floor((n1(nmax)+n2(nmax))/2);
            iadd=iadd+1; kadd(iadd)=Jest;
        else
            nc=nc+Deltan;
        end
    else % the anterior extremum is an maximum
        if flagmin==1 % the current extremum is a minimum
            tmaxante=abs(n_j(Jest)); xmaxante=x(tmaxante);
            if xmaxante>xmin
                nc=n1(nmin); flagante=-1; Jest=Jest+1;
                n_j(Jest)=-floor((n1(nmin)+n2(nmin))/2);
            else % the current minimum is larger than the anterior maximum
                 % an additional minimum is added ([ATE] p. 82 and 83)
                [xminn nminn]=min(x(tmaxante:nmin));
                nminn=nminn+tmaxante-1;
                nc=n1(nminn); flagante=-1; Jest=Jest+1;
                n_j(Jest)=-floor((n1(nminn)+n2(nminn))/2);
                iadd=iadd+1; kadd(iadd)=Jest;
            end
        elseif flagmax==1% the current extremum is also an maximum
                         %  an additional minimum is added ([ATE] p. 82)
            nc=n1(nmax)-1; flagante=-1;
            tmaxante=abs(n_j(Jest));
            [xmin nmin]=min(x(tmaxante:nc));
            nmin=nmin+tmaxante-1;
            Jest=Jest+1; n_j(Jest)=-floor((n1(nmin)+n2(nmin))/2);
            iadd=iadd+1; kadd(iadd)=Jest;
        else
            nc=nc+Deltan;
        end
    end
end

% x(ni) is not included in the partition of scale Deltan 
nj1=abs(n_j(1));
if nj1>ni
    if n1(nj1)>ni %the boundary ni is not included in the plateau containing
                  %the first local extremum at n_j(1) and it is
                  %added as an additional local extremum ([ATE] p. 83)
        n_j(2:Jest+1)=n_j(1:Jest);
        n_j(1)=-sign(n_j(1))*ni;
        Jest=Jest+1;
        kadd(2:iadd+1)=kadd(1:iadd)+1;
        iadd=iadd+1;
        kadd(1)=1;
    else %the boundary ni is included in the plateau containing
         %the first local extremum at n_j(1) and then the first local
         %extremum is moved at the boundary of the plateau
        n_j(1)=sign(n_j(1))*ni;
    end
end

% the same situation as before but for the other boundary nf
njJ=abs(n_j(Jest));
if njJ<nf
    if n2(njJ)<nf
        Jest=Jest+1;
        n_j(Jest)=-sign(n_j(Jest-1))*nf;
        iadd=iadd+1;
        kadd(iadd)=Jest;
    else
        n_j(Jest)=sign(n_j(Jest))*nf;
    end
end

end