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 have a large list (about 2.2 million) of pairs of integers. Each pair defines an Interval[{}]. Given some integer x, I need to find the position of the interval in my list that bounds x (or alternatively, return the element of the list).

None of the intervals overlap, thus x can belong to only one of the intervals. You can assume the large list is sorted. The challenge is that I need to run this "look-up" millions of times, so even a small gain in speed would be very beneficial

I tried a few obvious brute force approaches. First I created a list of intervals from the list of integer pairs, and then tried

Pick[myBigList, 
     Map[IntervalMemberQ[#, x] &, myBigList]]; //AbsoluteTiming
4.196432

and then

Select[myBigList, IntervalMemberQ[#,x]& ]; // AbsoluteTiming
8.814063

Assuming the list is ordered, then

LengthWhile[mySortedBigList, #[[1]] <= x &]
3.556826

All of these methods seem too slow. Any suggestions would be much appreciated.

(For those of you interested, here's the background on this problem. I need a fast way of translating an IP address to an approximate set of latitude-longitude coordinates. Each pair of integers in the list above corresponds to a range of IP addresses that have been converted from the zz.zz.zz.zz format to an integer. The interval defined by each pair of integers is associated with a set of lat-long coordinates. Again, this needs to be run millions or tens of millions of time on a regular basis. When this look-up was done on a 64-bit Windows desktop system running standard SQL, it took about 8 hours to geocode ten million IPs)

Thanks,

Mark

share|improve this question
    
Since you can sufficiently narrow down the lat-long with just the first three octets (sometimes, even with two), can you not create a hash table (downvalues/dispatch) for ~16.7 million IPs (first three octets) to the lat-long coordinates? –  rm -rf Mar 12 at 18:44
    
Will a bounding interval exist for every input? –  belisarius Mar 12 at 18:53
    
Are the intervals "touching", as in $[0,1)$, $[1,4)$, $[4,6)$, etc., or are there gaps inbetween? It sounds like good application for a binary search. –  Szabolcs Mar 12 at 19:37

5 Answers 5

Here is a solution based on binary search (compiled).

Implementation

First, this is a version of a binary search, which would return the position of a maximal number in a list, smaller or equal to yours, and -1 if no such is found:

ClearAll[bsminComp];
bsminComp=
    Compile[            
        {{lst,_Integer,1},{elem,_Integer}},
        Module[{n0=1,n1=Length[lst],m=0,pos=-1},                
            While[                    
                n0<=n1,                    
                m=Floor[(n0+n1)/2];
                If[lst[[m]]==elem,pos=m;Break[];];
                If[lst[[m]]<elem,n0=m+1,n1=m-1]
            ];
            If[pos==-1,pos=If[lst[[m]]<elem,m,m-1]];
            pos
        ]
    ]

Here is then the main function:

ClearAll[makePositionFunction];
makePositionFunction[ints_List]:=
    Module[{starts=ints[[All,1]],ends=ints[[All,2]]},
        Function[                
            num,
            Module[{},
                With[{pos=bsminComp[starts,num]},                        
                    If[pos==-1||pos==0,Return[-1,Module]];
                    If[ends[[pos]]>=num,pos,-1]
                ]
            ]
        ]
    ];

How it works

Basically, the strategy is simple: find the lower end of the interval possibly containing your number, using binary search on the lower ends of your intervals. Then, find the upper end of the interval and check if it is larger or equal to the element in question. If yes, return the position found by a binary search. If no, return -1, meaning that the element is not within any of the intervals.

Tests and benchmarks

ints = Partition[Sort@RandomSample[Range[20], 10], 2]

(* {{1, 3}, {4, 5}, {9, 13}, {14, 15}, {17, 18}}

pfun = makePositionFunction[ints];

pfun /@ Range[15]

(*  {1, 1, 1, 2, 2, -1, -1, -1, 3, 3, 3, 3, 3, 4, 4} *)

Larger test:

largeInts = Partition[Sort@RandomSample[Range[10000000],2000000],2];
pfLarge = makePositionFunction[largeInts];
pfLarge/@Range[1000]//Short//AbsoluteTiming

(* {0.018322,{-1,1,1,1,1,1,2,2,2,3,<<980>>,104,104,-1,105,105,105,105,105,105,-1}} *)

Module vs With - a cautionary tale

You could've noticed that I used Module to localize some of the variables, as:

Module[{starts = ints[[All, 1]], ends = ints[[All, 2]]}, ...]

Question is, why I did not use With (since nothing changes here), and what would happen if I did. You can try yourself, replacing Module -> With, and you will find that the function will still work, but be thousands times slower. You will also find that it would have a huge ByteCount, if you apply ByteCount to the resulting pure function (for large lists).

Now, here is a general caution: when you embed a large enough list (or other expression) into your pure function, using With, you are asking for trouble. The way Mathematica symbolic engine works is that when you apply such a function, it has to apply this huge expression (it is huge because Function makes its internals fully transparent to the evaluator, and if the body happens to contain huge expressions, Function will be a huge expression). By using Module, I delay the actual computation of the list until run-time, and make the body of my Function tiny - which makes all the difference here. This trick saved me more than once, and is IMO essential to know for anyone who wants to routinely use closures in Mathematica.

share|improve this answer
    
Heh, I was too lazy to implement a binary search function, but I should've remembered that I probably could've found it in quite a few of your answers... would've been simpler to steal it ;) –  rm -rf Mar 12 at 19:28
    
@rm-rf Yeah, right - I was luckier since I also happen to use these for work, so all I had to do was to go to the right source file, which was quick :). –  Leonid Shifrin Mar 12 at 19:30
2  
@rm-rf Actually, I personally view the last remark on Module vs With as the most valuable part of this answer, since I doubt that many people were aware of this subtlety. And there are real-world cases when this is critical. At the same time, this isn't very easy to discover on your own - I've spent a lot of time more than once figuring this out. This can be particularly bad when you have a large and complex code, because it is very hard to track the slowdowns to pure functions (they don't have a name after all, and don't exist as a part of global state). –  Leonid Shifrin Mar 12 at 19:31
    
I agree. I didn't know that either. I don't think I've run into such a use case, but I know that if I came across a need, I would've definitely used With there. Good to know! –  rm -rf Mar 12 at 19:34
    
@rm-rf The second part of this story, which I did not tell above, is that for such pure functions with captured Module - variables (mutable), these variables will be garbage-collected all right, but only if $HistoryLength was set to 0 before those were produced (in the interactive session), otherwise those variables won't be garbage-collected. IIRC, they are garbage-collected just fine when they are not part of any Out result, in any case. But this is important to know, of course, and we are lucky with this - otherwise, the practical uses for such closures would've been very limited. –  Leonid Shifrin Mar 12 at 19:38

An alternative approach, which uses binary search under the hood, is to have a zero-order interpolating function that maps from the jth value in the list to j. I'll illustrate using the same example as from my other response.

SeedRandom[1111];                                                      
ll = DeleteDuplicates[RandomInteger[10^9, 3*10^6]];             
          If[Mod[Length[ll], 2] == 1, ll = Most[ll]];                        
myBigList = Partition[Sort[ll], 2];

vals = Map[Round[Mean[#]] &, myBigList];                               
vals2 = RandomSample[vals];

That sets up the data for the example. Next is the off-line preprocessing.

Timing[interp =                                                        
            Interpolation[Transpose[{Flatten[myBigList], Range[Length[ll]]}],   
             InterpolationOrder -> 0];]                                         

(* Out[12]= {3.756000, Null} *)

Here we time how long to locate all query points.

Timing[res =                                                           
            Table[nbr = interp[vals2[[j]]];                                     
             If[OddQ[nbr], myBigList[[(nbr + 1)/2]], myBigList[[nbr/2]]], {j,   
              Length[vals2]}];]                                                 

(* Out[13]= {14.288000, Null} *)

Now check that the sorted result is the same as myBigList.

Timing[resb = Sort[res];]                                              

(* Out[14]= {3.644000, Null} *)

resb===myBigList                                                       

(* Out[16]= True *)

Somewhat faster than the other method I showed (factor of 2-3 or so). Not sure how it compares to Leonid's compiled explicit binary search.

share|improve this answer
    
+1. I did not fully optimize my function, though: first, binary search can be compiled to C, and then, if we have many elements in a list (rather than, say, separate function calls for each element), it is possible to write a compiled vectorized binary search (which will also compile the loop over the searched elements) - I've done that before - and also vectorize the top-level closure. If someone else with (too much) free time on their hands would come up with a detailed timing comparisons between all posted solutions, I will make and post the vectorized version :). –  Leonid Shifrin Mar 12 at 20:56
    
@Leonid Shifrin What about setting the compiled search to be listable via RuntimeAttributes? That should handle the vectorizing part more or less for free. –  Daniel Lichtblau Mar 12 at 21:00
    
Unfortunately, upon a superficial look at least, I think this won't work, in such a direct form at least, because we need it Listable in only the second argument, but not the first (which also is a List), while setting RuntimeAttributes -> Listable would affect both arguments. I could compile the binary search with embedded list, first, and then inline / use that in another compiled function which I can make runtime - Listable - this might work. I'll test this when I have more time. –  Leonid Shifrin Mar 12 at 21:09
    
@Leonid Shifrin Right, you'd feed in the pairs list and return a listable compiled function that does binary search of that list on input query points. –  Daniel Lichtblau Mar 12 at 21:24

You can use Nearest to find possible pair endpoints, then do post-processing to get the actual intervals.

Here is an example using around 1.5 million intervals.

SeedRandom[1111];                                                      
Timing[ll = DeleteDuplicates[RandomInteger[10^9, 3*10^6]];             
          If[Mod[Length[ll], 2] == 1, ll = Most[ll]];]                          

(* Out[25]= {0.356000, Null} *)

Length[ll]                                                             

(* Out[26]= 2995472 *)

Timing[myBigList = Partition[Sort[ll], 2];]                            

(* Out[27]= {0.612000, Null} *)

Timing[nf = Nearest[ll];]                                              

(* Out[28]= {10.096000, Null} *)

We'll use the interval midpoints for our sample queries, and shuffle them so that we know order is not playing a role in the speed.

vals = Map[Round[Mean[#]] &, myBigList];                               
vals2 = RandomSample[vals]; 

Now for each query we find the three closest neighbors. We could find the two closest but in rare cases where intervals are separated by one unit we might bet the wron result.

Timing[res = Map[nf[#, 3] &, vals2];]                                  

(* Out[45]= {26.920000, Null} *)

Timing[res2 = Sort[Map[Sort, res]];]                                   

(* Out[46]= {1.804000, Null} *)

Now figure out whether the correct interval is from the first two or last two points. This next line is really part of the preprocessing.

Timing[Map[(interval[#] = True)&, myBigList];]                         

(* Out[90]= {4.664000, Null} *)

Timing[res3 = Map[If[TrueQ[interval[Most[#]]], Most[#], Rest[#]] &, res
2];]                                                                            

(* Out[91]= {6.476000, Null} *)

Check that we get back that original pair list; the example is rigged to do exactly that.

res3===myBigList                                                       

(* Out[93]= True *)
share|improve this answer

Since the integers are all ordered, one approach would be to put them all in a big list and then use a NearestFunction to find the closest one. As a simple example, here's a sorted list of integers and a nearest function `nf' based on that list.

sorted = Sort[RandomInteger[{0, 100000}, 500]];
nf = Nearest[sorted];

To find the closest integer to 20000 in the list, use nf[20000], which returns the nearest value. You would then need to do some post-processing to figure out whether the actual match was below, in, or above the relevant interval.

share|improve this answer

Because the value to be looked up can only belong to one interval, a trivial optimization is to use the Select variant which picks only the first element for which IntervalMemberQ is True:

Select[myBigList, IntervalMemberQ[#,x]&, 1]
share|improve this answer
    
This won't change the linear complexity though, and for lists of this size (millions of elements), this is a big deal. –  Leonid Shifrin Mar 12 at 20:18

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.