Efficiently slice windows from a 1D numpy array, around indices given by second 2D array

Question

I want to extract multiple slices from the same 1D numpy array, where the slice indices are drawn from a random distribution. Basically, I want to achieve the following:

import numpy as np
import numpy.random

# generate some 1D data
data = np.random.randn(500)

# window size (slices are 2*winsize long)
winsize = 60

# number of slices to take from the data
inds_size = (100, 200)

# get random integers that function as indices into the data
inds = np.random.randint(low=winsize, high=len(data)-winsize, size=inds_size)

# now I want to extract slices of data, running from inds[0,0]-60 to inds[0,0]+60
sliced_data = np.zeros( (winsize*2,) + inds_size )
for k in range(inds_size[0]):
    for l in range(inds_size[1]):
        sliced_data[:,k,l] = data[inds[k,l]-winsize:inds[k,l]+winsize]

# sliced_data.shape is now (120, 100, 200)

The above nested loop works fine, but is very slow. In my real code, I will need to do this thousands of times, for data arrays a lot bigger than these. Is there any way to do this more efficiently?

Note that inds will always be 2D in my case, but after getting the slices I will always be summing over one of these two dimensions, so an approach that only accumulates the sum across the one dimension would be fine.

I found this question and this answer which seem almost the same. However, the question is only about a 1D indexing vector (as opposed to my 2D). Also, the answer lacks a bit of context, as I don't really understand how the suggested as_strided works. Since my problem does not seem uncommon, I thought I'd ask again in the hope of a more explanatory answer rather than just code.

Sammelsurium · Accepted Answer · 2015-12-10 00:28:56Z

6

Using as_strided in this way appears to be somewhat faster than Divakar's approach (20 ms vs 35 ms here), although memory usage might be an issue.

data_wins = as_strided(data, shape=(data.size - 2*winsize + 1, 2*winsize), strides=(8, 8))
inds = np.random.randint(low=0, high=data.size - 2*winsize, size=inds_size)
sliced = data_wins[inds]
sliced = sliced.transpose((2, 0, 1))    # to use the same index order as before

Strides are the steps in bytes for the index in each dimension. For example, with an array of shape (x, y, z) and a data type of size d (8 for float64), the strides will ordinarily be (y*z*d, z*d, d), so that the second index steps over whole rows of z items. Setting both values to 8, data_wins[i, j] and data_wins[j, i] will refer to the same memory location.

>>> import numpy as np
>>> from numpy.lib.stride_tricks import as_strided
>>> a = np.arange(10, dtype=np.int8)
>>> as_strided(a, shape=(3, 10 - 2), strides=(1, 1))
array([[0, 1, 2, 3, 4, 5, 6, 7],
       [1, 2, 3, 4, 5, 6, 7, 8],
       [2, 3, 4, 5, 6, 7, 8, 9]], dtype=int8)

answered Dec 10, 2015 at 0:28

Sammelsurium

5163 silver badges6 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

EelkeSpaak Over a year ago

Thanks for this! I now understand as_strided very well, love it. Also, your approach is faster than Divakar's approach (thanks to him too though!) by about a factor ~3 on my tests. Finally, memory usage is fine using as_strided, since it does not make copies of the data. Memory usage is identical as my own looped approach (I profiled using memory_profiler), whereas Divakar's broadcasting approach uses double memory (it stores all indices + the data).

Divakar · Accepted Answer · 2015-12-09 21:12:17Z

4

Here's a vectorized approach using broadcasting -

# Get 3D offsetting array and add to inds for all indices
allinds = inds + np.arange(-60,60)[:,None,None]

# Index into data with all indices for desired output
sliced_dataout = data[allinds]

Runtime test -

In [20]: # generate some 1D data
    ...: data = np.random.randn(500)
    ...: 
    ...: # window size (slices are 2*winsize long)
    ...: winsize = 60
    ...: 
    ...: # number of slices to take from the data
    ...: inds_size = (100, 200)
    ...: 
    ...: # get random integers that function as indices into the data
    ...: inds=np.random.randint(low=winsize,high=len(data)-winsize, size=inds_size)
    ...: 

In [21]: %%timeit 
    ...: sliced_data = np.zeros( (winsize*2,) + inds_size )
    ...: for k in range(inds_size[0]):
    ...:     for l in range(inds_size[1]):
    ...:         sliced_data[:,k,l] = data[inds[k,l]-winsize:inds[k,l]+winsize]
    ...: 
10 loops, best of 3: 66.9 ms per loop

In [22]: %%timeit 
    ...: allinds = inds + np.arange(-60,60)[:,None,None]
    ...: sliced_dataout = data[allinds]
    ...: 
10 loops, best of 3: 24.1 ms per loop

Memory consumption : Compromise solution

If memory consumption is an issue, here's a compromise solution with one loop -

sliced_dataout = np.zeros( (winsize*2,) + inds_size )
for k in range(sliced_data.shape[0]):
    sliced_dataout[k] = data[inds-winsize+k]

edited Dec 9, 2015 at 21:12

answered Dec 9, 2015 at 17:54

Divakar

222k19 gold badges273 silver badges374 bronze badges

4 Comments

EelkeSpaak Over a year ago

Thanks, that's in principle a good one and time-efficient. However, unfortunately it uses a lot of memory...

Divakar Over a year ago

@EelkeSpaak Hmm I guess you could alternatively inline it like so - sliced_dataout = data[inds + np.arange(-60,60)[:,None,None]].

Divakar Over a year ago

@EelkeSpaak Check out the just added Compromise solution.

EelkeSpaak Over a year ago

Thanks for your answer, and the compromise solution. I've accepted Sammelsurium's answer though, as it turns out that solution is the better option in all aspects (time and memory).

Collectives™ on Stack Overflow

Efficiently slice windows from a 1D numpy array, around indices given by second 2D array

2 Answers 2

1 Comment

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related