5
$\begingroup$

Warning: run the following code in a fresh Mma session, as some symbols could be shadowed (depending on your Mma version)

While trying to answer this question, I fell into the following:

 (* Let's load a large Directed Graph and convert it to Combinatorica *)
g = Graph@Union@Flatten[
     Thread[DirectedEdge @@ ##] & /@ Select[{#, IsotopeData[#, "DaughterNuclides"]} & /@ 
        IsotopeData[], #[[2]] != {} &]];
Needs["GraphUtilities`"]
<< Combinatorica`
cg = ToCombinatoricaGraph[g];

Girth[cg] gives the length of a shortest cycle in a simple graph g.

So, let's check if cg is Simple and calculate its Girth:

{SimpleQ@cg, Girth@cg}
(*
-> {True, 3}
*)

So there is at least one Cycle in cg of length 3.

But look what happens when we try to find it by the two available methods in Combinatorica:

{ExtractCycles@cg, FindCycle@cg}
(*
-> {{},{}}
*)

So, two questions:

  1. Is this a bug?
  2. What is the easiest way to find all cycles in g without using Combinatorica?

Edit

BTW, the (now) standard Graph functionality also detects cycles:

AcyclicGraphQ[g]
(* -> False *)
$\endgroup$
6
  • $\begingroup$ You've seen this, no? $\endgroup$ Commented Aug 24, 2012 at 17:35
  • $\begingroup$ @J.M. Nope, thanks. If the answers there work for a directed graph I'll delete this q $\endgroup$ Commented Aug 24, 2012 at 17:37
  • $\begingroup$ @J.M. I don't think those solutions scale well for this graph (VertexCount > 3000) $\endgroup$ Commented Aug 24, 2012 at 17:46
  • $\begingroup$ Apparently so, since I just tried them out. An efficient (in space, time, or preferably both) cycle finder would be a very nice thing. $\endgroup$ Commented Aug 24, 2012 at 17:48
  • 1
    $\begingroup$ @TomD I used g1 = FromOrderedPairs[ EdgeList@g /. Thread[Rule[VertexList[g], Range@VertexCount@g]] /. DirectedEdge -> List, Type -> Directed]; Print@{Girth@g1,FindCycle@g1}; and the result is the same. Thanks for the pointer! $\endgroup$ Commented Aug 24, 2012 at 18:46

2 Answers 2

4
$\begingroup$

I don't think it's a bug. Directed radioactive decay graphs shouldn't have cycles by definition and AcyclicGraphQ doesn't see them either:

AcyclicGraphQ@g

True

You don't seem to have heeded your own warning to start with a fresh kernel, as playing around with the code for a while gave me False too.

The output of ExtractCycles and FindCycle is therefore correct.

Girth doesn't seem to take directionality into account when determining cycles:

Cyclic directed graph:

Girth@ ToCombinatoricaGraph@
     System`Graph[{1 \[DirectedEdge] 2, 2 \[DirectedEdge] 3, 3 \[DirectedEdge] 1}]

3

Acyclic directed graph (note the reversed direction of the last edge):

Girth@ToCombinatoricaGraph@
     System`Graph[{1 \[DirectedEdge] 2, 2 \[DirectedEdge] 3, 1 \[DirectedEdge] 3}]

3

$\endgroup$
4
  • $\begingroup$ :) Of course I expected no cycles, that was what bothered me. I think WRI has a great opportunity to improve in the graph's features robustness $\endgroup$ Commented Aug 24, 2012 at 20:45
  • 1
    $\begingroup$ @verde The combination of the built-in Graph and Combinatorica is just poison. $\endgroup$ Commented Aug 24, 2012 at 20:47
  • $\begingroup$ I only hope v9 will bring all Combinatorica functionality into System. $\endgroup$ Commented Aug 25, 2012 at 5:46
  • $\begingroup$ @verde Same here. BTW I just loph the graph. It's strange, monstrous, caterpillar-like; almost like it conveys a message. $\endgroup$ Commented Aug 25, 2012 at 6:37
7
$\begingroup$

Edited for correctness:

I use the variant for directed graphs from here. I take your graph as above, extract edges, rename so vertices are integers from 1 to #vertices. After finishing we revert to the original names.

ee = EdgeList[g];
vv = VertexList[g];
reprule = Thread[vv -> Range[Length[vv]]];
revrule = Map[Reverse, reprule];
pairs = ee /. reprule /. DirectedEdge -> List;

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

cycles[omat_, k_] := Module[
  {n = Length[Union[Flatten@omat]], m2, cyc, cyclist, mat},
  mat = Join[omat, Thread[{Range[n], 0}]];
  m2 = Map[Last, SplitBy[Sort[mat], First], {2}];
  m2 = m2 /. 0 :> Sequence[];
  cyclist = 
   Flatten[Drop[MapIndexed[{#2[[1]], #1} &, m2, {2}], -k + 1], 1];
  cyclist = Select[cyclist, #[[2]] > #[[1]] &];
  Do[cyclist = 
    Flatten[Map[extendCycle[#, m2] &, cyclist], 1], {k - 2}];
  Map[If[MemberQ[m2[[Last[#]]], First[#]], Append[#, First[#]], 
      Null] &, cyclist] /. Null :> Sequence[]]

I get no cycles in any length between 3 and 20. Have not tried further. I now believe the code will behave correctly. Let me know if not.

$\endgroup$
6
  • $\begingroup$ Those cyc3 don't look good (from a physical viewpoint). For example, in the first one check IsotopeData["Iron62", "DaughterNuclides"] ... its only daughter is {"Cobalt62"} $\endgroup$ Commented Aug 24, 2012 at 19:53
  • $\begingroup$ @Verde Thanks. Probably programmer error. Though I suspect the underlying method is sound. $\endgroup$ Commented Aug 24, 2012 at 20:00
  • $\begingroup$ I corrected the code. It was not taking into account vertices with no outgoing edges. Should behave better now. $\endgroup$ Commented Aug 24, 2012 at 23:22
  • $\begingroup$ cycles[pairs, 3] gives me {{1025, 2035, 1823, 1025}, {1027, 2036, 1825, 1027}, {1027, 2037, 1827, 1027}}. Am I missing something? $\endgroup$ Commented Aug 25, 2012 at 0:13
  • $\begingroup$ @Verde Sorry, I left out an important part of the preprocessing. It now needs: pairs = ee /. reprule /. DirectedEdge -> List; Should be better with that. $\endgroup$ Commented Aug 25, 2012 at 15:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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