wolffd@0: function [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG)
wolffd@0: % MK_SRONG_JTREE Make a strong junction tree.
wolffd@0: % [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG)
wolffd@0: %
wolffd@0: % Here is a definition of a strong jtree from Jensen et al. 1994:
wolffd@0: % "A junction tree is said to be strong if it has at least one distinguished clique R,
wolffd@0: % called a strong root, s.t. for each pair (C1,C2) of adjacent cliques in the tree,
wolffd@0: % with C1 closer to R than C2, there exists and ordering of [the nodes below] C2
wolffd@0: % that respects [the partial order] and with the vertices of the separator C1 intersect C2
wolffd@0: % preceeding the vertices [below C2] of C2 \ C1."
wolffd@0: %          
wolffd@0: % For details, see
wolffd@0: % - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
wolffd@0: %
wolffd@0: % MTG is the moralized, triangulated graph.
wolffd@0: % elim_order is the elimination ordering used to compute MTG.
wolffd@0: 
wolffd@0: 
wolffd@0: % Warning: this is a very naive implementation of the algorithm in Jensen et al.
wolffd@0: 
wolffd@0: n = length(elim_order);
wolffd@0: alpha(elim_order) = n:-1:1;
wolffd@0: % alpha(u) = i if we eliminate u at step n-i+1
wolffd@0: % i.e., vertices with higher alpha numbers are eliminated before vertices with lower numbers.
wolffd@0: % e.g., from the Jensen et al paper
wolffd@0: % node a=1 eliminated at step 6, so alpha(a)=16-6+1=11.
wolffd@0: % alpha = [11 1 2 10 9 3 4 7 5 8 13 12 6 16 15 14]
wolffd@0: 
wolffd@0: 
wolffd@0: % We sort the cliques in order of increasing index. The index of a clique C is defined as follows.
wolffd@0: % Let lower = {u | alpha(u) < alpha(v)}, and
wolffd@0: % let v in C be the highest-numbered vertex s.t. the vertices in W = lower intersect C
wolffd@0: % have a common neighbor u in U, where U = lower \ C.
wolffd@0: % If such a v exists, define index(C) = alpha(v), otherwise, index(C) = 1.
wolffd@0: % Intuitively, index(C) is the step in the elimination process at which C disappears.
wolffd@0: 
wolffd@0: num_cliques = length(cliques);
wolffd@0: index = zeros(1, num_cliques);
wolffd@0: for c = 1:num_cliques
wolffd@0:   C = cliques{c};
wolffd@0:   highest_num = -inf;
wolffd@0:   for vi = 1:length(C)
wolffd@0:     v = C(vi);
wolffd@0:     lower = find(alpha < alpha(v));
wolffd@0:     W = myintersect(lower, C);
wolffd@0:     U = mysetdiff(lower, C);
wolffd@0:     found = 0;
wolffd@0:     for ui=1:length(U)
wolffd@0:       u = U(ui);
wolffd@0:       if mysubset(W, neighbors(MTG, u))
wolffd@0: 	found = 1;
wolffd@0: 	break;
wolffd@0:       end
wolffd@0:     end
wolffd@0:     if found
wolffd@0:       if alpha(v) > highest_num
wolffd@0: 	highest_num = alpha(v);
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if highest_num == -inf
wolffd@0:     index(c) = 1;
wolffd@0:   else
wolffd@0:     index(c) = highest_num;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Permute the cliques so that they are ordered according to index
wolffd@0: [dummy, clique_order] = sort(index);
wolffd@0: cliques = cliques(clique_order);
wolffd@0: 
wolffd@0: w = zeros(num_cliques, 1); 
wolffd@0: B = sparse(num_cliques, 1);
wolffd@0: for i=1:num_cliques
wolffd@0:   B(i, cliques{i}) = 1;
wolffd@0:   w(i) = prod(ns(cliques{i}));
wolffd@0: end
wolffd@0: 
wolffd@0: % Pearl p113 suggests ordering the cliques by rank of the highest vertex in each clique.
wolffd@0: % However, this will only work if we use maximum cardinality search.
wolffd@0: 
wolffd@0: 
wolffd@0: % Join up the cliques so that they satisfy the Running Intersection Property.
wolffd@0: % This states that, for all k > 1, S(k) subseteq C(j) for some j < k, where
wolffd@0: % S(k) = C(k) intersect (union_{i=1}^{k-1} C(i))
wolffd@0: jtree = sparse(num_cliques, num_cliques);
wolffd@0: for k=2:num_cliques
wolffd@0:   S = [];
wolffd@0:   for i=1:k-1
wolffd@0:     S = myunion(S, cliques{i});
wolffd@0:   end
wolffd@0:   S = myintersect(S, cliques{k});
wolffd@0:   found = 0;
wolffd@0:   for j=1:k-1
wolffd@0:     if mysubset(S, cliques{j})
wolffd@0:       found = 1;
wolffd@0:       break;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if ~found
wolffd@0:     disp(['RIP is violated for clique ' num2str(k)]);
wolffd@0:   end
wolffd@0:   jtree(k,j)=1;
wolffd@0:   jtree(j,k)=1;
wolffd@0: end
wolffd@0: 
wolffd@0: % Pearl p113 suggests connecting Ci to a predecessor Cj (j < i) sharing
wolffd@0: % the highest number of vertices with Ci (i.e., the heaviest i-j edge
wolffd@0: % in the jgraph). However, this will only work if we use maximum cardinality search.
wolffd@0: 
wolffd@0: root = 1;
wolffd@0: 
wolffd@0: