Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/static/@stab_cond_gauss_inf_engine/stab_cond_gauss_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 = 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