Оптимизация фрагмента кода, включая несколько операторов np.multipy - PullRequest
Новая
       8

Оптимизация фрагмента кода, включая несколько операторов np.multipy

3 голосов
/ 10 января 2020

Я работаю над оптимизацией некоторого кода в модуле python. Я установил узкое место и являюсь фрагментом кода, который выполняет некоторые вычисления в numpy. А именно следующий код: 

    xh = np.multiply(K_Rinv[0, 0], x )
    xh += np.multiply(K_Rinv[0, 1],  y) 
    xh += np.multiply(K_Rinv[0, 2],  h)
    yh = np.multiply(K_Rinv[1, 0],  x) 
    yh += np.multiply(K_Rinv[1, 1],  y) 
    yh += np.multiply(K_Rinv[1, 2],  h)
    q = np.multiply(K_Rinv[2, 0],  x) 
    q += np.multiply(K_Rinv[2, 1],  y) 
    q += np.multiply(K_Rinv[2, 2],  h)

, где x, y и h - np.ndarray с формой (4206,5749), а K_Rinv - np.ndarray с формой (3,3). Этот фрагмент кода вызывается несколько раз и занимает более 50% времени всего кода. Есть ли способ ускорить это? Или это так, как есть и не может быть ускорено.

Edit1:
Спасибо за ответы. После проблем с Numba (см. Сообщение об ошибке в конце) я попробовал предложение с Numberxpr. Однако мой код сломался при использовании этого решения. Так что я проверил, если результаты, где то же самое, и они не. Вот код, который я использую: 

    xh_1 = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[0, 0], 'b1': x,
                                                'a2': K_Rinv[0, 1], 'b2': y,
                                                'a3': K_Rinv[0, 2], 'b3': h})
    yh_1 = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[1, 0], 'b1': x,
                                                'a2': K_Rinv[1, 1], 'b2': y,
                                                'a3': K_Rinv[1, 2], 'b3': h})
    q_1 = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[2, 0], 'b1': x,
                                            'a2': K_Rinv[2, 1], 'b2': y,
                                            'a3': K_Rinv[2, 2], 'b3': h})
    xh_2 = np.multiply(K_Rinv[0, 0], x )
    xh_2 += np.multiply(K_Rinv[0, 1],  y) 
    xh_2 += np.multiply(K_Rinv[0, 2],  h)
    yh_2 = np.multiply(K_Rinv[1, 0],  x) 
    yh_2 += np.multiply(K_Rinv[1, 1],  y) 
    yh_2 += np.multiply(K_Rinv[1, 2],  h)
    q_2 = np.multiply(K_Rinv[2, 0],  x) 
    q_2 += np.multiply(K_Rinv[2, 1],  y) 
    q_2 += np.multiply(K_Rinv[2, 2],  h)
    check1 = xh_1.all() == xh_2.all() 
    check2 = yh_1.all() == yh_2.all() 
    check3 = q_2.all() == q_1.all()
    print ( " Check1 :{} , Check2: {} , Check3:{}" .format (check1,check2,check3))

У меня нет опыта работы с Numberxpr. Обычно это не одно и то же?

Ошибка от numba: 

 File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 420, in _compile_for_args
    raise e
  File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 353, in _compile_for_args
    return self.compile(tuple(argtypes))
  File "/usr/local/lib/python3.6/dist-packages/numba/compiler_lock.py", line 32, in _acquire_compile_lock
    return func(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 768, in compile
    cres = self._compiler.compile(args, return_type)
  File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 77, in compile
    status, retval = self._compile_cached(args, return_type)
  File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 91, in _compile_cached
    retval = self._compile_core(args, return_type)
  File "/usr/local/lib/python3.6/dist-packages/numba/dispatcher.py", line 109, in _compile_core
    pipeline_class=self.pipeline_class)
  File "/usr/local/lib/python3.6/dist-packages/numba/compiler.py", line 551, in compile_extra
    return pipeline.compile_extra(func)
  File "/usr/local/lib/python3.6/dist-packages/numba/compiler.py", line 327, in compile_extra
    raise e
  File "/usr/local/lib/python3.6/dist-packages/numba/compiler.py", line 321, in compile_extra
    ExtractByteCode().run_pass(self.state)
  File "/usr/local/lib/python3.6/dist-packages/numba/untyped_passes.py", line 67, in run_pass
    bc = bytecode.ByteCode(func_id)
  File "/usr/local/lib/python3.6/dist-packages/numba/bytecode.py", line 215, in __init__
    self._compute_lineno(table, code)
  File "/usr/local/lib/python3.6/dist-packages/numba/bytecode.py", line 237, in _compute_lineno
    known = table[_FIXED_OFFSET].lineno
KeyError: 2

Edit2 Резервуары для комментариев и ответов. Так что после того, как вы ознакомитесь с кодом, решение nuxx работает. Большое спасибо. Я все еще проводил некоторое тестирование в отдельном файле python, чтобы проверить, работает ли там код numba, но работает ли он очень медленно. Посмотрите код, который я использую ниже: 

import numpy as np
import numba as nb
import numexpr
from datetime import datetime

def calc(x,y,h,K_Rinv):
        xh_2 = np.multiply(K_Rinv[0, 0], x )
        xh_2 += np.multiply(K_Rinv[0, 1],  y) 
        xh_2 += np.multiply(K_Rinv[0, 2],  h)
        yh_2 = np.multiply(K_Rinv[1, 0],  x) 
        yh_2 += np.multiply(K_Rinv[1, 1],  y) 
        yh_2 += np.multiply(K_Rinv[1, 2],  h)
        q_2 = np.multiply(K_Rinv[2, 0],  x) 
        q_2 += np.multiply(K_Rinv[2, 1],  y) 
        q_2 += np.multiply(K_Rinv[2, 2],  h)
        return xh_2, yh_2, q_2

def calc_numexpr(x,y,h,K_Rinv):
    xh = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[0, 0], 'b1': x,
                                            'a2': K_Rinv[0, 1], 'b2': y,
                                            'a3': K_Rinv[0, 2], 'b3': h})
    yh = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[1, 0], 'b1': x,
                                            'a2': K_Rinv[1, 1], 'b2': y,
                                            'a3': K_Rinv[1, 2], 'b3': h})
    q = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[2, 0], 'b1': x,
                                           'a2': K_Rinv[2, 1], 'b2': y,
                                           'a3': K_Rinv[2, 2], 'b3': h})
    return xh,yh,q


@nb.njit(fastmath=True,parallel=True)
def calc_nb(x,y,h,K_Rinv):
    xh=np.empty_like(x)
    yh=np.empty_like(x)
    q=np.empty_like(x)

    for i in nb.prange(x.shape[0]):
        for j in range(x.shape[1]):
            xh[i,j]=K_Rinv[0, 0]*x[i,j]+K_Rinv[0, 1]*  y[i,j]+K_Rinv[0, 2]*h[i,j]
            yh[i,j]=K_Rinv[1, 0]*x[i,j]+K_Rinv[1, 1]*  y[i,j]+K_Rinv[1, 2]*h[i,j]
            q[i,j] =K_Rinv[2, 0]*x[i,j]+K_Rinv[2, 1]*  y[i,j]+K_Rinv[2, 2]*h[i,j]
    return xh,yh,q


x = np.random.random((4206, 5749))
y = np.random.random((4206, 5749))
h = np.random.random((4206, 5749))
K_Rinv = np.random.random((3, 3))

start = datetime.now()
x_calc,y_calc,q_calc = calc(x,y,h,K_Rinv)
end = datetime.now()
print("Calc took:           {} ".format(end - start))

start = datetime.now()
x_numexpr,y_numexpr,q_numexpr = calc_numexpr(x,y,h,K_Rinv)
end = datetime.now()
print("Calc_numexpr took:   {} ".format(end - start))

start = datetime.now()
x_nb,y_nb,q_nb = calc_nb(x,y,h,K_Rinv)
end = datetime.now()
print("Calc nb took:        {} ".format(end - start))

check_nb_q = (q_calc==q_nb).all()
check_nb_y = (y_calc==y_nb).all()
check_nb_x = (x_calc==x_nb).all()

check_numexpr_q = (q_calc==q_numexpr).all()
check_numexpr_y = (y_calc==y_numexpr).all()
check_numexpr_x = (x_calc==x_numexpr).all()

print("Checks for numexpr: {} , {} ,{} \nChecks for nb: {} ,{}, {}" .format(check_numexpr_x,check_numexpr_y,check_numexpr_q,check_nb_x,check_nb_y,check_nb_q))

, который превосходит следующее: 

Calc took:           0:00:01.944150 
Calc_numexpr took:   0:00:00.616224 
Calc nb took:        0:00:01.553058 
Checks for numexpr: True , True ,True 
Checks for nb: False ,False, False

, поэтому версия numba не работает должным образом. Любая идея, что я делаю не так? Хотелось бы, чтобы решение numba также работало.

Ps. nb. версия равна '0.47.0'

Ответы [ 3 ]

2 голосов
/ 10 января 2020

Другой возможностью было бы использование Numba.

Пример 

import numpy as np
import numba as nb

@nb.njit(fastmath=True,parallel=True)
def calc_nb(x,y,h,K_Rinv):
    xh=np.empty_like(x)
    yh=np.empty_like(x)
    q=np.empty_like(x)

    for i in nb.prange(x.shape[0]):
        for j in range(x.shape[1]):
            xh[i,j]=K_Rinv[0, 0]*x[i,j]+K_Rinv[0, 1]*  y[i,j]+K_Rinv[0, 2]*h[i,j]
            yh[i,j]=K_Rinv[1, 0]*x[i,j]+K_Rinv[1, 1]*  y[i,j]+K_Rinv[1, 2]*h[i,j]
            q[i,j] =K_Rinv[2, 0]*x[i,j]+K_Rinv[2, 1]*  y[i,j]+K_Rinv[2, 2]*h[i,j]
    return xh,yh,q

Ограничена ли эта пропускная способность памяти вычислений? 

def copy(x,y,h,K_Rinv):
    return np.copy(x),np.copy(y),np.copy(h)

%timeit copy(x,y,h,K_Rinv)
#147 ms ± 4.98 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

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

Времена 

x = np.random.random((4206, 5749))
y = np.random.random((4206, 5749))
h = np.random.random((4206, 5749))
K_Rinv = np.random.random((3, 3))

%timeit calc(x,y,h,K_Rinv) #Your implementation
#581 ms ± 8.05 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit calc_nb(x,y,h,K_Rinv)
#145 ms ± 3.81 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit calc_numexpr_scleronomic(x,y,h,K_Rinv)
#175 ms ± 1.83 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit calc_Daniel_F(x,y,h,K_Rinv)
#589 ms ± 24.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Возможна дальнейшая оптимизация: повторное использование уже выделенной памяти 

@nb.njit(fastmath=True,parallel=True)
def calc_nb_2(x,y,h,K_Rinv,xh,yh,q):
    for i in nb.prange(x.shape[0]):
        for j in range(x.shape[1]):
            xh[i,j]=K_Rinv[0, 0]*x[i,j]+K_Rinv[0, 1]*  y[i,j]+K_Rinv[0, 2]*h[i,j]
            yh[i,j]=K_Rinv[1, 0]*x[i,j]+K_Rinv[1, 1]*  y[i,j]+K_Rinv[1, 2]*h[i,j]
            q[i,j] =K_Rinv[2, 0]*x[i,j]+K_Rinv[2, 1]*  y[i,j]+K_Rinv[2, 2]*h[i,j]
    return xh,yh,q

#allocate memory only once if you call this function repeatedly
xh=np.empty_like(x)
yh=np.empty_like(x)
q=np.empty_like(x)

%timeit calc_nb_2(x,y,h,K_Rinv,xh,yh,q)
69.2 ms ± 194 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 голосов
/ 10 января 2020

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

import numpy as np
import numexpr

x = np.random.random((4206, 5749))
y = np.random.random((4206, 5749))
h = np.random.random((4206, 5749))
K_Rinv = np.random.random((3, 3))

xh = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[0, 0], 'b1': x,
                                            'a2': K_Rinv[0, 1], 'b2': y,
                                            'a3': K_Rinv[0, 2], 'b3': h})
yh = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[1, 0], 'b1': x,
                                            'a2': K_Rinv[1, 1], 'b2': y,
                                            'a3': K_Rinv[1, 2], 'b3': h})
q = numexpr.evaluate('a1*b1+a2*b2+a3*b3', {'a1': K_Rinv[2, 0], 'b1': x,
                                           'a2': K_Rinv[2, 1], 'b2': y,
                                           'a3': K_Rinv[2, 2], 'b3': h})

На моей машине это примерно в 5 раз быстрее, чем без numexpr.

Другое дело, что я бы предпочитайте использовать умножение матриц и широковещание с numpy, если вы имеете дело с матрицами, вместо того, чтобы разбивать умножения и сложения: 

xyh_mat = np.concatenate([x[:, :, np.newaxis],
                          y[:, :, np.newaxis],
                          h[:, :, np.newaxis]], axis=-1)[:, :, :, np.newaxis]  
# (4206, 5749, 3, 1)
K_Rinv_mat = K_Rinv[np.newaxis, np.newaxis, :, :]  
# (1, 1, 3, 3)


xyh_mat_2 = np.einsum("ijkl, ijlk->ijk", K_Rinv_mat, xyh_mat)
# 1.25x faster

xyh_mat_2 = K_Rinv_mat @ xyh_mat
# 3x slower

# xh = xyh_mat_2[:, :, 0]
# yh = xyh_mat_2[:, :, 1]
# q = xyh_mat_2[:, :, 2]

Однако, похоже, что при использовании * 1012 выигрыша в скорости нет. в этом случае, что немного удивляет меня.

РЕДАКТИРОВАТЬ

относительно комментария к дальнейшим вычислениям: 

np.divide(xh, q, x)
np.divide(yh, q, y)
# should translate to:
x = numexpr.evaluate('a/b', {'a': xh , 'b': q })
y = numexpr.evaluate('a/b', {'a': yh , 'b': q })
1 голос
/ 10 января 2020

Уверен, что это просто и расширенный точечный продукт: 

x_y_h = np.stack([x, y, h], axis = 0)
xh_yh_q = np.einsum('ij, jkl -> ikl', K_Rinv, x_y_h)
[xh, yh, q] = list(xh_yh_q)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...