I have a matrix (2d numpy ndarray, to be precise):
A = np.array([[4, 0, 0], [1, 2, 3], [0, 0, 5]])
And I want to roll each row of
A independently, according to roll values in another array:
r = np.array([2, 0, -1])
That is, I want to do this:
print np.array([np.roll(row, x) for row,x in zip(A, r)]) [[0 0 4] [1 2 3] [0 5 0]]
Is there a way to do this efficiently? Perhaps using fancy indexing tricks?
Sure you can do it using advanced indexing, whether it is the fastest way probably depends on your array size (if your rows are large it may not be):
rows, column_indices = np.ogrid[:A.shape, :A.shape] # Use always a negative shift, so that column_indices are valid. # (could also use module operation) r[r < 0] += A.shape column_indices = column_indices - r[:, np.newaxis] result = A[rows, column_indices]
numpy.lib.stride_tricks.as_strided stricks (abbrev pun intended) again!
Speaking of fancy indexing tricks, there’s the infamous –
np.lib.stride_tricks.as_strided. The idea/trick would be to get a sliced portion starting from the first column until the second last one and concatenate at the end. This ensures that we can stride in the forward direction as needed to leverage
np.lib.stride_tricks.as_strided and thus avoid the need of actually rolling back. That’s the whole idea!
Now, in terms of actual implementation we would use
scikit-image's view_as_windows to elegantly use
np.lib.stride_tricks.as_strided under the hoods. Thus, the final implementation would be –
from skimage.util.shape import view_as_windows as viewW def strided_indexing_roll(a, r): # Concatenate with sliced to cover all rolls a_ext = np.concatenate((a,a[:,:-1]),axis=1) # Get sliding windows; use advanced-indexing to select appropriate ones n = a.shape return viewW(a_ext,(1,n))[np.arange(len(r)), (n-r)%n,0]
Here’s a sample run –
In : A = np.array([[4, 0, 0], ...: [1, 2, 3], ...: [0, 0, 5]]) In : r = np.array([2, 0, -1]) In : strided_indexing_roll(A, r) Out: array([[0, 0, 4], [1, 2, 3], [0, 5, 0]])
# @seberg's solution def advindexing_roll(A, r): rows, column_indices = np.ogrid[:A.shape, :A.shape] r[r < 0] += A.shape column_indices = column_indices - r[:,np.newaxis] return A[rows, column_indices]
Let’s do some benchmarking on an array with large number of rows and columns –
In : np.random.seed(0) ...: a = np.random.rand(10000,1000) ...: r = np.random.randint(-1000,1000,(10000)) # @seberg's solution In : %timeit advindexing_roll(a, r) 10 loops, best of 3: 71.3 ms per loop # Solution from this post In : %timeit strided_indexing_roll(a, r) 10 loops, best of 3: 44 ms per loop
In case you want more general solution (dealing with any shape and with any axis), I modified @seberg’s solution:
def indep_roll(arr, shifts, axis=1): """Apply an independent roll for each dimensions of a single axis. Parameters ---------- arr : np.ndarray Array of any shape. shifts : np.ndarray How many shifting to use for each dimension. Shape: `(arr.shape[axis],)`. axis : int Axis along which elements are shifted. """ arr = np.swapaxes(arr,axis,-1) all_idcs = np.ogrid[[slice(0,n) for n in arr.shape]] # Convert to a positive shift shifts[shifts < 0] += arr.shape[-1] all_idcs[-1] = all_idcs[-1] - shifts[:, np.newaxis] result = arr[tuple(all_idcs)] arr = np.swapaxes(result,-1,axis) return arr
I implement a pure
numpy.lib.stride_tricks.as_strided solution as follows
from numpy.lib.stride_tricks import as_strided def custom_roll(arr, r_tup): m = np.asarray(r_tup) arr_roll = arr[:, [*range(arr.shape),*range(arr.shape-1)]].copy() #need `copy` strd_0, strd_1 = arr_roll.strides n = arr.shape result = as_strided(arr_roll, (*arr.shape, n), (strd_0 ,strd_1, strd_1)) return result[np.arange(arr.shape), (n-m)%n] A = np.array([[4, 0, 0], [1, 2, 3], [0, 0, 5]]) r = np.array([2, 0, -1]) out = custom_roll(A, r) Out: array([[0, 0, 4], [1, 2, 3], [0, 5, 0]])
Building on divakar’s excellent answer, you can apply this logic to 3D array easily (which was the problematic that brought me here in the first place). Here’s an example – basically flatten your data, roll it & reshape it after::
def applyroll_30(cube, threshold=25, offset=500): flattened_cube = cube.copy().reshape(cube.shape*cube.shape, cube.shape) roll_matrix = calc_roll_matrix_flattened(flattened_cube, threshold, offset) rolled_cube = strided_indexing_roll(flattened_cube, roll_matrix, cube_shape=cube.shape) rolled_cube = triggered_cube.reshape(cube.shape, cube.shape, cube.shape) return rolled_cube def calc_roll_matrix_flattened(cube_flattened, threshold, offset): """ Calculates the number of position along time axis we need to shift elements in order to trig the data. We return a 1D numpy array of shape (X*Y, time) elements """ # armax(...) finds the position in the cube (3d) where we are above threshold roll_matrix = np.argmax(cube_flattened > threshold, axis=1) + offset # ensure we don't have index out of bound roll_matrix[roll_matrix>cube_flattened.shape] = cube_flattened.shape return roll_matrix def strided_indexing_roll(cube_flattened, roll_matrix_flattened, cube_shape): # Concatenate with sliced to cover all rolls # otherwise we shift in the wrong direction for my application roll_matrix_flattened = -1 * roll_matrix_flattened a_ext = np.concatenate((cube_flattened, cube_flattened[:, :-1]), axis=1) # Get sliding windows; use advanced-indexing to select appropriate ones n = cube_flattened.shape result = viewW(a_ext,(1,n))[np.arange(len(roll_matrix_flattened)), (n - roll_matrix_flattened) % n, 0] result = result.reshape(cube_shape) return result
Divakar’s answer doesn’t do justice to how much more efficient this is on large cube of data. I’ve timed it on a 400x400x2000 data formatted as int8. An equivalent for-loop does ~5.5seconds, Seberg’s answer ~3.0seconds and strided_indexing…. ~0.5second.