Match pick between intervals and data points

Question

I have two sets of data as follows:

1-Data1 is a series of intervals. These intervals rarely overlap but they could.

2-Data2 consists of pairs of data where the first element of each pair could be or could not be a member of the intervals in Data1.

What I am looking for is a fast way to extract the data so that I get list which each element contains the interval from Data1 and all points from Data2 that their first elements lie within this interval.

As an example of what I am looking for is as follows:

n = 200;
intervals = 
  Sort[Partition[RandomReal /@ Table[{i, i + 5}, {i, 0, n, 6}], 2]];
v = RandomReal[n, Length[intervals] + RandomInteger[40]];
v2 = {#, RandomReal[{0.5, 0.8}]} & /@ v;
Plot[UnitStep[(x - #1) (#2 - x)] & @@@ intervals, {x, 0, 201}, 
 Epilog -> {PointSize[0.005], Red, Point[v2]}, Filling -> Bottom, 
 GridLines -> None, Exclusions -> None, PlotPoints -> 100]

As you can see from this plot that some of the data points are within the intervals, but some are not.

I want to get pairs of data like following:

intervals2 = 
  Pick[intervals, 
   Or @@@ Transpose[
     IntervalMemberQ[Interval /@ intervals, #] & /@ v2[[;; , 1]]]];
result = {#, Pick[v2, IntervalMemberQ[Interval@#, v2[[;; , 1]]]]} & /@
    intervals2;

And now if i plot that i got:

Plot[UnitStep[(x - #1) (#2 - x)] & @@@ result[[;; , 1]], {x, 0, 201}, 
 Epilog -> {{PointSize[0.005], Black, 
    Point[Flatten[result[[;; , 2]], 1]]}}, Filling -> Bottom, 
 GridLines -> None, PlotPoints -> 100]

This is really not an efficient way. My real Data1 has hundreds of thousands intervals and Data2 also contains hundreds of thousands pairs.

If i use the method mentioned above it will go forever.

I would appreciate it if someone can suggest a better and faster way.

Thank you

Answer 1

Here's a brute force alternative that's 5 times faster than @halirutan 's (subject to verification by others).

AbsoluteTiming[
 (* Set things up as before *)
 n = 200000;
 intervals = Sort[Partition[RandomReal /@ Table[{i, i + 5}, {i, 0, n, 6}], 2]];
 v = RandomReal[n, Length[intervals] + RandomInteger[40]];
 v2 = {#, RandomReal[{0.5, 0.8}]} & /@ v;

 (* Create starting and ending points for boundaries of intervals *)
 (* {2.07388,"S"},{6.19068,"E"},{12.1195,"S"},{19.04680,"E"},... *)
 i0 = Flatten[Table[{"S", "E"}, Length[intervals]]];
 boundaries = Sort[Transpose[{Flatten[intervals], i0}], #1[[1]] < #2[[1]] &];

 (* Merge with data *)
 all = Sort[Join[boundaries, v2], #1[[1]] < #2[[1]] &];

 (* Find the positions of the starting and ending points *)
 s = Flatten[Position[all[[All, 2]], "S"]];
 e = Flatten[Position[all[[All, 2]], "E"]];

 (* Points to keep are those between the starting and ending points *)
 i1 = Select[Transpose[{s + 1, e - 1}], #[[2]] >= (#[[1]]) &];
 i2 = Flatten[Range[#[[1]], #[[2]]] & /@ i1];
 result = all[[i2, All]];
 ]

This takes about a half a second with n = 200000.

Now this won't work if the intervals overlap.

Answer 2

You could break things down. Assume you have your data but only one interval and you just collect all elements that fall into this single interval. The naive approach is to go through all your entries and compare the interval-boundaries with the position of the data-point.

In compiled code this is almost straight-forward with the exception that I'm using a Bag because you cannot know upfront how many points will fall into this interval. A Bag is nice for continuously collecting data-points.

As soon as you make this compiled function parallel, you can apply it to a list of intervals without changing a single line of code. The getIntervals function is a wrapper for the compiled function that will finally build the same structure as you have. I will keep intervals with no points, but the can be deleted in almost no-time.

getIntervalsC = Compile[{{intv, _Real, 1}, {data, _Real, 2}},
   Module[{min = intv[[1]], max = intv[[2]], 
     b = Internal`Bag[Most[{0.}]], res = {0.}},
    Do[
     If[min <= val[[1]] <= max,
      Internal`StuffBag[b, val, 2];
      ], {val, data}
     ];
    res = Internal`BagPart[b, All];
    If[Length[res] > 1,
     Partition[res, 2],
     {res}
     ]
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ];
getIntervals[int_, data_] := 
 Transpose[{int, getIntervalsC[int, data]}]

Your example with n=200

result = DeleteCases[getIntervals[intervals, v2], {{_, _}, {{}}}];

Plot[UnitStep[(x - #1) (#2 - x)] & @@@ result[[;; , 1]], {x, 0, 201}, 
 Epilog -> {{PointSize[0.005], Black, 
    Point[Flatten[result[[;; , 2]], 1]]}}, Filling -> Bottom, 
 GridLines -> None, PlotPoints -> 100]

Runtime depends on how many processors you have. If I run your example with 16 threads and n=200000 it takes about 2.5 seconds. Your function takes about 200x this time.