wolffd@0: % make an unrolled HMM, convert to factor graph, and check that 
wolffd@0: % loopy propagation on the fgraph gives the exact answers.
wolffd@0: 
wolffd@0: seed = 1;
wolffd@0: rand('state', seed);
wolffd@0: randn('state', seed);
wolffd@0: 
wolffd@0: T = 3;
wolffd@0: Q = 3;
wolffd@0: O = 3;
wolffd@0: cts_obs = 0;
wolffd@0: param_tying = 1;
wolffd@0: bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
wolffd@0: 
wolffd@0: data = sample_bnet(bnet);
wolffd@0: 
wolffd@0: fgraph = bnet_to_fgraph(bnet);
wolffd@0: big_bnet = fgraph_to_bnet(fgraph);
wolffd@0: % converting factor graph back does not recover the structure of the original bnet
wolffd@0: 
wolffd@0: max_iter = 2*T;
wolffd@0: 
wolffd@0: engine = {};
wolffd@0: engine{1} = jtree_inf_engine(bnet);
wolffd@0: engine{2} = belprop_inf_engine(bnet, 'max_iter', max_iter);
wolffd@0: engine{3} = belprop_fg_inf_engine(fgraph, 'max_iter', max_iter);
wolffd@0: engine{4} = jtree_inf_engine(big_bnet);
wolffd@0: nengines = length(engine);
wolffd@0: 
wolffd@0: big_engine = 4;
wolffd@0: fgraph_engine = 3;
wolffd@0: 
wolffd@0: 
wolffd@0: N = 2*T;
wolffd@0: evidence = cell(1,N);
wolffd@0: onodes = bnet.observed;
wolffd@0: evidence(onodes) = data(onodes);
wolffd@0: hnodes = mysetdiff(1:N, onodes);
wolffd@0: 
wolffd@0: bigN = length(big_bnet.dag);
wolffd@0: big_evidence = cell(1, bigN);
wolffd@0: big_evidence(onodes) = data(onodes);
wolffd@0: big_evidence(N+1:end) = {1}; % factors are observed to be 1
wolffd@0: 
wolffd@0: ll = zeros(1, nengines);
wolffd@0: for i=1:nengines
wolffd@0:   if i==big_engine
wolffd@0:     tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, big_evidence); toc
wolffd@0:   else
wolffd@0:     tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence); toc
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % compare all engines to engine{1}
wolffd@0: 
wolffd@0: % the log likelihood values may be bogus...
wolffd@0: for i=2:nengines
wolffd@0:   %assert(approxeq(ll(1), ll(i)));
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: marg = zeros(T, nengines, Q); % marg(t,e,:)
wolffd@0: for t=1:T
wolffd@0:   for e=1:nengines
wolffd@0:     m = marginal_nodes(engine{e}, t);
wolffd@0:     marg(t,e,:) = m.T;
wolffd@0:   end
wolffd@0: end
wolffd@0: marg
wolffd@0: 
wolffd@0: 
wolffd@0: m = cell(nengines, T);
wolffd@0: for i=1:T
wolffd@0:   for e=1:nengines
wolffd@0:     m{e,i} = marginal_nodes(engine{e}, hnodes(i));
wolffd@0:   end
wolffd@0:   for e=2:nengines
wolffd@0:     assert(approxeq(m{e,i}.T, m{1,i}.T));
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: mpe = {};
wolffd@0: ll = zeros(1, nengines);
wolffd@0: for e=1:nengines
wolffd@0:   if e==big_engine
wolffd@0:     mpe{e} = find_mpe(engine{e}, big_evidence);
wolffd@0:     mpe{e} = mpe{e}(1:N); % chop off dummy nodes
wolffd@0:   else
wolffd@0:     mpe{e} = find_mpe(engine{e}, evidence);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % fgraph can't compute loglikelihood for software reasons
wolffd@0: % jtree on the big_bnet gives the wrong ll
wolffd@0: for e=2:nengines
wolffd@0:   %assert(approxeq(ll(1), ll(e)));
wolffd@0:   assert(approxeq(cell2num(mpe{1}), cell2num(mpe{e})))
wolffd@0: end