wolffd@0: function engine = frontier_inf_engine(bnet)
wolffd@0: % FRONTIER_INF_ENGINE Inference engine for DBNs which which uses the frontier algorithm.
wolffd@0: % engine = frontier_inf_engine(bnet)
wolffd@0: %
wolffd@0: % The frontier algorithm extends the forwards-backwards algorithm to DBNs in the obvious way,
wolffd@0: % maintaining a joint distribution (frontier) over all the nodes in a time slice.
wolffd@0: % When all the hidden nodes in the DBN are persistent (have children in the next time slice),
wolffd@0: % its theoretical running time is often similar to that of the junction tree algorithm,
wolffd@0: % although in practice, this algorithm seems to very slow (at least in matlab).
wolffd@0: % However, it is extremely simple to describe and implement.
wolffd@0: %
wolffd@0: % Suppose there are n binary nodes per slice, so the frontier takes O(2^n) space.
wolffd@0: % Each time step takes between O(n 2^{n+1}) and O(n 2^{2n}) operations, depending on the graph structure.
wolffd@0: % The lower bound is achieved by a set of n independent chains, as in a factorial HMM.
wolffd@0: % The upper bound is achieved by a set of n fully interconnected chains, as in an HMM.
wolffd@0: %
wolffd@0: % The factor of n arises because we need to multiply in each CPD from slice t+1.
wolffd@0: % The second factor depends on the size of the frontier to which we add the new node.
wolffd@0: % In an FHMM, once we have added X(i,t+1), we can marginalize out X(i,t) from the frontier, since
wolffd@0: % no other nodes depend on it; hence the frontier never contains more than n+1 nodes.
wolffd@0: % In a fully coupled HMM, we must leave X(i,t) in the frontier until all X(j,t+1) have been
wolffd@0: % added; hence the frontier will contain 2*n nodes at its peak.
wolffd@0: %
wolffd@0: % For details, see
wolffd@0: %   "The Factored Frontier Algorithm for Approximate Inference in DBNs",
wolffd@0: %   Kevin Murphy and Yair Weiss, UAI 01.
wolffd@0: 
wolffd@0: ns = bnet.node_sizes_slice;
wolffd@0: onodes = bnet.observed;
wolffd@0: ns(onodes) = 1;
wolffd@0: ss = length(bnet.intra);
wolffd@0: 
wolffd@0: [engine.ops, engine.fdom] = best_first_frontier_seq(ns, bnet.dag);
wolffd@0: engine.ops1 = 1:ss;
wolffd@0: 
wolffd@0: engine.fwdback = [];
wolffd@0: engine.fwd_frontier = [];
wolffd@0: engine.back_frontier = [];
wolffd@0: 
wolffd@0: engine.fdom1 = cell(1,ss);
wolffd@0: for s=1:ss
wolffd@0:   engine.fdom1{s} = 1:s;
wolffd@0: end
wolffd@0: 
wolffd@0: engine = class(engine, 'frontier_inf_engine', inf_engine(bnet));
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%
wolffd@0: 
wolffd@0: function [ops, frontier_set] = best_first_frontier_seq(ns, dag)
wolffd@0: % BEST_FIRST_FRONTIER_SEQ Do a greedy search for the sequence of additions/removals to the frontier.
wolffd@0: % [ops, frontier_set] = best_first_frontier_seq(ns, dag)
wolffd@0: %
wolffd@0: % We maintain 3 sets: the frontier (F), the right set (R), and the left set (L).
wolffd@0: % The invariant is that the nodes in R are d-separated from L given F.
wolffd@0: % We start with slice 1 in F and slice 2 in R.
wolffd@0: % The goal is to move slice 1 from F to L, and slice 2 from R to F, so as to minimize the size
wolffd@0: % of the frontier at each step, where the size(F) = product of the node-sizes of nodes in F.
wolffd@0: % A node may be removed (from F to L) if it has no children in R.
wolffd@0: % A node may be added (from R to F) if its parents are in F.
wolffd@0: %
wolffd@0: % ns(i) = num. discrete values node i can take on (i=1..ss, where ss = slice size)
wolffd@0: % dag is the (2*ss) x (2*ss) adjacency matrix for the 2-slice DBN.
wolffd@0: 
wolffd@0: % Example:
wolffd@0: %
wolffd@0: % 4    9
wolffd@0: % ^    ^
wolffd@0: % |    |
wolffd@0: % 2 -> 7
wolffd@0: % ^    ^
wolffd@0: % |    |
wolffd@0: % 1 -> 6
wolffd@0: % |    |
wolffd@0: % v    v
wolffd@0: % 3 -> 8
wolffd@0: % |    |
wolffd@0: % v    V
wolffd@0: % 5    10
wolffd@0: %
wolffd@0: % ops = -4, -5, 6, -1, 7, -2, 8, -3, 9, 10
wolffd@0: 
wolffd@0: ss = length(ns);
wolffd@0: ns = [ns(:)' ns(:)'];
wolffd@0: ops = zeros(1,ss);
wolffd@0: L = []; F = 1:ss; R = (1:ss)+ss;
wolffd@0: frontier_set = cell(1,2*ss);
wolffd@0: for s=1:2*ss
wolffd@0:   remcost = inf*ones(1,2*ss);
wolffd@0:   %disp(['L: ' num2str(L) ', F: ' num2str(F) ', R: ' num2str(R)]);
wolffd@0:   maybe_removable = myintersect(F, 1:ss);
wolffd@0:   for n=maybe_removable(:)'
wolffd@0:     cs = children(dag, n);
wolffd@0:     if isempty(myintersect(cs, R))
wolffd@0:       remcost(n) = prod(ns(mysetdiff(F, n)));
wolffd@0:     end
wolffd@0:   end
wolffd@0:   %remcost
wolffd@0:   if any(remcost < inf)
wolffd@0:     n = argmin(remcost);
wolffd@0:     ops(s) = -n;
wolffd@0:     L = myunion(L, n);
wolffd@0:     F = mysetdiff(F, n);
wolffd@0:   else
wolffd@0:     addcost = inf*ones(1,2*ss);
wolffd@0:     for n=R(:)'
wolffd@0:       ps = parents(dag, n);
wolffd@0:       if mysubset(ps, F)
wolffd@0: 	addcost(n) = prod(ns(myunion(F, [ps n])));
wolffd@0:       end
wolffd@0:     end
wolffd@0:     %addcost
wolffd@0:     assert(any(addcost < inf));
wolffd@0:     n = argmin(addcost);
wolffd@0:     ops(s) = n;
wolffd@0:     R  = mysetdiff(R, n);
wolffd@0:     F = myunion(F, n);
wolffd@0:   end
wolffd@0:   %fprintf('op at step %d = %d\n\n', s, ops(s));
wolffd@0:   frontier_set{s} = F;
wolffd@0: end