wolffd@0: function order = best_first_elim_order(G, node_sizes, stage)
wolffd@0: % BEST_FIRST_ELIM_ORDER Greedily search for an optimal elimination order.
wolffd@0: % order = best_first_elim_order(moral_graph, node_sizes)
wolffd@0: %
wolffd@0: % Find an order in which to eliminate nodes from the graph in such a way as to try and minimize the
wolffd@0: % weight of the resulting triangulated graph.  The weight of a graph is the sum of the weights of each
wolffd@0: % of its cliques; the weight of a clique is the product of the weights of each of its members; the
wolffd@0: % weight of a node is the number of values it can take on.
wolffd@0: %
wolffd@0: % Since this is an NP-hard problem, we use the following greedy heuristic:
wolffd@0: % at each step, eliminate that node which will result in the addition of the least
wolffd@0: % number of fill-in edges, breaking ties by choosing the node that induces the lighest clique.
wolffd@0: % For details, see
wolffd@0: % - Kjaerulff, "Triangulation of graphs -- algorithms giving small total state space",
wolffd@0: %      Univ. Aalborg tech report, 1990 (www.cs.auc.dk/~uk)
wolffd@0: % - C. Huang and A. Darwiche, "Inference in Belief Networks: A procedural guide",
wolffd@0: %      Intl. J. Approx. Reasoning, 11, 1994
wolffd@0: %
wolffd@0: 
wolffd@0: % Warning: This code is pretty old and could probably be made faster.
wolffd@0: 
wolffd@0: n = length(G);
wolffd@0: if nargin < 3, stage = { 1:n }; end % no constraints
wolffd@0: 
wolffd@0: % For long DBNs, it may be useful to eliminate all the nodes in slice t before slice t+1.
wolffd@0: % This will ensure that the jtree has a repeating structure (at least away from both edges).
wolffd@0: % This is why we have stages.
wolffd@0: % See the discussion of splicing jtrees on p68 of
wolffd@0: % Geoff Zweig's PhD thesis, Dept. Comp. Sci., UC Berkeley, 1998.
wolffd@0: % This constraint can increase the clique size significantly.
wolffd@0: 
wolffd@0: MG = G; % copy the original graph
wolffd@0: uneliminated = ones(1,n);
wolffd@0: order = zeros(1,n);
wolffd@0: t = 1;  % Counts which time slice we are on        
wolffd@0: for i=1:n
wolffd@0:   U = find(uneliminated);
wolffd@0:   valid = myintersect(U, stage{t});
wolffd@0:   % Choose the best node from the set of valid candidates
wolffd@0:   score1 = zeros(1,length(valid));
wolffd@0:   score2 = zeros(1,length(valid));
wolffd@0:   for j=1:length(valid)
wolffd@0:     k = valid(j);
wolffd@0:     ns = myintersect(neighbors(G, k), U);
wolffd@0:     l = length(ns);
wolffd@0:     M = MG(ns,ns);
wolffd@0:     score1(j) = l^2 - sum(M(:)); % num. added edges
wolffd@0:     score2(j) = prod(node_sizes([k ns])); % weight of clique
wolffd@0:   end
wolffd@0:   j1s = find(score1==min(score1));
wolffd@0:   j = j1s(argmin(score2(j1s)));
wolffd@0:   k = valid(j);
wolffd@0:   uneliminated(k) = 0;
wolffd@0:   order(i) = k;
wolffd@0:   ns = myintersect(neighbors(G, k), U);
wolffd@0:   if ~isempty(ns)
wolffd@0:     G(ns,ns) = 1;
wolffd@0:     G = setdiag(G,0);
wolffd@0:   end
wolffd@0:   if ~any(logical(uneliminated(stage{t}))) % are we allowed to the next slice?
wolffd@0:     t = t + 1;
wolffd@0:   end   
wolffd@0: end
wolffd@0: