function [LOGLIK,VAR,PE,A0,A1,A2] = kalman(THIS,A,X,sy,O,OPT)
% kalman  [Not a public function] Kalman filter for rhsmodel objects.
%
% Backend IRIS function.
% No help provided.

retlikonly = nargout <= 3;
dosmooth = nargout > 3;

% TODO: Return k-step-ahead predictions.

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

ndata = size(X,3);
if OPT.ahead > ndata
    X = cat(3,X,nan(size(X,1),size(X,2),OPT.ahead-ndata));
end

maxlag = size(THIS.occur,3) - 1;
nper = size(X,2);
ny = sum(THIS.nametype == 1);
ne = sum(THIS.nametype == 3);
na = length(THIS.systemid);
first = maxlag + 1;
ylog = THIS.log(1:ny);

% Update parameters.
p = THIS.Assign(1,THIS.nametype == 4,:).';
X(THIS.nametype == 4,:,:) = p(:,ones(1,nper),ones(1,ndata));
Omg = covfun.stdcorr2cov(THIS.stdcorr,ne);

Y1 = X(1:ny,:,:);

A1 = A(:,:,1);
A0 = nan(na,nper);

P1 = zeros(na,na,nper);
P0 = nan(na,na,nper);
if dosmooth
    P2 = nan(na,na,nper);
    A2 = nan(na,nper);
end

PE = nan(ny,nper,size(Y1,3));
Z = eye(ny,na);
Zt = Z.';

logdetF = 0;
peFipe = 0;
nobs = 0;

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

% Pre-allocate matrices re-used in smoothing.
C = cell(1,nper);
ZtFi = cell(1,nper);
Zk = cell(1,nper);
F = cell(1,nper);
pe = cell(1,nper);
T = cell(1,nper);
R = cell(1,nper);

for t = first : nper

    y1 = Y1(:,t,:);
    y1 = y1(:,:);
    sy1 = sy(:,t,:);
    sy1 = sy1(:,:);
    lastpred = find(any(~isnan(y1),1),1,'last');
    ahead = max([OPT.ahead,lastpred]);
    y1 = y1(:,1:ahead);
    sy1 = sy1(:,1:ahead);
    ylog1 = ylog(:);
    ylog1 = ylog1(:,ones(1,ahead));
    y1(ylog1) = log(y1(ylog1));
    sy1(ylog1 & sy1 ~= 0) = log(sy1(ylog1 & sy1 ~= 0));

    % Prediction step.

    a0 = [A1(:,1:t-1),nan(na,ahead)];
    for k = 1 : ahead
        [a0,res,flag] = prediction(THIS,a0,X,t-1+k,O);
        if flag < 1
            do_crash();
            return
        end
    end
    
    A0(:,t) = a0(:,t);
    A1(:,t) = a0(:,t);
    
    % Re-arrange predictions for Kalman calculations.
    y0 = a0(1:ny,t-1+(1:ahead));
    y0 = y0(:,:);

    % Updating step.
    
    % Linearise the system around the current state.
    [T{t},R{t}] = system(THIS,A1,X,t);

    S = R{t}*Omg*R{t}.';

    P0(:,:,t) = T{t}*P1(:,:,t-1)*T{t}.' + S;
    P0(:,:,t) = symmetric(P0(:,:,t));

    % Create matrix `Zk` so that
    % y0(:) = Zk*a0(:,t) + independent terms.
    Zk{t} = [];
    % ZTp{i} := Z*T^(i-1).
    ZTp = cell(1,ahead);
    ZTp{1} = Z;
    for k = 1 : ahead
        Zk{t} = [Zk{t};ZTp{k}];
        if k < ahead
            ZTp{k+1} = ZTp{k}*T{t};
        end
    end

    % Create MSE matrix for y0(:).
    F{t} = [];
    F1 = cell(ahead,ahead);
    % Pk is MSE of a(t-1+k|t-1).
    Pk = P0(:,:,t);
    for k = 1 : ahead
        ZPk = Z*Pk;
        F1{k,k} = ZPk*Zt + diag(sy1(:,k).^2);
        for m = 1 : ahead-k
            F1{k,k+m} = ZPk*ZTp{m+1}.';
            F1{k+m,k} = F1{k,k+m}.';
        end
        F{t} = [F{t};F1{k,:}];
        if k < ahead
            Pk = T{t}*Pk*T{t}.' + S;
        end
    end
    F{t} = symmetric(F{t});

    % Reduce `Zk` and `F` down to available observations.
    jk = ~isnan(y1);
    j1 = jk(1:ny);
    Zk{t} = Zk{t}(jk(:),:);
    F{t} = F{t}(jk(:),jk(:));

    % Prediction error.
    pe{t} = y1(jk) - y0(jk);
    % PE(:,t,1:ahead) = reshape(y1k - y0k,[ny,1,ahead]);

    if ~any(jk(:))
        % Note that A1 has been already updated with a0 for us to be able
        % to compute time t system matrices.
        P1(:,:,t) = P0(:,:,t);
        ZtFi{t} = zeros(na,0);
        C{t} = eye(na);
        continue
    end

    % Updating step.
    ZtFi{t} = Zk{t}.' / F{t};
    K = P0(:,:,t)*ZtFi{t};

    C{t} = eye(na) - K*Zk{t};
    A1(:,t) = A0(:,t) + K*pe{t};
    P1(:,:,t) = C{t}*P0(:,:,t);
    P1(j1,:,t) = 0;
    P1(:,j1,t) = 0;

    logdetF = logdetF + log(det(F{t}));
    peFipe = peFipe + pe{t}.'*(F{t}\pe{t});

    nobs = nobs + sum(jk(:));
    
    if imag(logdetF) ~= 0 || ~isfinite(logdetF)
        do_crash();
        return
    end

end

VAR = 1;
if OPT.relative
    VAR = peFipe / nobs;
    peFipe = nobs;
    logdetF = logdetF + nobs*log(VAR);
end

LOGLIK = (logdetF + peFipe)/2;

if retlikonly
    return
end

if dosmooth
    % Run and return smoothing.
    P2(:,:,nper) = P1(:,:,nper);
    A2(:,nper) = A1(:,nper);
    M = zeros(na);
    m = zeros(na,1);
    for t = nper : -1 : first
        if t < nper
            P2(:,:,t) = P1(:,:,t) - P1(:,:,t)*M*P1(:,:,t);
            A2(:,t) = A1(:,t) - P1(:,:,t)*m;
        end
        if t > first
            % N{t}, n{t}.
            N = ZtFi{t}*Zk{t} + C{t}.'*M*C{t};
            n = -ZtFi{t}*pe{t} + C{t}.'*m;
            % M{t-1}, m{t-1}.
            M = T{t}.'*N*T{t};
            m = T{t}.'*n;
        end
    end
end

    function do_crash()
        LOGLIK = Inf;
        VAR = NaN;        
    end

end