% LA lawn example: see www.cs.berkeley.edu/~murphyk/Bayes/usage.html for a more detailed explanation.

if 1
  N = 4; 
  dag = zeros(N,N);
  dag(1,2) = 1; dag(1,3) = 1; dag(2,4) = 1; dag(3,4) = 1;
else
  names = {'WetGrass', 'Sprinkler', 'Cloudy', 'Rain'}; 
  topology = {'Cloudy', 'Sprinkler'; 'Cloudy', 'Rain'; 'Sprinkler', 'WetGrass'; 'Rain', 'WetGrass'}; 
  dag = mk_adj_mat(topology, names, 1); 
end

node_sizes = twos(1,N);
bnet = mk_bnet(dag, node_sizes);

for i=1:N
  bnet.bnodes{i} = mk_tabular_node(bnet, i);
end

bnet.bnodes{1}.CPT(:) = [0.5 0.5];
bnet.bnodes{2}.CPT(:) = [0.5 0.9 0.5 0.1];
bnet.bnodes{3}.CPT(:) = [0.8 0.2 0.2 0.8];
bnet.bnodes{4}.CPT(:) = [1 0.1 0.1 0.01 0 0.9 0.9 0.99];

bnet = compile_params(bnet);

obs_nodes = [4];
vals = [2];
[bnet, ll] = enter_evidence(bnet, vals, obs_nodes);
m = marginal_nodes(bnet, 3);
p = m.T(2)
assert(approxeq(p, 0.7079));

obs_nodes = [2 4];
vals = [2 2];
[bnet, ll] = enter_evidence(bnet, vals, obs_nodes);
m = marginal_nodes(bnet, 3);
p = m.T(2) 
assert(approxeq(p, 0.3204));


% Try some learning

nsamples = 50;
samples = sample_bnet(bnet, nsamples);
seed = 0;
bnet2 = randomize_params(bnet, seed);

% hide 20% of the nodes
nnodes = 4;
hide = rand(nsamples, nnodes) > 0.8;
vals = cell(1, nsamples);
obs = cell(1, nsamples);
for i=1:nsamples
  obs{i} = find(~hide(i,:));
  vals{i} = samples(i, obs{i});
end

max_iter = 5;
[bnet3, LL] = learn_params(bnet2, vals, obs, max_iter);

%plot(LL / nsamples,'o-')
%xlabel('EM iteration')
%ylabel('loglik. per sample')
%title('learning curve')
