You don't provide much information, which is unfortunate, since this kind of optimization problem relies on knowing how to "skew" the code in order to get better performance.
For example, are your points constant? If you are always searching the same set of 27k points - for example, cities or road termini or railroad stations or whatever - then you can pre-compute lots of data about those points to make the code faster.
Are the points accumulated gradually over time? Could you perform an expensive computation on each point safely, knowing that the various points are collected at some interval that allows for that?
Are the points close together or far apart, relative to the likely neighborhood distances? That is, can you estimate the number of points that will be in the likely neighborhood queries?
Is your program memory-constrained? Can we trade memory allocation for performance?
eps values - the neighborhood distances - truly variable? Or is there a static list (1km, 10km, 25km, 50km) that the user chooses from?
All these questions, and more, could be used to tune the code for speed. When you're asking a speed question, there's no such thing as too much data!
Finally, what flavor of python are you using? PyPy? Cpython? Jython? Cython?
With that out of the way, let me make a few "low hanging fruit" suggestions. I'll assume you're using "the python that came with my system" - meaning Cpython 2.7 or 3.4+.
The cpython compiler isn't nearly as smart as you might like. Let's look at the generated code from your function (I put it in
test.py in a class named
Python 2.7.13 (default, Apr 23 2017, 16:50:35)
[GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import test
>>> import dis
3 0 BUILD_LIST 0
3 STORE_FAST 2 (nearby)
4 6 SETUP_LOOP 230 (to 239)
9 LOAD_FAST 0 (self)
12 LOAD_ATTR 0 (data)
>> 16 FOR_ITER 219 (to 238)
19 STORE_FAST 3 (pt)
6 22 LOAD_GLOBAL 1 (map)
25 LOAD_GLOBAL 2 (math)
28 LOAD_ATTR 3 (radians)
31 LOAD_FAST 1 (p)
34 LOAD_ATTR 4 (x)
37 LOAD_FAST 3 (pt)
40 LOAD_ATTR 4 (x)
43 LOAD_FAST 1 (p)
46 LOAD_ATTR 5 (y)
49 LOAD_FAST 3 (pt)
52 LOAD_ATTR 5 (y)
55 BUILD_LIST 4
58 CALL_FUNCTION 2
61 UNPACK_SEQUENCE 4
64 STORE_FAST 4 (px)
67 STORE_FAST 5 (ptx)
70 STORE_FAST 6 (py)
73 STORE_FAST 7 (pty)
7 76 LOAD_FAST 4 (px)
79 LOAD_FAST 5 (ptx)
83 STORE_FAST 8 (dx)
8 86 LOAD_FAST 6 (py)
89 LOAD_FAST 7 (pty)
93 STORE_FAST 9 (dy)
9 96 LOAD_GLOBAL 2 (math)
99 LOAD_ATTR 6 (sin)
102 LOAD_FAST 9 (dy)
105 LOAD_CONST 1 (2)
109 CALL_FUNCTION 1
112 LOAD_CONST 1 (2)
116 LOAD_GLOBAL 2 (math)
119 LOAD_ATTR 7 (cos)
122 LOAD_FAST 6 (py)
125 CALL_FUNCTION 1
128 LOAD_GLOBAL 2 (math)
131 LOAD_ATTR 7 (cos)
134 LOAD_FAST 7 (pty)
137 CALL_FUNCTION 1
141 LOAD_GLOBAL 2 (math)
144 LOAD_ATTR 6 (sin)
147 LOAD_FAST 8 (dx)
150 LOAD_CONST 1 (2)
154 CALL_FUNCTION 1
157 LOAD_CONST 1 (2)
163 STORE_FAST 10 (a)
10 166 LOAD_CONST 1 (2)
169 LOAD_GLOBAL 2 (math)
172 LOAD_ATTR 8 (asin)
175 LOAD_GLOBAL 2 (math)
178 LOAD_ATTR 9 (sqrt)
181 LOAD_FAST 10 (a)
184 CALL_FUNCTION 1
187 CALL_FUNCTION 1
191 STORE_FAST 11 (c)
11 194 LOAD_CONST 2 (6371.0)
197 LOAD_FAST 11 (c)
201 STORE_FAST 12 (d)
12 204 LOAD_FAST 12 (d)
207 LOAD_FAST 0 (self)
210 LOAD_ATTR 10 (eps)
213 COMPARE_OP 0 (<)
216 POP_JUMP_IF_FALSE 16
13 219 LOAD_FAST 2 (nearby)
222 LOAD_ATTR 11 (append)
225 LOAD_FAST 3 (pt)
228 CALL_FUNCTION 1
232 JUMP_ABSOLUTE 16
235 JUMP_ABSOLUTE 16
>> 238 POP_BLOCK
14 >> 239 LOAD_FAST 2 (nearby)
First, have a look at "line 10" - the block of code that starts with a "10" in the left-most column. That code corresponds to this python:
c = 2*math.asin(math.sqrt(a))
The Python VM is a stack machine, so every opcode either takes something off the stack, puts something on the stack, or both. For example, the various
LOAD_ ops put things onto the stack. The
CALL_FUNCTION op calls a function (which is stored on the stack) with a number of arguments. And
BINARY_MULTIPLY takes two arguments off the stack, multiplies them, and pushes the result back onto the stack.
There's a couple of optimization opportunities here that are good "rules of thumb." You can take this knowledge and apply it to any Cpython code you need to speed up (note: Cpython).
First, there's that
2 * at the front of the expression. That's fairly harmless, except that on the very next line you multiply by 6371.0, and you never do anything else with the
c variable - it's just a temp. So to squeeze out some performance, go ahead and move the
2* to the next line. It turns out that the compiler will go ahead and do constant math in the same expression, so you can say
6371.0 * 2 * c and it will collapse that to a single constant.
Second, note that there's this sequence of lookups all over the place:
169 LOAD_GLOBAL 2 (math)
172 LOAD_ATTR 8 (asin)
In general, Python believes that any operation might cause the result of a name lookup to change. Because you can return any type, and any type can overload pretty much any operation, so why shouldn't it be possible for
2 * math.asin(x) to change the function pointed to by
Because of this, and because name lookups like this really can be expensive (remember
property methods and descriptors), you generally want to cache dotted-name lookups in a local variable if you're going to be looping with them:
asin = math.asin
sqrt = math.sqrt
Please note: This is not the same as doing
from math import asin because that will still create a "global" name, and require a
LOAD_GLOBAL operation, instead of a
LOAD_FAST. Notice the
FAST in the opcode name. ;-)
The same thing is true for almost every dotted name - if you are going to use it more than once, put it in a local variable.
One thing that might not be obvious is that you can cache "bound methods" - that is, method lookups like
nearby.append. So you can also do:
nearby_append = nearby.append
And squeeze out another lookup.
But wait, there's more!
Because before you start doing that "peephole" optimization stuff, you need to look at improving your code in some other ways. Take a look at this:
def find(self, p):
nearby = 
for pt in self.data:
# Haversine algorithm
px, ptx, py, pty = map(math.radians, [p.x, pt.x, p.y, pt.y])
For the purposes of this function, the
p object is a constant. So why do you keep doing the same conversion on it? Remember that
map just applies the given function to each argument in turn. So you're calling
p.y 27000 times. Guess what? The answer isn't going to change!
Of course, I suspect that none of your points are going to change. That is, I suspect that once you create a point object, it doesn't get updated. So I think you would benefit from performing this conversion in the class initializer and storing the radian values in the object:
def __init__(self, lon, lat):
self.x = lon
self.x_radians = math.radians(lon) # <- HERE
self.y = lat
self.y_radians = math.radians(lat) # <- HERE
self.cluster = None
self.neighbors = 0
If you do that, you can eliminate that entire
map call, and just retrieve the values. Yes, you pay an up-front cost of 27000 extra function calls, but it beats doing it in your find loop. (If you can't pay that cost, you could certainly do it as a separate loop at the front of your find function. You're going to do it anyway, and this lets you store the result.)
Also of note:
math.cos(pty). You can certainly compute the cosine of
py once at the top of the function. And if you're going to call this function a lot, it's probably worth storing the cosine(y) for all the points in the initializer.
Another suggestion is simply to guess what's a neighbor and what's not. You've got a nice calculation there, with lots of trigonometry. But keep in mind what it's doing: it's trying to compute distance on a sphere, which is basically distance along a circle. The key word is circle.
At θ = 45°, cosθ is √2/2 as is sinθ. The sum of the two is √2, which is the maximum sum for a circle of radius 1. Thus, you know that for whatever distance N defines your neighborhood, the sum of cosine and sine will be less than or equal to √2 * N.
If your points are widely separated - that is, if you're dealing with cities and not the entrances to buildings in the same town - you can implement a heuristic: change √2 into 2.5 to allow for spherical error, and compute the corners of a "square" centered at
p with sides 2 * 2.5 * N. Then just compare pt.x and pt.y with the min/max x and y of this square. Only if the point is inside the square do you perform the "expensive" computations. This adds a set of 4 comparisons to your loop, but eliminates the trig calls for all the non-matching candidates. Hopefully, that's most of the points.
This is important! If your points are heavily clustered, there's no sense having a heuristic to guess if one point is a neighbor of another - because you know the answer is yes! In that case, it would make more sense to develop a reverse heuristic, that says "here's a square contained in the circle. If the point is inside this square you know it's a neighbor, otherwise do the computation." The objective with heuristics it to produce a fast computation that handles 75% or more of cases. If your data is biased, make the heuristic align with the biases!
Finally, another simple optimization would be to take your
self.eps value and divide it by the various factors you are presently using to multiply. This can be done outside the loop, once, and saves you some cycles.
First, don't forget that if A is a neighbor of B, then B is a neighbor of A. Make sure you are only looping over the "higher" values in your N-squared neighbor search:
for i, p1 in enumerate(points):
for p2 in points[i+1:]:
You just have to be sure and add A to B's neighborhood and B to A's neighborhood at the same time.
Depending on what you're doing, there are a couple of geospatial data structures that you might profitably use. Two obvious ones are a hash of some sort and a K-d tree. The hash approach will work if your points are clustered and if you can determine a key generation algorithm that will reflect this. For example, if you know that your points are GPS coordinates in major metropolitan areas (New York, Los Angeles, Dallas, Chicago) then a hash key that "rounded" the coordinates to the nearest city center might be adequate. (Or truncated the lat/lon to, say, the nearest 5 degrees.) In this case, you'd want to have a hash value with lots of collisions. The collisions would be the likely neighbors. So you would iterate over the hash keys (equal to, say, 10% of the points) and identify which buckets you would search in. On the other hand, if the points are more distributed, or the clustering isn't obvious, then go with the K-d tree.
The K-d tree is a binary tree (or binary-search-tree if you're a child of the 21st century) that can be used for locating any point. To find neighboring points using your distance function, you would simply do the distance computation at each node of the tree. The special case would be when a node's key value is within the heuristic neighborhood of your target point, you'll have to traverse both children of the node, to find any other nodes in the other half of the tree that might be in the neighborhood. For example, if a node is a lat node, with lat=A, and your target point has lat=A-2km, then for neighborhood=5km you'll have to search on either side of that node.
Also, note that you'll need to develop a strategy for nodes with equal values. For example, if two points have the same longitude, would you store a list-of-points in the node, or make a rule that says "sort by lowest latitude and put the other points in the right hand subtree"? It's up to you, but you need to handle it. (I suggest the latter.)