Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/dynamic/@frontier_inf_engine/frontier_inf_engine.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function engine = frontier_inf_engine(bnet)
wolffd@0 2 % FRONTIER_INF_ENGINE Inference engine for DBNs which which uses the frontier algorithm.
wolffd@0 3 % engine = frontier_inf_engine(bnet)
wolffd@0 4 %
wolffd@0 5 % The frontier algorithm extends the forwards-backwards algorithm to DBNs in the obvious way,
wolffd@0 6 % maintaining a joint distribution (frontier) over all the nodes in a time slice.
wolffd@0 7 % When all the hidden nodes in the DBN are persistent (have children in the next time slice),
wolffd@0 8 % its theoretical running time is often similar to that of the junction tree algorithm,
wolffd@0 9 % although in practice, this algorithm seems to very slow (at least in matlab).
wolffd@0 10 % However, it is extremely simple to describe and implement.
wolffd@0 11 %
wolffd@0 12 % Suppose there are n binary nodes per slice, so the frontier takes O(2^n) space.
wolffd@0 13 % Each time step takes between O(n 2^{n+1}) and O(n 2^{2n}) operations, depending on the graph structure.
wolffd@0 14 % The lower bound is achieved by a set of n independent chains, as in a factorial HMM.
wolffd@0 15 % The upper bound is achieved by a set of n fully interconnected chains, as in an HMM.
wolffd@0 16 %
wolffd@0 17 % The factor of n arises because we need to multiply in each CPD from slice t+1.
wolffd@0 18 % The second factor depends on the size of the frontier to which we add the new node.
wolffd@0 19 % In an FHMM, once we have added X(i,t+1), we can marginalize out X(i,t) from the frontier, since
wolffd@0 20 % no other nodes depend on it; hence the frontier never contains more than n+1 nodes.
wolffd@0 21 % In a fully coupled HMM, we must leave X(i,t) in the frontier until all X(j,t+1) have been
wolffd@0 22 % added; hence the frontier will contain 2*n nodes at its peak.
wolffd@0 23 %
wolffd@0 24 % For details, see
wolffd@0 25 % "The Factored Frontier Algorithm for Approximate Inference in DBNs",
wolffd@0 26 % Kevin Murphy and Yair Weiss, UAI 01.
wolffd@0 27
wolffd@0 28 ns = bnet.node_sizes_slice;
wolffd@0 29 onodes = bnet.observed;
wolffd@0 30 ns(onodes) = 1;
wolffd@0 31 ss = length(bnet.intra);
wolffd@0 32
wolffd@0 33 [engine.ops, engine.fdom] = best_first_frontier_seq(ns, bnet.dag);
wolffd@0 34 engine.ops1 = 1:ss;
wolffd@0 35
wolffd@0 36 engine.fwdback = [];
wolffd@0 37 engine.fwd_frontier = [];
wolffd@0 38 engine.back_frontier = [];
wolffd@0 39
wolffd@0 40 engine.fdom1 = cell(1,ss);
wolffd@0 41 for s=1:ss
wolffd@0 42 engine.fdom1{s} = 1:s;
wolffd@0 43 end
wolffd@0 44
wolffd@0 45 engine = class(engine, 'frontier_inf_engine', inf_engine(bnet));
wolffd@0 46
wolffd@0 47
wolffd@0 48 %%%%%%%%%
wolffd@0 49
wolffd@0 50 function [ops, frontier_set] = best_first_frontier_seq(ns, dag)
wolffd@0 51 % BEST_FIRST_FRONTIER_SEQ Do a greedy search for the sequence of additions/removals to the frontier.
wolffd@0 52 % [ops, frontier_set] = best_first_frontier_seq(ns, dag)
wolffd@0 53 %
wolffd@0 54 % We maintain 3 sets: the frontier (F), the right set (R), and the left set (L).
wolffd@0 55 % The invariant is that the nodes in R are d-separated from L given F.
wolffd@0 56 % We start with slice 1 in F and slice 2 in R.
wolffd@0 57 % 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 58 % of the frontier at each step, where the size(F) = product of the node-sizes of nodes in F.
wolffd@0 59 % A node may be removed (from F to L) if it has no children in R.
wolffd@0 60 % A node may be added (from R to F) if its parents are in F.
wolffd@0 61 %
wolffd@0 62 % ns(i) = num. discrete values node i can take on (i=1..ss, where ss = slice size)
wolffd@0 63 % dag is the (2*ss) x (2*ss) adjacency matrix for the 2-slice DBN.
wolffd@0 64
wolffd@0 65 % Example:
wolffd@0 66 %
wolffd@0 67 % 4 9
wolffd@0 68 % ^ ^
wolffd@0 69 % | |
wolffd@0 70 % 2 -> 7
wolffd@0 71 % ^ ^
wolffd@0 72 % | |
wolffd@0 73 % 1 -> 6
wolffd@0 74 % | |
wolffd@0 75 % v v
wolffd@0 76 % 3 -> 8
wolffd@0 77 % | |
wolffd@0 78 % v V
wolffd@0 79 % 5 10
wolffd@0 80 %
wolffd@0 81 % ops = -4, -5, 6, -1, 7, -2, 8, -3, 9, 10
wolffd@0 82
wolffd@0 83 ss = length(ns);
wolffd@0 84 ns = [ns(:)' ns(:)'];
wolffd@0 85 ops = zeros(1,ss);
wolffd@0 86 L = []; F = 1:ss; R = (1:ss)+ss;
wolffd@0 87 frontier_set = cell(1,2*ss);
wolffd@0 88 for s=1:2*ss
wolffd@0 89 remcost = inf*ones(1,2*ss);
wolffd@0 90 %disp(['L: ' num2str(L) ', F: ' num2str(F) ', R: ' num2str(R)]);
wolffd@0 91 maybe_removable = myintersect(F, 1:ss);
wolffd@0 92 for n=maybe_removable(:)'
wolffd@0 93 cs = children(dag, n);
wolffd@0 94 if isempty(myintersect(cs, R))
wolffd@0 95 remcost(n) = prod(ns(mysetdiff(F, n)));
wolffd@0 96 end
wolffd@0 97 end
wolffd@0 98 %remcost
wolffd@0 99 if any(remcost < inf)
wolffd@0 100 n = argmin(remcost);
wolffd@0 101 ops(s) = -n;
wolffd@0 102 L = myunion(L, n);
wolffd@0 103 F = mysetdiff(F, n);
wolffd@0 104 else
wolffd@0 105 addcost = inf*ones(1,2*ss);
wolffd@0 106 for n=R(:)'
wolffd@0 107 ps = parents(dag, n);
wolffd@0 108 if mysubset(ps, F)
wolffd@0 109 addcost(n) = prod(ns(myunion(F, [ps n])));
wolffd@0 110 end
wolffd@0 111 end
wolffd@0 112 %addcost
wolffd@0 113 assert(any(addcost < inf));
wolffd@0 114 n = argmin(addcost);
wolffd@0 115 ops(s) = n;
wolffd@0 116 R = mysetdiff(R, n);
wolffd@0 117 F = myunion(F, n);
wolffd@0 118 end
wolffd@0 119 %fprintf('op at step %d = %d\n\n', s, ops(s));
wolffd@0 120 frontier_set{s} = F;
wolffd@0 121 end