wolffd@0
|
1 function engine = stab_cond_gauss_inf_engine(bnet)
|
wolffd@0
|
2 % STAB_COND_GAUSS_INF_ENGINE Junction tree using stable CG potentials
|
wolffd@0
|
3 % engine = cond_gauss_inf_engine(bnet)
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % This class was written by Shan Huang (shan.huang@intel.com) 2001
|
wolffd@0
|
6 % and fixed by Rainer Deventer deventer@informatik.uni-erlangen.de March 2003
|
wolffd@0
|
7 N = length(bnet.dag);
|
wolffd@0
|
8 clusters = {};
|
wolffd@0
|
9 root = N;
|
wolffd@0
|
10 stages = { 1:N };
|
wolffd@0
|
11 onodes = [];
|
wolffd@0
|
12 engine = init_fields;
|
wolffd@0
|
13 engine.evidence = [];
|
wolffd@0
|
14 engine = class(engine, 'stab_cond_gauss_inf_engine', inf_engine(bnet));
|
wolffd@0
|
15
|
wolffd@0
|
16 ns = bnet.node_sizes(:);
|
wolffd@0
|
17 ns(onodes) = 1; % observed nodes have only 1 possible value
|
wolffd@0
|
18
|
wolffd@0
|
19 %[engine.jtree, dummy, engine.cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ...
|
wolffd@0
|
20 % dag_to_jtree(bnet, onodes, stages, clusters);
|
wolffd@0
|
21
|
wolffd@0
|
22
|
wolffd@0
|
23 partial_order = determine_elim_constraints(bnet, onodes);
|
wolffd@0
|
24 strong = ~isempty(partial_order);
|
wolffd@0
|
25 stages = {};
|
wolffd@0
|
26 clusters = {};
|
wolffd@0
|
27 [engine.jtree, dummy_root, engine.cliques, B, w, elim_order] =
|
wolffd@0
|
28 graph_to_jtree(moralize(bnet.dag), ns, partial_order, stages, clusters);
|
wolffd@0
|
29
|
wolffd@0
|
30
|
wolffd@0
|
31 engine.cliques_bitv = B;
|
wolffd@0
|
32 engine.clique_weight = w;
|
wolffd@0
|
33 C = length(engine.cliques);
|
wolffd@0
|
34 engine.clpot = cell(1,C);
|
wolffd@0
|
35
|
wolffd@0
|
36 % A node can be a member of many cliques, but is assigned to exactly one, to avoid
|
wolffd@0
|
37 % double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that
|
wolffd@0
|
38 % contains i's family, so it can accomodate its CPD.
|
wolffd@0
|
39
|
wolffd@0
|
40 engine.clq_ass_to_node = zeros(1, N);
|
wolffd@0
|
41 num_cliques = length(engine.cliques);
|
wolffd@0
|
42 for i=1:N
|
wolffd@0
|
43 clqs_containing_family = find(all(B(:,family(bnet.dag, i)), 2)); % all selected columns must be 1
|
wolffd@0
|
44 c = clqs_containing_family(argmin(w(clqs_containing_family)));
|
wolffd@0
|
45 engine.clq_ass_to_node(i) = c;
|
wolffd@0
|
46 end
|
wolffd@0
|
47
|
wolffd@0
|
48 % Compute the separators between connected cliques.
|
wolffd@0
|
49 [is,js] = find(engine.jtree > 0);
|
wolffd@0
|
50 engine.separator = cell(num_cliques, num_cliques);
|
wolffd@0
|
51 for k=1:length(is)
|
wolffd@0
|
52 i = is(k); j = js(k);
|
wolffd@0
|
53 engine.separator{i,j} = find(B(i,:) & B(j,:)); % intersect(cliques{i}, cliques{j});
|
wolffd@0
|
54 end
|
wolffd@0
|
55 %keyboard;
|
wolffd@0
|
56 engine.seppot = cell(C,C);
|
wolffd@0
|
57
|
wolffd@0
|
58 pot_type = 'scg';
|
wolffd@0
|
59 check_for_cd_arcs([], bnet.cnodes, bnet.dag);
|
wolffd@0
|
60
|
wolffd@0
|
61 % Make the jtree rooted, so there is a fixed message passing order.
|
wolffd@0
|
62 if strong
|
wolffd@0
|
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
64 % Start the search for the strong root at the clique with the %
|
wolffd@0
|
65 % highest number. %
|
wolffd@0
|
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
67 root = length(engine.cliques);
|
wolffd@0
|
68 root_found = 0;
|
wolffd@0
|
69
|
wolffd@0
|
70 while ((~root_found) & (root >= 1))
|
wolffd@0
|
71 root_found = test_strong_root(engine.jtree,engine.cliques,bnet.dnodes,root);
|
wolffd@0
|
72 if ~root_found
|
wolffd@0
|
73 root = root - 1;
|
wolffd@0
|
74 end
|
wolffd@0
|
75 end
|
wolffd@0
|
76 assert(root > 0)
|
wolffd@0
|
77 engine.root = root;
|
wolffd@0
|
78 % the last clique is guaranteed to be a strong root
|
wolffd@0
|
79 %engine.root = length(engine.cliques);
|
wolffd@0
|
80 else
|
wolffd@0
|
81 % jtree_dbn_inf_engine requires the root to contain the interface.
|
wolffd@0
|
82 % This may conflict with the strong root requirement! *********** BUG *************
|
wolffd@0
|
83 engine.root = clq_containing_nodes(engine, root);
|
wolffd@0
|
84 if engine.root <= 0
|
wolffd@0
|
85 error(['no clique contains ' num2str(root)]);
|
wolffd@0
|
86 end
|
wolffd@0
|
87 end
|
wolffd@0
|
88
|
wolffd@0
|
89 [engine.jtree, engine.preorder, engine.postorder] = mk_rooted_tree(engine.jtree, engine.root);
|
wolffd@0
|
90
|
wolffd@0
|
91 % Evaluate CPDs with evidence, and convert to potentials
|
wolffd@0
|
92 pot = cell(1, N);
|
wolffd@0
|
93 inited = zeros(1, C);
|
wolffd@0
|
94 clpot = cell(1, C);
|
wolffd@0
|
95 evidence = cell(1, N);
|
wolffd@0
|
96 for n=1:N
|
wolffd@0
|
97 fam = family(bnet.dag, n);
|
wolffd@0
|
98 e = bnet.equiv_class(n);
|
wolffd@0
|
99 %pot{n} = CPD_to_scgpot(bnet.CPD{e}, fam, ns, bnet.cnodes, evidence);
|
wolffd@0
|
100 pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
|
wolffd@0
|
101 cindex = engine.clq_ass_to_node(n);
|
wolffd@0
|
102 if inited(cindex)
|
wolffd@0
|
103 clpot{cindex} = direct_combine_pots(pot{n}, clpot{cindex});
|
wolffd@0
|
104 else
|
wolffd@0
|
105 clpot{cindex} = pot{n};
|
wolffd@0
|
106 inited(cindex) = 1;
|
wolffd@0
|
107 end
|
wolffd@0
|
108 end
|
wolffd@0
|
109
|
wolffd@0
|
110 for i=1:C
|
wolffd@0
|
111 if inited(i) == 0
|
wolffd@0
|
112 clpot{i} = scgpot([], [], [], []);
|
wolffd@0
|
113 end
|
wolffd@0
|
114 end
|
wolffd@0
|
115
|
wolffd@0
|
116 seppot = cell(C, C);
|
wolffd@0
|
117 % separators are is not need to initialize
|
wolffd@0
|
118
|
wolffd@0
|
119 % collect to root (node to parents)
|
wolffd@0
|
120 % Unlike the HUGIN architecture the complements are stored in the cliques during COLLECT
|
wolffd@0
|
121 % and the separators are not playing a specific role during this process
|
wolffd@0
|
122 for n=engine.postorder(1:end-1)
|
wolffd@0
|
123 for p=parents(engine.jtree, n)
|
wolffd@0
|
124 if ~isempty(engine.separator{p,n})
|
wolffd@0
|
125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
126 % The empty case might happen for unlinked nodes, i.e. the DAG is not %
|
wolffd@0
|
127 % a single tree, but a forest %
|
wolffd@0
|
128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
129 [margpot, comppot] = complement_pot(clpot{n}, engine.separator{p,n});
|
wolffd@0
|
130 clpot{n} = comppot;
|
wolffd@0
|
131 clpot{p} = combine_pots(clpot{p}, margpot);
|
wolffd@0
|
132 end
|
wolffd@0
|
133 end
|
wolffd@0
|
134 end
|
wolffd@0
|
135
|
wolffd@0
|
136 % distribute message from root
|
wolffd@0
|
137 % We have not to store the weak clique marginals and keep the original complement potentials.
|
wolffd@0
|
138 % This is a minor variation of HUGIN architecture.
|
wolffd@0
|
139 temppot = clpot;
|
wolffd@0
|
140 for n=engine.preorder
|
wolffd@0
|
141 for c=children(engine.jtree, n)
|
wolffd@0
|
142 seppot{n,c} = marginalize_pot(temppot{n}, engine.separator{n,c});
|
wolffd@0
|
143 temppot{c} = direct_combine_pots(temppot{c}, seppot{n,c});
|
wolffd@0
|
144 end
|
wolffd@0
|
145 end
|
wolffd@0
|
146
|
wolffd@0
|
147 engine.clpot = clpot;
|
wolffd@0
|
148 engine.seppot = seppot;
|
wolffd@0
|
149
|
wolffd@0
|
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
151 % init_fields() %
|
wolffd@0
|
152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
wolffd@0
|
153 function engine = init_fields()
|
wolffd@0
|
154
|
wolffd@0
|
155 engine.evidence = [];
|
wolffd@0
|
156 engine.jtree = [];
|
wolffd@0
|
157 engine.cliques = [];
|
wolffd@0
|
158 engine.cliques_bitv = [];
|
wolffd@0
|
159 engine.clique_weight = [];
|
wolffd@0
|
160 engine.preorder = [];
|
wolffd@0
|
161 engine.postorder = [];
|
wolffd@0
|
162 engine.root = [];
|
wolffd@0
|
163 engine.clq_ass_to_node = [];
|
wolffd@0
|
164 engine.separator = [];
|
wolffd@0
|
165 engine.clpot =[];
|
wolffd@0
|
166 engine.seppot = [];
|
wolffd@0
|
167
|
wolffd@0
|
168
|
wolffd@0
|
169
|
wolffd@0
|
170
|
wolffd@0
|
171
|
wolffd@0
|
172
|
wolffd@0
|
173
|
wolffd@0
|
174
|
wolffd@0
|
175
|
wolffd@0
|
176
|
wolffd@0
|
177
|
wolffd@0
|
178
|