function [trend_rcma,varargout] =trendrcma(x0,flagacd,dn0)
%TRENDRCMA Automatically estimate the nonmonotonic trend of a noisy time series
%trend_rcma=trendrcma(x0) estimates the nonmonotonic trend of the noisy
%       time series x0 using the smoothing algorithm AutRCMA presented in
%       Section 7.1 of the book [ATE] C. Vamos and M. Craciun, Automatic Trend
%       Estimation, Springer 2012.
%[trend_rcma,pest] =trendrcma(x0) in addition supplies the estimated
%       parameters of the AutRCMA and of the time series x0, the rcma trend and the residuals
%       pest=Structure with the following fields
%       pest.dn = the minimum distance between two succesive local extrema of
%            the estimated trend
%       pest.K = semilength of the averaging window
%       pest.nrep = number of repetitions of the CMA
%       pest.fi=the estimated value of noise correlation parameter 
%[trend_rcma,pest,trend_acd]=trendrcma(x0,1) in addition, by applying the ACD
%       algorithm (function 'trendacd'), it computes a trend composed only
%       by the monotonic parts of the trendrcma from which the unsignificant local
%       extrema are removed and graphically represents the iterative process
%       (see Section 7.2 in [ATE])
%[...]=trendrcma(x0,flag,dn0) dn0 is the minimum distance between two extrema of
%       the estimated trend in case where the characteristic time scale of 
%       the phenomenum is known (see [ATE], sect.6.4)
        

N=length(x0);
minx0=min(x0); maxx0=max(x0);

if nargin==1 
    flagacd=0;
end

%estimation of the serial correlation parameter \phi from Eq. (2.6) in [ATE]
normdif=norm(x0(2:N)-x0(1:N-1))^2;
normdif2=norm(x0(3:N)-x0(1:N-2))^2;
fiest=(N-1)/(N-2)*normdif2/normdif-1;
if fiest>0.99 fiest=0.99; end;
%estimation of the ratio r between trend and noise amplitudes from
%Eq. (F.1) in Appendix F of [ATE]
nj=abs(local_extrema(x0));
E=length(nj); %the number of extrema of the time series
dnj=nj(2:E)-nj(1:E-1); kdnj=find(dnj<11);
amplm=abs(x0(nj(kdnj+1))-x0(nj(kdnj)))/(maxx0-minx0);
rest=1/max(amplm)-1;
rest=max(rest-0.2573,0);
% determination of the optimal parameters of the RCMA from table 7.1
% ([ATE] p. 102) for the estimated values fiest and rest
    
   
if nargin<3
    K=10; dn=20;
    if fiest>0.6 K=20; dn=40; end
    if rest<1.5 K=K+10; dn=dn+40; end
    if N<500 dn=ceil(dn*N/1000); end%%% ar fi mai bine N<1000????
else
    dn=dn0;
    if dn<=20 K=10;
    elseif 20<dn && dn<40 K=20;
    else K=30;
    end
end
%%for short series the semilength of the averanging window must decrease
if N<500
        K=ceil(K*N/1000); end

% the trend is estimated by AutRCMA algorithm described in Sect. 7.1 [ATE]
dn_rcma=1; irep=0; trend_rcma=x0;
while dn_rcma<=dn
    irep=irep+1;
    if irep>1 trend_rcma=rcma(trend_rcma,K,1); end
    nj_rcma=abs(local_extrema(trend_rcma));
    E_rcma=length(nj_rcma); % the number of extrema of the rcma trend
    dn_rcma=min(nj_rcma(2:E_rcma)-nj_rcma(1:E_rcma-1));
end

z_rcma=x0-trend_rcma;
%computation of the correlation parameter for the estimated noise
zmod=z_rcma-mean(z_rcma);
sigz=std(zmod); acvf=zmod(1:N-1)*zmod(2:N)';
fi_rcma=acvf/(sigz^2*(N-1));
if nargout>=2
  pest.dn = dn;
  pest.K = K;
  pest.nrep=irep-1;
  pest.fi=fiest;
  pest.r=rest;
  varargout={pest};
 else
   pest=[];
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimation of the significant local extrema of the rcma trend which limit
%the monotonic segments of the estimated trend (see Section 7.2 in [ATE])
if flagacd==1
    
    % we compute the average on 30 surrogate time series of the fluctuations
    % magnitude of the smoothed estimated noise (p. 95 [ATE])
    % the AR(1) surrogates are generated by the algorithm described in Sect. 1.3 of [ATE]
    sig_s=std(z_rcma,1)*sqrt(1-fi_rcma^2); % Eq. (1.22) in [ATE]
    zeta=zeros(1,N);
    for R=1:30
        Z=sig_s*randn(1,N);
        zeta(1)=Z(1)/sig_s;
        for n=2:N, zeta(n)=fi_rcma*zeta(n-1)+Z(n); end
        zeta=zeta-mean(zeta);
        zeta_rcma=rcma(zeta,K,irep);
        Azeta(R)=max(zeta_rcma)-min(zeta_rcma); % p. 95 in [ATE]
    end
    meanAzeta=mean(Azeta);


    J_acd=1; %initially there is a single monotonic segment of the rcma trend
    n1acd=1; n2acd=N; % initially the entire rcma trend is considered monotonic
    Nmin=14; %the minimum number of values for which the acd algorithm can be applied
    flagm=0; % index for the monotony of the rcma fragments; if flagm(j)=0 then the
             % segment does not contain significant local extrema
    trend_acd=trendacd(trend_rcma);
    Aj=(max(trend_rcma-trend_acd)-min(trend_rcma-trend_acd)); % the amplitude of the
                  % deviations from monotony of the trend_rcma
    theta=zeros(1,E_rcma); % the significance of the local extrema of the rcma trend
    % defined at p. 106 of [ATE] equal with the nonmonotony index given by eq. (6.3) in [ATE]
    i=0; %the number of iterations
    while 1
    % the rcma trend is iteratively partitioned using its local extrema
    % which are not related to smoothed noise fluctuations and which limit
    % the monotonic segments
        JJ_acd=0; nn1acd=[]; nn2acd=[]; %initiation of the parameters describing
        AAj=[]; flagmm=[];              %the next partition of the rcma trend
        for jacd=1:J_acd
            if flagm(jacd)==1
                % the jacd segment was previously recognized as monotonic
                % and it is included in the new partition without modification
                JJ_acd=JJ_acd+1;
                flagmm(JJ_acd)=1;
                nn1acd(JJ_acd)=n1acd(jacd);
                nn2acd(JJ_acd)=n2acd(jacd);
                AAj(JJ_acd)=Aj(jacd);
            else
                %determination of the local extrema with the largest time
                %scales within the jacd segment
                Deltan=n2acd(jacd)-n1acd(jacd)+1; Jfrag=2;
                while Jfrag==2 && Deltan>dn;
                    [nfrag,nadd]=local_extrema(trend_rcma,Deltan,n1acd(jacd),n2acd(jacd));
                    Jfrag=length(nfrag);
                    nfrag=abs(nfrag);
                    Deltan=floor(Deltan/2);
                end
                if Jfrag==2
                    %the jacd segment has no local extrema with time scale
                    %smaller than dn and it is added to the acd trend
                    JJ_acd=JJ_acd+1;
                    flagmm(JJ_acd)=1;
                    nn1acd(JJ_acd)=n1acd(jacd);
                    nn2acd(JJ_acd)=n2acd(jacd);
                    AAj(JJ_acd)=Aj(jacd);
                    xfrag=trend_rcma(n1acd(jacd):n2acd(jacd));
                    trend_acd(n1acd(jacd):n2acd(jacd))=trendacd(xfrag);
                    kredus1=find(nj_rcma==n1acd(jacd));
                    kredus2=find(nj_rcma==n2acd(jacd));
                    if theta(kredus1)==0
                        theta(kredus1)=theta(kredus2); end
                    if theta(kredus2)==0
                        theta(kredus2)=theta(kredus1); end
                else
                    %the Jfrag-1 fragments of the jacd segment are analyzed
                    %to see if they contain other significant local
                    %extrema, otherwise they are added to the partition of the
                    %rcma trend
                    for jfrag=1:Jfrag-1
                        NN=nfrag(jfrag+1)-nfrag(jfrag)+1;
                        xfrag=trend_rcma(nfrag(jfrag):nfrag(jfrag+1));
                        tend=trendacd(xfrag);
                        trend_acd(nfrag(jfrag):nfrag(jfrag+1))=tend;
                        JJ_acd=JJ_acd+1;
                        nn1acd(JJ_acd)=nfrag(jfrag);
                        nn2acd(JJ_acd)=nfrag(jfrag+1);
                        AAj(JJ_acd)=(max(xfrag-tend)-min(xfrag-tend));
                        if NN<Nmin  || AAj(JJ_acd)/meanAzeta<1
                            flagmm(JJ_acd)=1;
                        else flagmm(JJ_acd)=0;
                        end
                    end
                    %the values theta are not atributed to the additional
                    %local extrema of the partition
                    abstextr_redus=nfrag(setdiff([1:Jfrag],nadd));
                    Nredus=length(abstextr_redus);
                    for nredus=1:Nredus
                        kredus=find(nj_rcma==abstextr_redus(nredus));
                        if theta(kredus)==0
                            theta(kredus)=Aj(jacd)/meanAzeta; end
                    end
                end
            end
        end
        i=i+1;
        trend_graph{i}=trend_acd;
        n1graph{i}=nn1acd; n2graph{i}=nn2acd;
        theta_graph{i}=theta;
        AAgraph{i}=AAj;
        JJgraph{i}=JJ_acd;
        if length(find(flagmm==0))==0 
            break;
        else
            J_acd=JJ_acd; flagm=flagmm; Aj=AAj;
            n1acd=nn1acd; n2acd=nn2acd;
        end
    end

    %%%%graphical representation of the iterative construction of the acd
    %%%%trend; the last panel in each figure is repeated as the first panel
    %%%%of the next figure
    nsp=3; %the number of subplots in a figure
    dx=(maxx0-minx0)/20;
    figure;
    set(gcf, 'Units', 'centimeters','PaperOrientation', 'landscape','Position', [1,1,25,15])
    subplot(nsp,1,1);hold all;
    plot(1:N,x0,'.g');plot(1:N,trend_rcma,'r');
    ylabel('x_n');
    xlim([-N/50,51*N/50]); ylim([minx0-dx,maxx0+dx])
    title('estimated AutRCMA trend')
    iteration=1; %the iteration number of the partition of the rcma trend
    ipanel=1; %the panel of the figure under construction
    flagf=0; %flag for the final panel
    while iteration<=i  && flagf==0
        nsr=mod(ipanel,nsp)+1;
        if nsr-1==0
            figure;
            set(gcf, 'Units', 'centimeters','PaperOrientation', 'landscape','Position', [1,1,25,15])
        end
        subplot(nsp,1,nsr); hold
        plot(trend_rcma,'r');
        plot(trend_graph{iteration},'--b');
        ming=min([trend_rcma,trend_graph{iteration}]);
        maxg=max([trend_rcma,trend_graph{iteration}]);
        dg=(maxg-ming)/20;
        xlim([-N/50,51*N/50]); ylim([ming-dg,maxg+dg])
        for jacd=1:JJgraph{iteration}
            x1=n1graph{iteration}(jacd);
            x2=n2graph{iteration}(jacd);
            plot([x1 x1],[ming+3*dg maxg],'k');
            plot([x2 x2],[ming+3*dg maxg],'k');
            delta=(x2-x1)/5;
            chi=AAgraph{iteration}(jacd)/meanAzeta;
            text((x1+x2)/2-delta,maxg-dg,['\chi=',num2str(chi,'%3.1f')]);
            kfl=find(nj_rcma==x1);
            teta=theta_graph{iteration}(kfl);
            if teta~=0
                if kfl==1
                    text(x1,ming+dg,['\theta=',num2str(teta,'%3.1f')]);
                else
                    text(x1,ming+dg,['\theta=',num2str(teta,'%3.1f')]);
                end
            end
            if jacd==JJgraph{iteration}
                tetaf=theta_graph{iteration}(E_rcma);
                if tetaf~=0
                    text(N,ming+dg,['\theta=',num2str(tetaf,'%3.1f')]);
                end
            end
        end
        ylabel('f_n'); if iteration==nsp xlabel('n'); end
        box on
        if nsr<nsp
            iteration=iteration+1;
        end
        ipanel=ipanel+1;
        if iteration>i || (iteration==i && nsr==nsp)
            flagf=1;end
    end
end

if nargout==3
   varargout={pest,trend_acd};  
end


end

%%% the repeated central moving average (RCMA) with averaging window
%%% semilength K and irep repetitions (Section 4.1 of [ATE])
function xm=rcma(x,K,irep)
N=length(x);
T=2*K+1;
for inrep=1:irep
    tablou=zeros(T,N); numitor=zeros(1,N);
    numitor(1:K)=[K+1:2*K];
    numitor(N-K+1:N)=[2*K:-1:K+1];
    numitor(K+1:N-K)=T*ones(1,N-2*K);
    for t=1:K+1 tablou(t,1:N+t-K-1)=x(K+2-t:N); end
    for t=1:K tablou(K+1+t,t+1:N)=x(1:N-t); end
    xm=sum(tablou)./numitor;
    x=xm;
end

end