Как сделать быстрый цикл для расчета матрицы в python - PullRequest
Новая
       101

Как сделать быстрый цикл для расчета матрицы в python

1 голос
/ 21 марта 2020

У меня проблема в том, что для l oop требуется так много времени для завершения. Я хочу более быстрый способ завершить это. мой код: 

dx = 20
dy = 20
dz = 20

x = np.arange(0, 1201, dx)
y = np.arange(0, 1001, dy)
z = np.arange(20, 501, dz)
drho = 3000  # Delta Rho (Density Contrast) kg/m^3

# Input Rho to Model
M = np.zeros((len(z), len(x), len(y)))
M[6:16, 26:36, 15:25] = drho
m = np.array(M.flat)
# p (61, 1525)
# M(25, 61)
#  m(1525,)
# Station Position
stx, sty = np.meshgrid(x, y)
stx = np.array(stx.flat)
sty = np.array(sty.flat)
stz = np.zeros(len(stx))

# Make meshgrid
X, Y, Z = np.meshgrid(x, y, z)
X = np.array(X.flat)
Y = np.array(Y.flat)
Z = np.array(Z.flat)

p = np.zeros((len(stx), len(X)))

# p(3111, 77775)
for i in range(len(X)):
    for j in range(len(stx)):
        p[j, i] = (Z[i] - stz[j]) / ((Z[i] - stz[j]) ** 2 + (X[i] - stx[j]) ** 2 + (Y[i] - sty[j]) ** 2) ** (3/2)

Переменная итераций иногда превышает 241 миллион, и это происходит как всегда.

Ответы [ 2 ]

2 голосов
/ 21 марта 2020

Вы можете улучшить производительность здесь, используя широковещательную рассылку с: 

p = (Z - stz[:,None]) / ((Z - stz[:,None])**2  + (X - stx[:,None])**2 + (Y - sty[:,None])**2) ** (3/2)

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

Проверьте и время, используя вместо этого: 

x = np.arange(0, 121, dx)
y = np.arange(0, 101, dy)
z = np.arange(20, 51, dz)

def op():
    p = np.zeros((len(stx), len(X)))
    for i in range(len(X)):
        for j in range(len(stx)):
            p[j, i] = (Z[i] - stz[j]) / ((Z[i] - stz[j]) ** 2 + (X[i] - stx[j]) ** 2 + (Y[i] - sty[j]) ** 2) ** (3/2)
    return p

def ap_1():
    return (Z - stz[:,None]) / ((Z - stz[:,None])**2  + (X - stx[:,None])**2 + (Y - sty[:,None])**2) ** (3/2)

%timeit p = op()
# 44.5 ms ± 3.58 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit p_ = ap_1()
# 169 µs ± 2.27 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

np.allclose(p_, p)
# True

Мы можем еще больше повысить производительность и повысить эффективность памяти, позволив numexpr занять забота об арифметике c: 

import numexpr as ne

def ap_2():
    return ne.evaluate('(Z - stz2D) / ((Z - stz2D)**2  + (X - stx2D)**2 + (Y - sty2D)**2) ** (3/2)',
           {'stz2D':stz[:,None], 'stx2D':stx[:,None], 'sty2D':sty[:,None]})

%timeit ap_2()
# 106 µs ± 6.34 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Итак, при втором подходе мы получаем 420x ускорение

0 голосов
/ 21 марта 2020

Действительно классная вещь в python заключается в том, что функции map (), filter () и redu () сильно оптимизированы для больших наборов данных.

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

...