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 [http://en.wikipedia.org/wiki/Feedback_arc_set “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)
;;