function jtree = mk_strong_jtree_topology(MG, ns, partial_order)
% MK_STRONG_JTREE_TOPOLOGY Make a strong junction tree.
% jtree = mk_strong_jtree_topology(MG, ns, partial_order)
%
% partial_order(i,j)=1 if we must marginalize j *before* i
% (so i will be nearer the strong root).
% e.g., if j is a decision node and i is its information set:
%   we cannot maximize j if we have marginalized out some of i
% e.g., if j is a continuous child and i is its discrete parent:
%   we want to integrate out the cts nodes before the discrete ones
%
% For details, see
% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
% - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical
%   association models", JASA 87(420):1098--1108, 1992.

% Here is a definition from Jensen et al.
% "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."
%
% A graph with different kinds of nodes (e.g., discrete and cts) is called marked.
% A strong root is guaranteed to exist if the marked graph is triangulated and does not have any paths of
% the form discrete -> cts -> discrete. In general we need to add extra edges to
% the moral graph to ensure this (see example in Lauritzen fig 3b).
% However, a simpler sufficient condition is to eliminate all the cts nodes before the discrete ones,
% because then, as we move from the leaves to the root, the cts nodes get marginalized away
% and we are left with purely discrete cliques.

total_order = topsort(partial_order);
elim_order = reverse(total_order); % no attempt to find an optimal constrained ordering!
jtree.elim_order = elim_order;
% e.g., from the Jensen et al paper
% elim_order = [14 15 16 11 12 1 4 5 10 8 13 9 7 6 3 2];

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
% alpha = [11 1 2 10 9 3 4 7 5 8 13 12 6 16 15 14]

[MTG, jtree.cliques, fill_in_edges]  = triangulate(MG, elim_order); % MTG = triangulate(MG)
num_cliques = length(jtree.cliques);

index = zeros(1, num_cliques);
% We sort the cliques in order of increasing index. The index of a clique C is defined as follows:
% let v in C be the highest-numbered vertex s.t. the vertices in W = lower intersect C,
% where lower = {u | alpha(u) < alpha(v)},
% 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.

for c = 1:num_cliques
  C = jtree.cliques{c};
  highest_num = -inf;
  for vi = 1:length(C)
    v = C(vi);
    lower = find(alpha < alpha(v));
    W = intersect(lower, C);
    U = setdiff(lower, C);
    found = 0;
    for ui=1:length(U)
      u = U(ui);
      if subset(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

% 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.
if 0
  for i=1:num_cliques
    index(i) = max(order(jtree.cliques{i}));
  end
end

% Permute the cliques so that they are ordered according to index
[index2, clique_order] = sort(index);
cliques = jtree.cliques;
for c=1:num_cliques
  jtree.cliques{c} = cliques{clique_order(c)};
end


jtree.adj_mat = sparse(num_cliques, num_cliques);

num_bn_nodes = length(MG);
[jgraph, jgraph2, jtree.clique_sizes] = jgraph_costs(jtree.cliques, num_bn_nodes, ns);

% 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 join graph). However, this will only work if we use maximum cardinality search.
if 0
  for i=2:num_cliques
    c = clique_order(i);
    lns = intersect(neighbors(jgraph, c), clique_order(1:i-1)); 
    cc = lns(1);
    jtree.adj_mat(cc,c) = 1;
  end
end

% 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)
for k=2:num_cliques
  S = [];
  for i=1:k-1
    S = union(S, jtree.cliques{i});
  end
  S = intersect(S, jtree.cliques{k});
  found = 0;
  for j=1:k-1
    if subset(S, jtree.cliques{j})
      found = 1;
      break;
    end
  end
  if ~found
    disp(['RIP is violated for clique ' num2str(k)]);
  end
  jtree.adj_mat(k,j)=1;
end

for i=1:num_cliques
  for j=1:num_cliques
    if jtree.adj_mat(i,j)==1
      jtree.adj_mat(j,i)=1;
    end
  end
end

jtree.cliques_containing_node = inverse_relation(jtree.cliques, num_bn_nodes);
check_jtree_property(jtree.cliques_containing_node, jtree.adj_mat);
root = 1; % This will be a strong root by construction
directed = 0;
% Although the jtree is directed, to do a dfs down from the root we must pretend it is undirected.
[d, jtree.preorder, jtree.postorder] = dfs(jtree.adj_mat, root, directed);
