Daniel@0: function [time, engine] = cmp_inference_dbn(bnet, engine, T, varargin)
Daniel@0: % CMP_INFERENCE_DBN Compare several inference engines on a DBN
Daniel@0: % function [time, engine] = cmp_inference_dbn(bnet, engine, T, ...)
Daniel@0: %
Daniel@0: % engine{i} is the i'th inference engine.
Daniel@0: % time(e) = elapsed time for doing inference with engine e
Daniel@0: %
Daniel@0: % The list below gives optional arguments [default value in brackets].
Daniel@0: %
Daniel@0: % exact - specifies which engines do exact inference [ 1:length(engine) ]
Daniel@0: % singletons_only - if 1, we only call marginal_nodes, else this  and marginal_family [0]
Daniel@0: % check_ll - 1 means we check that the log-likelihoods are correct [1]
Daniel@0: 
Daniel@0: % set default params
Daniel@0: exact = 1:length(engine);
Daniel@0: singletons_only = 0;
Daniel@0: check_ll = 1;
Daniel@0: onodes = bnet.observed;
Daniel@0: 
Daniel@0: args = varargin;
Daniel@0: nargs = length(args);
Daniel@0: for i=1:2:nargs
Daniel@0:   switch args{i},
Daniel@0:    case 'exact', exact = args{i+1};
Daniel@0:    case 'singletons_only', singletons_only = args{i+1};
Daniel@0:    case 'check_ll', check_ll = args{i+1};
Daniel@0:    case 'observed', onodes = args{i+1};
Daniel@0:    otherwise,
Daniel@0:     error(['unrecognized argument ' args{i}])
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: E = length(engine);
Daniel@0: ref = exact(1); % reference
Daniel@0: 
Daniel@0: ss = length(bnet.intra);
Daniel@0: ev = sample_dbn(bnet, 'length', T);
Daniel@0: evidence = cell(ss,T);
Daniel@0: evidence(onodes,:) = ev(onodes, :);
Daniel@0: 
Daniel@0: for i=1:E
Daniel@0:   tic;
Daniel@0:   [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence);
Daniel@0:   time(i)=toc;
Daniel@0:   fprintf('engine %d took %6.4f seconds\n', i, time(i));
Daniel@0: end
Daniel@0: 
Daniel@0: cmp = mysetdiff(exact, ref);
Daniel@0: if check_ll
Daniel@0:   for i=cmp(:)'
Daniel@0:     if ~approxeq(ll(ref), ll(i))
Daniel@0:       error(['engine ' num2str(i) ' has wrong ll'])
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0: ll
Daniel@0: 
Daniel@0: hnodes = mysetdiff(1:ss, onodes);
Daniel@0: 
Daniel@0: if ~singletons_only
Daniel@0:   get_marginals(engine, hnodes, exact, 0, T);
Daniel@0: end
Daniel@0: get_marginals(engine, hnodes, exact, 1, T);
Daniel@0: 
Daniel@0: %%%%%%%%%%
Daniel@0: 
Daniel@0: function get_marginals(engine, hnodes, exact, singletons, T)
Daniel@0: 
Daniel@0: bnet = bnet_from_engine(engine{1});
Daniel@0: N = length(bnet.intra);
Daniel@0: cnodes_bitv = zeros(1,N);
Daniel@0: cnodes_bitv(bnet.cnodes) = 1;
Daniel@0: ref = exact(1); % reference
Daniel@0: cmp = exact(2:end);
Daniel@0: E = length(engine);
Daniel@0: m = cell(1,E);
Daniel@0: 
Daniel@0: for t=1:T
Daniel@0:   for n=1:N
Daniel@0:   %for n=hnodes(:)'
Daniel@0:     for e=1:E
Daniel@0:       if singletons
Daniel@0: 	m{e} = marginal_nodes(engine{e}, n, t);
Daniel@0:       else
Daniel@0: 	m{e} = marginal_family(engine{e}, n, t);
Daniel@0:       end
Daniel@0:     end
Daniel@0:     for e=cmp(:)'
Daniel@0:       assert(isequal(m{e}.domain, m{ref}.domain));
Daniel@0:       if cnodes_bitv(n) & isfield(m{e}, 'mu') & isfield(m{ref}, 'mu')
Daniel@0: 	wrong = ~approxeq(m{ref}.mu, m{e}.mu) | ~approxeq(m{ref}.Sigma, m{e}.Sigma);
Daniel@0:       else
Daniel@0: 	wrong = ~approxeq(m{ref}.T(:), m{e}.T(:));
Daniel@0:       end
Daniel@0:       if wrong
Daniel@0: 	error(sprintf('engine %d is wrong; n=%d, t=%d, fam=%d', e, n, t, ~singletons))
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: end