wolffd@0: function [err, time, engine] = cmp_inference(bnet, engine, exact, T, filter, singletons, maximize)
wolffd@0: % CMP_INFERENCE Compare several inference engines on a DBN
wolffd@0: % [err, time, engine] = cmp_inference(bnet, engine, exact, T, filter, singletons, maximize)
wolffd@0: %
wolffd@0: % engine{i} is the i'th inference engine.
wolffd@0: % 'exact' specifies which engines do exact inference - 
wolffd@0: %   we check that these all give the same results.
wolffd@0: % 'T' is the length of the random sequence we generate.
wolffd@0: % If filter=1, we do filtering, else smoothing (default: smoothing)
wolffd@0: % If singletons=1, we compare marginal_nodes, else marginal_family (default: family)
wolffd@0: %
wolffd@0: % err(e,n,t) = sum_i | Pr_exact(X(n,t)=i) - Pr_e(X(n,t)=i) |
wolffd@0: %   where Pr_e = prob. according to engine e
wolffd@0: % time(e) = elapsed time for doing inference with engine e
wolffd@0: 
wolffd@0: err = [];
wolffd@0: 
wolffd@0: if nargin < 5, filter = 0; end
wolffd@0: if nargin < 6, singletons = 0; end
wolffd@0: if nargin < 7, maximize = 0; end
wolffd@0: 
wolffd@0: check_ll = 1;
wolffd@0: 
wolffd@0: assert(~maximize);
wolffd@0: 
wolffd@0: E = length(engine);
wolffd@0: ref = exact(1); % reference
wolffd@0: 
wolffd@0: ss = length(bnet.intra);
wolffd@0: ev = sample_dbn(bnet, 'length', T);
wolffd@0: evidence = cell(ss,T);
wolffd@0: onodes = bnet.observed;
wolffd@0: evidence(onodes,:) = ev(onodes, :);
wolffd@0: 
wolffd@0: assert(~filter);
wolffd@0: for i=1:E
wolffd@0:   tic;
wolffd@0:   %[engine{i}, ll(i)] = enter_evidence(engine{i}, evidence, 'maximize', maximize);
wolffd@0:   [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence);
wolffd@0:   time(i)=toc;
wolffd@0:   fprintf('engine %d took %6.4f seconds\n', i, time(i));
wolffd@0: end
wolffd@0: 
wolffd@0: cmp = mysetdiff(exact, ref);
wolffd@0: if check_ll
wolffd@0: for i=cmp(:)'
wolffd@0:   if ~approxeq(ll(ref), ll(i))
wolffd@0:     error(['engine ' num2str(i) ' has wrong ll'])
wolffd@0:   end
wolffd@0: end
wolffd@0: end
wolffd@0: ll
wolffd@0: 
wolffd@0: hnodes = mysetdiff(1:ss, onodes);
wolffd@0: m = cell(1,E);
wolffd@0: for t=1:T
wolffd@0:   for n=hnodes(:)'
wolffd@0:     for e=1:E
wolffd@0:       if singletons
wolffd@0: 	m{e} = marginal_nodes(engine{e}, n, t);
wolffd@0:       else
wolffd@0: 	m{e} = marginal_family(engine{e}, n, t);
wolffd@0:       end
wolffd@0:     end
wolffd@0:     for e=1:E
wolffd@0:       assert(isequal(m{e}.domain, m{ref}.domain));
wolffd@0:     end
wolffd@0:     for e=cmp(:)'
wolffd@0:       if ~approxeq(m{ref}.T(:), m{e}.T(:))
wolffd@0: 	str= sprintf('engine %d is wrong; n=%d, t=%d', e, n, t);
wolffd@0: 	error(str)
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end