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.

To build a graph, I need to apply a function f[a_, b_] to all pairs of a list (3500 elements). The function itself returns a link {a <-> b} if a particular relation holds - I collect all the results into a list and use it as input to Graph[].

The question is: is there an efficient and elegant way to do this? I've tried two ways: a recursive method and a (similar) iterative method. The former went over the usual recursion limits, the latter was slow and (I believe) not the optimal way to do it.

Both of these performed by applying f[] to the First[] element vs. the 2nd, 3rd, ...Last[] elements, and collecting the results. Then I'd remove the First[] from the list, and repeat - doing n*(n-1)/2 evaluations of f[]. This many evaluations is required, but I definitely do not have a clean, functional implementation.

So, in short, if this can be turned into an efficient 1-liner, instead of a loop, please do let me know! Thanks in advance.

share|improve this question
    
Thanks all for the great, instructive answers. –  michar Aug 18 '12 at 2:07

3 Answers 3

up vote 9 down vote accepted

There are two built-in functions to generate pairs, either with (Tuples) or without (Subsets) duplication. Since your question states the number of iterations as $n*(n-1)/2$ I believe you want the latter:

set = {1, 2, 3, 4};
Subsets[set, {2}]
{{1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}}

The short notation for Apply at level 1 is @@@, so this gives f for each pair:

f @@@ Subsets[set, {2}]
{f[1, 2], f[1, 3], f[1, 4], f[2, 3], f[2, 4], f[3, 4]}

This is in my opinion the most elegant code to produce this specific result, and it is quite fast, but it is not memory efficient if you only need to collect a result for a low percentage of the pairs. Let's define f as follows (I use {} here generically):

f[a_, b_] /; b~Divisible~a := {a, b}

f[___] = Sequence[];

If we now compute all the pairs for Range[5000] it takes one gigabyte of memory:

pairs = Range@5000 ~Subsets~ {2};
pairs // ByteCount

1099780032

And applying f we see that of the nearly 12.5 million pairs we need a return for only 38376 of them which takes only 3MB of storage:

r1 = f @@@ pairs;

ByteCount[r1]

Length[r1]

3377120

38376

Yet, the maximum memory used is 1.6GB:

 MaxMemoryUsed[]

1629793536

A simple method to reduce memory consumption is to process the subsets in blocks, rather than all at once, as follows:

set = Range@5000;
n = Length@set;
max = n (n - 1)/2;
block = 10000;

Timing[
 r2 =
   Join @@ Array[
     f @@@ Subsets[set, {2}, block {# - 1, #} + {1, 0}] &,
     ⌈max/block⌉
   ];
]

Length[r2]

MaxMemoryUsed[]

{8.299, Null}

38376

19769800

This only uses a maximum of ~20MB of memory, only a few MB over the baseline on my system.
(It issues a Subsets::take message but there is no error.)

My preferred method

Another method, and the one that I prefer, is to compute the pairs more manually allowing f to be embedded in the process so as to not generate all pairs beforehand. This uses Outer to effect the pair generation for each element separately (that is, all pairs starting with a certain element).

pairMap[f_, s_] := Module[{ss = s}, 
  Flatten[Outer[f, {#}, ss = Rest@ss, 1] & /@ Most@s, 2] ]

pairMap[f, Range@5000] // Length // Timing

MaxMemoryUsed[]

{7.816, 38376}

19430512

This also uses only a small amount of memory, and testing should bear out that it is faster as well. A variation of this method that may be even faster is to not build an output expression at all, relying instead on Sow and Reap to gather your results:

g[a_, b_] /; b ~Divisible~ a := Sow[ {a, b} ]

pairScan[f_, s_] := Module[{ss = s}, Outer[f, {#}, ss = Rest@ss, 1] & ~Scan~ Most@s ]

Reap[ pairScan[g, Range@5000] ][[2, 1]] // Length // Timing

MaxMemoryUsed[]

{6.583, 38376}

18757552

An argument for Do loops

While Outer is somewhat faster, after further consideration and conferring with jVincent I think perhaps after all a Do loop is as good as anything. One could write pairScan in this way:

pairScan2[f_, s_] := Module[{ss = s}, Do[f[i, j], {i, ss}, {j, ss = Rest@ss}] ]

Reap[ pairScan2[g, Range@5000] ][[2, 1]] // Length // Timing

MaxMemoryUsed[]

{7.613, 38376}

18711080

share|improve this answer
    
Very helpful, I can solve my problem and learned a lot at the same time. –  michar Aug 18 '12 at 2:16
vertices = Range[10];
pairs = Tuples[vertices, 2];
func[x__] := First@x <= Last@x;
edges = Pick[pairs, func[#] & /@ pairs];
(* or *)
edges = Pick[pairs, Boole[func[#]] & /@ pairs, 1];
Graph[DirectedEdge @@@ edges, VertexLabels -> "Name", ImagePadding -> 20]

enter image description here

 Graph[DirectedEdge @@@ edges, VertexLabels -> "Name",  ImagePadding -> 20]

enter image description here

vertices = Range[100];
pairs = Tuples[vertices, 2];
func2[x__] := PrimeQ[(First@x)^(Last@x) + 3];
edges = Pick[pairs, func2[#] & /@ pairs];
Graph[DirectedEdge @@@ edges, ImagePadding -> 20,
  VertexLabels -> "Name",
  VertexLabelStyle -> Red,
  VertexStyle -> White,
  VertexShapeFunction -> "Point",
  EdgeStyle -> Directive[Opacity[.5], Hue[.15, .5, .8]],
  Background -> Black,
  EdgeShapeFunction -> ({Arrowheads[.01], Arrow[#1, 0.01]} &), 
  ImageSize -> 800]

enter image description here

Update: Using SparseArray with AdjacencyGraph and Graph:

func2[x__] := PrimeQ[(First@x)^(Last@x) + 3]; (* or your choice of test function *)
adjmtrx = SparseArray[{{i_, j_} /; func2[{i, j}] -> 1}, {20, 20}];
styling = {ImagePadding -> 20,
     VertexLabels -> "Name",
     VertexLabelStyle -> Red,
     VertexStyle -> White,
     VertexShapeFunction -> "Point",
     EdgeStyle -> Directive[Opacity[.5], Hue[.15, .5, .8]],
     Background -> Black,
     EdgeShapeFunction -> ({Arrowheads[.02], Arrow[#1, 0.01]} &), 
     ImageSize -> 400};

Row[{AdjacencyGraph[adjmtrx, styling], Spacer[5], 
      Graph[Range@First@Dimensions[adjmtrx], 
           DirectedEdge @@@ adjmtrx["NonzeroPositions"], styling]}]

enter image description here

share|improve this answer

Yes, there are definitely shorter, non-loop ways to do this. Define your link-defining function like this (with whatever test you need as the first argument to If:

islinked[a_Integer?Positive, b_Integer?Positive] := 
 If[Mod[a, 3] == 0 && Mod[b, 2] == 1, a -> b, {}]

You can then Apply this at the level of each row using the @@@ shorthand.
I would use Tuples rather than Subsets to allow for both 1->2 and 2->1. Eliminate the empty cases with a replacement rule as shown.

gg = islinked @@@ Tuples[Range[10], 2] /. {} -> Sequence[]

(* {3 -> 1, 3 -> 3, 3 -> 5, 3 -> 7, 3 -> 9, 6 -> 1, 6 -> 3, 6 -> 5, 
    6 -> 7, 6 -> 9, 9 -> 1, 9 -> 3, 9 -> 5, 9 -> 7, 9 -> 9} *)

Graph[gg]

enter image description here

share|improve this answer
    
Thanks for this and the next answer. One question though, and it cuts to the center of my post: does generating the list of subsets (or tuples) beforehand cause a difficulty? If the original list has n elements, there are n(n-1) pairs. Maybe mathematica knows to generate them iteratively, as it is handling the Apply? –  michar Aug 17 '12 at 4:02
    
@michar - I find Tuples to be pretty quick: but note, it has n*n pairs - things like 2 -> 2 are possible, thus the self-loops in the picture above. –  Verbeia Aug 17 '12 at 4:14

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.