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$ May 31 '14 at 3:52
  • 1
    $\begingroup$ What result do you expect from {{1,2,3},{2,3,4},{3,4,5}}? $\endgroup$ May 31 '14 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$ May 31 '14 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
    May 31 '14 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
    May 31 '14 at 22:09

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$ Jun 4 '14 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
    Jun 2 '14 at 2:45
  • $\begingroup$ the " {q_,2}-> q];Transpose[ " got massacred in copy-paste. Oopsie! $\endgroup$
    – Wouter
    Jun 2 '14 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
    Jun 2 '14 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
    Jun 3 '14 at 21:50
  • $\begingroup$ @george2079: I have shared test data here: drive.google.com/file/d/0BwgIi7LSP4luR2xlU2Y5UnhoQ0E/… $\endgroup$
    – Mahdi
    Jun 4 '14 at 8:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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