Делая подобные вещи часто, я нашел пару лучших практик:
1) Старайтесь использовать как можно больше numpy и numba
2) Попробуйте использовать рычагимаксимально возможное распараллеливание
3) Пропустите циклы для векторизованного кода (здесь мы используем цикл с numba для усиления распараллеливания).
В данном конкретном случае я хочу указать на замедление, введенноепо геопии.Хотя это отличный пакет и дает довольно точные расстояния (по сравнению с методом Haversine), он намного медленнее (не рассматривал реализацию почему).
import numpy as np
from geopy import distance
origin = (np.random.uniform(-90,90), np.random.uniform(-180,180))
dest = (np.random.uniform(-90,90), np.random.uniform(-180,180))
%timeit distance.distance(origin, dest)
216 мкс ± 363 нс на цикл(среднее ± стандартное отклонение из 7 запусков, по 1000 циклов в каждом)
Это означает, что на этом временном интервале для вычисления расстояний 10 миллионов x 1 миллион потребуется приблизительно 2160000000 секунд или 600 тысяч часов.Даже параллелизм только очень поможет.
Поскольку вам интересно, когда точки очень близки, я бы предложил использовать Расстояние Хаверсайна (которое менее точно на больших расстояниях).
from numba import jit, prange, vectorize
@vectorize
def haversine(s_lat,s_lng,e_lat,e_lng):
# approximate radius of earth in km
R = 6373.0
s_lat = s_lat*np.pi/180.0
s_lng = np.deg2rad(s_lng)
e_lat = np.deg2rad(e_lat)
e_lng = np.deg2rad(e_lng)
d = np.sin((e_lat - s_lat)/2)**2 + np.cos(s_lat)*np.cos(e_lat) * np.sin((e_lng - s_lng)/2)**2
return 2 * R * np.arcsin(np.sqrt(d))
%timeit haversine(origin[0], origin[0], dest[1], dest[1])
1,85 мкс ± 53,9 нс на цикл (среднее ± стандартное отклонение из 7 циклов, по 100000 циклов в каждом)
Это уже 100-кратное улучшение.Но мы можем сделать лучше.Возможно, вы заметили @vectorize декоратор, который я добавил из numba.Это позволяет ранее скалярной функции Хаверсайна векторизоваться и принимать векторы в качестве входных данных.Мы воспользуемся этим на следующем шаге:
@jit(nopython=True, parallel=True)
def get_nearby_count(coords, coords2, max_dist):
'''
Input: `coords`: List of coordinates, lat-lngs in an n x 2 array
`coords2`: Second list of coordinates, lat-lngs in an k x 2 array
`max_dist`: Max distance to be considered nearby
Output: Array of length n with a count of coords nearby coords2
'''
# initialize
n = coords.shape[0]
k = coords2.shape[0]
output = np.zeros(n)
# prange is a parallel loop when operations are independent
for i in prange(n):
# comparing a point in coords to the arrays in coords2
x, y = coords[i]
# returns an array of length k
dist = haversine(x, y, coords2[:,0], coords2[:,1])
# sum the boolean of distances less than the max allowable
output[i] = np.sum(dist < max_dist)
return output
Надеемся, теперь у вас будет массив, равный длине первого набора координат (10 миллионов в вашем случае).Затем вы можете просто назначить это для вашего фрейма данных в качестве количества!
Время теста 100 000 x 10 000:
n = 100_000
k = 10_000
coords1 = np.zeros((n, 2))
coords2 = np.zeros((k, 2))
coords1[:,0] = np.random.uniform(-90, 90, n)
coords1[:,1] = np.random.uniform(-180, 180, n)
coords2[:,0] = np.random.uniform(-90, 90, k)
coords2[:,1] = np.random.uniform(-180, 180, k)
%timeit get_nearby_count(coords1, coords2, 1.0)
2,45 с ± 73,2 мс на цикл (среднее ± стандартное отклонение от7 прогонов, 1 цикл каждый)
К сожалению, это все еще означает, что вы будете смотреть на что-то около 20 000 с лишним секунд.И это было на машине с 80 ядрами (с использованием 76ish, на основе top использования).
Это лучшее, что я могу сделать на данный момент, удачи (также, первое сообщение, спасибо за то, что вдохновили меня внести свой вклад!)
PS: Вы можете также изучить массивы Dask и функцию map_block (), чтобы распараллелить эту функцию (вместо того, чтобы полагаться на prange).То, как вы разбиваете данные, может повлиять на общее время выполнения.
PPS: 1 000 000 x 100 000 (в 100 раз меньше, чем у вашего полного набора): 3 мин 27 с (207 секунд), поэтому масштабирование выглядит линейным и немного прощающим.
PPPS: реализовано спростой фильтр разности широт:
@jit(nopython=True, parallel=True)
def get_nearby_count_vlat(coords, coords2, max_dist):
'''
Input: `coords`: List of coordinates, lat-lngs in an n x 2 array
`coords2`: List of port coordinates, lat-lngs in an k x 2 array
`max_dist`: Max distance to be considered nearby
Output: Array of length n with a count of coords nearby coords2
'''
# initialize
n = coords.shape[0]
k = coords2.shape[0]
coords2_abs = np.abs(coords2)
output = np.zeros(n)
# prange is a parallel loop when operations are independent
for i in prange(n):
# comparing a point in coords to the arrays in coords2
point = coords[i]
# subsetting coords2 to reduce haversine calc time. Value .02 is from playing with Gmaps and will need to change for max_dist > 1.0
coords2_filtered = coords2[np.abs(point[0] - coords2[:,0]) < .02]
# in case of no matches
if coords2_filtered.shape[0] == 0: continue
# returns an array of length k
dist = haversine(point[0], point[1], coords2_filtered[:,0], coords2_filtered[:,1])
# sum the boolean of distances less than the max allowable
output[i] = np.sum(dist < max_dist)
return output