Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}}
share|improve this question
    
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? –  Jimeree Aug 15 '12 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. –  Mr.Wizard Aug 15 '12 at 8:12
    
Sure. Just added them at the bottom. –  Jimeree Aug 15 '12 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. –  Mr.Wizard Aug 16 '12 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). –  Jimeree Aug 19 '12 at 14:20
show 2 more comments

4 Answers 4

up vote 8 down vote accepted

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}}
share|improve this answer
1  
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! –  Jimeree Aug 15 '12 at 9:39
    
Oh, I just tried it for S4 and it seemed to be very quick. Great! –  Jimeree Aug 15 '12 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. –  Mr.Wizard Aug 15 '12 at 9:59
    
Sorry, I meant {1,2,3,1}. That and {3,1,2,3} are the same cycle up to cyclic reordering. –  Jimeree Aug 15 '12 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). –  Mr.Wizard Aug 15 '12 at 10:29
show 4 more comments

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]

enter image description here

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]
   ];
share|improve this answer
add comment

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.

share|improve this answer
add comment

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]}], _?(# === {} &)]]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.