Effective numpy.roll to numpy.sum () or mean ()

Question

I have several (about 1000) 3D form arrays (1000, 800, 1024) that I want to study. I need to calculate the average value along the axis = 0, but before I can do this, I need to collapse the data along axis 2 until it is “in the right place”.

That sounds weird, so I'll try to explain. The 1D sub-array of form (1024) is data from the physical ring buffer. The ring buffer is read in different ways that I know. Therefore, I have several pos arrays of the form (1000, 800). Telling me where the ring buffer stuck. And my 3D arrays of data form (1000, 800, 1024) that I need to roll over pos .

Only after rolling the 3D arrays make sense to me, and I can begin to analyze them. In C, you can write code that does this quite simply, so I wonder if I can somehow "talk" about the numpy mean () or sum () methods, which they should start with different indexes and roll around at the end 1D subarray.

I am currently doing the following:

 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]) 

It takes ~ 60 seconds. And then I do, for example:

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

It only takes 15 seconds.

My point is that creating a thumbnail takes a lot of space and time (well, I could save space by writing the thumbnail back to data ), although there is definitely a way (in C) to implement this averaging and rolling in a single loop , which saves a lot of time. My question is ... if there is a way to do something like this with numpy?

Answer 1

I got bored and wrote your function in Cython. It runs about 10 times faster than the code you sent without allocating an intermediate array.

 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 

Note that this uses the naive mean + stdv algorithm , so you may encounter floating point precision errors. However, my tests did not show a significant effect.