Как векторизовать эту разность циклов в numpy?

Question

Я чувствую, что должен быть быстрый способ ускорить этот код. Я думаю, что ответ here, но я не могу найти свою проблему в этом формате. Основная проблема, которую я пытаюсь решить, - найти точную разницу в терминах параллельных и перпендикулярных компонентов и создать 2d-гистограмму этих различий.

out = np.zeros((len(rpbins)-1,len(pibins)-1)) 
tmp = np.zeros((len(x),2)) 
for i in xrange(len(x)): 
    tmp[:,0] = x - x[i] 
    tmp[:,1] = y - y[i] 

    para = np.sum(tmp**2,axis=-1)**(1./2) 
    perp = np.abs(z - z[i]) 

    H, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins]) 
    out += H 

Answer 1

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

mask = np.triu_indices(x.shape[0], 1) 
para = np.sqrt((x[:, None] - x)**2 + (y[:, None] - y)**2) 
perp = np.abs(z[:, None] - z) 
hist, _, _ = np.histogram2d(para[mask], perp[mask], bins=[rpbins, pibins]) 

mask, чтобы избежать подсчета каждого расстояния в два раза. Я также установил смещение по диагонали до 1, чтобы избежать включения в гистограмму расстояний каждой точки к себе. Но если вы не указали para и perp, вы получите тот же результат, что и ваш код.

С помощью этой выборки данных:

items = 100 
rpbins, pibins = np.linspace(0, 1, 3), np.linspace(0, 1, 3) 
x = np.random.rand(items) 
y = np.random.rand(items) 
z = np.random.rand(items) 

Я получаю это для моего hist и ваш out:

>>> hist 
array([[ 1795., 651.], 
     [ 1632., 740.]]) 
>>> out 
array([[ 3690., 1302.], 
     [ 3264., 1480.]]) 

и out[i, j] = 2 * hist[i, j] для i = j = 0, за исключением случаев out[0, 0] = 2 * hist[0, 0] + items из-за 0 расстояния от каждого элемента к сам.

---

EDIT Пробовал следующее после комментария tcaswell в:

items = 1000 
rpbins, pibins = np.linspace(0, 1, 3), np.linspace(0, 1, 3) 
x, y, z = np.random.rand(3, items) 

def hist1(x, y, z, rpbins, pibins) : 
    mask = np.triu_indices(x.shape[0], 1) 
    para = np.sqrt((x[:, None] - x)**2 + (y[:, None] - y)**2) 
    perp = np.abs(z[:, None] - z) 
    hist, _, _ = np.histogram2d(para[mask], perp[mask], bins=[rpbins, pibins]) 
    return hist 

def hist2(x, y, z, rpbins, pibins) : 
    mask = np.triu_indices(x.shape[0], 1) 
    para = np.sqrt((x[:, None] - x)[mask]**2 + (y[:, None] - y)[mask]**2) 
    perp = np.abs((z[:, None] - z)[mask]) 
    hist, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins]) 
    return hist 

def hist3(x, y, z, rpbins, pibins) : 
    mask = np.triu_indices(x.shape[0], 1) 
    para = np.sqrt(((x[:, None] - x)**2 + (y[:, None] - y)**2)[mask]) 
    perp = np.abs((z[:, None] - z)[mask]) 
    hist, _, _ = np.histogram2d(para, perp, bins=[rpbins, pibins]) 
    return hist 

In [10]: %timeit -n1 -r10 hist1(x, y, z, rpbins, pibins) 
1 loops, best of 10: 289 ms per loop 

In [11]: %timeit -n1 -r10 hist2(x, y, z, rpbins, pibins) 
1 loops, best of 10: 294 ms per loop 

In [12]: %timeit -n1 -r10 hist3(x, y, z, rpbins, pibins) 
1 loops, best of 10: 278 ms per loop 

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