I have a very large network of atoms ($\approx 10^6$ atoms) with fixed positions that resemble a cellular structure:

enter image description here

I have two files including:

  • fixed locations of each atom
  • set of points which makes each cell.

The first one is like:

 1 1.72907 3.50783
 2 3.89771 506.561
 3 514.767 4.35252

The second dataset has $\approx 5 \times 10^5$ rows (but size of each row is different). For example:

{1, 485, 486, 487, 966, 968}
{2, 99222, 99223, 99224, 99225, 99706}

I need to find out all sets which have exactly two elements in common ($=$ two cells share an edge or in other words, are neighbors). I already have a code but it's inefficient because it compares the sets line by line. Here's my code ($n$ is number of rows):

 ParallelEvaluate[file = OpenWrite["RN" <> ToString[$KernelID] <> ".dat"]]    

 WriteString[file, i, " ",
 Flatten[Last[Reap[Do[If[Length[Intersection[ring[[i]], ring[[j]]]] == 2, 
 Sow[j]], {j, 1, n}]]]], "\n"]
 , {i, 1, n}];


I find intersection length of a specific set $i$ (ring[[i]]) with all other sets and if it is equal to two, I write the set number in a file. Is there anyway to improve efficiency of this code?


I have an alternative solution without using Intersection and with only one loop, as follows:

 ring = ReadList["rings.dat", Number, RecordLists -> True];
 ParallelEvaluate[file = OpenWrite["RN" <> ToString[$KernelID] <> ".dat"]]    

 RN = Complement[First/@Tally[Flatten[First /@ Position[ring, #] & /@ ring[[i]]]],   
 WriteString[file, i, " ", RN ,"\n"]
 , {i, 1, n}];


But it seems it is not that much better than previous one.

  • 1
    $\begingroup$ {i,1,n} loop is thing missing in you code. $\endgroup$ Commented May 31, 2014 at 3:52
  • 1
    $\begingroup$ What result do you expect from {{1,2,3},{2,3,4},{3,4,5}}? $\endgroup$ Commented May 31, 2014 at 4:31
  • 1
    $\begingroup$ You could let the comparison not go twice, i.e., let the iterators go from {j, 1, n}, {i, j + 1, n} (check out Block[{n = 4}, Do[Print[{i, j}], {j, 1, n}, {i, j + 1, n}]] ) $\endgroup$ Commented May 31, 2014 at 5:40
  • 1
    $\begingroup$ I think you'll need to further specify details to get any useful answer: What is the data (any possibility of doing work as the data is generated?). Is there structure to the data (it appears to be sorted in rows, includes consecutive sequences, and perhaps first element is only increasing or staying the same through rows?). How many unique elements? Can there be duplicate elements in a given row? Do you need only indices of rows, or also what the match pairs are? Short of info leading to shortcuts, you're at 1.3X10^11 conparisons, 30+ hours even if you could do each in a microsecond. $\endgroup$
    – ciao
    Commented May 31, 2014 at 6:56
  • 2
    $\begingroup$ @Mehdi: Ah, "...on a plane with fixed positions...", that's potential for a huge boost: can you scan the plane with an nXn plane that contains just the local atoms that could share a loop? In other words, do you have the needed correspondence between atoms and positions to do that? $\endgroup$
    – ciao
    Commented May 31, 2014 at 22:09

2 Answers 2


apologies for typos I had to retype this. (edit there was one now fixed)

 amax = Max@Flatten@idx;

construct complementary connectivity list.

 atomc = Flatten@# &/@ Last@Reap[Do[Sow[i,#]&/@idx[[i]],{i,Length[idx]}],Range[amax]];

extract neighbors (should be fast):

 celln = Flatten@(First/@Select[Tally@Flatten[atomc[[#]]&/@#],#[[2]]==2 &]) &/@idx ;

timing results for the example set w/ 823 cells:


for the two steps.

example result: celln[[400]]


This takes a minute for ~10^6 cells.

  • 1
    $\begingroup$ +1 Nice use of tags in Sow $\endgroup$ Commented Jun 4, 2014 at 20:34
idx = Import["c:/bazaar/rings.dat"];Max[Flatten[idx]]

gives 1799 distinct points;

First /@ Position[idx, 7]

shows that point 7 occurs in 'cycles' 5,7 and 8

v = Cases[Tally[Flatten@idx[[First /@ Position[idx, 7]]]], {q_, 2} -> q] 

tells us that, apart from point 7 (occurring 3 times), points 5, 27 and 28 occur twice. Now, we need to list the results as pairs of pairs: cycles {i,j} share points (p,q) or, {{i,j},{p,q}}, ...

We can do that by:

  • $\begingroup$ Thank you very much! That helped a lot but final line doesn't return correct answer. $\endgroup$
    – Mahdi
    Commented Jun 2, 2014 at 2:45
  • $\begingroup$ the " {q_,2}-> q];Transpose[ " got massacred in copy-paste. Oopsie! $\endgroup$
    – Wouter
    Commented Jun 2, 2014 at 7:17
  • $\begingroup$ Oh I used the correct form, but the result is not what I want. But the question has an update now based on what you suggested. It gives intended results and I am happy with that. On my current machine (A core 2 Duo with 4GB ram) for first 100 rings, the original one takes 15 seconds but updated solution takes 9 seconds to do the task. $\endgroup$
    – Mahdi
    Commented Jun 2, 2014 at 7:30
  • $\begingroup$ must have missed it somewhere, is the test data available? How big is it? Obviously not 10^6 cells..? $\endgroup$
    – george2079
    Commented Jun 3, 2014 at 21:50
  • $\begingroup$ @george2079: I have shared test data here: drive.google.com/file/d/0BwgIi7LSP4luR2xlU2Y5UnhoQ0E/… $\endgroup$
    – Mahdi
    Commented Jun 4, 2014 at 8:49

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.