function [X,Px,e,u,Y,Py,yindex] = ...
    varsmoother(A,B,K,Z,D,Omg,Sgm,y,e,x0,P0,invfunc,allobserved, ...
    tolerance,reuse)
% smoothervar  [Not a public function] Kalman smoother for VAR-based systems.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright 2007-2012 Jaromir Benes.

% The VAR-based state-space system is given by
%
% y = Z x + D + u,   Euu' = Sgm
% x = A(L) x(-1) + K + B e,   Eee' = Omg
%
% If allobserved == true, then a full set of observations y(t) exactly
% determines x(t).

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

persistent REUSE;

if ~exist('tolerance','var')
    tolerance = 0;
end

if ~exist('reuse','var')
    reuse = false;
end

nx = size(A,1);
[ny,nper] = size(y);
p = size(A,2)/nx;

if isempty(D)
    D = zeros(ny,1);
end

symm = @(X) (X+X.')/2;

T = [A;eye((p-1)*nx,p*nx)];
X = nan(nx*p,nper);
Xf = nan(nx*p,nper);
Px = zeros(nx*p,nx*p,nper);
Y = y;
Py = nan(ny,ny,nper);
F = nan(ny,ny,nper);
Fi = nan(ny,ny,nper);
pe = nan(ny,nper);
Fipe = nan(ny,nper);
L = nan(nx*p,nx*p,nper);
G = nan(nx*p,ny,nper);

if isempty(e)
    if isempty(B)
        ne = ny;
    else
        ne = size(B,2);
    end
    e = zeros(ne,nper);
end

if isempty(B)
    Be = e;
    BOmg = Omg;
    BOmgBt = symm(Omg);
else
    Be = B*e;
    BOmg = B*Omg;
    BOmgBt = symm(BOmg*B');
end

deviation = isempty(K);

issgm = any(any(Sgm ~= 0));
if allobserved
    PP = BOmgBt;
    FF = Z*PP*Z';
    if issgm
        FF = FF + Sgm;
    end
    FFi = invfunc(FF);
    FF = symm(FF);
    FFi = symm(FFi);
end

% First prediction step.
nonzero = any(A ~= 0,1);
if any(nonzero)
    X(1:nx,1) = A(:,nonzero)*x0(nonzero);
else
    X(1:nx,1) = 0;
end
X(nx+1:end,1) = x0(1:end-nx);
if ~deviation
    X(1:nx,1) = X(1:nx,1) + K;
end
if all(P0 == 0)
    Px(1:nx,1:nx,1) = BOmgBt;
else
    Px(:,:,1) = T*P0*T';
    Px(1:nx,1:nx,1) = Px(1:nx,1:nx,1) + BOmgBt;
end
Px(:,:,1) = symm(Px(:,:,1));

j = false(ny,1);
yindex = ~isnan(y);
nz2 = size(Z,2);
for t = 1 : nper
    j0 = j;
    j = yindex(:,t);
    pe(j,t) = y(j,t) - Z(j,:)*X(1:nz2,t) - D(j,:);
    if reuse && ~isempty(REUSE)
        F(:,:,t) = REUSE.F(:,:,t);
        Fi(:,:,t) = REUSE.Fi(:,:,t);
    else
        % allobserved == true indicates that all states are perfectly observed and
        % FMSE is static whenever all observations are available in this and the
        % previous period.
        if allobserved && all(j0) && all(j)
            F(:,:,t) = FF;
            Fi(:,:,t) = FFi;
        else
            F(:,:,t) = Z*Px(1:nz2,1:nz2,t)*Z.';
            if issgm
                F(:,:,t) = F(:,:,t) + Sgm;
            end
            F(:,:,t) = symm(F(:,:,t));
            recompute = true;
            if t > 1 && all(j0 == j)
                f = F(j,j,t);
                f0 = F(j,j,t-1);
                if  max(abs(f(:) - f0(:))) < tolerance
                    Fi(j,j,t) = Fi(j,j,t-1);
                    recompute = false;
                end
            end
            if recompute
                Fi(j,j,t) = invfunc(F(j,j,t));
                Fi(j,j,t) = symm(Fi(j,j,t));
            end
        end
    end
    Fipe(j,t) = F(j,j,t)\pe(j,t);
    G(:,j,t) = T*Px(:,1:nz2,t)*Z(j,:).'/F(j,j,t);
    Xf(:,t) = X(:,t) + Px(:,1:nz2,t)*Z(j,:).'*Fipe(j,t);
    % L = T - G*[Z,0]
    L(:,:,t) = T;
    L(:,1:nz2,t) = L(:,1:nz2,t) - G(:,j,t)*Z(j,:);
    if t < nper
        X(1:nx,t+1) = A(:,nonzero)*X(nonzero,t) + Be(:,t+1);
        if ~deviation
            X(1:nx,t+1) = X(1:nx,t+1) + K;
        end
        X(nx+1:end,t+1) = X(1:end-nx,t);
        if any(j)
            X(:,t+1) = X(:,t+1) + G(:,j,t)*pe(j,t);
        end
        % Px = T*Px*L' + B*Omg*B';
        Px(:,:,t+1) = T*Px(:,:,t)*L(:,:,t)';
        if allobserved && all(j)
            % All states are perfectly observed, hence t+1 forecast MSE is
            % just PP = B*Omg*B'.
            Px(1:nx,1:nx,t+1) = PP;
        else
            Px(1:nx,1:nx,t+1) = Px(1:nx,1:nx,t+1) + BOmgBt;
        end
        Px(:,:,t+1) = symm(Px(:,:,t+1));
    end
end

if reuse
    if isempty(REUSE)
        REUSE = struct();
        REUSE.F = F;
        REUSE.Fi = Fi;
    end
else
    REUSE = [];
end

lastobs = max([0,find(any(yindex,1),1,'last')]);
if lastobs < nper
    Y(:,lastobs+1:end) = Z*X(1:nz2,lastobs+1:end);
    Py(:,:,lastobs+1:end) = F(:,:,lastobs+1:end);
end

r = zeros(p*nx,1);
N = 0;

% Any measurement errors?
isu = issgm && size(Sgm,1) == size(y,1);
if isu
    u = zeros(ny,nper);
end

for t = lastobs : -1 : 1
    j = yindex(:,t);
    if isu
        u(:,t) = Sgm(j,:)'*(Fipe(j,t) - G(:,j,t)'*r);
    end
    % r = [Z,0]'*Fi*pe + L'*r;
    r = L(:,:,t)'*r;
    r(1:nz2) = r(1:nz2) + Z(j,:)'*Fipe(j,t);
    X(:,t) = X(:,t) + Px(:,:,t)*r;
    e(:,t) = e(:,t) + BOmg'*r(1:nx);
    if any(~j)
        Y(~j,t) = Z(~j,:)*X(1:nz2,t);
        if isu
            Y(~j,t) = Y(~j,t) + u(~j,t);
        end
    end
    % N = [Z,0]'*Fi*[Z,0] + L'*N*L;
    N = L(:,:,t)'*N*L(:,:,t);
    N(1:nz2,1:nz2) = N(1:nz2,1:nz2) + Z(j,:).'*Fi(j,j,t)*Z(j,:);
    PxNPx = Px(:,:,t)*N*Px(:,:,t);
    Px(:,:,t) = symm(Px(:,:,t) - PxNPx);
    Py(:,:,t) = symm(F(:,:,t) - Z*PxNPx(1:nz2,1:nz2)*Z');
    Py(j,:,t) = 0;
    Py(:,j,t) = 0;
    if allobserved && all(j)
        Px(1:nx,1:nx,t) = 0;
    end
end

end