Fast way of finding position of Interval containing an Integer

Question

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

Answer 1

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.

Answer 2

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.

Answer 3

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 *)

Answer 4

Here is a function based on my answer to Given an ordered set S, find the points in S corresponding to another list:

IntervalFunction[ints_?MatrixQ]:= With[{S=Developer`ToPackedArray[ints]},
    IntervalFunction[
        S,
        Nearest[S[[All, 1]] -> "Index"],
        Nearest[S[[All, 2]] -> "Index"]]
]

IntervalFunction[s_, nl_, nr_][list_] := With[
    {L = nl[list][[All, 1]], R = nr[list][[All, 1]]},
    With[{ldiff = UnitStep[list-s[[L,1]]], rdiff = UnitStep[s[[R,2]]-list]},
        (L - 1 + ldiff)(L - R - 1 + ldiff + rdiff)
    ]
]

MakeBoxes[i:IntervalFunction[s_, nf_, _],StandardForm] ^:= Module[
    {
    len = Length[s]/2,
    g = FirstCase[ToBoxes[nf], _GraphicsBox, GraphicsBox[Point[{0,0}]], Infinity]
    },
    BoxForm`ArrangeSummaryBox[
        IntervalFunction,
        i,
        RawBoxes@g,
        {
        BoxForm`MakeSummaryItem[{"Intervals: ", Length[s]/2}, StandardForm],
        BoxForm`MakeSummaryItem[{"Range: ", MinMax[s]}, StandardForm]},
        {},
        StandardForm,
        "Interpretable"->True
    ]
]

Example usage (comparing to @LeonidShifrin's answer):

largeInts = Partition[Sort@RandomSample[Range[10000000],2000000],2];

f = IntervalFunction[largeInts];
r1 = f[Range[1000]]; //AbsoluteTiming

pfLarge = makePositionFunction[largeInts];
r2 = pfLarge /@ Range[1000]; //AbsoluteTiming

(r1 /. 0->-1) === r2

{0.001145, Null}

{0.006707, Null}

True

Note that one can use @jkuczm'2 answer instead and achieve a further speed up.

Answer 5

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.

Answer 6

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]