# Finding elementary cycles of (directed) graphs

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}}

• I've added the ordered pairs for the Cayley graph of S3, as well as the code to produce the Cayley graphs. Is this what's needed? Aug 15, 2012 at 8:10
• Yes, thanks. If you don't mind including the S4 data too that would be helpful for version 7 users like me, who don't have CayleyGraph. Aug 15, 2012 at 8:12
• Sure. Just added them at the bottom. Aug 15, 2012 at 8:19
• Would you tell me if Daniel's method is giving the right output? It does not agree with the output from mine. Aug 16, 2012 at 13:28
• I think his method gives the right output, but there is some redundancy (for example, {1,4,1} and {4,1,4} are both included for the graph of S3). Aug 19, 2012 at 14:20

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}}

• That's great, thanks! I have two questions though: 1) That list is complete, but it now gives multiple entries for the same cycle, in that it gives all the possible cyclic permutations of each cycle. For example, {1,2,3,2} and {3,1,2,3} are both included, but they are, in fact, the same cycle. Is there any way of only including one entry for each cycle? 2) Out of interest, what is the complexity of this algorithm? How long would you expect it to take to compute all the cycles for S4? Anyway, thanks very much - that's really useful! Aug 15, 2012 at 9:39
• Oh, I just tried it for S4 and it seemed to be very quick. Great! Aug 15, 2012 at 9:49
• @Jimeree (1) Glad I can help; I got lucky as I wasn't sure this was what you needed. (2) Don't be too quick to Accept an answer, as it may discourage someone from posting a better one. (3) {1,2,3,2} should not be included and I do not see it in my output, nevertheless I shall think about that. (4) I have no idea what the computational complexity of this method is (yet) but off hand I would expect it to be pretty good with short cycles and struggle on really long ones. Aug 15, 2012 at 9:59
• Sorry, I meant {1,2,3,1}. That and {3,1,2,3} are the same cycle up to cyclic reordering. Aug 15, 2012 at 10:04
• @Jimeree how large are the graphs you are working with? It should be practical to remove that kind of duplicate at the end of the calculation, but I'm struggling to use that duplication to make the calculation more efficient (though it should be possible). Aug 15, 2012 at 10:29

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;
While[
n = RandomInteger@nodeRange;
adj = RandomChoice[{1-edgeProb, edgeProb} -> {0, 1}, {n, n}];
(* only weakly conncected graphs should be checked *)
(* create adjacency matrix and edge list separately for the two functions *)
];

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

{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},
];
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];

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]
];


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]]


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}}}} *)


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.

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]}], _?(# === {} &)]]
`