function [R,V,sE_R,trendidx,tfilt]=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
%
% 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)
%
% (c) Aslak Grinsted 2004 
%

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

if nargin<3, sE_x=[]; end


x=formatts(x);
t=x(:,1);
x=x(:,2);
meanx=mean(x);
sE_x=sE_x(:);
if length(sE_x)==1
    sE_x=ones(size(x))*sE_x;
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!!!!!!!



%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;
switch Args.ExtrapMethod
    case {'minimumroughness','minrough','mr'}
        idx=(1:mp)';
        pleft=polyfit(idx,x(idx),1);
        epleft=sqrt(max(sE_x(idx))^2+std(detrend(x(idx)))^2);
        
        idx=(n-mp+(1:mp)');
        pright=polyfit(idx,x(idx),1);
        epright=sqrt(max(sE_x(idx))^2+std(detrend(x(idx)))^2);
        
        paddedX=[polyval(pleft,(-(M-1):0)');x;polyval(pright,n+(1:M)')];
    
    case {'minimumslope','minslope','ms'}
        idx=(1:mp)';
        pleft=mean(x(idx));
        epleft=sqrt(max(sE_x(idx))^2+std(x(idx))^2);
        idx=(n-mp+(1:mp)');
        pright=mean(x(idx));
        epright=sqrt(max(sE_x(idx))^2+std(x(idx))^2);
        paddedX=[pleft+zeros(size(idx));x;pright+zeros(size(idx))];
end


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

%What happens to whitenoise if it goes through this tfilt?
if noerrorbar
    sE_R=[];
else
    %pad sE_x
    sE_x=[zeros(M,1)+epleft;sE_x;zeros(M,1)+epright];
    if 1
        %MONTE CARLO WAY
        mcreal=zeros(length(x),1000);
        for ii=1:size(mcreal,2)
            temp=conv(paddedX+randn(size(paddedX)).*sE_x,tfilt);
            mcreal(:,ii)=temp(floor((length(temp)-n)/2)+(1:n));
        end
        sE_R=std(mcreal')';
    else
        et=randn(100000,1); %100000 might be overkill but conv is so fast anyway.
        et=conv(et,tfilt);
        sE_R=std(et);

        sE_R=conv(sE_x,tfilt/sum(tfilt))*sE_R;

        sE_R=sE_R(floor((length(sE_R)-n)/2)+(1:n));
    end
end

if (nargout==0)
    if ~isempty(sE_R)
        h=fill([t;flipud(t)],[R+sE_R*2;flipud(R-sE_R*2)],[.85 .85 .85]);
        set(h,'linestyle','none')
        hold on
    end
    plot(t,R,'k',t,x,':')
    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
