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