function [bnet, ll] = propagate_lds(bnet, data)
% PROPOGATE_LDS Call the Kalman smoother to do inference in the Linear Dynamical System.
% [bnet, ll] = propogate_lds(bnet, vals)
% 
% We compute the posterior distribution on the hidden variables in each slice and each pair of
% slices, and store in bnet.hmm.one_slice_marginal{t}.moment and bnet.hmm.two_slice_marginal{t}.moment respectively.

T = size(data,2);

if bnet.filter_only
  [ll, xhat, Vhat, VVhat] = kalman_filter(data, bnet.hmm.trans_mat, bnet.hmm.obs_mat, ...
    bnet.hmm.trans_cov, bnet.hmm.obs_cov, bnet.hmm.init_state, bnet.hmm.init_cov);
else
  [xhat, Vhat, VVhat, ll] = kalman_smoother(data, bnet.hmm.trans_mat, bnet.hmm.obs_mat, ...
      bnet.hmm.trans_cov, bnet.hmm.obs_cov, bnet.hmm.init_state, bnet.hmm.init_cov);
end

A = bnet.hmm.trans_mat;
Q = bnet.hmm.trans_cov;
ns = bnet.node_sizes_slice;
hs = sum(ns(bnet.hidden_slice));
os = sum(ns(bnet.observed_slice));
ss = bnet.slice_size;

if length(bnet.hmm.one_slice_marginal) ~= T
  bnet = mk_space_for_hmm_marginals(bnet, 0);
end

for t=1:T
  % The following is when we include the observable in one_slice_counts
  %mom.mu = [xhat(:,t) ; bnet.values(block(2+(t-1)*ss, bs))];
  %mom.Sigma = [Vhat(:,:,t) zeros(hs,os); zeros(os,hs) zeros(os,os)];
  % The following is when we don't include the observable in one_slice_counts
  bnet.hmm.one_slice_marginal{t}.mu = xhat(:,t);
  bnet.hmm.one_slice_marginal{t}.Sigma = Vhat(:,:,t);
  bnet.hmm.one_slice_marginal{t}.logp = ll;
  
  if t>1
    x1 = xhat(:,t-1);
    x2 = xhat(:,t);
    mu = [x1; x2];
    V1 = Vhat(:,:,t-1);
    V2 = Vhat(:,:,t);
    V12 = VVhat(:,:,t);
    %Sigma = [V1  V1*A'+Q; A*V1+Q  A*V1*A'+V2+Q];
    Sigma = [V1 V12'; V12 V2];
    logp = ll;
    bnet.hmm.two_slice_marginal{t-1}.mu = mu;
    bnet.hmm.two_slice_marginal{t-1}.Sigma = Sigma;
    bnet.hmm.two_slice_marginal{t-1}.logp = ll;
  end
end

