% Compare results with those on p46 of
% "A Bayesian Decision-Theoretic Framework for Real-Time Monitoring and Diagnosis of Complex
% Systems: Theory and Application", S. Alag, 1996. U.C. Berkeley, Dept. Mech. Eng. PhD thesis.
%
%
%  x1   x2
%    \ /
%     x3
%   /   \
%  x4    x5


N=5;
A=zeros(N,N);
A(1,3) = 1; A(2,3) = 1; A(3,4) = 1; A(3, 5) = 1;
node_sizes = [2 1 2 1 2];
bnet = mk_bnet(A, node_sizes, []);
for i=1:N
  bnet.bnodes{i} = mk_gaussian_node(bnet, i, 'full');
end

bnet.bnodes{1}.mean = [1 0]';
bnet.bnodes{1}.cov = [4 1; 1 4];

bnet.bnodes{2}.mean = 1;
bnet.bnodes{2}.cov = 1;

bnet.bnodes{3}.mean = [0 0]';
bnet.bnodes{3}.cov = [2 1; 1 1];
bnet.bnodes{3}.weights = [ [1 2; 1 0] [2 ; 1] ]; % block matrix

bnet.bnodes{4}.mean = 0;
bnet.bnodes{4}.cov = 1;
bnet.bnodes{4}.weights = [1 1];

bnet.bnodes{5}.mean = [0 0]';
bnet.bnodes{5}.cov = [1 0; 0 1];
bnet.bnodes{5}.weights = [1 0; 1 1];

bnet = compile_params(bnet);


% Now check the BN gives the same results as his thesis and a brute force method.

% No evidence
vals = [];
bnet = enter_evidence(bnet, vals);

s1 = [31 39; 39 55];
m2 = marginal_nodes(bnet, 5);
[m3.mu, m3.Sigma] = naive_gaussian_marginal(bnet, 5);
assert(approxeq(s1, m2.Sigma));
assert(approxeq(s1, m3.Sigma));


% Diagnostic estimation
obs_nodes = 5;
vals = [5 5];
bnet = enter_evidence(bnet, vals, obs_nodes);

s1 = [1.8261 -0.1957; -0.1957 1.0924];
m2 = marginal_nodes(bnet, 1);
[m3.mu, m3.Sigma] = naive_gaussian_marginal(bnet, 1);
assert(approxeq(s1, m2.Sigma));
assert(approxeq(s1, m3.Sigma));
s1 = [0.3478 1.1413]';
assert(approxeq(s1, m2.mu));
assert(approxeq(s1, m3.mu));


% Predictive and diagnostic estimation
obs_nodes = [2 4 5];
vals = [0 10 5 5];
bnet = enter_evidence(bnet, vals, obs_nodes);
m2 = marginal_nodes(bnet, 1);
[m3.mu, m3.Sigma] = naive_gaussian_marginal(bnet, 1);
assert(approxeq(m2.mu, m3.mu));
assert(approxeq(m2.Sigma, m3.Sigma));






