wolffd@0: % Make a QMR-like network 
wolffd@0: % This is a bipartite graph, where the top layer contains hidden disease nodes,
wolffd@0: % and the bottom later contains observed finding nodes.
wolffd@0: % The diseases have Bernoulli CPDs, the findings noisy-or CPDs.
wolffd@0: % See quickscore_inf_engine for references.
wolffd@0: 
wolffd@0: pMax = 0.01;
wolffd@0: Nfindings = 10;
wolffd@0: Ndiseases = 5;
wolffd@0: %Nfindings = 20;
wolffd@0: %Ndiseases = 10;
wolffd@0: 
wolffd@0: N=Nfindings+Ndiseases;
wolffd@0: findings = Ndiseases+1:N;
wolffd@0: diseases = 1:Ndiseases;
wolffd@0: 
wolffd@0: G = zeros(Ndiseases, Nfindings);
wolffd@0: for i=1:Nfindings
wolffd@0:   v= rand(1,Ndiseases);
wolffd@0:   rents = find(v<0.8);
wolffd@0:   if (length(rents)==0)
wolffd@0:     rents=ceil(rand(1)*Ndiseases);
wolffd@0:   end
wolffd@0:   G(rents,i)=1;
wolffd@0: end       
wolffd@0: 
wolffd@0: prior = pMax*rand(1,Ndiseases);
wolffd@0: leak = 0.5*rand(1,Nfindings); % in real QMR, leak approx exp(-0.02) = 0.98     
wolffd@0: %leak = ones(1,Nfindings); % turns off leaks, which makes inference much harder
wolffd@0: inhibit = rand(Ndiseases, Nfindings);
wolffd@0: inhibit(not(G)) = 1;
wolffd@0: 
wolffd@0: 
wolffd@0: % first half of findings are +ve, second half -ve
wolffd@0: % The very first and last findings are hidden
wolffd@0: pos = 2:floor(Nfindings/2);
wolffd@0: neg = (pos(end)+1):(Nfindings-1);
wolffd@0: 
wolffd@0: % Make the bnet in the straightforward way
wolffd@0: tabular_leaves = 0;
wolffd@0: obs_nodes = myunion(pos, neg) + Ndiseases;
wolffd@0: big_bnet = mk_qmr_bnet(G, inhibit, leak, prior, tabular_leaves, obs_nodes);
wolffd@0: big_evidence = cell(1, N);
wolffd@0: big_evidence(findings(pos)) = num2cell(repmat(2, 1, length(pos)));
wolffd@0: big_evidence(findings(neg)) = num2cell(repmat(1, 1, length(neg)));
wolffd@0: 
wolffd@0: %clf;draw_layout(big_bnet.dag);
wolffd@0: %filename = '../public_html/Bayes/Figures/qmr.rnd.jpg';
wolffd@0: %% 3x3 inches
wolffd@0: %set(gcf,'units','inches');
wolffd@0: %set(gcf,'PaperPosition',[0 0 3 3])  
wolffd@0: %print(gcf,'-djpeg','-r100',filename);
wolffd@0: 
wolffd@0: 
wolffd@0: % Marginalize out hidden leaves apriori
wolffd@0: positive_leaves_only = 1;
wolffd@0: [bnet, vals] = mk_minimal_qmr_bnet(G, inhibit, leak, prior, pos, neg, positive_leaves_only);
wolffd@0: obs_nodes = bnet.observed;
wolffd@0: evidence = cell(1, Ndiseases + length(obs_nodes));
wolffd@0: evidence(obs_nodes) = num2cell(vals);
wolffd@0: 
wolffd@0: 
wolffd@0: clear engine;
wolffd@0: engine{1} = quickscore_inf_engine(inhibit, leak, prior);
wolffd@0: engine{2} = jtree_inf_engine(big_bnet);
wolffd@0: engine{3} = jtree_inf_engine(bnet);
wolffd@0: 
wolffd@0: %fname = '/home/cs/murphyk/matlab/Misc/loopybel.txt';
wolffd@0: global BNT_HOME
wolffd@0: fname = sprintf('%s/loopybel.txt', BNT_HOME);
wolffd@0: 
wolffd@0: 
wolffd@0: max_iter = 6;
wolffd@0: engine{4} = pearl_inf_engine(bnet, 'protocol', 'parallel', 'max_iter', max_iter);
wolffd@0: %engine{5} = belprop_inf_engine(bnet, 'max_iter', max_iter, 'filename', fname);
wolffd@0: engine{5} = belprop_inf_engine(bnet, 'max_iter', max_iter);
wolffd@0: 
wolffd@0: E = length(engine);
wolffd@0: exact = 1:3;
wolffd@0: loopy = [4 5];
wolffd@0: 
wolffd@0: ll = zeros(1,E);
wolffd@0: tic; engine{1} = enter_evidence(engine{1}, pos, neg); toc
wolffd@0: tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, big_evidence); toc
wolffd@0: tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
wolffd@0: tic; [engine{4}, ll(4), niter(4)] = enter_evidence(engine{4}, evidence); toc
wolffd@0: tic; [engine{5}, niter(5)] = enter_evidence(engine{5}, evidence); toc
wolffd@0: 
wolffd@0: ll
wolffd@0: 
wolffd@0: post = zeros(E, Ndiseases);
wolffd@0: for e=1:E
wolffd@0:   for i=diseases(:)'
wolffd@0:     m = marginal_nodes(engine{e}, i);
wolffd@0:     post(e, i) = m.T(2);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: for e=exact(:)'
wolffd@0:   for i=diseases(:)'
wolffd@0:     assert(approxeq(post(1, i), post(e, i)));
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: a = zeros(Ndiseases, 2);
wolffd@0: for ei=1:length(loopy)
wolffd@0:   for i=diseases(:)'
wolffd@0:     a(i,ei) = approxeq(post(1, i), post(loopy(ei), i));
wolffd@0:   end
wolffd@0: end
wolffd@0: disp('is the loopy posterior correct?');
wolffd@0: disp(a)