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.

I'm trying to optimize code that uses Position MANY times.

The following toy code works:

list1 = RandomInteger[{1, 10000}, 500000];
list2 = RandomInteger[{1, 10000}, {3, 500}];

test = ParallelTable[
        Flatten[
          Map[
            Position[list1, #] &, list2[[i]]
          ]
         ], {i, 1,  Length[list2]}];

My trouble begins because I have so much data to process. When I use my actual data in the code above, Mathematica takes 34 minutes to complete the computation. This is just too long to be practical, so I began exploring Compile.

When I try this:

  test2 = 
    Compile[{{list1, _Integer, 1}, {list2, _Integer, 2}}, 
      Table[
        Flatten[
           Map[
              Position[list1, #] &, list2[[i]]
           ]
        ], {i, 1, Length[list2]}]];

   wow = test2[list1, list2]

 CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. 

You can see that Mathematica complains.

My questions are:

  1. How can I debug what the problem is in the Compiled statement above? I really don't know how to start figuring out the problem.

  2. Is there an obvious alternative to using Position that itself would be faster than Position?

I am not new to Mathematica, but I admit to being completely new to compiling code.

Any thoughts?

share|improve this question
add comment

3 Answers

up vote 7 down vote accepted

I will try to address the second question you asked, since the first one has been answered already. The following code will not win a beauty contest, but will speed up your code 800 times for your test list sizes (note that I will use Table in place of ParellelTable in your code, since ParallelTable gives some really wrong results for your tests on my machine).

This is a massive binary search function, which takes a sorted list and a list of elements, and returns a list of position pairs for minimal and maximal position of an element in a list:

ClearAll[bsearchMassiveC];
bsearchMassiveC =
  Compile[{{list, _Real, 1}, {elems, _Real, 1}},
   Module[{n0 = 1, n1 = Length[list], m = 0, 
      res = Table[{-1,-1}, {Length[elems]}], ctr = 0, found = False
     },
     Do[
       found = False;
       n0 = 1;
       n1 = Length[list];
       m = 0;
       While[n0 <= n1, m = Floor[(n0 + n1)/2];
         If[list[[m]] == elem,
            n0 = n1 = m;
            While[n0 > 0 && list[[n0]] == elem , n0--];
            While[n1 <= Length[list] && list[[n1]] == elem , n1++];
            res[[++ctr]] = {n0 + 1, n1 - 1};
            found = True;
            Break[]
         ];
         If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]
       ];
       If[! found, res[[++ctr]] = {-1,-1};], 
       {elem, elems}
     ];
     Take[res, ctr]
   ],
   CompilationTarget -> "C"];

For example,

bsearchMassiveC[{1,1,1,2,2,3,4,4,6,6,6,7,7,7,7},{1,3,5,7}]

(*  {{1,3},{6,6},{-1,-1},{12,15}}  *)

This is the main function. We first sort the list, then find positions in a sorted list, and finally reconstruct the positions in the original list:

ClearAll[getPositions];
getPositions[list_, elems_] :=
  With[{ord = Ordering[list]},
    With[{sorted = list[[ord]]},
      Table[
         ord[[Flatten[Range @@@ Select[
              bsearchMassiveC [sorted, elems[[i]]],
              # != {-1, -1} &]]
         ]],
         {i, Length[elems]}]]];

Here are the benchmarks:

list1=RandomInteger[{1,10000},500000];
list2=RandomInteger[{1,10000},{3,500}];

(res1= getPositions[list1,list2]);//AbsoluteTiming

(*   {0.1035157,Null}   *)

(res2 = 
    Table[Flatten[Map[Position[list1,#]&,list2[[i]]]],i,1,Length[list2]}]
 );//AbsoluteTiming

 (*  {82.6718750,Null}  *)

res1==res2

(*  True  *)

I did not attempt to parallelize my code, but I think it should be possible. Again, my attempts to use your code verbatim with ParallelTable resulted in wrong results. It would be good if others could confirm it, and even better to understand the reason behind it.

share|improve this answer
    
+1, binary search is the obvious solution but I just got time now and apparently it's too late :) –  acl Aug 13 '12 at 22:52
    
@acl Thanks:). I didn't have much time for SE recently, but today I had some, and the temptation was too high. I voted for your answer too, right after I posted mine. –  Leonid Shifrin Aug 13 '12 at 23:00
add comment

This is not as fast as Leonid's code, but runs in about 2.5 seconds without using Compile. It uses an indexing scheme based on the fact that list1 is composed of integers within a relatively small range. A new list L1 is created such that the element L1[[n]] contains a list of the positions of the value n in list1.

test2 = Module[{L1},
   L1 = ConstantArray[{}, {10000}];
   MapIndexed[(L1[[#]] = {L1[[#]], #2}) &, list1];
   L1 = Flatten /@ L1;
   Flatten /@ Map[L1[[#]] &, list2, {2}]];
share|improve this answer
    
A nice one, +1. I also often use such "hashing" schemes, but for some reasons this did not occur to me here. –  Leonid Shifrin Aug 13 '12 at 21:03
    
while generally I'd go for binary search, this is clever –  acl Aug 13 '12 at 22:53
add comment

The problem is that the returned object would have been ragged, and Compile can't deal with that.

To see it's ragged, you try running it for smaller inputs and look at the output (even though Compile complains, you do get the output obtained by calling the uncompiled code):

list1 = RandomInteger[{1, 10000}, 500]; 
list2 = RandomInteger[{1, 10000}, {3, 50}];
f[list1, list2]

Mathematica graphics

And to see that Compile can't deal with ragged arrays, try to construct a list of (lists with different lengths):

Compile[{}, 
  Table[RandomInteger[{0, 1}, RandomInteger[{0, 10}]], {i, 1, 10}]][]

Mathematica graphics

while if the sublists are of the same length, no complaints:

Compile[{}, Table[RandomInteger[{0, 1}, 2], {i, 1, 10}]][]

Mathematica graphics

I'll not get into speeding this up.

share|improve this answer
add 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.