function [THIS,OUTP,V,DELTA,PE] = filter(THIS,INP,RANGE,varargin)
% filter  Kalman smoother and estimator of out-of-likelihood parameters.
%
% Syntax
% =======
%
%     [M,OUTP,V,DELTA,PE] = filter(M,INP,RANGE,...)
%     [M,OUTP,V,DELTA,PE] = filter(M,INP,RANGE,J,...)
%
% Input arguments
% ================
%
% * `M` [ model ] - Solved model object.
%
% * `INP` [ struct | cell ] - Input database or datapack from which the
% measurement variables will be taken.
%
% * `RANGE` [ numeric ] - Filter date range.
%
% * `J` [ struct ] - Database with tunes on the mean of shocks and/or
% time-varying std devs of shocks.
%
% Output arguments
% =================
%
% * `M` [ model ] - Model object with updates of std devs (if `'relative='`
% is true) and/or updates of out-of-likelihood parameters (if `'outoflik='`
% is non-empty).
%
% * `OUTP` [ struct | cell ] - Output struct with smoother or prediction
% data.
%
% * `V` [ numeric ] - Estimated variance scale factor if the `'relative='`
% options is true; otherwise `V` is 1.
%
% * `DELTA` [ struct ] - Database with estimates of out-of-likelihood
% parameters.
%
% * `PE` [ struct ] - Database with prediction errors for measurement
% variables.
%
% Options
% ========
%
% * `'ahead='` [ numeric | *`1`* ] - Predictions will be computed this number
% of period ahead.
%
% * `'deviation='` [ `true` | *`false`* ] - Treat input and output data as
% deviations from balanced-growth path.
%
% * `'dtrends='` [ *'auto'* | `true` | `false` ] - Measurement data contain
% deterministic trends.
%
% * `'data='` [ `'predict'` | *`'smooth'`* | `'predict,smooth'` ] - Return
% smoother data or prediction data or both.
%
% * `'exclude='` [ cellstr | *empty* ] - List of measurement variables that
% will be excluded from the likelihood function.
%
% * `'initCond='` [ `'fixed'` | `'optimal'` | *`'stochastic'`* | struct ] - Method
% or data to initialise the Kalman filter.
%
% * `'lastSmooth='` [ numeric | *`Inf`* ] - Last date up to which to smooth
% data backward from the end of the range; if `Inf` smoother will run on the
% entire range.
%
% * `'meanOnly='` [ `true` | *`false`* ] - Return a plain database with mean
% data only.
%
% * `'outOfLik='` [ cellstr | empty ] - List of parameters in deterministic
% trends that will be estimated by concentrating them out of the likelihood
% function.
%
% * `'objective='` [ *`'-loglik'`* | `'prederr'` ] - Objective function
% computed; can be either minus the log likelihood function or weighted sum
% of prediction errors.
%
% * `'objRange='` [ numeric | *`Inf`* ] - The objective function will
% be computed on this subrange only; `Inf` means the entire filter range.
%
% * `'precision='` [ *`'double'`* | `'single'` ] - Numeric precision to which
% output data will be stored; all calculations themselves always run to
% double precision.
%
% * `'rollback='` [ numeric | *empty* ] - Date up to which to roll back
% individual observations on measurement variables from the end of the
% sample.
%
% * `'relative='` [ *`true`* | `false` ] - Std devs of shocks assigned in the
% model object will be treated as relative std devs, and a common variance
% scale factor will be estimated.
%
% * `'returnMse='` [ *`true`* | `false` ] - Return MSE matrices for
% predetermined state variables; these can be used for settin up initial
% condition in subsequent call to another `filter` or `jforecast`.
%
% * `'returnStd='` [ *`true`* | `false` ] - Return database with std devs
% of model variables.
%
% * `'tolMse='` [ numeric | *`0`* ] - Tolerance under which two MSE
% matrices in two consecutive periods will be treated as equal, and the
% Kalman gain system will be re-used, not re-computed.
%
% * `'weighting='` [ numeric | *empty* ] - Weighting vector or matrix for
% prediction errors when `'objective=' 'prederr'`; empty means prediction
% errors are weighted equally.
%
% Options for models with non-linearised equations
% =================================================
%
% * `'nonlinearise='` [ numeric | *`0`* ] - If non-zero the prediction step
% in the Kalman filter will be run in an exact non-linear mode using the
% same technique as [`model/simulate`](model/simulate).
%
% * `'simulate='` [ cell | empty ] - Options passed in to `simulate` when
% invoking the non-linear simulation in the prediction step; only used when
% `nonlinearise=` is greater than `0`.
%
% Description
% ============
%
% The `'ahead='` and `'rollback='` options cannot be combined with one
% another, or with multiple data sets, or with multiple parameterisations.
%
% Example
% ========
%

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

isstrfind = @(a,b) ~isempty(strfind(a,b));
NARGOUT = nargout;

% Database with tunes.
j = [];
if ~isempty(varargin) && (isstruct(varargin{1}) || isempty(varargin{1}))
    j = varargin{1};
    varargin(1) = [];
end

P = inputParser();
P.addRequired('model',@ismodel);
P.addRequired('data',@(x) isstruct(x) || iscell(x) || isempty(x));
P.addRequired('range',@isnumeric);
P.addRequired('tune',@(x) isempty(x) || isstruct(x) || iscell(x));
P.parse(THIS,INP,RANGE,j);

% This FILTER function options.
[opt,varargin] = passvalopt('model.filter',varargin{:});

% Process Kalman filter options; `mypreloglik` also expands solution
% forward if needed for tunes on the mean of shocks.
[kalmanopt,THIS] = mypreloglik(THIS,RANGE,'t',j,varargin{:});
kalmanopt.ahead = opt.ahead;

% Get observables and determine the format of output data.
INP = datarequest('y',THIS,INP,RANGE);
ndata = size(INP,3);
nalt = size(THIS.Assign,3);

% Check option conflicts.
dochkconflicts();

% Create additional data sets for rollback.
dorollback();

% Throw a warning if some of the data sets have no observations.
nandata = all(all(isnan(INP),1),2);
if any(nandata)
    utils.warning('model', ...
        'No observations available in input data set(s)%s.', ...
        sprintf(' #%g',(nandata)));
end

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

ny = sum(THIS.nametype == 1);

% Requested output arguments.
Y = struct();
dorequestoutp();

% Run the Kalman filter.
[obj,V,delta,Y] = mykalman(THIS,INP,Y,kalmanopt); %#ok<ASGLU>

% Assign out-of-lik params in model object, and create out-of-lik parameter
% database.
DELTA = struct();
if ~isempty(kalmanopt.outoflik)
    for i = 1 : length(kalmanopt.outoflik)
        namepos = kalmanopt.outoflik(i);
        name = THIS.name{namepos};
        THIS.Assign(1,namepos,:) = delta(1,i,:);
        DELTA.(name) = permute(delta(1,i,:),[1,3,2]);
    end
end

% Scale std devs by estimated factor in model object.
if kalmanopt.relative
    ne = sum(THIS.nametype == 3);
    se = sqrt(V);
    for ialt = 1 : nalt
        THIS.stdcorr(1,1:ne,ialt) = ...
            THIS.stdcorr(1,1:ne,ialt)*se(ialt);
    end
end

% No need to refresh links after assigning out-of-liks and re-scaling
% stdevs because neither of them are allowed in dynamic links.

xrange = RANGE(1)-1 : RANGE(end);
[Y,~,PE,~,pred,smooth] = mykalmanoutput(THIS,Y,xrange);

if isfield(Y,'predmean') && isfield(Y,'smoothmean')
    OUTP.pred = pred;
    OUTP.smooth = smooth;
elseif isfield(Y,'predmean')
    OUTP = pred;
elseif isfield(Y,'smoothmean')
    OUTP = smooth;
end

% Nested functions.

%**************************************************************************
    function dochkconflicts()
        if kalmanopt.ahead > 1 && (ndata > 1 || nalt > 1)
            utils.error('model', ...
                ['Cannot run the option ''ahead'' with ', ...
                'multiple data sets or multiple parameterisations.']);
        end
        if ~isempty(opt.rollback) ...
                && (ndata > 1 || nalt > 1 || kalmanopt.ahead > 1)
            utils.error('model', ...
                ['Cannot run the option ''rollback='' with ', ...
                'multiple data sets, multiple parameterisations, or ''ahead=''.']);
        end
        if ~isempty(kalmanopt.outoflik) && (ndata ~= nalt)
            utils.error('model', ...
                ['When the option ''outoflik='' is used, ', ...
                'the number of parameterisations must match ', ...
                'the number of data sets.']);
        end
    end
% dochkconflicts().

%**************************************************************************
    function dorollback()
        if ~isempty(opt.rollback)
            index = round(opt.rollback(:).' - RANGE(1)) + 1;
            opt.rollback = index(index >= 1 & index <= nper);
            opt.rollback = sort(opt.rollback,'descend');
            % Create additional sets of observables for rollbacks.
            if ~isempty(opt.rollback)
                for ii = opt.rollback
                    adddata = Y(:,:,end*ones(1,ny));
                    for jj = 1 : ny
                        adddata(ny+1-j:ny,ii,jj) = NaN;
                    end
                    Y = cat(3,Y,adddata);
                end
                ndata = size(Y,3);
            end
        end
    end
% dorollback().

%**************************************************************************
    function dorequestoutp()
        Y.delta = [];
        if NARGOUT >= 5
            Y.pe = [];
        end
        if NARGOUT >= 2
            if isstrfind(opt.data,'pred')
                Y.predmean = [];
                if ~kalmanopt.meanonly
                    if kalmanopt.returnstd
                        Y.predstd = [];
                    end
                    if kalmanopt.returnmse
                        Y.predmse = [];
                    end
                end
            end
            if isstrfind(opt.data,'smooth')
                Y.smoothmean = [];
                if ~kalmanopt.meanonly
                    if kalmanopt.returnstd
                        Y.smoothstd = [];
                    end
                    if kalmanopt.returnmse
                        Y.smoothmse = [];
                    end
                end                
            end
        end
    end
% dorequestoutp().

end