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