wolffd@0: function [time, engine] = cmp_inference_static(bnet, engine, varargin)
wolffd@0: % CMP_INFERENCE Compare several inference engines on a BN
wolffd@0: % function [time, engine] = cmp_inference_static(bnet, engine, ...)
wolffd@0: %
wolffd@0: % engine{i} is the i'th inference engine.
wolffd@0: % time(e) = elapsed time for doing inference with engine e
wolffd@0: %
wolffd@0: % The list below gives optional arguments [default value in brackets].
wolffd@0: %
wolffd@0: % exact - specifies which engines do exact inference [ 1:length(engine) ]
wolffd@0: % singletons_only - if 1, we only call marginal_nodes, else this  and marginal_family [0]
wolffd@0: % maximize - 1 means we do max-propagation, 0 means sum-propagation [0]
wolffd@0: % check_ll - 1 means we check that the log-likelihoods are correct [1]
wolffd@0: % observed - list of the observed ndoes [ bnet.observed ]
wolffd@0: % check_converged - list of loopy engines that should be checked for convergence [ [] ]
wolffd@0: %    If an engine has converged, it is added to the exact list.
wolffd@0: 
wolffd@0: 
wolffd@0: % set default params
wolffd@0: exact = 1:length(engine);
wolffd@0: singletons_only = 0;
wolffd@0: maximize = 0;
wolffd@0: check_ll = 1;
wolffd@0: observed = bnet.observed;
wolffd@0: check_converged = [];
wolffd@0: 
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'exact', exact = args{i+1};
wolffd@0:    case 'singletons_only', singletons_only = args{i+1};
wolffd@0:    case 'maximize', maximize = args{i+1};
wolffd@0:    case 'check_ll', check_ll = args{i+1};
wolffd@0:    case 'observed', observed = args{i+1};
wolffd@0:    case 'check_converged', check_converged = args{i+1};
wolffd@0:    otherwise,
wolffd@0:     error(['unrecognized argument ' args{i}])
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: E = length(engine);
wolffd@0: ref = exact(1); % reference
wolffd@0: 
wolffd@0: N = length(bnet.dag);
wolffd@0: ev = sample_bnet(bnet);
wolffd@0: evidence = cell(1,N);
wolffd@0: evidence(observed) = ev(observed);
wolffd@0: %celldisp(evidence(observed))
wolffd@0: 
wolffd@0: for i=1:E
wolffd@0:   tic;
wolffd@0:   if check_ll
wolffd@0:     [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence, 'maximize', maximize);
wolffd@0:   else
wolffd@0:     engine{i} = enter_evidence(engine{i}, evidence, 'maximize', maximize);
wolffd@0:   end
wolffd@0:   time(i)=toc;
wolffd@0: end
wolffd@0: 
wolffd@0: for i=check_converged(:)'
wolffd@0:   niter = loopy_converged(engine{i});
wolffd@0:   if niter > 0
wolffd@0:     fprintf('loopy engine %d  converged in %d iterations\n', i, niter);
wolffd@0: %    exact = myunion(exact, i);
wolffd@0:   else
wolffd@0:     fprintf('loopy engine %d has not converged\n', i);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: cmp = exact(2:end);
wolffd@0: if check_ll
wolffd@0:   for i=cmp(:)'
wolffd@0:     assert(approxeq(ll(ref), ll(i)));
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: hnodes = mysetdiff(1:N, observed);
wolffd@0: 
wolffd@0: if ~singletons_only
wolffd@0:   get_marginals(engine, hnodes, exact, 0);
wolffd@0: end
wolffd@0: get_marginals(engine, hnodes, exact, 1);
wolffd@0: 
wolffd@0: %%%%%%%%%%
wolffd@0: 
wolffd@0: function get_marginals(engine, hnodes, exact, singletons)
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine{1});
wolffd@0: N = length(bnet.dag);
wolffd@0: cnodes_bitv = zeros(1,N);
wolffd@0: cnodes_bitv(bnet.cnodes) = 1;
wolffd@0: ref = exact(1); % reference
wolffd@0: cmp = exact(2:end);
wolffd@0: E = length(engine);
wolffd@0: 
wolffd@0: for n=hnodes(:)'
wolffd@0:   for e=1:E
wolffd@0:     if singletons
wolffd@0:       m{e} = marginal_nodes(engine{e}, n);
wolffd@0:     else
wolffd@0:       m{e} = marginal_family(engine{e}, n);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   for e=cmp(:)'
wolffd@0:     if cnodes_bitv(n)
wolffd@0:       assert(approxeq(m{ref}.mu, m{e}.mu))
wolffd@0:       assert(approxeq(m{ref}.Sigma, m{e}.Sigma))
wolffd@0:     else
wolffd@0:       assert(approxeq(m{ref}.T, m{e}.T))
wolffd@0:     end
wolffd@0:     assert(isequal(m{e}.domain, m{ref}.domain));
wolffd@0:   end
wolffd@0: end