Вопрос по python, list-comprehension, optimization – Оптимизация расчета расстояния Python с учетом периодических граничных условий

10

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

def PBCdist(coord1,coord2,UC):
    dx = coord1[0] - coord2[0]
    if (abs(dx) > UC[0]*0.5):
       dx = UC[0] - dx
    dy = coord1[1] - coord2[1]
    if (abs(dy) > UC[1]*0.5):
       dy = UC[1] - dy
    dz = coord1[2] - coord2[2]
    if (abs(dz) > UC[2]*0.5):
       dz = UC[2] - dz
    dist = np.sqrt(dx**2 + dy**2 + dz**2)
    return dist

Затем я вызываю функцию так

for i, coord2 in enumerate(coordlist):
  if (PBCdist(coord1,coord2,UC) < radius):
      do something with i

Недавно я прочитал, что могу значительно повысить производительность, используя понимание списка. Следующее работает для случая без PBC, но не для случая PBC

coord_indices = [i for i, y in enumerate([np.sqrt(np.sum((coord2-coord1)**2)) for coord2 in coordlist]) if y < radius]
for i in coord_indices:
   do something

Есть ли какой-нибудь способ сделать эквивалент этого для дела PBC? Есть ли альтернатива, которая будет работать лучше?

Координат представляет собой массив с формой ок (5711,3). Координат сам по себе получается из большого списка, поэтому я эффективно зацикливаю список координат 20000 раз, и этот список координат перебирается примерно 50 раз ... Вы получаете картину. johnjax
& Quot; Vectorising & Quot; в контексте NumPy означает, что вы используете Ufuncs NumPy для перемещения циклов, которые в противном случае вам пришлось бы делать в коде Python для кода C (как я делал в своем ответе). Я рекомендую немного почитать о NumPy, если вы беспокоитесь о производительности. Sven Marnach
Я посмотрел функцию векторизации в NumPy. В документации сказано:"The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop." johnjax
Вы используете NumPy, поэтому вы должны векторизовать цикл для повышения производительности. Какова именно структураcoordlist? Это должен быть двумерный массив NumPy, чтобы можно было оптимизировать цикл с помощью ufuncs NumPy. Sven Marnach

Ваш Ответ

4   ответа
0

высокопроизводительное руководство по питону, Он содержит множество предложений о том, куда вы можете пойти дальше.

В том числе:

vectorization cython pypy numexpr pycuda multiprocesing parallel python
0

meshgrid очень полезно для генерации расстояний. Например:

import numpy as np
row_diff, col_diff = np.meshgrid(range(7), range(8))
radius_squared = (row_diff - x_coord)**2 + (col_diff - y_coord)**2

У меня теперь есть массив (radius_squared) где каждая запись указывает квадрат расстояния от позиции массива[x_coord, y_coord].

Чтобы округлить массив, я могу сделать следующее:

row_diff, col_diff = np.meshgrid(range(7), range(8))
row_diff = np.abs(row_diff - x_coord)
row_circ_idx = np.where(row_diff > row_diff.shape[1] / 2)
row_diff[row_circ_idx] = (row_diff[row_circ_idx] - 
                         2 * (row_circ_idx + x_coord) + 
                         row_diff.shape[1])
row_diff = np.abs(row_diff)
col_diff = np.abs(col_diff - y_coord)
col_circ_idx = np.where(col_diff > col_diff.shape[0] / 2)
col_diff[row_circ_idx] = (row_diff[col_circ_idx] - 
                         2 * (col_circ,_idx + y_coord) + 
                         col_diff.shape[0])
col_diff = np.abs(row_diff)
circular_radius_squared = (row_diff - x_coord)**2 + (col_diff - y_coord)**2

Теперь у меня есть все расстояния массива, округленные с векторной математикой.

12

distance() функционировать таким образом, что вы можете векторизовать цикл над 5711 точками. Следующая реализация принимает массив точек какx0 или жеx1 параметр:

def distance(x0, x1, dimensions):
    delta = numpy.abs(x0 - x1)
    delta = numpy.where(delta > 0.5 * dimensions, delta - dimensions, delta)
    return numpy.sqrt((delta ** 2).sum(axis=-1))

Пример:

>>> dimensions = numpy.array([3.0, 4.0, 5.0])
>>> points = numpy.array([[2.7, 1.5, 4.3], [1.2, 0.3, 4.2]])
>>> distance(points, [1.5, 2.0, 2.5], dimensions)
array([ 2.22036033,  2.42280829])

В результате получается массив расстояний между точками, передаваемый в качестве второго параметраdistance() и каждая точка вpoints.

@Патрикnumpy.hypot() работает только для двух измерений, тогда как OP нужен код, который работает для трехмерных точек. И относительно знакаdeltaНу, если вам все равно, не стесняйтесь редактировать пост. :)
Нет никакой разницы после принятия нормы, но мне нравится получать правильный знак дельты, заменяя 3-ю строку наdelta = numpy.where(", delta - dimensions, "), Также обратите внимание, чтоnp.hypot лучше избежать переполнения, чем sqrt (сумма (x ** 2))
Это привело к увеличению скорости моего кода примерно в 5 раз. Спасибо! Для тех, кто ищет альтернативу, ответ от lazyr также дал хороший результат. johnjax
4
import numpy as np

bounds = np.array([10, 10, 10])
a = np.array([[0, 3, 9], [1, 1, 1]])
b = np.array([[2, 9, 1], [5, 6, 7]])

min_dists = np.min(np.dstack(((a - b) % bounds, (b - a) % bounds)), axis = 2)
dists = np.sqrt(np.sum(min_dists ** 2, axis = 1))

Вотa а такжеb списки векторов, которые вы хотите рассчитать расстояние между иbounds являются границами пространства (поэтому здесь все три измерения идут от 0 до 10, а затем переносятся). Он рассчитывает расстояния междуa[0] а такжеb[0], a[1] а такжеb[1], и так далее.

Я уверен, что опытные эксперты могли бы добиться большего успеха, но это, вероятно, будет на порядок быстрее, чем то, что вы делаете, поскольку большая часть работы в настоящее время выполняется в C.

@johnjax Для чего бы я ни стоил, я бы тоже принял ответ Свена Марнача, будь я на твоем месте. Он более применим, чем мой.
Я тоже попробовал этот подход. Это также привело к примерно 5-кратному улучшению кода. К сожалению, я могу проверить только 1 ответ как правильный: / johnjax

Похожие вопросы