wolffd@0: function [G, cliques, fill_ins] = triangulate(G, order)
wolffd@0: % TRIANGULATE Ensure G is triangulated (chordal), i.e., every cycle of length > 3 has a chord.
wolffd@0: % [G, cliques, fill_ins, cliques_containing_node] = triangulate(G, order)
wolffd@0: % 
wolffd@0: % cliques{i} is the i'th maximal complete subgraph of the triangulated graph.
wolffd@0: % fill_ins(i,j) = 1 iff we add a fill-in arc between i and j.
wolffd@0: %
wolffd@0: % To find the maximal cliques, we save each induced cluster (created by adding connecting
wolffd@0: % neighbors) that is not a subset of any previously saved cluster. (A cluster is a complete,
wolffd@0: % but not necessarily maximal, set of nodes.)
wolffd@0: 
wolffd@0: MG = G;
wolffd@0: n = length(G);
wolffd@0: eliminated = zeros(1,n);
wolffd@0: cliques = {};
wolffd@0: for i=1:n
wolffd@0:   u = order(i);
wolffd@0:   U = find(~eliminated); % uneliminated
wolffd@0:   nodes = myintersect(neighbors(G,u), U); % look up neighbors in the partially filled-in graph
wolffd@0:   nodes = myunion(nodes, u); % the clique will always contain at least u
wolffd@0:   G(nodes,nodes) = 1; % make them all connected to each other
wolffd@0:   G = setdiag(G,0);  
wolffd@0:   eliminated(u) = 1;
wolffd@0:   
wolffd@0:   exclude = 0;
wolffd@0:   for c=1:length(cliques)
wolffd@0:     if mysubset(nodes,cliques{c}) % not maximal
wolffd@0:       exclude = 1;
wolffd@0:       break;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if ~exclude
wolffd@0:     cnum = length(cliques)+1;
wolffd@0:     cliques{cnum} = nodes;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: fill_ins = sparse(triu(max(0, G - MG), 1));
wolffd@0: 
wolffd@0: %assert(check_triangulated(G)); % takes 72% of the time!
wolffd@0: 
wolffd@0: