Blitz код производит различные выходные

Question

Я хочу использовать weave.blitz, чтобы улучшить производительность следующей Numpy коды:

def fastIteration(self): 
    g = self.grid 
    nx,ny = g.ux.shape 

    uxold = g.old_ux 
    ux = g.ux 
    ux[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2]) 

    g.setBC() 
    g.old_ux = ux.copy() 

В этом коде г является расчетной сеткой. Которые состоят из двух разных полей ux и uxold. Старое просто используется для временного хранения переменных. В полном коде приблизительно 95% времени выполнения тратится на метод fastIteration, поэтому даже простое увеличение производительности уменьшит количество часов, затрачивающих значительное количество времени на выполнение этого кода.

Выхода методы Numpy выглядит так, как будто:

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

def blitzIteration(self): 
    ### does not work correct so far 
    g = self.grid 
    nx,ny = g.ux.shape 

    uxold = g.old_ux 
    ux = g.ux 
    expr = "ux[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2])" 
    weave.blitz(expr, check_size=0) 
    g.setBC() 
    g.old_ux = ux.copy() 

Однако это не дает правильный результат: (. Воспроизведенным, поданной и fixed Там больше информации о фактической ошибке там)

Answer 1

Это похоже на ошибку в weave.blitz ,

Мне показалось странным написать 0: вместо более короткого :, чтобы получить полный фрагмент, чтобы я заменил все эти фрагменты и вуаля, это сработало.

Я действительно не знаю, где лежит ошибка, но expr_code генерируется weave.blitz немного отличается:

• При использовании 0:

  ipdb> expr_code 
  'ux_blitz_buggy(blitz::Range(0,_end),blitz::Range(1,Nux_blitz_buggy(1)-1-1))=uxold(blitz::Range(0,_end),blitz::Range(1,Nuxold(1)-1-1))+ReI*(uxold(blitz::Range(0,_end),blitz::Range(2,_end))-2*uxold(blitz::Range(0,_end),blitz::Range(1,Nuxold(1)-1-1))+uxold(blitz::Range(0,_end),blitz::Range(0,Nuxold(1)-2-1)));\n' 
  

• При использовании :

  ipdb> expr_code 
  'ux_blitz_not_buggy(_all,blitz::Range(1,Nux_blitz_not_buggy(1)-1-1))=uxold(_all,blitz::Range(1,Nuxold(1)-1-1))+ReI*(uxold(_all,blitz::Range(2,_end))-2*uxold(_all,blitz::Range(1,Nuxold(1)-1-1))+uxold(_all,blitz::Range(0,Nuxold(1)-2-1)));\n' 
  

Итак, blitz::Range(0,_end) становится _all, и они ведут себя по-другому.

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

import numpy as np 
from scipy.weave import blitz 


def test_blitz_bug(N=4): 
    ReI = 1.2 
    ux_blitz_buggy, ux_blitz_not_buggy, ux_np = np.zeros((N, N)), np.zeros((N, N)), np.zeros((N, N)) 
    uxold = np.random.randn(N, N) 
    ux_np[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2]) 
    expr_buggy = 'ux_blitz_buggy[0:,1:-1] = uxold[0:,1:-1] + ReI* (uxold[0:,2:] - 2*uxold[0:,1:-1] + uxold[0:,0:-2])' 
    expr_not_buggy = 'ux_blitz_not_buggy[:,1:-1] = uxold[:,1:-1] + ReI* (uxold[:,2:] - 2*uxold[:,1:-1] + uxold[:,0:-2])' 
    blitz(expr_buggy) 
    blitz(expr_not_buggy) 
    assert not np.allclose(ux_blitz_buggy, ux_np) 
    assert np.allclose(ux_blitz_not_buggy, ux_np) 

if __name__ == '__main__': 
    test_blitz_bug()