[Numpy-discussion] Efficient numpy slicing for a "sliding window approach".

[Nicolas Pinto]

Thanks a lot for the pointer to segmentaxis. I'm trying to use it "as is"
and it seems that I need to a big reshape before the matrix multiplication.
Am I missing something ?


========================================

import numpy as np
from numpy import dot, transpose

arrh, arrw, arrd = 480,640,96
arr = np.random.randn(arrh, arrw, arrd).astype("float32")

stride = 16
winh, winw, wind = 128,64,96

limit = 100

clas_w = np.random.randn(8,4,96).astype("float32").ravel()

from segmentaxis import segment_axis

nh, nw = arrh-winh+1, arrw-winw+1

@profile
def func_loop(arr):

    resps = np.empty((nh,nw), dtype="float32")
    for j in xrange(nh):
         for i in xrange(nw):
            win = arr[j:j+winh:stride, i:i+winw:stride]
            win = win.ravel()
            resp = dot(win, clas_w)
            resps[j,i] = resp

    resps = resps.ravel()
    print resps.mean()


@profile
def func_segment(arr):

    arr = segment_axis(arr, winh, winh-1, axis=0)
    arr = arr[:, ::stride]

    arr = segment_axis(arr, winw, winw-1, axis=2)
    arr = arr[:, :, :, ::stride]

    arr = transpose(arr, [0, 2, 1, 3, 4])
    arr = arr.reshape(-1, clas_w.size)

    resps2 = dot(arr, clas_w)

    resps2 = resps2.ravel()
    print resps2.mean()

    # ...

func_loop(arr)

func_segment(arr)


========================================


% python kernprof.py -l -v sliding_win_all.py
w-0.0500485824034
segmentaxis.py:94: UserWarning: Problem with ndarray creation forces copy.
  warnings.warn("Problem with ndarray creation forces copy.")
-0.0500485824034
Wrote profile results to sliding_win_all.py.lprof
Timer unit: 1e-06 s

File: sliding_win_all.py
Function: func_loop at line 18
Total time: 3.9785 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    18                                           @profile
    19                                           def func_loop(arr):
    20
    21         1           28     28.0      0.0      resps =
np.empty((nh,nw), dtype="float32")
    22       354          308      0.9      0.0      for j in xrange(nh):
    23    204034       199753      1.0      5.0           for i in
xrange(nw):
    24    203681       691915      3.4     17.4              win =
arr[j:j+winh:stride, i:i+winw:stride]
    25    203681      1341174      6.6     33.7              win =
win.ravel()
    26    203681      1417998      7.0     35.6              resp = dot(win,
clas_w)
    27    203681       326520      1.6      8.2              resps[j,i] =
resp
    28
    29         1            2      2.0      0.0      resps = resps.ravel()
    30         1          805    805.0      0.0      print resps.mean()

File: sliding_win_all.py
Function: func_segment at line 33
Total time: 3.82026 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    33                                           @profile
    34                                           def func_segment(arr):
    35
    36         1           43     43.0      0.0      arr = segment_axis(arr,
winh, winh-1, axis=0)
    37         1            4      4.0      0.0      arr = arr[:, ::stride]
    38
    39         1       546505 546505.0     14.3      arr = segment_axis(arr,
winw, winw-1, axis=2)
    40         1           12     12.0      0.0      arr = arr[:, :, :,
::stride]
    41
    42         1           12     12.0      0.0      arr = transpose(arr,
[0, 2, 1, 3, 4])
    43         1      2435047 2435047.0     63.7      arr = arr.reshape(-1,
clas_w.size)
    44
    45         1       837700 837700.0     21.9      resps2 = dot(arr,
clas_w)
    46
    47         1           41     41.0      0.0      resps2 = resps2.ravel()
    48         1          892    892.0      0.0      print resps2.mean()


On Fri, Feb 20, 2009 at 11:55 PM, David Cournapeau <cournape@gmail.com>wrote:

> On Sat, Feb 21, 2009 at 1:46 PM, Nicolas Pinto <pinto@mit.edu> wrote:
> > Dear all,
> >
> > I'm trying to optimize the code below and I was wondering if there is an
> > efficient method that could reduce the numpy slicing overheard without
> going
> > with cython. Is there anyway I could use mgrid to get a matrix with all
> my
> > "windows" and then do a large matrix multiply instead?
>
> If you only care about removing the two loops for the per-window
> processing, Anne Archibald and Robert Kern wrote a very useful
> function, segment_axis, which is like the matlab buffer function on
> steroids, using numpy stride tricks (to avoid copies in many cases). I
> think that would do everything you want, right ? I use it a lot in my
> own code, it may be worth being included in numpy or scipy proper.
>
>
> http://projects.scipy.org/scipy/scikits/browser/trunk/talkbox/scikits/talkbox/tools/segmentaxis.py
>
> David
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>


Thanks again!

-- 
Nicolas Pinto
Ph.D. Candidate, Brain & Computer Sciences
Massachusetts Institute of Technology, USA
http://web.mit.edu/pinto
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20090221/c03623b8/attachment-0001.html