function [mll,C,L] = kalman(this,X,C,options)
% kalman  Kalman filter for bkwmodel objects.
%
% Backend IRIS function.
% No help provided.

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

maxlag = size(this.occur,3) - 1;
nper = size(X,2);
n = sum(this.nametype == 1);
ne = sum(this.nametype == 3);
nn = length(this.systemid);
first = maxlag + 1;
islog = this.log(1:n);

F = nan([n,n,nper]);
Ts = nan([nn,nn,nper]);
Rs = nan([nn,ne,nper]);

P0 = zeros([nn,nn,nper]);
P1 = zeros([nn,nn,nper]);

X0 = X;
X(2,first:nper) = NaN;
x0 = X;
x1 = X;
y1 = X(this.nametype == 1,:);
x0(this.nametype == 1,first:nper) = NaN;
pe = nan([n,nper]);
pe1 = pe;
Z = eye([n,nn]);
a0 = nan([nn,nper]);
a1 = nan([nn,nper]);

%{
if any(this.log)
   x0(this.log,:) = log(x0(this.log,:));
   x1(this.log,:) = log(x1(this.log,:));
   y1(islog,:) = log(y1(islog,:));
end
%}

logdetF = 0;
peFipe = 0;
nobs = 0;
exitflag = NaN;
L = [];
mll = 1e10;

for t = first : nper

   % Prediction step.
   if false % t > first+1
      x0(1:n,t) = a0(1:n,t-1) + T(1:n,:)*(a1(:,t-1) - a1(:,t-2));
      %if ~options.approxpred
      %     x1(1:n,t) = x0(1:n,t);
      %     x0(1:n,t) = predict_(x1,t);
      %end
   else
      x0(1:n,t) = predict_(x1,t);
      if exitflag <= 0
         mll = 1e10;
         return
      end
   end

   temp = x1;
   temp(:,t) = x0(:,t);
   [A,B,T,R] = system(this,temp,t,false);
   if rank(A) < size(A,1)
      mll = 2e10;
      return
   end
   
   Ts(:,:,t) = T;
   Rs(:,:,t) = R;
   Tt = T.';
   Rt = R.';
   P0(:,:,t) = T*P1(:,:,t-1)*Tt + R*C*Rt;
   j = ~isnan(y1(:,t)).';
   F(j,j,t) = P0(j,j,t);
   if rank(F(j,j,t)) < sum(j)
      mll = 3e10;
      return
   end
   
   % Prediction error.
   pe(j,t) = y1(j,t) - x0(j,t);
   pe1(j,t) = y1(j,t) - x1(j,t);
   
   % Updating step.
   ZtFi = Z(j,:).'/F(j,j,t);
   Fipe = F(j,j,t) \ pe(j,t);
   K = P0(:,:,t)*ZtFi;
   
   % Window of x0 from t-maxlag to t.
   xw0 = x0(this.nametype == 1,t-maxlag:t);
   a0(:,t) = xw0(this.systemxi);
   a1(:,t) = a0(:,t) + K*pe(j,t);
   xw1 = xw0;
   xw1(this.systemxi) = a1(:,t);
   x1(this.nametype == 1,t-maxlag:t) = xw1;
   P1(:,:,t) = (eye(nn) - K*Z(j,:))*P0(:,:,t);
   P1(j,j,t) = 0;
   
   if any(j)
      logdetF = logdetF + log(det(F(j,j,t)));
      peFipe = pe(j,t).'*Fipe;
      nobs = nobs + sum(j);
   end
   
end

%{
%if options.relative
   v = peFipe / nobs;
   peFipe = nobs;
   logdetF = logdetF + nobs*log(v);
%end
%}

mll = (logdetF + peFipe)/2;
%{
      x = X0;
      res = this.eqtnEvalAll(x,first:nper);
      W = sum(diag(C\res(j,:)*res(j,:).'));
      v1 = W / nobs;
      W = nobs;
      logdetC = (nper-maxlag)*log(det(C));
      logdetC = logdetC + nobs*log(v1);
      L = (logdetC + W)/2;
      
      for j = first : nper
         A = system(this,x,j,false);
         A = A(1:4,1:4);
         A(3,:) = [];
         A(:,2) = [];
         % A = A(~exclude,~exclude);
         logdetA = log(det(A));
         L = L - logdetA;
      end

%}

%if nargout > 2
%     C = C*v;
%end

% Nested functions follow.

   %***********************************************************************
   % Nested function.
   function u1 = predict_(x,t)
      x(this.log,:) = log(x(this.log,:));
      if t > first + 1
         u0 = x(this.nametype == 1,t-1); 
      else
         u0 = x(this.nametype == 1,t-1);
      end
      
      % u0(isnan(u0)) = 0;
      %if isequal(options.solver,'lsqnonlin')
      %     [u1,resnorm,ans,exitflag] = ...
      %            lsqnonlin(@objfunc_,u0,[],[],options.optimset);
      %else
         [u1,fval,exitflag] = ...
                  fsolve(@objfunc_,u0,options.optimset);
      u1(islog,:) = exp(u1(islog,:));
      %end

      %***********************************************************************
      % Nested^2 function.
      function obj = objfunc_(u)
         xx = x;
         xx(this.nametype == 1,t) = u;
         if any(this.log)
            xx(this.log,:) = exp(xx(this.log,:));
         end
         obj = this.eqtnEvalAll(xx,t);
         obj = obj(:);
      end
      % End of nested function objfunc_().

   end
   % End of nested function predict_().


end