effiziente numpy.roll vor numpy.sum() oder mean()

Question

Ich habe mehrere (Bestellung von 1000) 3D-Arrays von Form (1000, 800, 1024) Ich möchte studieren. Ich muss den Mittelwert entlang der Achse = 0 berechnen, aber bevor ich das tun kann, muss ich die Daten entlang der Achse 2 rollen, bis sie "an der richtigen Stelle" ist.

Das klingt seltsam, also werde ich versuchen, es zu erklären. Das 1D-Sub-Array von Shape (1024,) ist Daten von einem physischen Ring-Puffer. Der Ringpuffer wird in verschiedenen Positionen ausgelesen, die ich kenne. Also habe ich mehrere Arrays pos von Form (1000, 800). Er sagte mir, in welcher Position der Ringpuffer ausgelesen wurde. Und meine 3D-Arrays data von Form (1000, 800, 1024), die ich nach pos rollen muss.

Erst nach dem Rollen .. sind die 3D-Arrays für mich sinnvoll und ich kann sie analysieren. In C kann man Code schreiben, der das ziemlich einfach macht, also frage ich mich, ob ich die numpy mean() - oder sum() - Routinen "erzählen" kann, sie sollten bei verschiedenen Indizes anfangen und am Ende von 'rollen' das 1D-Subarray.

Was ich derzeit tun, ist dies:

rolled = np.zeros_like(data) # shape (1000, 800, 1024) 
for a in range(rolled.shape[0]): 
    for b in range(rolled.shape[1]): 
     rolled[a,b] = np.roll(data[a,b], pos[a,b]) 

Dieser Vorgang dauert ~ 60sec Und dann mache ich z:

m = rolled.mean(axis=0) 
s = rolled.std(axis=0) 

, die nur 15 Sekunden oder so dauert.

Mein Punkt ist, dass das Erstellen der gerollten Kopie viel Platz und Zeit in Anspruch nimmt (okay, ich könnte Platz sparen, indem ich die gerollten Sachen zurück in data schreibe), während es definitiv einen Weg (in C) gibt, diese Mittelung zu implementieren und rollen in einer Schleife, wodurch viel Zeit gespart wird. Meine Frage ist ... wenn es eine Möglichkeit gibt, etwas Ähnliches mit Numpy zu machen?

Answer 1

Ich wurde gelangweilt und schrieb Ihre Funktion in Cython. Es läuft etwa 10x schneller als der Code, den Sie gepostet haben, ohne ein intermediäres Array zuzuweisen.

import numpy as np 
cimport numpy as np 
cimport cython 
from libc.math cimport sqrt 

@cython.boundscheck(False) 
@cython.wraparound(False) 
@cython.nonecheck(False) 
@cython.cdivision(True) 
def rolled_mean_std(double[:,:,::1] data, int[:,::1] pos): 
    cdef Py_ssize_t s1,s2,s3,a,b,c 
    s1 = data.shape[0] 
    s2 = data.shape[1] 
    s3 = data.shape[2] 
    cdef double[:,::1] sums = np.zeros((s2,s3)) 
    cdef double[:,::1] sumsq = np.zeros((s2,s3)) 
    cdef double d 
    cdef int p 
    # Compute sums and sum-of-squares. 
    for a in range(s1): 
    for b in range(s2): 
     p = pos[a,b] 
     for c in range(s3): 
     d = data[a,b,(c+s3-p)%s3] 
     sums[b,c] += d 
     sumsq[b,c] += d * d 
    # Calculate mean + std in place. 
    for b in range(s2): 
    for c in range(s3): 
     d = sums[b,c] 
     sums[b,c] /= s1 
     sumsq[b,c] = sqrt((s1*sumsq[b,c] - (d*d)))/s1 
    return sums, sumsq 

Beachten Sie, dass dies die Naive mean+stdv algorithm verwendet, so dass Sie in Gleitkommagenauigkeit Fehler führen könnten. Meine Tests zeigten jedoch keine große Wirkung.