Daniel@0: function [jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ...
Daniel@0:     dag_to_jtree(dag, node_sizes, partial_order, stages, clusters)
Daniel@0: % DAG_TO_JTREE Moralize and triangulate a DAG, and make a junction tree from its cliques.
Daniel@0: % [jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ...
Daniel@0: %    dag_to_jtree(dag, node_sizes, partial_order, stages, clusters)
Daniel@0: %
Daniel@0: % Input:
Daniel@0: % dag(i,j) 
Daniel@0: % jtree(i,j) = 1 iff there is an arc between clique i and clique j 
Daniel@0: % root = the root clique
Daniel@0: % cliques{i} = the nodes in clique i
Daniel@0: % B(i,j) = 1 iff node j occurs in clique i
Daniel@0: % w(i) = weight of clique i
Daniel@0: 
Daniel@0: N = length(bnet.dag);
Daniel@0: if nargin < 2, obs_nodes = []; end
Daniel@0: if nargin < 3, stages = { 1:N }; end
Daniel@0: if nargin < 4, clusters = {}; end
Daniel@0: 
Daniel@0: [MG, moral_edges]  = moralize(bnet.dag);
Daniel@0: 
Daniel@0: % Add extra arcs between nodes in each cluster to ensure they occur in the same clique
Daniel@0: for i=1:length(clusters)
Daniel@0:   c = clusters{i};
Daniel@0:   MG(c,c) = 1;
Daniel@0: end
Daniel@0: MG = setdiag(MG, 0);
Daniel@0: 
Daniel@0: % Find an optimal elimination ordering (NP-hard problem!)
Daniel@0: ns = bnet.node_sizes(:);
Daniel@0: ns(obs_nodes) = 1; % observed nodes have only 1 possible value
Daniel@0: partial_order = determine_elim_constraints(bnet, obs_nodes);
Daniel@0: 
Daniel@0: if isempty(partial_order)
Daniel@0:   strong = 0;
Daniel@0:   elim_order = best_first_elim_order(MG, ns, stages);
Daniel@0: else
Daniel@0:   strong = 1;
Daniel@0:   elim_order = strong_elim_order(MG, ns, partial_order);
Daniel@0: end
Daniel@0: 
Daniel@0: [MTG, cliques, fill_in_edges]  = triangulate(MG, elim_order);
Daniel@0: 
Daniel@0: % Connect the cliques up into a jtree,
Daniel@0: [jtree, root, B, w] = cliques_to_jtree(cliques, ns);
Daniel@0: 
Daniel@0: if 0
Daniel@0:   disp('testing dag to jtree');
Daniel@0:   % Find the cliques containing each node, and check they form a connected subtree
Daniel@0:   clqs_con_node = cell(1,N);
Daniel@0:   for i=1:N
Daniel@0:     clqs_con_node{i} = find(B(:,i))';
Daniel@0:   end
Daniel@0:   check_jtree_property(clqs_con_node, jtree);
Daniel@0: end