view toolboxes/FullBNT-1.0.7/graph/cliques_to_strong_jtree.m @ 0:e9a9cd732c1e tip

first hg version after svn

function [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG)
% MK_SRONG_JTREE Make a strong junction tree.
% [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG)
%
% Here is a definition of a strong jtree from Jensen et al. 1994:
% "A junction tree is said to be strong if it has at least one distinguished clique R,
% called a strong root, s.t. for each pair (C1,C2) of adjacent cliques in the tree,
% with C1 closer to R than C2, there exists and ordering of [the nodes below] C2
% that respects [the partial order] and with the vertices of the separator C1 intersect C2
% preceeding the vertices [below C2] of C2 \ C1."
%          
% For details, see
% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
%
% MTG is the moralized, triangulated graph.
% elim_order is the elimination ordering used to compute MTG.


% Warning: this is a very naive implementation of the algorithm in Jensen et al.

n = length(elim_order);
alpha(elim_order) = n:-1:1;
% alpha(u) = i if we eliminate u at step n-i+1
% i.e., vertices with higher alpha numbers are eliminated before vertices with lower numbers.
% e.g., from the Jensen et al paper
% node a=1 eliminated at step 6, so alpha(a)=16-6+1=11.
% alpha = [11 1 2 10 9 3 4 7 5 8 13 12 6 16 15 14]


% We sort the cliques in order of increasing index. The index of a clique C is defined as follows.
% Let lower = {u | alpha(u) < alpha(v)}, and
% let v in C be the highest-numbered vertex s.t. the vertices in W = lower intersect C
% have a common neighbor u in U, where U = lower \ C.
% If such a v exists, define index(C) = alpha(v), otherwise, index(C) = 1.
% Intuitively, index(C) is the step in the elimination process at which C disappears.

num_cliques = length(cliques);
index = zeros(1, num_cliques);
for c = 1:num_cliques
  C = cliques{c};
  highest_num = -inf;
  for vi = 1:length(C)
    v = C(vi);
    lower = find(alpha < alpha(v));
    W = myintersect(lower, C);
    U = mysetdiff(lower, C);
    found = 0;
    for ui=1:length(U)
      u = U(ui);
      if mysubset(W, neighbors(MTG, u))
	found = 1;
	break;
      end
    end
    if found
      if alpha(v) > highest_num
	highest_num = alpha(v);
      end
    end
  end
  if highest_num == -inf
    index(c) = 1;
  else
    index(c) = highest_num;
  end
end


% Permute the cliques so that they are ordered according to index
[dummy, clique_order] = sort(index);
cliques = cliques(clique_order);

w = zeros(num_cliques, 1); 
B = sparse(num_cliques, 1);
for i=1:num_cliques
  B(i, cliques{i}) = 1;
  w(i) = prod(ns(cliques{i}));
end

% Pearl p113 suggests ordering the cliques by rank of the highest vertex in each clique.
% However, this will only work if we use maximum cardinality search.


% Join up the cliques so that they satisfy the Running Intersection Property.
% This states that, for all k > 1, S(k) subseteq C(j) for some j < k, where
% S(k) = C(k) intersect (union_{i=1}^{k-1} C(i))
jtree = sparse(num_cliques, num_cliques);
for k=2:num_cliques
  S = [];
  for i=1:k-1
    S = myunion(S, cliques{i});
  end
  S = myintersect(S, cliques{k});
  found = 0;
  for j=1:k-1
    if mysubset(S, cliques{j})
      found = 1;
      break;
    end
  end
  if ~found
    disp(['RIP is violated for clique ' num2str(k)]);
  end
  jtree(k,j)=1;
  jtree(j,k)=1;
end

% Pearl p113 suggests connecting Ci to a predecessor Cj (j < i) sharing
% the highest number of vertices with Ci (i.e., the heaviest i-j edge
% in the jgraph). However, this will only work if we use maximum cardinality search.

root = 1;