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