Mise en œuvre vectorisée pour un script Numpy alambiqué

Question

Je travaille sur un projet personnel qui consiste à prédire les mouvements des modèles météorologiques à partir du radar. J'ai trois n tableaux de numpy m; une avec des valeurs d'intensité de précipitation, une avec le mouvement (en pixels) dans la direction X de cette précipitation et l'autre avec le mouvement (en pixels) dans la direction Y de cette précipitation. Je veux utiliser ces trois tableaux pour déterminer l'emplacement des pixels de précipitation en utilisant les décalages dans les deux autres réseaux.

xMax = currentReflectivity.shape[0] 
yMax = currentReflectivity.shape[1]        
for x in xrange(currentReflectivity.shape[0]): 
    for y in xrange(currentReflectivity.shape[1]): 
     targetPixelX = xOffsetArray[x,y] + x 
     targetPixelY = yOffsetArray[x,y] + y 
     targetPixelX = int(targetPixelX) 
     targetPixelY = int(targetPixelY) 
     if targetPixelX < xMax and targetPixelY < yMax: 
      interpolatedReflectivity[targetPixelX,targetPixelY] = currentReflectivity[x,y] 

Je ne peux pas imaginer un moyen de vectoriser ceci; des idées?

Answer 1

est ici une approche vectorisé faisant usage de broadcasting -

x_arr = np.arange(currentReflectivity.shape[0])[:,None] 
y_arr = np.arange(currentReflectivity.shape[1]) 

targetPixelX_arr = (xOffsetArray[x_arr, y_arr] + x_arr).astype(int) 
targetPixelY_arr = (yOffsetArray[x_arr, y_arr] + y_arr).astype(int) 

valid_mask = (targetPixelX_arr < xMax) & (targetPixelY_arr < yMax) 

R = targetPixelX_arr[valid_mask] 
C = targetPixelY_arr[valid_mask] 
interpolatedReflectivity[R,C] = currentReflectivity[valid_mask] 

Runtime tests

approches -

def org_app(currentReflectivity, xOffsetArray, yOffsetArray): 
    m,n = currentReflectivity.shape 
    interpolatedReflectivity = np.zeros((m,n)) 

    xMax = currentReflectivity.shape[0] 
    yMax = currentReflectivity.shape[1]        
    for x in xrange(currentReflectivity.shape[0]): 
     for y in xrange(currentReflectivity.shape[1]): 
      targetPixelX = xOffsetArray[x,y] + x 
      targetPixelY = yOffsetArray[x,y] + y 
      targetPixelX = int(targetPixelX) 
      targetPixelY = int(targetPixelY) 

      if targetPixelX < xMax and targetPixelY < yMax: 
       interpolatedReflectivity[targetPixelX,targetPixelY] = \ 
       currentReflectivity[x,y] 
    return interpolatedReflectivity 

def broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray): 
    m,n = currentReflectivity.shape 
    interpolatedReflectivity = np.zeros((m,n)) 

    xMax, yMax = m,n   
    x_arr = np.arange(currentReflectivity.shape[0])[:,None] 
    y_arr = np.arange(currentReflectivity.shape[1]) 

    targetPixelX_arr = (xOffsetArray[x_arr, y_arr] + x_arr).astype(int) 
    targetPixelY_arr = (yOffsetArray[x_arr, y_arr] + y_arr).astype(int) 

    valid_mask = (targetPixelX_arr < xMax) & (targetPixelY_arr < yMax)  
    R = targetPixelX_arr[valid_mask] 
    C = targetPixelY_arr[valid_mask] 
    interpolatedReflectivity[R,C] = currentReflectivity[valid_mask] 
    return interpolatedReflectivity 

minutage et vérification -

In [276]: # Setup inputs 
    ...: m,n = 100,110 # currentReflectivity.shape 
    ...: max_r = 120 # xOffsetArray's extent 
    ...: max_c = 130 # yOffsetArray's extent 
    ...: 
    ...: currentReflectivity = np.random.rand(m, n) 
    ...: xOffsetArray = np.random.randint(0,max_r,(m, n)) 
    ...: yOffsetArray = np.random.randint(0,max_c,(m, n)) 
    ...: 

In [277]: out1 = org_app(currentReflectivity, xOffsetArray, yOffsetArray) 
    ...: out2 = broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray) 
    ...: print np.allclose(out1, out2) 
    ...: 
True 

In [278]: %timeit org_app(currentReflectivity, xOffsetArray, yOffsetArray) 
100 loops, best of 3: 6.86 ms per loop 

In [279]: %timeit broadcasting_app(currentReflectivity, xOffsetArray, yOffsetArray) 
1000 loops, best of 3: 212 µs per loop 

In [280]: 6860.0/212  # Speedup number 
Out[280]: 32.35849056603774 

Answer 2

Je suis assez sûr que vous pouvez vectoriser ce simplement en prenant tout de la boucle:

targetPixelX = (xOffsetArray + np.arange(xMax).reshape(xMax, 1)).astype(np.int) 
targetPixelY = (yOffsetArray + np.arange(yMax)).astype(np.int) 
mask = ((targetPixelX < xMax) & (targetPixelY < yMax)) 
interpolatedReflectivity[mask] = currentReflectivity[mask] 

Ce sera beaucoup plus rapide mais plus gourmande en mémoire. Fondamentalement, targetPixelX et targetPixelY sont maintenant des tableaux contenant les valeurs pour chaque pixel qui ont été préalablement calculées sur une base par itération.

Seules les valeurs masquées sont définies dans interpolatedReflectivity, de façon similaire à ce que la déclaration if faisait dans la boucle.