% This file is the Gaussian analog of discrete2.m

% Make the following network (from Jensen (1996) p84 fig 4.17)
%    1
%  / | \
% 2  3  4
% |  |  |
% 5  6  7
%  \/ \/
%  8   9
% where all arcs point downwards

N = 9;
dag = zeros(N,N);
dag(1,2)=1; dag(1,3)=1; dag(1,4)=1;
dag(2,5)=1; dag(3,6)=1; dag(4,7)=1;
dag(5,8)=1; dag(6,8)=1; dag(6,9)=1; dag(7,9) = 1;

% Make a linear Gaussian network with random params
node_sizes = [5 4 3 2 2 1 2 2 2]; % vector-valued nodes
%node_sizes = ones(1,9); % scalar nodes
discrete_nodes = [];
bnet = mk_bnet(dag, node_sizes, discrete_nodes);
for i=1:N
  bnet.bnodes{i} = mk_gaussian_node(bnet, i, 'diag');
end
bnet = compile_params(bnet);


% Check correctness of inference by comparing with the results of using the equivalent
% global multivariate Gaussian (mvg)

onodes = [5 6];
vals = rand(sum(node_sizes(onodes)),1);
[bnet, ll1] = enter_evidence(bnet, vals, onodes);

m1 = marginal_family(bnet, 3);
clear m2
[m2.mu, m2.Sigma, ll2] = naive_gaussian_marginal(bnet, family(bnet.adj_mat, 3));
assert(approxeq(ll1, ll2));
assert(approxeq(m1.mu, m2.mu));
assert(approxeq(m1.Sigma, m2.Sigma));


% Instantiate all the nodes

mvg = gbn_to_mvg(bnet);
vals = mvg.mu; % + randn(sum(node_sizes),1); % prevent underflow
onodes = 1:N;
[bnet, ll1] = enter_evidence(bnet, vals, onodes);
[m2.mu, m2.Sigma, ll2] = naive_gaussian_marginal(bnet, []);
ll3 = gaussian_prob(vals, mvg.mu, mvg.Sigma, 1);
assert(approxeq(ll1, ll2));
assert(approxeq(ll1, ll3));

i = 7;
fam = family(bnet.adj_mat, i);
m1 = marginal_family(bnet, i);
[m2.mu, m2.Sigma, ll2] = naive_gaussian_marginal(bnet, fam);

assert(approxeq(ll1, ll2));
assert(approxeq(m1.mu, m2.mu));
assert(approxeq(m1.Sigma, m2.Sigma));


