Python – как ускорить расчет расстояний между городами

У меня 55249 городов в моей базе данных. Каждый из них имеет значения долготы широты. Для каждого города я хочу рассчитать расстояния до любого другого города и хранить те, которые находятся не дальше 30 км. Вот мой алгоритм:

# distance function from math import sin, cos, sqrt, atan2, radians def distance(obj1, obj2): lat1 = radians(obj1.latitude) lon1 = radians(obj1.longitude) lat2 = radians(obj2.latitude) lon2 = radians(obj2.longitude) dlon = lon2 - lon1 dlat = lat2 - lat1 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 c = 2 * atan2(sqrt(a), sqrt(1-a)) return round(6373.0 * c, 2) def distances(): cities = City.objects.all() # I am using Django ORM for city in cities: closest = list() for tested_city in cities: distance = distance(city, tested_city) if distance <= 30. and distance != 0.: closest.append(tested_city) city.closest_cities.add(*closest) # again, Django thing city.save() # Django 

Это работает, но занимает очень много времени. Понадобится несколько недель. В любом случае я мог бы ускорить его?

5 Solutions collect form web for “Python – как ускорить расчет расстояний между городами”

Вы не можете позволить себе вычислить расстояние между каждой парой городов. Вместо этого вам нужно поместить ваши города в структуру данных разделения пространства, для которой вы можете выполнять быстрые запросы ближайшего соседа. SciPy поставляется с реализацией kd- scipy.spatial.KDTree , scipy.spatial.KDTree , подходящей для этого приложения.

Здесь есть две трудности. Во-первых, scipy.spatial.KDTree использует евклидово расстояние между точками, но вы хотите использовать большое расстояние по поверхности Земли. Во-вторых, долгота обходит вокруг, так что ближайшие соседи могут иметь долготы, которые отличаются на 360 °. Обе проблемы могут быть решены, если вы примете следующий подход:

  1. Преобразуйте свои местоположения из геодезических координат ( широта , долгота ) в координаты ECEF (с землей, с землей) ( x , y , z ).

  2. Поместите эти координаты scipy.spatial.KDTree в scipy.spatial.KDTree .

  3. Преобразуйте расстояние вашего большого круга (например, 30 км) в евклидово расстояние.

  4. Вызовите scipy.spatial.KDTree.query_ball_point чтобы получить города в радиусе действия.

Вот пример кода, иллюстрирующий этот подход. Функция geodetic2ecef принадлежит PySatel Дэвидом Парунакяном и лицензируется под лицензией GPL.

 from math import radians, cos, sin, sqrt # Constants defined by the World Geodetic System 1984 (WGS84) A = 6378.137 B = 6356.7523142 ESQ = 6.69437999014 * 0.001 def geodetic2ecef(lat, lon, alt=0): """Convert geodetic coordinates to ECEF.""" lat, lon = radians(lat), radians(lon) xi = sqrt(1 - ESQ * sin(lat)) x = (A / xi + alt) * cos(lat) * cos(lon) y = (A / xi + alt) * cos(lat) * sin(lon) z = (A / xi * (1 - ESQ) + alt) * sin(lat) return x, y, z def euclidean_distance(distance): """Return the approximate Euclidean distance corresponding to the given great circle distance (in km). """ return 2 * A * sin(distance / (2 * B)) 

Давайте составим пятьдесят тысяч случайных городов и преобразуем их в координаты ECEF:

 >>> from random import uniform >>> cities = [(uniform(-90, 90), uniform(0, 360)) for _ in range(50000)] >>> ecef_cities = [geodetic2ecef(lat, lon) for lat, lon in cities] 

Поместите их в scipy.spatial.KDTree :

 >>> import numpy >>> from scipy.spatial import KDTree >>> tree = KDTree(numpy.array(ecef_cities)) 

Найти все города в пределах 100 км от Лондона:

 >>> london = geodetic2ecef(51, 0) >>> tree.query_ball_point([london], r=euclidean_distance(100)) array([[37810, 15755, 16276]], dtype=object) 

Этот массив содержит для каждой точки, которую вы запросили, массив городов на расстоянии r . Каждому соседу присваивается его индекс в исходном массиве, который вы передали в KDTree . Таким образом, есть три города в пределах примерно 100 км от Лондона, а именно города с индексами 37810, 15755 и 16276 в первоначальном списке:

 >>> from pprint import pprint >>> pprint([cities[i] for i in [37810, 15755, 16276]]) [(51.7186871990946, 359.8043453670437), (50.82734317063884, 1.1422052710187103), (50.95466110717763, 0.8956257749604779)] 

Заметки:

  1. Вы можете видеть из примера вывода, что соседи с долготами, которые отличаются примерно на 360 °, правильно обнаружены.

  2. Такой подход кажется достаточно быстрым. Здесь мы находим соседей в пределах 30 км для первой тысячи городов, занимающих около 5 секунд:

     >>> from timeit import timeit >>> timeit(lambda:tree.query_ball_point(ecef_cities[:1000], r=euclidean_distance(30)), number=1) 5.013611573027447 

    Экстраполируя, мы ожидаем найти соседей в пределах 30 км для всех 50 000 городов за четыре минуты.

  3. Функция euclidean_distance переоценивает евклидово расстояние, соответствующее заданному расстоянию большого круга (чтобы не пропустить ни одного города). Это может быть достаточно для некоторых приложений, ведь города – это не точечные объекты, но если вам нужна более высокая точность, то вы можете фильтровать полученные точки, используя, скажем, одну из величайших функций расстояния круга от геофизики .

Вы можете ускорить вычисление расстояния, не входя в сложные тригонометрические формулы, если знаете, что города расположены на расстоянии более 30 км друг от друга, поскольку их различие в широте соответствует более чем 30-километровой дуге. Дуге длины a = 30 км соответствует угол a / r = 0,00470736, следовательно:

 def distance(obj1, obj2): lat1 = radians(obj1.latitude) lon1 = radians(obj1.longitude) lat2 = radians(obj2.latitude) lon2 = radians(obj2.longitude) dlon = lon2 - lon1 dlat = lat2 - lat1 if dlat > 0.00471: return 32 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 c = 2 * atan2(sqrt(a), sqrt(1-a)) return round(6373.0 * c, 2) 

Радиус 32 является просто фиктивным значением, указывающим, что города находятся дальше, чем на расстоянии 30 км друг от друга. Вы должны применить аналогичную логику для долготы, для которой вы должны учитывать наибольшую абсолютную широту:

  if cos(lat1) * dlon > 0.00471 and cos(lat2) * dlon > 0.00471: return 32 

Если вы знаете, что ваши города находятся в фиксированном диапазоне широт, вы можете настроить постоянный предел на наихудший случай. Например, если все ваши города находятся в смежных Соединенных Штатах, они должны быть ниже 49 ° северной широты, а ваш предел составляет 0,00471 / cos (49 °) = 0,00718.

  if dlon > 0.00718: return 32 

Этот более простой критерий означает, что вы вводите точный расчет для слишком многих городов в Техасе или Флориде. Вы также можете связать эти критерии. сначала используйте приблизительный предел, затем точный предел, основанный на максимальной абсолютной широте, затем вычислите точное расстояние для всех остальных кандидатов.

Вы можете рассчитать этот предел заранее с максимальной абсолютной широтой. Эта эвристика также должна помочь вам поместить города в ведра фиксированной долготы и широты, как предположил РемкоГерлих. Его метод должен значительно ускорить ваш процесс, предварительно рассмотрев только разумные пары городов.

Редактировать Мне немного стыдно видеть, что мой код выше не проверяет абсолютное значение предела. Во всяком случае, настоящий урок здесь заключается в том, что независимо от того, насколько вы ускоряете расчет расстояний, реальная выгода для больших наборов данных исходит от выбора интеллектуального механизма поиска, такого как поиск в bucket или kd деревьях, которые другие комментаторы предложили, возможно, вместе с некоторыми комментариями к отсекают двойные проверки.

Сначала я создавал бы «секторы», каждый из которых был бы ограничен двумя широтами в Х км и двумя долготами на расстоянии в 6 км друг от друга. X должен быть как можно большим, с одним ограничением: все города в пределах сектора находятся не более чем в 30 км.

Секторы могут храниться в массиве:

 Sector[][] sectors; 

Внутри этого массива легко идентифицировать сектор, содержащий конкретные координаты. Также легко идентифицировать смежные сектора для определенного сектора.

Затем:

(1) Каждому городу присваивается свой сектор. В каждом секторе есть список городов, лежащих в нем.

(2) Для каждого города найдите все города в своем секторе. Те немедленно соответствуют критерию 30 км.

(3) Для каждого города C найдите все города C 'во всех 8 смежных секторах. Для каждого C 'проверьте расстояние CC' и выведите CC ', если оно составляет <30 км.

Этот алгоритм по-прежнему равен O (n ^ 2), но он должен быть быстрее, поскольку для каждого города вы проверяете только небольшое подмножество всего набора.

  1. Попробуйте использовать различные алгоритмы, чтобы ускорить вычисление на единицу расстояния, как кто-то заявил.
  2. Используйте только перестановки городов, чтобы не пересчитывать повторения.
  3. Используйте multiprocessing модуль для распределения работы на нескольких ядрах.

1 и 2 являются простыми. Для третьего пункта я предлагаю использовать imap_unordered() для достижения максимальной скорости с рабочим процессом, подобным этому:

  • получить все перестановки
  • в основной петле загружаются все модели городов в одном вызове хранилища данных
  • распределить расчет расстояния до рабочего пула
  • попробуйте сохранить результаты в одиночном работнике или в основном потоке. Я подозреваю, что было бы лучше в рабочем, но я не знаю, как django будет обрабатывать параллелизм в автономном скрипте (если кто-нибудь знает, пожалуйста, добавьте комментарии, чтобы мы могли интегрироваться).
  • либо, если вы сохраняете основной поток или отдельных работников, попробуйте сохранить большие куски моделей с использованием транзакций.

Но

Вам также нужно немного изменить свои модели. Для распределенной обработки вам нужно отделить переменную closest_cities . Поскольку различные процессы изменят его. Вы можете использовать словарь списков на основном уровне процесса, сохраняя все ближайшие города для любого города в качестве ключа, а затем сохраняя его для каждой модели, t конец цикла или между тем.

Вы делаете много ненужной работы.

Как и другие, вы можете ограничить количество вычислений, изменив структуру циклов. У тебя есть:

 for city in cities: for tested_city in cities: 

Таким образом, вы не только сравните каждый город с собой, вы также сравните city1 с city2 а затем сравните city2 с city1 .

Я не программист на Python, поэтому я не могу сказать вам, какой синтаксис использовать здесь, но то, что вы хотите, это структура вложенного цикла, подобная:

 for (i = 0; i < cities.Length-1; ++i) { for (j = i+1; j < cities.Length; ++j) { compare_cities(cities[i], cities[j]); } } 

Это позволит сократить вдвое количество сравниваемых городов. Это уменьшает его с примерно 3 миллиардов расстояний до примерно 1,5 миллиарда.

Другие также упоминали ранний потенциал, сравнив dlat и dlong прежде чем попасть в дорогостоящие триггерные функции.

Вы также можете сэкономить некоторое время, переведя lat1 и lon1 в радианы один раз, а также вычислив cos(lat1) один раз и передав эти значения на расчёт расстояния, а не вычисляя их каждый раз. Например:

 for (i = 0; i < cities.Length-1; ++i) { lat1 = radians(cities[i].latitude lon1 = radians(cities[i].longitude cos1 = cos(lat1) for (j = i+1; j < cities.Length; ++j) { compare_cities(lat1, lon1, cos1, cities[j]); } } 

И вам не нужно конвертировать c в километры. Например, у вас есть:

 return round(6373.0 * c, 2) 

Результат должен быть <= 30.0 . Почему умножение и округление? Вы можете просто return c , а в вашем коде сравнить возвращаемое значение с 0.0047 (что равно 30.0/6373 ).

Python - лучший язык программирования в мире.