Finding elementary cycles of (directed) graphs

Question

I need to enumerate all the simple cycles (i.e. elementary cycles where no vertex is repeated other than the starting) of a graph which has both directed and undirected edges, where we can treat the undirected edges as doubly directed. (Specifically, I am looking at the Cayley graphs of S3 and S4, which can be produced using CayleyGraph[SymmetricGroup[3]] and CayleyGraph[SymmetricGroup[4]] respectively.)

I have tried two ways of doing this so far. First, I have tried using ExtractCycles in the Combinatorica package, as detailed in this answer by TomD. For example, entering the ordered pairs for S3 (as "el"):

el = {{1, 2}, {1, 4}, {2, 3}, {2, 5}, {3, 1}, {3, 6}, {4, 1}, {4, 6}, {5, 4}, {5, 2}, 
{6, 5}, {6, 3}}

And then using:

Needs["Combinatorica`"]
ExtractCycles@FromOrderedPairs@el

returns:

{{5, 4, 6, 5}, {6, 3, 6}, {5, 2, 5}, {4, 1, 4}, {3, 1, 2, 3}}

However, that set is incomplete - what about e.g. {{5, 4, 1, 2, 5}, {5, 4, 1, 2, 3, 6, 5}}? These are simple cycles, so why are they not included in the list? The output for S4 is also much truncated (obviously we would expect a lot of cycles there, but the output gives only eighteen).

The second thing I tried was this answer by kguler. Taking the example of S3, using CycleGraph[3, DirectedEdges -> True] gave the right results, but CycleGraph[4, DirectedEdges -> True] did not - i.e. the cycles {5, 4, 1, 2, 5} etc. were not picked up again. Presumably this is something to do with the doubly directed edges?

Any help with this would be much appreciated!

Edit: As requested, the ordered pairs for the Cayley graph of S4 are:

el2 = {{1, 2}, {1, 9}, {2, 3}, {2, 17}, {3, 4}, {3, 13}, {4, 1}, {4, 5}, {5, 4}, {5, 6},
{6, 16}, {6, 7}, {7, 8}, {7, 22}, {8, 5}, {8, 10}, {9, 1}, {9, 10}, {10, 8}, {10, 11},
{11, 12}, {11, 21}, {12, 9}, {12, 18}, {13, 14}, {13, 3}, {14, 15}, {14, 20}, {15, 16},
{15, 23}, {16, 13}, {16, 6}, {17, 2}, {17, 18}, {18, 19}, {18, 12}, {19, 24}, {19, 20},
{20, 17}, {20, 14}, {21, 11}, {21, 22}, {22, 7}, {22, 23}, {23, 15}, {23, 24}, {24, 19},
{24, 21}}

Answer 1

Here is a brute force method:

cycles[el_] :=
 Module[{f, edges = Rule @@@ el // Dispatch},
  f[x_, b___, x_] := {{x, b, x}};
  f[___, x_, ___, x_] = {};
  f[c___, v_] := Join @@ (f[c, v, #] & /@ ReplaceList[v, edges]);
  Join @@ f /@ Union @@ el
 ]

In the code above the line f[___, x_, ___, x_] = {}; was used for clarity, but faster duplicate tests exist. For short cycles f[c__] /; ! UnsameQ@c = {}; should be fast(est), but if long cycles may be present you should use f[c__] /; Signature@{c} === 0 = {}; or f[c__] /; {c} =!= DeleteDuplicates@{c} = {};

Test:

el = {{1, 2}, {1, 4}, {2, 3}, {2, 5}, {3, 1}, {3, 6}, {4, 1}, {4, 6},
      {5, 4}, {5, 2}, {6, 5}, {6, 3}};

cycles[el]

{{1, 2, 3, 1}, {1, 2, 3, 6, 5, 4, 1}, {1, 2, 5, 4, 1}, {1, 2, 5, 4, 6,
   3, 1}, {1, 4, 1}, {1, 4, 6, 5, 2, 3, 1}, {1, 4, 6, 3, 1}, {2, 3, 1,
   2}, {2, 3, 1, 4, 6, 5, 2}, {2, 3, 6, 5, 4, 1, 2}, {2, 3, 6, 5, 
  2}, {2, 5, 4, 1, 2}, {2, 5, 4, 6, 3, 1, 2}, {2, 5, 2}, {3, 1, 2, 
  3}, {3, 1, 2, 5, 4, 6, 3}, {3, 1, 4, 6, 5, 2, 3}, {3, 1, 4, 6, 
  3}, {3, 6, 5, 4, 1, 2, 3}, {3, 6, 5, 2, 3}, {3, 6, 3}, {4, 1, 2, 3, 
  6, 5, 4}, {4, 1, 2, 5, 4}, {4, 1, 4}, {4, 6, 5, 4}, {4, 6, 5, 2, 3, 
  1, 4}, {4, 6, 3, 1, 2, 5, 4}, {4, 6, 3, 1, 4}, {5, 4, 1, 2, 3, 6, 
  5}, {5, 4, 1, 2, 5}, {5, 4, 6, 5}, {5, 4, 6, 3, 1, 2, 5}, {5, 2, 3, 
  1, 4, 6, 5}, {5, 2, 3, 6, 5}, {5, 2, 5}, {6, 5, 4, 1, 2, 3, 6}, {6, 
  5, 4, 6}, {6, 5, 2, 3, 1, 4, 6}, {6, 5, 2, 3, 6}, {6, 3, 1, 2, 5, 4,
   6}, {6, 3, 1, 4, 6}, {6, 3, 6}}

To remove the duplicate cycles from this output one can use:

DeleteDuplicates[RotateLeft[#, Ordering[#, 1] - 1] & /@ Most /@ #]&

{{1, 2, 3}, {1, 2, 3, 6, 5, 4}, {1, 2, 5, 4}, {1, 2, 5, 4, 6, 3}, {1, 4},
 {1, 4, 6, 5, 2, 3}, {1, 4, 6, 3}, {2, 3, 6, 5}, {2, 5}, {3, 6}, {4, 6, 5}}

Here is an alternative approach to remove the duplicates on-the-fly with a form of memoization. It seems it may be more or less efficient depending on the graph; I have not tested it well yet.

cycles2[el_] :=
 Module[{f, edges = Rule @@@ el // Dispatch},
  f[x_, b___, x_] := (
    (f[##] = {}) & @@@ NestList[RotateLeft, {x, b}, Length@{x, b} - 1];
    {{x, b}}
   );
  f[c__] /; Signature@{c} === 0 = {};
  f[c___, v_] := Join @@ (f[c, v, #] & /@ ReplaceList[v, edges]);
  Join @@ f /@ Union @@ el
  ]

cycles2[el]

{{1, 2, 3}, {1, 2, 3, 6, 5, 4}, {1, 2, 5, 4}, {1, 2, 5, 4, 6, 3}, {1, 4},
 {1, 4, 6, 5, 2, 3}, {1, 4, 6, 3}, {2, 3, 6, 5}, {2, 5}, {3, 6}, {4, 6, 5}}

Answer 2

Below is an implementation of Johnson's 1975 exhaustive algorithm (see PDF, AFAIK the fastest exhaustive algorithm), improved upon the rather procedural version of Daniel Skates (see Mathematica demonstration). A hand-crafted C-version of the code is also available (if you mail me), which adds a further tenfold increase of speed compared to the Mathematica version. Code is at end of post (findAllCycles).

Note:

• the input should be either a graph or a binary adjacency matrix (SparseArray is ok);
• self-loops are allowed (Johnson's algorithm is unable to correctly find all, so they are assessed beforehand and the diagonal of the adjacency matrix is zeroed for the algorithm);
• multiple edges (other than doubly directed ones) or wheights are not allowed.

Result comparison:

el = {{1, 2}, {1, 4}, {2, 3}, {2, 5}, {3, 1}, {3, 6},
      {4, 1}, {4, 6}, {5, 4}, {5, 2}, {6, 5}, {6, 3}};
g = Graph[DirectedEdge @@@ el, DirectedEdges -> True, 
  VertexLabels -> "Name", ImagePadding -> 10, ImageSize -> 300]

mrW = DeleteDuplicates[RotateLeft[#, Ordering[#, 1] - 1] &@Most@# & /@ (cycles@el)];
z = findAllCycles@g

Sort@mrW === Sort@z

{{1, 2, 3}, {1, 2, 3, 6, 5, 4}, {1, 2, 5, 4}, {1, 2, 5, 4, 6, 3}, {1, 4},
 {1, 4, 6, 3}, {1, 4, 6, 5, 2, 3}, {2, 3, 6, 5}, {2, 5}, {4, 6, 5}, {3, 6}}

True

Time comparison, by generating 1000 weakly connected directed graphs (since MrWizard's method is redundant, Johnson's method is expected to be faster):

iteration = 1000;
nodeRange = {9, 12};
edgeProb = .3;
adjList = elList = {};
While[
  n = RandomInteger@nodeRange;
  adj = RandomChoice[{1-edgeProb, edgeProb} -> {0, 1}, {n, n}];
  (* only weakly conncected graphs should be checked *)
  Length@adjList < iteration && ConnectedGraphQ@AdjacencyGraph[adj, DirectedEdges -> False],
  (* create adjacency matrix and edge list separately for the two functions *)
  AppendTo[adjList, adj];
  AppendTo[elList, Position[Normal@adj, 1, {2}]];
  ];

AbsoluteTiming[
 mrW = (DeleteDuplicates[RotateLeft[#, Ordering[#, 1] - 1] &@
          Most@# & /@ (cycles@#)]) & /@ elList;]
AbsoluteTiming[z = findAllCycles /@ adjList;]

{23.590349, Null}   (* MrW *)
{3.780216, Null}    (* Z *)

Note that cycles and findAllCycles return results in different orders, so I omitted the direct comparison and only the number of found cycles are checked:

Length /@ mrW === Length /@ z

True

Code:

(* wrappers for Graph and SparseArray input *)
findAllCycles[g_Graph] := Module[{nodes = VertexList@g},
   Replace[findAllCycles@Normal@AdjacencyMatrix@g, 
    Thread[Range@Length@nodes -> nodes], {2}]
   ];
findAllCycles[s_SparseArray] := findAllCycles@Normal@s;

(* Johnson's algorithm *)
findAllCycles[a_?MatrixQ] := 
  Module[{unblock, circuit, n = Length@a, adj, loops, AK, s = 1, 
    stack = {}, circuitsFound = {}, B, blocked},

   (* remove self-loops, as Johnson's algorithm is unable to deal with them. *)
   loops = Position[Diagonal@a, 1];
   adj = a*(1 - IdentityMatrix@n);

   B = Table[{}, {n}];
   blocked = Table[False, {n}];

   unblock[u_Integer] := Module[{w},
     blocked[[u]] = False;
     While[Length@B[[u]] > 0,
      w = B[[u, 1]];
      B[[u]] = Drop[B[[u]], 1];
      If[blocked[[w]], unblock@w]];
     ];

   circuit[v_Integer] := Module[{f = False},
     stack = Append[stack, v];
     blocked[[v]] = True;
     Do[
      If[w == s,
       circuitsFound = Append[circuitsFound, stack]; f = True,
       If[! blocked[[w]] && circuit@w, f = True]
       ], {w, AK[[v]]}];
     If[f,
      unblock@v
      ,
      Do[If[FreeQ[B[[w]], v], B[[w]] = Join[B[[w]], {v}]], {w, AK[[v]]}];
      ];
     stack = Drop[stack, -1];
     f
     ];

   While[s < n,
    AK = Take[adj, {s, n}, {s, n}];
    AK = Flatten@Position[#, 1] & /@ AK;
    AK = AK /. {x_ :> (x + s - 1)};
    AK = Join[Table[{}, {s - 1}], AK];
    If[AK === {},
     s = n,
     Do[blocked[[i]] = False; B[[i]] = {}, {i, s, n}];
     circuit@s;
     s = s + 1;
     ]
    ];
   Join[loops, circuitsFound]
   ];

Answer 3

Mostly I cribbed this from here. I simply changed the line that converts input to the particular sparse representation used by the main part of the code.

Important caveat: I do not know for a fact that this code is correct for the directed case.

extendCycle[cyc_List, edges_List] := 
 Map[If[# > First[cyc] && ! MemberQ[cyc, #], Append[cyc, #], 
    Null /. Null :> Sequence[]] &, edges[[Last[cyc]]]]

cycles[mat_, k_] := Module[{n = Length[mat], m2, cyc, cyclist},
  m2 = Map[Last, Split[Sort[mat], First[#1] == First[#2] &], {2}] ;
  cyclist = 
   Flatten[Drop[MapIndexed[{#2[[1]], #1} &, m2, {2}], -k + 1], 1];
  Do[cyclist = 
    Flatten[Map[extendCycle[#, m2] &, cyclist], 1], {k - 2}];
  Map[If[MemberQ[m2[[Last[#]]], First[#]], Append[#, First[#]], 
     Null /. Null :> Sequence[]] &, cyclist]]

Your simpler example:

el = {{1, 2}, {1, 4}, {2, 3}, {2, 5}, {3, 1}, {3, 6}, {4, 1}, {4, 
    6}, {5, 4}, {5, 2}, {6, 5}, {6, 3}};

Timing[
 Table[cycles[el, j], {j, 2, Length[Union[Flatten[el]]]}] /. {} :> 
   Sequence[]]

(* Out[436]= {0., {{{1, 4, 1}, {2, 5, 2}, {3, 6, 3}, {4, 1, 4}, {5, 2, 
    5}}, {{1, 2, 3, 1}, {4, 6, 5, 4}}, {{1, 2, 5, 4, 1}, {1, 4, 6, 3, 
    1}, {2, 3, 6, 5, 2}, {3, 1, 4, 6, 3}}, {{1, 2, 3, 6, 5, 4, 1}, {1,
     2, 5, 4, 6, 3, 1}, {1, 4, 6, 5, 2, 3, 1}}}} *)

Your bigger case:

el2 = {{1, 2}, {1, 9}, {2, 3}, {2, 17}, {3, 4}, {3, 13}, {4, 1}, {4, 
    5}, {5, 4}, {5, 6}, {6, 16}, {6, 7}, {7, 8}, {7, 22}, {8, 5}, {8, 
    10}, {9, 1}, {9, 10}, {10, 8}, {10, 11}, {11, 12}, {11, 21}, {12, 
    9}, {12, 18}, {13, 14}, {13, 3}, {14, 15}, {14, 20}, {15, 
    16}, {15, 23}, {16, 13}, {16, 6}, {17, 2}, {17, 18}, {18, 
    19}, {18, 12}, {19, 24}, {19, 20}, {20, 17}, {20, 14}, {21, 
    11}, {21, 22}, {22, 7}, {22, 23}, {23, 15}, {23, 24}, {24, 
    19}, {24, 21}};

Timing[
 cycs2 = Table[
     cycles[el2, j], {j, 2, Length[Union[Flatten[el2]]]}] /. {} :> 
     Sequence[];]

(* Out[439]= {0.340000, Null} *)

Now check sizes.

Length[cycs2]

(* Out[440]= 10 *)

Length[Flatten[cycs2, 1]]

(* Out[443]= 199 *)

Upshot: Not too bad in performance, if it happens to be correct.

Answer 4

One approach, born out of curiosity more than anything, is to take the determinant, extract the permutations from the indices, and then cycle decompose. I don't think this is a very efficient method as Det is slow, but it seems to work.

To cycle decompose, I use combinatoricatocycles[] with Mathematica 7 (see here), rather than the new Mathematica 8 function PermutationCycles

detPerms[ei_] := List @@@ (List @@ (Det@SparseArray[# -> 
            a[#] & /@ ei]) /. {-x_ -> x}) /. a[{x_, y_}] -> y

Applying to el:

Union@Flatten[combinatoricatocycles /@ #, 1] &@detPerms[el]

(* {{4, 1}, {5, 2}, {6, 3}, {2, 3, 1}, {6, 5, 4}, {2, 5, 4, 1}, {3, 6, 5, 2}, {4, 6, 3, 1}, {2, 3, 6, 5, 4, 1}, {2, 5, 4, 6, 3, 1}, {4, 6, 5, 2, 3, 1}} *)

And to el2:

Length@Union@Flatten[combinatoricatocycles /@ #, 1] &@detPerms[el2]

(* 180 *)

combinatoricatocycles (see here)

combinatoricatocycles[p_] := 
 Module[{k, j, first, np = p, q = Table[0, {Length[p]}], i}, 
  DeleteCases[
   Table[If[np[[i]] == 0, {}, j = 1; first = np[[i]]; np[[i]] = 0;
     k = q[[j++]] = first;
     While[np[[k]] != 0, q[[j++]] = np[[k]]; np[[k]] = 0;
      k = q[[j - 1]]];
     Take[q, j - 1]], {i, Length[p]}], _?(# === {} &)]]