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.

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 *)
share|improve this question
    
You've seen this, no? –  J. M. Aug 24 '12 at 17:35
    
@J.M. Nope, thanks. If the answers there work for a directed graph I'll delete this q –  belisarius Aug 24 '12 at 17:37
    
@J.M. I don't think those solutions scale well for this graph (VertexCount > 3000) –  belisarius Aug 24 '12 at 17:46
    
Apparently so, since I just tried them out. An efficient (in space, time, or preferably both) cycle finder would be a very nice thing. –  J. M. Aug 24 '12 at 17:48
1  
@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! –  belisarius Aug 24 '12 at 18:46
show 1 more comment

2 Answers

up vote 3 down vote accepted

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

share|improve this answer
    
:) 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 –  belisarius Aug 24 '12 at 20:45
1  
@verde The combination of the built-in Graph and Combinatorica is just poison. –  Sjoerd C. de Vries Aug 24 '12 at 20:47
    
I only hope v9 will bring all Combinatorica functionality into System. –  belisarius Aug 25 '12 at 5:46
    
@verde Same here. BTW I just loph the graph. It's strange, monstrous, caterpillar-like; almost like it conveys a message. –  Sjoerd C. de Vries Aug 25 '12 at 6:37
add comment

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.

share|improve this answer
    
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"} –  belisarius Aug 24 '12 at 19:53
    
@Verde Thanks. Probably programmer error. Though I suspect the underlying method is sound. –  Daniel Lichtblau Aug 24 '12 at 20:00
    
I corrected the code. It was not taking into account vertices with no outgoing edges. Should behave better now. –  Daniel Lichtblau Aug 24 '12 at 23:22
    
cycles[pairs, 3] gives me {{1025, 2035, 1823, 1025}, {1027, 2036, 1825, 1027}, {1027, 2037, 1827, 1027}}. Am I missing something? –  belisarius Aug 25 '12 at 0:13
    
@Verde Sorry, I left out an important part of the preprocessing. It now needs: pairs = ee /. reprule /. DirectedEdge -> List; Should be better with that. –  Daniel Lichtblau Aug 25 '12 at 15:59
show 1 more 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.