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];
radii=Table[r, {r, 1., 40., 0.5}];
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
Thanks for the comments everyone. Here's some more info about what I'm trying to do that will hopefully clarify my question.
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.
radii
? $\endgroup$