[Numpy-discussion] Example Usage of Neighborhood Iterator in Cython

[eat]

Hi,

On Mon, Oct 17, 2011 at 9:19 PM, T J <tjhnson at gmail.com> wrote:

> I recently put together a Cython example which uses the neighborhood
> iterator.  It was trickier than I thought it would be, so I thought to
> share it with the community.  The function takes a 1-dimensional array
> and returns a 2-dimensional array of neighborhoods in the original
> area. This is somewhat similar to the functionality provided by
> segment_axis (http://projects.scipy.org/numpy/ticket/901), but I
> believe this slightly different in that neighborhood can extend to the
> left of the current index (assuming circular boundaries).  Keep in
> mind that this is just an example, and normal usage probably is not
> concerned with creating a new array.
>
> External link:  http://codepad.org/czRIzXQl
>
> --------------
>
> import numpy as np
> cimport numpy as np
>
> cdef extern from "numpy/arrayobject.h":
>
>    ctypedef extern class numpy.flatiter [object PyArrayIterObject]:
>        cdef int nd_m1
>        cdef np.npy_intp index, size
>        cdef np.ndarray ao
>        cdef char *dataptr
>
>    # This isn't exposed to the Python API.
>    # So we can't use the same approach we used to define flatiter
>    ctypedef struct PyArrayNeighborhoodIterObject:
>        int nd_m1
>        np.npy_intp index, size
>        np.PyArrayObject *ao # note the change from np.ndarray
>        char *dataptr
>
>    object PyArray_NeighborhoodIterNew(flatiter it, np.npy_intp* bounds,
>                                       int mode, np.ndarray fill_value)
>    int PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject *it)
>    int PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject *it)
>
>    object PyArray_IterNew(object arr)
>    void PyArray_ITER_NEXT(flatiter it)
>    np.npy_intp PyArray_SIZE(np.ndarray arr)
>
>    cdef enum:
>        NPY_NEIGHBORHOOD_ITER_ZERO_PADDING,
>        NPY_NEIGHBORHOOD_ITER_ONE_PADDING,
>        NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING,
>        NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING,
>        NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING
>
> np.import_array()
>
> def windowed(np.ndarray[np.int_t, ndim=1] arr, bounds):
>
>    cdef flatiter iterx = <flatiter>PyArray_IterNew(<object>arr)
>    cdef np.npy_intp size = PyArray_SIZE(arr)
>    cdef np.npy_intp* boundsPtr = [bounds[0], bounds[1]]
>    cdef int hoodSize = bounds[1] - bounds[0] + 1
>
>    # Create the Python object and keep a reference to it
>    cdef object niterx_ = PyArray_NeighborhoodIterNew(iterx,
>        boundsPtr, NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING, None)
>    cdef PyArrayNeighborhoodIterObject *niterx = \
>        <PyArrayNeighborhoodIterObject *>niterx_
>
>    cdef int i,j
>    cdef np.ndarray[np.int_t, ndim=2] hoods
>
>    hoods = np.empty((arr.shape[0], hoodSize), dtype=np.int)
>    for i in range(iterx.size):
>        for j in range(niterx.size):
>            hoods[i,j] = (niterx.dataptr)[0]
>            PyArrayNeighborhoodIter_Next(niterx)
>        PyArray_ITER_NEXT(iterx)
>        PyArrayNeighborhoodIter_Reset(niterx)
>    return hoods
>
> def test():
>    x = np.arange(10)
>    print x
>    print
>    print windowed(x, [-1, 3])
>    print
>    print windowed(x, [-2, 2])
>
>
> ----------
>
> If you run test(), this is what you should see:
>
> [0 1 2 3 4 5 6 7 8 9]
>
> [[9 0 1 2 3]
>  [0 1 2 3 4]
>  [1 2 3 4 5]
>  [2 3 4 5 6]
>  [3 4 5 6 7]
>  [4 5 6 7 8]
>  [5 6 7 8 9]
>  [6 7 8 9 0]
>  [7 8 9 0 1]
>  [8 9 0 1 2]]
>
> [[8 9 0 1 2]
>  [9 0 1 2 3]
>  [0 1 2 3 4]
>  [1 2 3 4 5]
>  [2 3 4 5 6]
>  [3 4 5 6 7]
>  [4 5 6 7 8]
>  [5 6 7 8 9]
>  [6 7 8 9 0]
>  [7 8 9 0 1]]
>
> windowed(x, [0, 2]) is almost like segment_axis(x, 3, 2, end='wrap').
>
Just wondering what are the main benefits, of your approach, comparing to
simple:
In []: a= arange(5)
In []: n= 10
In []: b= arange(n)[:, None]
In []: mod(a+ roll(b, 1), n)
Out[]:
array([[9, 0, 1, 2, 3],
       [0, 1, 2, 3, 4],
       [1, 2, 3, 4, 5],
       [2, 3, 4, 5, 6],
       [3, 4, 5, 6, 7],
       [4, 5, 6, 7, 8],
       [5, 6, 7, 8, 9],
       [6, 7, 8, 9, 0],
       [7, 8, 9, 0, 1],
       [8, 9, 0, 1, 2]])
In []: mod(a+ roll(b, 2), n)
Out[]:
array([[8, 9, 0, 1, 2],
       [9, 0, 1, 2, 3],
       [0, 1, 2, 3, 4],
       [1, 2, 3, 4, 5],
       [2, 3, 4, 5, 6],
       [3, 4, 5, 6, 7],
       [4, 5, 6, 7, 8],
       [5, 6, 7, 8, 9],
       [6, 7, 8, 9, 0],
       [7, 8, 9, 0, 1]])

Regards,
eat

> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20111017/01d2c44e/attachment.html>