Finding all the elementary circuits of a directed graph

Below you can find an implementation in ocaml of the algorithm described in this paper by D. B. Johnson

  Finding all the elementary circuits of a directed graph.
  D. B. Johnson, SIAM Journal on Computing 4, no. 1, 77-84, 1975.
  http://dx.doi.org/10.1137/0204007
the idea of this algorithm is to enumerate all cycles in a strongly connected component. Next step will be to implement the "feedback arc set" of this connected component so to find the optimal way to break these loops in a way that the connected component can be partitioned in smaller DAGs.

This is only the main function, you can find everything else in the git repo given below. If I've time I'll also give a shot to implement tarjan algorithm and to compare the two... Johnson is supposed to be better complexity wise, but I'm curious to check the difference on my domain specific input.

let find_all_cycles_johnson g =
  if not G.is_directed then
    assert false;

  (*  stack of nodes in current path *)
  let path = Stack.create () in

  (* main function. return (true, acc) if in the last iteration we find a new cycle *)
  let rec circuit t result thisnode startnode component =

   (* add the node the candate cycles and block it *)
    Stack.push thisnode path;
    block t thisnode;

    let (closed,result) =
      G.fold_succ (fun nextnode (c,r) ->
        if G.V.equal nextnode startnode then
         (* this is a cycle ! *)
          (true,(Stack.copy path)::r)
        else begin
          if not(is_bloked t nextnode) then
            circuit t r nextnode startnode component
          else
            (c,r)
        end
      ) component thisnode (false,result)
    in
    if closed then
      (* this node is part of a cycles, but we don't want to exclude it from subsequent
          iterations. Longer cycles can contain smaller cycles *)

      unblock t thisnode
    else
      (* since it is not part of a cycles, then it must be part of a portion of the graph that
          does not contain any elementary circuits *)

      G.iter_succ (fun nextnode ->
        let l = get_notelem t nextnode in
        if List.mem thisnode l then
          Hashtbl.replace t.notelem nextnode (thisnode::l)
      ) component thisnode;

    ignore(Stack.pop path);
    (closed,result)
  in

  (* Johnson's algorithm requires some ordering of the nodes.
      We use the order induced by the compare function associated to the
      vertex type. See ocamlgraph Pack.Diagraph  *)

  let vertex_set = G.fold_vertex SV.add g SV.empty in
 
  SV.fold (fun s result ->
    (* Build the subgraph induced by s and following nodes in the ordering *)
    let subset = SV.add s (partition vertex_set s) in
    let subgraph = extract_subgraph g subset in

    if G.nb_edges subgraph > 0 then begin
      (* Find the strongly connected component in the subgraph
       * that contains the least node according to the ordering *)

      let scc = G.Components.scc_list subgraph in
      let minnode = SV.min_elt subset in
      let mincomp = List.find (fun l -> List.mem minnode l) scc in
     
      (* smallest node in the component according to the ordering *)
      let startnode = minnode in
      let component = extract_subgraph subgraph (to_set mincomp) in

      (* init the block table for this component *)
      let t = init_block component in

      snd(circuit t result startnode startnode component);
    end else
      result
  ) vertex_set []
;;

You can get the code here : git clone http://mancoosi.org/~abate/repos/cycles.git . There are two versions of the algorithm. One is a translation in ocaml of the imperative algorithm, the other, presented above, adds a bit of functional flavor ...

Update 1

the code repository now contains also two implementations of approximate algorithms the feedback arc set problem

Update 2

Since the number of cycles in a graph can be very large, I've added a lazy version of the algorithm using ExtLib enums (but can be easily done without). It would be cool to avoid using a global data structure to keep track of which nodes I've already visited like the persistent array of filliatr ... the code could be further optimized of course as many operations on lists and set can be avoided with a smarter DS to build the connected component of the subgraph at each iteration ... Complete code in the git repo above.

let find_all_cycles_enum g =
  let rec circuit path t thisnode startnode component =                                        
     let rec aux = function                                                                    
       |[] -> Enum.empty ()                                                                    
       |nextnode :: rest ->
         if G.V.equal nextnode startnode then begin                                            
           let e = aux rest in                                                                
           unblock t thisnode;
           Enum.push e (List.rev (thisnode::path));                                            
           e
         end else
           if not(is_bloked t nextnode) then begin
             let e = circuit (thisnode::path) t nextnode startnode component in                
             Enum.append e (aux rest)                                                          
           end else begin                                                                      
             aux rest                                                                          
           end                                                                                
     in
     block t thisnode;
     let e = aux (G.succ component thisnode) in                                                
     G.iter_succ (fun nextnode ->
       let l = get_notelem t nextnode in                                                      
       if List.mem thisnode l then
         Hashtbl.replace t.notelem nextnode (thisnode::l)                                      
     ) component thisnode;                                                                    
     e                                                                                        
  in                                                                                          
 
  let vertex_set = G.fold_vertex SV.add g SV.empty in                                          
  let rec aux = function                                                                      
    |[] -> Enum.empty ()                                                                      
    |s :: rest ->  
      let subset = SV.add s (partition vertex_set s) in                                        
      let subgraph = extract_subgraph g subset in                                              
      (* I need only one... not all scc *)
      (* actually I need only the scc that contains the min *)                                
      let scc = G.Components.scc_list subgraph in                                              
      let minnode = SV.min_elt subset in
      let mincomp = List.find (fun l -> List.mem minnode l) scc in                            
      let startnode = minnode in
      let component = extract_subgraph subgraph (to_set mincomp) in                            
      let t = init_block component in
      let partial = circuit [] t startnode startnode component in                              
      Enum.append partial (aux rest)                                                          
  in
  aux (SV.elements vertex_set)                                                                
;;

Average: 2.3 (4 votes)