function [R,V,sE_R,trendidx,tfilt,mcreal,paddedX,noise]=ssatrend(x,M,sE_x,varargin)
% Finds the ssa non-linear trend from a time series.
% 
%
% USAGE:
%     [R,V,sE_R,trendidx,tfilt]=ssatrend(x,M[,sE_x,parameterlist])
%
% INPUT:
%   x: the time series
%   M: embedding dimension (=maximum lag), (or a symmetric filter in the time domain, it will not be using ssa then)
%   SE_x: standard error of a measurement in x (assumend to be independent)
%         - can also be replaced by a function handle to a function that
%           generates noise-surrogates of the same length as x. E.g.
%           E.g: @()randn(size(x))*.6 is equivalent to .6
%
% PARAMETERS:
%   SSAMethod: see help on ssa
%   ExtrapMethod: MinimumRoughness(=Default) or MinimumSlope
%
% OUTPUT:
%   R: the nonlinear trend (a specific reconstructed component where
%      x has succesively padded with the linear trend of the last M points. For
%      details check source).
%      note that the mean(x) has been added to the RC.
%      (so that it is easy to plot on the same plot as the original
%      timeseries)
%   V: the variance explained by R (as a fraction of the total variance)
%   SE_R: standard error of the R (how much of sE_x is remaining in R).
%   trendidx: the index of the rc corresponding to the trend
%   tfilt: symmetric reconstructive SSA filter. (what is convoluted on x)
%
% if no output is specified then the function will make a plot where 2
% standard errors (2*sE_R) is shaded grey.
%
% (Requires ssa.m by Eric Breitenberger and formatts.m from my WTC toolbox)
%
% Note: the ssatrend filter very often results in virtually the same output
% as a simple triangular FIR filter of length 2*M-1.
%
% (c) Aslak Grinsted 2009
% 2009: added boostrapping estimates of padding errors. - please note the
% algorithm at present assumes white noise uncertainties of the trend.
%

Args=struct('SSAMethod','unbiased',...
    'ExtrapMethod','minimumroughness');
Args=parseArgs(varargin,Args);



if nargin<3, sE_x=[]; end

if nargin==0
    warning('No data specified, using debug data...')
    x=cumsum(randn(300,1)*.2)+randn(300,1);
    M=100;
    sE_x=1;
    Args.ExtrapMethod='mr';
end

x=formatts(x);
t=x(:,1);
x=x(:,2);
meanx=mean(x);
if ~isa(sE_x,'function_handle') %then make it into a function-handle
    sE_x=sE_x(:);
    %     if length(sE_x)==1
    %         sE_x=ones(size(x))*sE_x;
    %     end
    
    if length(sE_x)==1
        sE_x=@()randn(size(x))*sE_x;
    else
        sE_x=@()randn(size(x)).*sE_x;
    end
end

noerrorbar=isempty(sE_x);
if noerrorbar
    sE_x=zeros(size(x,1),1);
end


if mean(sE_x())>std(x)*1.01
    warning('sE_x>std(x)!!!')
end

n=length(x);

if numel(M)==1
    
    [E,V,A,R,C]=ssa(x,M,Args.SSAMethod);
    
    
    
    %find trend idx:
    [mx,trendidx]=max(abs(sum(E))./sum(abs(E))); %trend would be non oscillatory and therefore have a large sum in the index
    
    
    %throw away everything but the trend:
    E=E(:,trendidx);
    V=V(trendidx)/sum(V); %convert to percent
    
%     if V<.1
%         warning('ssa trend accounts for less than 10% of the variance.')
%     end
    R=R(:,trendidx)+meanx;
    tfilt=conv(E,flipud(E))/M;
else
    tfilt=M(:);
    V=nan;trendidx=nan;
    M=ceil(length(tfilt)/2); %used for padding.
end
%tfilt=tfilt/sum(tfilt); %NORMALIZE!!!!!!! (NO!)



%succesively pad the series with values from the trend and reconstruct
%again to get better end estimates.
% paddedX=x;
% for ii=1:M
%     leftR=mean(paddedX(1:M));
%     rightR=mean(paddedX(end+1-(1:M)));
%     paddedX=[leftR;paddedX;rightR];
% end
mp=M-1;
switch Args.ExtrapMethod
    case {'minimumroughness','minrough','mr'}
        idx=(1:mp)';
        
        %pre calculate all mp-windowed linear fits for the bootstrap. The
        %prediction residuals will be used to calculate the uncertainty due
        %to padding.
        pt=[idx ones(mp,1)];
        ptleft=[idx-mp ones(mp,1)];
        ptright=[idx+mp ones(mp,1)];
        residleft=nan(mp,n-mp*2);
        residright=nan(mp,n-mp*2);
        for ii=1:n-mp
            %pp(:,ii)=pt\x(idx+ii-1);
            pp=pt\x(idx+ii-1);
            if ii>mp
                residleft(:,ii-mp)=ptleft*pp-x(idx+ii-mp-1);
            end
            if ii<=n-mp*2
                residright(:,ii)=ptright*pp-x(idx+ii+mp-1);
            end
        end
        
        pleft=pt\x(idx);
        pright=pt\x(idx+n-mp);
        %
        paddedX=[ptleft*pleft;x;ptright*pright];
        %
        
        
    case {'minimumslope','minslope','ms'}
        idx=(1:mp)';
        
        residleft=nan(mp,n-mp*2);
        residright=nan(mp,n-mp*2);
        for ii=1:n-mp
            pp=mean(x(idx+ii-1));
            if ii>mp
                residleft(:,ii-mp)=pp-x(idx+ii-mp-1);
            end
            if ii<=n-mp*2
                residright(:,ii)=pp-x(idx+ii+mp-1);
            end
        end
        
        pleft=zeros(mp,1)+mean(x(idx));
        pright=zeros(mp,1)+mean(x(idx+n-mp));
        %
        paddedX=[pleft;x;pright];
        

end


if (mp*7>n)
%     warning('The lag is quite big. That makes it difficult for the boostrapping procedure. Assuming symmetry of residual errors to get more bootstrap data...')
%     
    residleft=[residleft flipud(residright)];
    residleft=[residleft -residleft];
    residright=flipud(residleft); 
end


% R=filter(tfilt,1,paddedX);
% R=R((1:n)+M*2-1);
R=filter2(tfilt,paddedX,'valid');

%What happens to whitenoise if it goes through this tfilt?
if noerrorbar
    sE_R=[];
else
    %pad sE_x
    %    sE_x=[zeros(mp,1)+epleft;sE_x;zeros(mp,1)+epright];
    
    %MONTE CARLO WAY
    mcreal=nan(length(x),2000);
    
    for ii=1:size(mcreal,2)
        noise=[residleft(:,ceil(rand*size(residleft,2)));sE_x();residright(:,ceil(rand*size(residright,2)))];
        mcreal(:,ii)=filter2(tfilt,paddedX+noise,'valid');
    end
    sE_R=std(mcreal')';
end

if (nargout==0)
    if ~isempty(sE_R)
        CI=prctile(mcreal',[5 95])';
        h=fill([t;flipud(t)],[CI(:,1);flipud(CI(:,2))],[.6 .6 .8]);
        set(h,'linestyle','none')
        hold on
    end
    hh=plot(t,x,'b',t,R,'k');
    set(hh(2),'linewidth',2);
    set(hh(1),'color',[.8 .8 1]);
    hold off
end


if nargout<1, clear R,end
if nargout<2, clear V,end
if nargout<3, clear sE_R,end
if nargout<4, clear trendidx,end
if nargout<5, clear tfilt,end
