wolffd@0: function order = strong_elim_order(G, node_sizes, partial_order)
wolffd@0: % STRONG_ELIM_ORDER Find an elimination order to produce a strongly triangulated graph.
wolffd@0: % order = strong_elim_order(moral_graph, node_sizes, partial_order)
wolffd@0: % 
wolffd@0: % partial_order(i,j)=1 if we must marginalize i *after* j
wolffd@0: % (so i will be nearer the strong root).
wolffd@0: % e.g., if j is a decision node and i is its information set:
wolffd@0: %   we cannot maximize j if we have marginalized out some of i
wolffd@0: % e.g., if j is a continuous child and i is its discrete parent:
wolffd@0: %   we want to integrate out the cts nodes before the discrete ones,
wolffd@0: %   so that the marginal is strong.
wolffd@0: %
wolffd@0: % For details, see
wolffd@0: % - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94.
wolffd@0: % - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical
wolffd@0: %   association models", JASA 87(420):1098--1108, 1992.
wolffd@0: %
wolffd@0: % On p369 of the Jensen paper, they state "the reverse of the elimination order must be some
wolffd@0: % extension of [the partial order] to a total order".
wolffd@0: % We make no attempt to find the best such total ordering, in the sense of minimizing the weight
wolffd@0: % of the resulting cliques.
wolffd@0: 
wolffd@0: % Example from the Jensen paper:
wolffd@0: % Let us number the nodes in Fig 1 from top to bottom, left to right,
wolffd@0: % so a=1,b=2,D1=3,c=4,...,l=14,j=15,k=16.
wolffd@0: % The elimination ordering they propose on p370 is [14 15 16 11 12 1 4 5 10 8 13 9 7 6 3 2];
wolffd@0: 
wolffd@0: if 0
wolffd@0:   total_order = topological_sort(partial_order);
wolffd@0:   order = total_order(end:-1:1); % no attempt to find an optimal constrained ordering!
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: % The following implementation is due to Ilya Shpitser and seems to give wrong
wolffd@0: % results on cg1
wolffd@0: 
wolffd@0: n = length(G);
wolffd@0: MG = G; % copy the original graph
wolffd@0: uneliminated = ones(1,n);
wolffd@0: order = zeros(1,n);
wolffd@0: 
wolffd@0: for i=1:n
wolffd@0:   roots = [];
wolffd@0:   k = 1;
wolffd@0:   for j=1:n
wolffd@0:     if sum(partial_order(j,:)) == 0
wolffd@0:       roots(k) = j;
wolffd@0:       k = k + 1;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   U = find(uneliminated);
wolffd@0:   valid = myintersect(U, roots);
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:   partial_order(:,k) = 0;
wolffd@0: end