annotate toolboxes/FullBNT-1.0.7/bnt/examples/static/cmp_inference_static.m @ 0:e9a9cd732c1e tip

first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [time, engine] = cmp_inference_static(bnet, engine, varargin)
wolffd@0 2 % CMP_INFERENCE Compare several inference engines on a BN
wolffd@0 3 % function [time, engine] = cmp_inference_static(bnet, engine, ...)
wolffd@0 4 %
wolffd@0 5 % engine{i} is the i'th inference engine.
wolffd@0 6 % time(e) = elapsed time for doing inference with engine e
wolffd@0 7 %
wolffd@0 8 % The list below gives optional arguments [default value in brackets].
wolffd@0 9 %
wolffd@0 10 % exact - specifies which engines do exact inference [ 1:length(engine) ]
wolffd@0 11 % singletons_only - if 1, we only call marginal_nodes, else this and marginal_family [0]
wolffd@0 12 % maximize - 1 means we do max-propagation, 0 means sum-propagation [0]
wolffd@0 13 % check_ll - 1 means we check that the log-likelihoods are correct [1]
wolffd@0 14 % observed - list of the observed ndoes [ bnet.observed ]
wolffd@0 15 % check_converged - list of loopy engines that should be checked for convergence [ [] ]
wolffd@0 16 % If an engine has converged, it is added to the exact list.
wolffd@0 17
wolffd@0 18
wolffd@0 19 % set default params
wolffd@0 20 exact = 1:length(engine);
wolffd@0 21 singletons_only = 0;
wolffd@0 22 maximize = 0;
wolffd@0 23 check_ll = 1;
wolffd@0 24 observed = bnet.observed;
wolffd@0 25 check_converged = [];
wolffd@0 26
wolffd@0 27 args = varargin;
wolffd@0 28 nargs = length(args);
wolffd@0 29 for i=1:2:nargs
wolffd@0 30 switch args{i},
wolffd@0 31 case 'exact', exact = args{i+1};
wolffd@0 32 case 'singletons_only', singletons_only = args{i+1};
wolffd@0 33 case 'maximize', maximize = args{i+1};
wolffd@0 34 case 'check_ll', check_ll = args{i+1};
wolffd@0 35 case 'observed', observed = args{i+1};
wolffd@0 36 case 'check_converged', check_converged = args{i+1};
wolffd@0 37 otherwise,
wolffd@0 38 error(['unrecognized argument ' args{i}])
wolffd@0 39 end
wolffd@0 40 end
wolffd@0 41
wolffd@0 42 E = length(engine);
wolffd@0 43 ref = exact(1); % reference
wolffd@0 44
wolffd@0 45 N = length(bnet.dag);
wolffd@0 46 ev = sample_bnet(bnet);
wolffd@0 47 evidence = cell(1,N);
wolffd@0 48 evidence(observed) = ev(observed);
wolffd@0 49 %celldisp(evidence(observed))
wolffd@0 50
wolffd@0 51 for i=1:E
wolffd@0 52 tic;
wolffd@0 53 if check_ll
wolffd@0 54 [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence, 'maximize', maximize);
wolffd@0 55 else
wolffd@0 56 engine{i} = enter_evidence(engine{i}, evidence, 'maximize', maximize);
wolffd@0 57 end
wolffd@0 58 time(i)=toc;
wolffd@0 59 end
wolffd@0 60
wolffd@0 61 for i=check_converged(:)'
wolffd@0 62 niter = loopy_converged(engine{i});
wolffd@0 63 if niter > 0
wolffd@0 64 fprintf('loopy engine %d converged in %d iterations\n', i, niter);
wolffd@0 65 % exact = myunion(exact, i);
wolffd@0 66 else
wolffd@0 67 fprintf('loopy engine %d has not converged\n', i);
wolffd@0 68 end
wolffd@0 69 end
wolffd@0 70
wolffd@0 71 cmp = exact(2:end);
wolffd@0 72 if check_ll
wolffd@0 73 for i=cmp(:)'
wolffd@0 74 assert(approxeq(ll(ref), ll(i)));
wolffd@0 75 end
wolffd@0 76 end
wolffd@0 77
wolffd@0 78 hnodes = mysetdiff(1:N, observed);
wolffd@0 79
wolffd@0 80 if ~singletons_only
wolffd@0 81 get_marginals(engine, hnodes, exact, 0);
wolffd@0 82 end
wolffd@0 83 get_marginals(engine, hnodes, exact, 1);
wolffd@0 84
wolffd@0 85 %%%%%%%%%%
wolffd@0 86
wolffd@0 87 function get_marginals(engine, hnodes, exact, singletons)
wolffd@0 88
wolffd@0 89 bnet = bnet_from_engine(engine{1});
wolffd@0 90 N = length(bnet.dag);
wolffd@0 91 cnodes_bitv = zeros(1,N);
wolffd@0 92 cnodes_bitv(bnet.cnodes) = 1;
wolffd@0 93 ref = exact(1); % reference
wolffd@0 94 cmp = exact(2:end);
wolffd@0 95 E = length(engine);
wolffd@0 96
wolffd@0 97 for n=hnodes(:)'
wolffd@0 98 for e=1:E
wolffd@0 99 if singletons
wolffd@0 100 m{e} = marginal_nodes(engine{e}, n);
wolffd@0 101 else
wolffd@0 102 m{e} = marginal_family(engine{e}, n);
wolffd@0 103 end
wolffd@0 104 end
wolffd@0 105 for e=cmp(:)'
wolffd@0 106 if cnodes_bitv(n)
wolffd@0 107 assert(approxeq(m{ref}.mu, m{e}.mu))
wolffd@0 108 assert(approxeq(m{ref}.Sigma, m{e}.Sigma))
wolffd@0 109 else
wolffd@0 110 assert(approxeq(m{ref}.T, m{e}.T))
wolffd@0 111 end
wolffd@0 112 assert(isequal(m{e}.domain, m{ref}.domain));
wolffd@0 113 end
wolffd@0 114 end