function jtree = mk_jtree(bnet, cluster)
% MK_JTREE  Make a rooted junction tree whose nodes are the maximal cliques of triangulate(moralize(dag)).
% jtree = mk_jtree(bnet, clusters)
%
% Input:
% clusters is an optional argument, needed for the BK algorithm.
%
% Output:
% jtree is a structure with the following fields
%
% adjmat(i,j) = 1 if there is an arc from clique i to j, = 0 otherwise
% cliques{i} = the nodes in clique i
% separator{s} = the nodes in separator s
% sep_index(i,j)=s means the separator between cliques i and j is in s (since we cant have sparse cell arrays.)
% preorder = the cliques in the order parents before children
% postorder = the cliques in the order children before parents
%
% cliques_containing_node{k}: we can get the marginal on k from any of these cliques
% cl_containing_family(k): this is the lightest clique containing k's family
%
% nodes_assigned_to_clique{i}:
% A node can be a member of many cliques, but is assigned to exactly
% one to avoid over counting. We assign node i to clique c if c is the lightest clique containing i's family.
% c has to contain i's family so that c's domain is big enough to support i's CPT.
%
% The clique and separator potentials will be created by mk_ev_specific_pot.
%
% If the dag is not fully connected, some of the separators will be empty, but the join tree
% will always be a tree, not a forest.

if nargin < 2, cluster = []; end

num_bn_nodes = length(bnet.adj_mat);
[MG, moral_edges]  = moralize(bnet.adj_mat); % MG = moralize(graph)

if ~isempty(bnet.uedges)
  MG = or(MG, bnet.uedges);
end

% Make sure there are cliques that contain each cluster in slice t and t+1 (for BK)
if bnet.is_dbn
  ss = bnet.slice_size;
  num_clusters = length(cluster);
  for i=1:num_clusters
    c = cluster{i};
    MG(c,c) = 1;
    if bnet.num_slices > 1
      c = cluster{i} + ss;
      MG(c,c) = 1;
    end
  end
else
  assert(isempty(cluster));
end
% Make sure there are cliques that contain each desired marginal
for i=1:length(bnet.desired_marginal)
  m = bnet.desired_marginal{i};
  if bnet.is_dbn & bnet.num_slices == 1
    m = intersect(m, 1:ss);
  end
  MG(m,m) = 1;
end
MG = setdiag(MG,0);

% Since the (effective) size of a node can only decrease once it is observed,
% we might get a more optimal jtree if we wait until we know what nodes are observed.
ns = bnet.node_sizes(:);
ns(intersect(bnet.observed_nodes, bnet.discrete_nodes)) = 1; % only 1 possible value
ns(intersect(bnet.observed_nodes, bnet.cts_nodes)) = 0; % obs. vectors get replaced by a point

if bnet.is_dbn & strcmp(bnet.algo, 'jtree')
  % Eliminate all the nodes in slice t before slice t+1.
  % This will ensure that the jtree has a repeating structure (at least away from
  % both edges). See the discussion of splicing jtrees on p68 of
  % Geoff Zweig's PhD thesis, Dept. Comp. Sci., UC Berkeley, 1998.
  % This constraint can increase the clique size significantly,
  % so we don't enforce it on the BK algorithm.
  ss = bnet.slice_size;
  T = bnet.num_slices;
  stages = num2cell(unroll_set(1:ss, ss, T), 1);
else
  stages = {1:(length(bnet.adj_mat))}; % no constraints
end


strong_triangulation = 0;
if strong_triangulation
  %partial_order = bnet.adj_mat'; % eliminate in reverse topological order
  n = length(bnet.adj_mat);
  partial_order = zeros(n, n);
  % discrete before cts
  partial_order(bnet.discrete_nodes, bnet.cts_nodes) = 1;
  jtree = mk_strong_jtree_topology(MG, ns, partial_order);
else
  elim_order = best_first_elim_order(MG, ns, stages);
  jtree = mk_jtree_topology(MG, ns, elim_order);
end
jtree.moral_edges = moral_edges;

num_cliques = length(jtree.cliques);

jtree.sep_index = sparse(num_cliques, num_cliques); % initializes to 0s
[is,js] = find(jtree.adj_mat > 0);
s = 1;
for k=1:length(is)
  i = is(k); j = js(k);
  jtree.sep_index(i,j) = s;
  jtree.separator{s} = intersect(jtree.cliques{i}, jtree.cliques{j});
  s = s + 1;
end

jtree.cl_containing_family = zeros(1, num_bn_nodes);
jtree.nodes_assigned_to_clique = cell(1,num_cliques);
for i=1:num_bn_nodes
  [c, cls] = find_clique(jtree, family(bnet.adj_mat, i));
  jtree.cl_containing_family(i) = c; 
  jtree.nodes_assigned_to_clique{c} = union(jtree.nodes_assigned_to_clique{c}, i);
end

