% Check correctness of inference in a Factorial HMM with Gaussian output.
% Requires the HMM toolbox.

M = 2; % num. chains
K = 2; % binary hidden variables
O = 4; % size of observed vector
ss = M+1; % slice size
obs = M+1; % the index of the observable node
[intra, inter] = mk_fhmm_topology(M);
node_sizes = [K*ones(1,M) O];
discrete_nodes = 1:M;
equiv = [1:M (M+1)];
equiv_class = [equiv equiv(1:end-1)+ss obs];
bnet = mk_dbn(intra, inter, node_sizes, discrete_nodes, equiv_class, obs, 'bk', 'exact');
for i=1:M
  bnet.bnodes{i} = mk_tabular_node(bnet, i);
end
bnet.bnodes{obs} = mk_gaussian_node(bnet, obs);
for i=[1:M] + ss
  bnet.bnodes{i} = mk_tabular_node(bnet, i);
end
bnet = randomize_params(bnet, 0);
% so randn doesn't produce wildly improbable sequences ...
bnet.bnodes{obs}.cov = repmat(eye(O), [1 1 K*ones(1,M)]);

bnet = compile_params(bnet);

% Make a random observation sequence
T = 5; % num slices
data = randn(O, T);


% Check inference

[bnet, llj] = enter_evidence(bnet, data);

Mu = reshape(bnet.bnodes{obs}.mean, [O M^K]);
Cov = reshape(bnet.bnodes{obs}.cov, [O O M^K]);
Pi = cell(1,M);
trans = cell(1,M);
for i=1:M
  Pi{i} = bnet.bnodes{i}.CPT;
  trans{i} = bnet.bnodes{i+ss}.CPT;
end
obsmat = mk_ghmm_obs_mat(data, Mu, Cov);

[ll2, gamma] = fhmm_infer(bnet.inter, Pi, trans, obsmat, node_sizes(1:end-1));

t = 2;
m = marginal_family(bnet, obs, t);
assert(approxeq(m.T(:), gamma(:,t)));
