Я пытаюсь оптимизировать фрагмент, который вызывается много (миллионы раз), поэтому любое улучшение скорости (надеюсь, удаление цикла for) было бы здорово.
Я вычисляю корреляционную функцию некоторой j-й частицы со всеми остальными
C_j(|r-r'|) = sqrt(E((s_j(r')-s_k(r))^2)) усредненный по k.
Моя идея состоит в том, чтобы иметь переменную corrfun, которая собирает данные в некоторые ячейки (г, определенный в другом месте). Я нахожу, к какой ячейке r принадлежит каждый s_k, и это сохраняется в ind. Таким образом, ind[0] — это индекс r (и, следовательно, corrfun), которому соответствует точка j=0. Несколько точек могут попасть в один и тот же бин (на самом деле я хочу, чтобы бины были достаточно большими, чтобы содержать несколько точек), поэтому я суммирую все (s_j(r')-s_k(r))^2, а затем делю на количество точки в этом бине (хранятся в переменной rw). Код, который я в итоге сделал для этого, следующий (np для numpy):
for k, v in enumerate(ind):
if j==k:
continue
corrfun[v] += (s[k]-s[j])**2
rw[v] += 1
rw2 = rw
rw2[rw < 1] = 1
corrfun = np.sqrt(np.divide(corrfun, rw2))
Обратите внимание, что дело с rw2 было связано с тем, что я хочу избежать проблем с делением на 0, но я возвращаю массив rw и хочу иметь возможность различать элементы rw=0 и rw=1. Возможно, для этого тоже есть более элегантное решение.
Есть ли способ ускорить цикл for? Хотя я не хотел бы добавлять взаимодействие с самим собой (j==k), я даже согласен с взаимодействием с самим собой, если это означает, что я могу получить значительно более быстрый расчет (длина ind ~ 1E6, поэтому взаимодействие с самим собой, вероятно, в любом случае незначительно).
Благодарю вас!
Илья
Редактировать:
Вот полный код. Обратите внимание, что в полном коде я также усредняю по j.
import numpy as np
def twopointcorr(x,y,s,dr):
width = np.max(x)-np.min(x)
height = np.max(y)-np.min(y)
n = len(x)
maxR = np.sqrt((width/2)**2 + (height/2)**2)
r = np.arange(0, maxR, dr)
print(r)
corrfun = r*0
rw = r*0
print(maxR)
''' go through all points'''
for j in range(0, n-1):
hypot = np.sqrt((x[j]-x)**2+(y[j]-y)**2)
ind = [np.abs(r-h).argmin() for h in hypot]
for k, v in enumerate(ind):
if j==k:
continue
corrfun[v] += (s[k]-s[j])**2
rw[v] += 1
rw2 = rw
rw2[rw < 1] = 1
corrfun = np.sqrt(np.divide(corrfun, rw2))
return r, corrfun, rw
Я отлаживаю тест следующим образом
from twopointcorr import twopointcorr
import numpy as np
import matplotlib.pyplot as plt
import time
n=1000
x = np.random.rand(n)
y = np.random.rand(n)
s = np.random.rand(n)
print('running two point corr functinon')
start_time = time.time()
r,corrfun,rw = twopointcorr(x,y,s,0.1)
print("--- Execution time is %s seconds ---" % (time.time() - start_time))
fig1=plt.figure()
plt.plot(r, corrfun,'-x')
fig2=plt.figure()
plt.plot(r, rw,'-x')
plt.show()
Опять же, основная проблема заключается в том, что в реальном наборе данных n~1E6. Конечно, я могу выполнить повторную выборку, чтобы сделать ее меньше, но мне бы хотелось на самом деле прокрутить набор данных.
range(0, n-1)вместоrange(0, n)? - person John Zwinck schedule 26.04.2016add.at- небуферизованной версии обычного дополнения. docs.scipy.org/doc/ numpy-1.10.1/ссылка/сгенерированный/ - person hpaulj schedule 26.04.2016