function varargout = myfilter(TYPE,INP,RANGE,ORDER,varargin)
% myfilter  [Not a public function] Generic low/high-pass filter with tunes.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

switch TYPE
    case 'hpf'
        system = @hpfsystem;
        lambda = @hpflambda;
    case 'llf'
        system = @llfsystem;
        lambda = @llflambda;
    case 'bwf'
        system = @bwfsystem;
        % BWF lambda is always supplied from the caller.
end

freq = datfreq(INP.start);
options = passvalopt('tseries.filter',varargin{:});

if ~isempty(options.growth)
    options.change = options.growth;
end

%**************************************************************************

% Default lambda.
if isempty(options.lambda) ...
        || (ischar(options.lambda) && strcmpi(options.lambda,'auto'))
    if freq == 0
        utils.error('tseries', ...
            ['No default lambda for tseries with ', ...
            'indeterminate or daily frequency.']);
    else
        options.lambda = lambda(freq);
    end
end

if any(options.lambda <= 0)
    error('Smoothing parameter must be a positive number.');
end

options.lambda = options.lambda(:).';
options.drift = options.drift(:).';

% Get the user range.
RANGE = specrange(INP,RANGE);
INP = resize(INP,RANGE);
xstart = RANGE(1);
xend = RANGE(end);

% Determine the filtering range.
if ~isempty(options.level) && isa(options.level,'tseries')
    lstart = options.level.start;
    lend = options.level.start + size(options.level.data,1) - 1;
else
    lstart = [];
    lend = [];
end
if ~isempty(options.change) && isa(options.change,'tseries')
    gstart = options.change.start - 1;
    gend = options.change.start + size(options.change.data,1) - 1;
else
    gstart = [];
    gend = [];
end
fstart = min([xstart,lstart,gstart]);
fend = max([xend,lend,gend]);
nper = round(fend - fstart + 1);

% Get input, level and growth data on the filtering range.
xdata = rangedata(INP,[fstart,fend]);
xsize = size(xdata);
xdata = xdata(:,:);
if ~isempty(lstart)
    ldata = rangedata(options.level,[fstart,fend]);
    ldata = ldata(:,:);
end
if ~isempty(gstart)
    gdata = rangedata(options.change,[fstart,fend]);
    gdata = gdata(:,:);
end

if options.log
    xdata = log(xdata);
    if ~isempty(lstart)
        ldata = log(ldata);
    end
    if ~isempty(gstart)
        gdata = log(gdata);
    end
end

nlambda = length(options.lambda);
ndrift = length(options.drift);
nx = size(xdata,2);
if ~isempty(lstart)
    nl = size(ldata,2);
else
    nl = [];
end
if ~isempty(gstart)
    ng = size(gdata,2);
else
    ng = [];
end
nloop = max([nlambda,ndrift,nx,nl,ng]);

tnd = nan(nper,nloop);
gap = nan(nper,nloop);

for iloop = 1 : nloop
    
    if iloop <= nlambda
        lambdai = options.lambda(iloop);
    end
    
    if iloop <= ndrift
        drifti = options.drift(iloop);
    end
    
    if iloop <= nx
        xi = xdata(:,iloop);
    end
    
    [X,B] = system(xi,lambdai,drifti,ORDER);
    
    % Add level constraints.
    if ~isempty(lstart)
        if iloop <= nl
            li = ldata(:,iloop);
        end
        [X,B] = addlevel(X,B,li);
    end
    
    % Add growth constraints.
    if ~isempty(gstart)
        if iloop <= ng
            gi = gdata(:,iloop);
        end
        [X,B] = addgrowth(X,B,gi);
    end
    
    XX = B \ X;
    tnd(:,iloop) = XX(1:nper);
    gap(:,iloop) = xi - tnd(:,iloop);
    
end

if options.log
    tnd = exp(tnd);
    gap = exp(gap);
end

varargout{1} = INP;
varargout{1}.start = fstart;
varargout{1}.data = reshape(tnd,[nper,xsize(2:end)]);
varargout{1} = mytrim(varargout{1});

varargout{2} = INP;
varargout{2}.start = fstart;
varargout{2}.data = reshape(gap,[nper,xsize(2:end)]);
varargout{2} = mytrim(varargout{2});

% Output argument swap requested.
if options.swap
    varargout([1,2]) = varargout([2,1]);
end

end

% Subfunctions.

%**************************************************************************
function [y,B] = addlevel(y,B,ldata)
ldata = ldata(:);
index = ~isnan(ldata).';
if any(index)
    y = [y;ldata(index)];
    for j = find(index)
        B(end+1,j) = 1; %#ok<*AGROW>
        B(j,end+1) = 1;
    end
end
end
% addlevel().

%**************************************************************************
function [y,B] = addgrowth(y,B,gdata)
gdata = gdata(:);
index = ~isnan(gdata).';
if any(index)
    y = [y;gdata(index)];
    for j = find(index)
        B(end+1,[j-1,j]) = [-1,1];
        B([j-1,j],end+1) = [-1;1];
    end
end
end
% addgrowth().

%**************************************************************************
function [y,B] = hpfsystem(y,lambda,drift,order) %#ok<INUSD>
n = size(y,1);
d = [1,-4,6,-4,1];
d = d(ones(1,n),:);
d(1,:) = [1,-2,1,NaN,NaN];
d(2,:) = [1,-4,5,-2,NaN];
d(end,:) = d(1,end:-1:1);
d(end-1,:) = d(2,end:-1:1);
d = d*lambda;
B = spdiags(d,-2:2,n,n);
index = isnan(y);
y(index) = 0;
e = speye(n);
e(index,index) = 0;
B = B + e;
end
% hpfsystem().

%**************************************************************************
function lambda = hpflambda(freq)
lambda = 100*freq^2;
end
% hpflambda().

%**************************************************************************
function [y,B] = llfsystem(y,lambda,drift,order) %#ok<INUSD>
n = size(y,1);
d = [-lambda,2*lambda,-lambda];
d = d(ones(1,n),:);
d(1,2) = lambda;
d(end,2) = lambda;
B = spdiags(d,-1:1,n,n);
index = isnan(y);
y(index) = 0;
y(1) = y(1) - lambda*drift;
y(end) = y(end) + lambda*drift;
e = speye(n);
e(index,index) = 0;
B = B + e;
end
% llfsystem().

%**************************************************************************
function lambda = llflambda(freq)
lambda = 10*freq;
end
% llflambda().

%*************************************************************************
function [y,B] = bwfsystem(y,lambda,drift,order) %#ok<INUSL>
n = size(y,1);
p = order;
x2 = pascalrow(order);
if rem(order,2) == 0
    % Even order.
    tmp = (-1).^(0 : 2*order);
else
    % Odd order.
    tmp = (-1).^(1 : 2*order+1);
end
x2 = x2 .* tmp(ones(1,order+1),:);
d = sum(x2,1);
d = d(ones(1,n),:);
tmp = cumsum(x2,1);
tmp(tmp == 0) = NaN;
d(1:p+1,1:2*p+1) = tmp;
d(end-p:end,end-2*p:end) = tmp(end:-1:1,end:-1:1);
d = d*lambda;
B = spdiags(d,-p:p,n,n);
index = isnan(y);
y(index) = 0;
e = speye(n);
e(index,index) = 0;
B = B + e;
end
% @ bwfsystem().

%**************************************************************************
function x2 = pascalrow(n)
% pascalrow  Get a particular row of the Pascal triangle.
if n == 0
    x2 = 1;
    return
end
% Pascal triangle.
x = [1,1];
for i = 2 : n
    x = sum([x(1:end-1);x(2:end)],1);
    x = [1,x,1];
end
% Row x row.
x2 = x;
for i = 2 : n+1
    x2(i,i:end+1) = x(i)*x;
end
end
% pascalrow().