# Counting number of points in fluctuating lines that fall within a circle

I have a list of points pts that fluctuate around some mean line. pts is pre-sorted by the projection of the points onto this mean line (for instance, in the special case where the points fluctuate about a horizontal line, then the points would be sorted by their x-values). I want to find the number of elements in pts that land within a circle centered at fixed point centre and for various radii.

My problem is that my current code is too slow, as I work with systems with a huge number of points, and I feel that there must be a way to further exploit the geometry of the problem to get significant speed-ups.

Here's some sample data for pts, which I'll have fluctuate about the mean line defined by meanY and unit vector vec:

(*The numbers I use in this block are not important*)

vec = {0.951057, 0.309017};
orthoVec = {-vec[[2]], vec[[1]]}; (*the normal vector*)
meanY = -10;

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec}; (*This line is the change of basis so that the data is along our mean line*)


and the circle parameters:

centre=RandomReal[{-1,1},2];


This gives a configuration that looks like this:

Here's my current attempt: I find the norm of all the points from the circle center, sort them, and them finding the index of the largest norm less than the circle radius. To find the indices, I use the function cf from Henrik Schumacher's answer here. The code for my attempt is:

norms = Sqrt[(pts[[All,1]] - centre[[1]])^2 + (pts[[All,2]] - centre[[2]])^2]; (*take the norm of the points*)
norms=Sort[norms]; (*sort the points*)
ctr = cf[norms, radii]; (*count the number of points inside the circle for each radius in radii*)


But given the structure of my point set (i.e., along a line), this feels like a waste, primarily because sorting the norms is very costly when I already have a point set that is sorted along a line. Any ideas on how I can exploit this fact to speed up my code?

Here is my own start of an idea on how this algorithm might be improved: because the points in pts are already sorted by their projections onto the mean line, we get the following plot for the norms by index:

This leads me to think that there should be a way to first rough sort the points according to just some function of their index, and then apply Sort to the resulting list, which I think should result in fewer computations, but I'm not sure how to go about coding this up.

Edit

The points fluctuating about a line is just one piece of my data. What I'm actually trying to do is to compute the number of points inside the circle from a 2D point system as a function of circle radius, where my point system can be broken into these fluctuating lines:

ptsArray = Table[
px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px];
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts.{vec, orthoVec}, {meanY, -20, 20, 2.5}];


A more complete picture of my point system would look like this:

I COULD have just used all of the data combined, computed the norms, and found the number of points in the circle that way, but it slows down the code significantly when we go to very large systems (on the order of millions of points in all) and I've found it much more efficient to break it into these 1D sequences. That is why the code needs to be so fast: I'm not just computing this for 1 line of points, I'm looping through thousands of such lines and then adding up the number of points in the circle from each.

• If many of the points aren't near the circle, you could throw away all the ones not in the circle's bounding square or octagon instead of calculating all of their norms.
Sep 7 at 2:54
• Why we need lots of radii ? Sep 7 at 7:27
• @Adam thanks, I'm working on adding this. But I think that the main problem is the complexity of the algorithm when the problem is highly geometric Sep 7 at 12:29
• @cvgmt the reason I need to get the number for many radii is that my code is meant to find how the number of points changes as a function of radius Sep 7 at 12:31
• @Adam see my edit. Yes, I am partitioning the data into lines to make the code more efficient. So I'm actually trying to compute the hyperuniformity coefficient of a 2D point system. Say I fix my radius at r. Once I've computed the number of points inside this circle, I have to move the circle to a different centre elsewhere and find the number of points falling in that one in turn, and again and again, and finally take the variance of the number of points found in my collection of circles. Sep 7 at 18:55

Perhaps the cleanest way to find this distribution is with Nearest. I.e.

pts=RandomReal[{-10,10},{10000000,2}];
center={1,1};
(* .08 seconds, 78624 *)
(* upwards of 40 seconds, 78624 *)


with the added bonus that if you omit the preceeding Length you have a list of points ordered by closeness to center. Not sure yet how to calculate the variance of the distribution Nearest[pts,center,10000], but I'm sure it can be done. Regardless, we can plot it (this time with a uniform grid, idk how that should affect the variance):

ListPlot[
Norm[#-{1,1}]&/@
Nearest[
Join@@CoordinateBoundsArray[{{-10,10},{-10,10}},0.0063210],
{1,1},{\[Infinity],1}]
]


• A different train of thought led me to stackoverflow.com/questions/6691491/… -- thanks Daniel Lichtbau!!! -- and I saw Nearest mentioned which prompted this answer.
Sep 8 at 5:12
• Interesting...I'll have to think about what this would imply for later steps of my code (like taking the variance, as you noted). Sep 9 at 19:14
• You'd have to invert the function (i.e. reverse each datapoint), but the problem then (besides needing to delete the nonincreasing entries) is that the radii aren't uniformly sampled, so you can't simply call Variance on the first n entries, you have to normalize somehow.
Sep 9 at 21:40

I might be very well not understanding the issue because with even 100,000 points obtaining the number of points within all possible radii doesn't seem to take much time:

n = 100000;
px = FoldList[Plus, RandomReal[{0.5, 1}, n]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, n];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec};

centre = {0, 0};
AbsoluteTiming[
distance = Sqrt[(Total[(# - centre)^2] & /@ pts)] // Sort;
index = Range[Length[distance]];]
(* {0.0292905, Null} *)

ListPlot[Transpose[{distance, index}], Frame -> True,


One can then turn that into a function:

points = Interpolation[Transpose[{distance, index}]]
Floor[points[5000]]
(* 13337 *)

• See my edit on speed. I agree that for one line it's not such an issue, but it becomes prohibitive because this one line of points is actually just one of thousands. Sep 7 at 18:57

EDIT: I've updated the answer to use BinCounts, which suits this problem in multiple ways.

Most of the old description below applies to this solution, but here the counts for intervals between (squared) radii is performed with BinCounts, and they're summed up to correspond with the UnitStep solution with Accumulate:

Needs["Developer"];
With[{length = 1000000},
With[
{pts = ToPackedArray[RandomReal[{-30, 30}, {length, 2}]],
centre = RandomReal[{-1, 1}, 2]},
Accumulate[
BinCounts[(pts + ConstantArray[-centre, length])^2 . {1, 1},
{Prepend[Range[1., 40., 0.5], 0.]^2}]]]]


The actual computation takes about 0.10 seconds for 79 radii and million points on my laptop, or 8.2 seconds for 100 million points(!). (Using packed arrays gives just a couple-percent benefit on calculation side but it's so easy to add that I did it.)

EDIT: This can be sped even further!

If we use regular binning, the binning of the large squared distance list is even faster. Thus we first perform bin counts for bins between zero and square of the largest radius, with a bin size of square of the step, and then re-bin this data using a second HistogramList run (which is a cousin of BinCounts, but also handles WeightedData).

fixedHistogramList corresponds to HistogramList, but fixes a slow implementation for specific inputs - maybe in the future it is unnecessary.

Needs["Developer"];

(* Fix an odd performance issue which is present at least up to v13.1.0. *)
ClearAll[fixedHistogramList];
fixedHistogramList[list_List, {xmin_?NumericQ, xmax_, dx_}] :=
{Range[xmin, xmax, dx], BinCounts[list, {xmin, xmax, dx}]};

With[{binrange = {1., 40., 0.5}},
With[{length = 1000000},
With[{
pts = ToPackedArray[RandomReal[{-30, 30}, {length, 2}]],
centre = RandomReal[{-1, 1}, 2]},
fixedHistogramList[(pts + ConstantArray[-centre, length])^2 . {1, 1},
{0, binrange[[2]], binrange[[3]]}^2] //
Accumulate[Round[Last[
HistogramList[WeightedData[Most@First@#, Last@#],
{Prepend[Range @@ binrange, 0.]}^2]]]] &]]]


The actual computation (that is, excluding generation of the pts list) takes 3.6 seconds on my laptop for 100 million points and 79 radii (end result bins).

This is very slightly faster than putting Sqrt around the vectored squared distance calculation and dropping ^2s from range definitions. The square-rootless version has the benefit of avoiding computing million square roots (although they can be quite fast in a vectored manner), while one with a square root has a benefit of having less "unnecessary" bins on regular binning step.

I simplify the problem and don't care about the lines - just those random points and circles of different radii sharing the same centre.

This computes counts for a single set of points with different radii:

With[{length = 1000000},
With[
{pts = RandomReal[{-30, 30}, {length, 2}],
centre = RandomReal[{-1, 1}, 2]},
With[{sqdistances = ((pts + ConstantArray[-centre, length])^2) . {1, 1}},
Table[Total[UnitStep[r^2 - sqdistances]], {r, 1., 40., 0.5}]]]]


Here centre is subtracted from all point coordinates, individual resulting $$x$$ and $$y$$ values squared with a vectored way, and sums computed using .. This results squared distances from centre. After this UnitStep is used to convert each distance inside r to $$1$$, outliers $$0$$ and a total of these is computed.

Computing these counts for 79 radii with million points takes about 0.6 seconds on my laptop.

When number of radii being tested is large, this computation is pretty much dominated by the counting part, in which the UnitStep approach seems to be unfortunately dramatically more efficient that Select or Count approaches for this sort of a problem. So, the following, easier-to-understand approach is only marginally slower in the end:

With[{length = 1000000},
With[
{pts = RandomReal[{-30, 30}, {length, 2}],
centre = RandomReal[{-1, 1}, 2]},
With[{distances = Norm[# - centre] & /@ pts},
Table[Total[UnitStep[r - distances]], {r, 1., 40., 0.5}]]]]

• This is really interesting. Do you know how BinCounts actually works algorithmically? Sep 9 at 19:17
• @az123p No, I don't. I actually experimented with arrays of different distributions, and ran into a conclusion that distribution can affect runtime quite a bit. My first assumption is that it uses some sort of a tree data structure, since for uniform distribution the growth of runtime seems roughly logarithmic regarding number of bins. It's reasonably obvious why it performs better in this case for large-ish number of bins in comparison to the UnitStep and Total approach - that approach has very clearly linear runtime growth, and performs a lot of repeated, unnecessary work. Sep 10 at 5:15
• For completely regular bin spacing, BinCounts can actually run at $O(1)$ time (theoretically speaking). There's the $\{x_{min}, x_{max}, dx\}$ binning form for this purpose! Sep 10 at 5:16
• @az123p I updated my answer with a form which first uses regular binning, and then re-bins the results. This is significantly faster still than the previous variation, but frankly quite a hack. Sep 10 at 7:47

To keep things simple I give an example for just one radius. You may adapt it for several radii yourself.

First we create the points and the circle according to your example:

SeedRandom[1];

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec};

circle = Circle[{0, 0}, 15];

Graphics[{Point[pts], circle}, Axes -> True]


Then we fit a line through the points:

line = LinearModelFit[pts, {x}, x]["BestFit"]


Then subtract from all points the y intercept of the line: pts1. As the points now go through the origin we can rotate them so that they lay along the x axis: pts1:

pts[[All, 2]] -= line[[1]];
pts1 = RotationMatrix[-line[[2, 1]]] . Transpose@pts // Transpose;


We also shift and rotate the circle: circle1:

circle1 = circle;
circle1[[1]] =
RotationMatrix[-line[[2, 1]]] . (circle1[[1]] - {0, line[[1]]});


This now looks like:

Graphics[{Point[pts1], circle1}, Axes -> True]


Next we determine the points where the transformed circle cuts the x axis. These points determine approx. the points to count. As the points a fuzzy, we add some margin for safety:

limit = Sqrt[circle1[[2]]^2 - circle1[[1, 2]]^2] + 1;


Possible points to count lay inside the limits:

possible = Select[pts1, (-limit <= #[[1]] <= limit) &];


Now because the points do not lay on a straight line we need to check and count more carefully:

Select[possible, Norm[circle1[[1]] - #] <= circle1[[2]] &] // Length

(* 27 *)


And for convenience, here is all together:

SeedRandom[1];

vec = {0.951057, 0.309017};
orthoVec = {-vec[[2]], vec[[1]]}; (*the normal vector*)
meanY = -10;

px = FoldList[Plus, RandomReal[{0.5, 1}, 100]];
px = px - Median[px]; (*just for some symmetry in the data*)
py = RandomReal[meanY + {-0.5, 0.5}, 100];
pts = Transpose@{px, py};
pts = pts . {vec, orthoVec};

circle = Circle[{0, 0}, 15];

line = LinearModelFit[pts, {x}, x]["BestFit"];
pts[[All, 2]] -= line[[1]];
pts1 = RotationMatrix[-line[[2, 1]]] . Transpose@pts // Transpose;
circle1 = circle;
circle1[[1]] =
RotationMatrix[-line[[2, 1]]] . (circle1[[1]] - {0, line[[1]]});
Graphics[{Point[pts1], circle1}, Axes -> True]
limit = Sqrt[circle1[[2]]^2 - circle1[[1, 2]]^2] + 1;
possible = Select[pts1, (-limit <= #[[1]] <= limit) &];
Select[possible, Norm[circle1[[1]] - #] <= circle1[[2]] &] // Length

• Thanks for the answer. I added an edit that should clarify my question. It also means that things like Select are simply too slow for my purposes. But I think the general idea is interesting. Sep 7 at 18:59
• I think for many points on different lines a pre-selection using an inscribed and a circumscribed square would be fitting. Points inside the inscribed square are accepted. Points between the 2 squares are checked by the distance from the center. Sep 8 at 14:00