1
$\begingroup$

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]

enter image description here

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]

enter image description here

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Brilliant. Thanks. $\endgroup$ – Algohi Mar 5 '18 at 21:29
  • $\begingroup$ You might want to re-think that. When there are no "ties" (i.e., no common values between the interval borders and the data), this works fine. But if some of the data values match the interval borders then some of those "ties" might get left off. That's fixable but I'll wait for you to tell me if there are ties or not. $\endgroup$ – JimB Mar 5 '18 at 22:42
  • $\begingroup$ I am aware of that. I have made slight modification to your idea which makes it more applicable to my problem. I was looking for fast method and this is super fast. Now my problem can be solved in a fraction of a second. $\endgroup$ – Algohi Mar 5 '18 at 22:46
2
$\begingroup$

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]

Mathematica graphics

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.

$\endgroup$

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.