Daniel@0: % Test jtree_compiled on a toy QMR network.
Daniel@0: 
Daniel@0: rand('state', 0);
Daniel@0: randn('state', 0);
Daniel@0: pMax = 0.01;
Daniel@0: Nfindings = 10;
Daniel@0: Ndiseases = 5;
Daniel@0: 
Daniel@0: N=Nfindings+Ndiseases;
Daniel@0: findings = Ndiseases+1:N;
Daniel@0: diseases = 1:Ndiseases;
Daniel@0: 
Daniel@0: G = zeros(Ndiseases, Nfindings);
Daniel@0: for i=1:Nfindings
Daniel@0:   v= rand(1,Ndiseases);
Daniel@0:   rents = find(v<0.8);
Daniel@0:   if (length(rents)==0)
Daniel@0:     rents=ceil(rand(1)*Ndiseases);
Daniel@0:   end
Daniel@0:   G(rents,i)=1;
Daniel@0: end       
Daniel@0: 
Daniel@0: prior = pMax*rand(1,Ndiseases);
Daniel@0: leak = 0.5*rand(1,Nfindings); % in real QMR, leak approx exp(-0.02) = 0.98     
Daniel@0: %leak = ones(1,Nfindings); % turns off leaks, which makes inference much harder
Daniel@0: inhibit = rand(Ndiseases, Nfindings);
Daniel@0: inhibit(not(G)) = 1;
Daniel@0: 
Daniel@0: % first half of findings are +ve, second half -ve
Daniel@0: % The very first and last findings are hidden
Daniel@0: pos = 2:floor(Nfindings/2);
Daniel@0: neg = (pos(end)+1):(Nfindings-1);
Daniel@0: 
Daniel@0: big = 1;
Daniel@0: 
Daniel@0: if big
Daniel@0:   % Make the bnet in the straightforward way
Daniel@0:   tabular_leaves = 1;
Daniel@0:   obs_nodes = myunion(pos, neg) + Ndiseases;
Daniel@0:   bnet = mk_qmr_bnet(G, inhibit, leak, prior, tabular_leaves, obs_nodes);
Daniel@0:   evidence = cell(1, N);
Daniel@0:   evidence(findings(pos)) = num2cell(repmat(2, 1, length(pos)));
Daniel@0:   evidence(findings(neg)) = num2cell(repmat(1, 1, length(neg)));
Daniel@0: else
Daniel@0:   % Marginalize out hidden leaves apriori
Daniel@0:   positive_leaves_only = 1;
Daniel@0:   [bnet, vals] = mk_minimal_qmr_bnet(G, inhibit, leak, prior, pos, neg, positive_leaves_only);
Daniel@0:   obs_nodes = bnet.observed;
Daniel@0:   evidence = cell(1, Ndiseases + length(obs_nodes));
Daniel@0:   evidence(obs_nodes) = num2cell(vals);
Daniel@0: end
Daniel@0: 
Daniel@0: engine = {};
Daniel@0: engine{end+1} = jtree_inf_engine(bnet);
Daniel@0: 
Daniel@0: E = length(engine);
Daniel@0: exact = 1:E;
Daniel@0: ll = zeros(1,E);
Daniel@0: for e=1:E
Daniel@0:   tic; [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence); toc
Daniel@0: end
Daniel@0: 
Daniel@0: assert(all(approxeq(ll(exact), ll(exact(1)))))
Daniel@0: 
Daniel@0: post = zeros(E, Ndiseases);
Daniel@0: for e=1:E
Daniel@0:   for i=diseases(:)'
Daniel@0:     m = marginal_nodes(engine{e}, i);
Daniel@0:     post(e, i) = m.T(2);
Daniel@0:   end
Daniel@0: end
Daniel@0: for e=exact(:)'
Daniel@0:   for i=diseases(:)'
Daniel@0:     assert(approxeq(post(1, i), post(e, i)));
Daniel@0:   end
Daniel@0: end
Daniel@0: