|
wolffd@0
|
1 function order = strong_elim_order(G, node_sizes, partial_order)
|
|
wolffd@0
|
2 % STRONG_ELIM_ORDER Find an elimination order to produce a strongly triangulated graph.
|
|
wolffd@0
|
3 % order = strong_elim_order(moral_graph, node_sizes, partial_order)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % partial_order(i,j)=1 if we must marginalize i *after* j
|
|
wolffd@0
|
6 % (so i will be nearer the strong root).
|
|
wolffd@0
|
7 % e.g., if j is a decision node and i is its information set:
|
|
wolffd@0
|
8 % we cannot maximize j if we have marginalized out some of i
|
|
wolffd@0
|
9 % e.g., if j is a continuous child and i is its discrete parent:
|
|
wolffd@0
|
10 % we want to integrate out the cts nodes before the discrete ones,
|
|
wolffd@0
|
11 % so that the marginal is strong.
|
|
wolffd@0
|
12 %
|
|
wolffd@0
|
13 % For details, see
|
|
wolffd@0
|
14 % - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
|
|
wolffd@0
|
15 % - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical
|
|
wolffd@0
|
16 % association models", JASA 87(420):1098--1108, 1992.
|
|
wolffd@0
|
17 %
|
|
wolffd@0
|
18 % On p369 of the Jensen paper, they state "the reverse of the elimination order must be some
|
|
wolffd@0
|
19 % extension of [the partial order] to a total order".
|
|
wolffd@0
|
20 % We make no attempt to find the best such total ordering, in the sense of minimizing the weight
|
|
wolffd@0
|
21 % of the resulting cliques.
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 % Example from the Jensen paper:
|
|
wolffd@0
|
24 % Let us number the nodes in Fig 1 from top to bottom, left to right,
|
|
wolffd@0
|
25 % so a=1,b=2,D1=3,c=4,...,l=14,j=15,k=16.
|
|
wolffd@0
|
26 % The elimination ordering they propose on p370 is [14 15 16 11 12 1 4 5 10 8 13 9 7 6 3 2];
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 if 0
|
|
wolffd@0
|
29 total_order = topological_sort(partial_order);
|
|
wolffd@0
|
30 order = total_order(end:-1:1); % no attempt to find an optimal constrained ordering!
|
|
wolffd@0
|
31 return;
|
|
wolffd@0
|
32 end
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % The following implementation is due to Ilya Shpitser and seems to give wrong
|
|
wolffd@0
|
35 % results on cg1
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 n = length(G);
|
|
wolffd@0
|
38 MG = G; % copy the original graph
|
|
wolffd@0
|
39 uneliminated = ones(1,n);
|
|
wolffd@0
|
40 order = zeros(1,n);
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 for i=1:n
|
|
wolffd@0
|
43 roots = [];
|
|
wolffd@0
|
44 k = 1;
|
|
wolffd@0
|
45 for j=1:n
|
|
wolffd@0
|
46 if sum(partial_order(j,:)) == 0
|
|
wolffd@0
|
47 roots(k) = j;
|
|
wolffd@0
|
48 k = k + 1;
|
|
wolffd@0
|
49 end
|
|
wolffd@0
|
50 end
|
|
wolffd@0
|
51 U = find(uneliminated);
|
|
wolffd@0
|
52 valid = myintersect(U, roots);
|
|
wolffd@0
|
53 % Choose the best node from the set of valid candidates
|
|
wolffd@0
|
54 score1 = zeros(1,length(valid));
|
|
wolffd@0
|
55 score2 = zeros(1,length(valid));
|
|
wolffd@0
|
56 for j=1:length(valid)
|
|
wolffd@0
|
57 k = valid(j);
|
|
wolffd@0
|
58 ns = myintersect(neighbors(G, k), U);
|
|
wolffd@0
|
59 l = length(ns);
|
|
wolffd@0
|
60 M = MG(ns,ns);
|
|
wolffd@0
|
61 score1(j) = l^2 - sum(M(:)); % num. added edges
|
|
wolffd@0
|
62 score2(j) = prod(node_sizes([k ns])); % weight of clique
|
|
wolffd@0
|
63 end
|
|
wolffd@0
|
64 j1s = find(score1==min(score1));
|
|
wolffd@0
|
65 j = j1s(argmin(score2(j1s)));
|
|
wolffd@0
|
66 k = valid(j);
|
|
wolffd@0
|
67 uneliminated(k) = 0;
|
|
wolffd@0
|
68 order(i) = k;
|
|
wolffd@0
|
69 ns = myintersect(neighbors(G, k), U);
|
|
wolffd@0
|
70 if ~isempty(ns)
|
|
wolffd@0
|
71 G(ns,ns) = 1;
|
|
wolffd@0
|
72 G = setdiag(G,0);
|
|
wolffd@0
|
73 end
|
|
wolffd@0
|
74 partial_order(:,k) = 0;
|
|
wolffd@0
|
75 end
|
