[Numpy-discussion] custom accumlators

[Jeff Whitaker]

Matt Knox wrote:
> I made a post about this a while ago on the scipy-user mailing list, but I didn't receive much of a response so I'm just throwing it out there again (with more detail) in case it got overlooked.
>
> Basically, I'd like to be able to do accumulate operations with custom functions. numpy.vectorize does not seem to provide an accumulate method with the functions it returns. I'm hoping I don't have to write ufuncs in C to accomplish this, but I fear that may the case. Either way, it would be nice to know if it can, or cannot be done in an easy manner.
>
> I have lots of examples of where this kind of thing is useful, but I'll just outline two for now.
>
> Assume the parameter x in all the functions below is a 1-d array
>
> ---------------------------------------------
> Example 1 - exponential moving average:
>
> # naive brute force method...
> def expmave(x, k):
>     result = numpy.array(x, copy=True)
>     for i in range(1, result.size):
>        result[i] = result[i-1] + k * (result[i] - result[i-1])
>     return result
>
> # slicker method (if it worked, which it doesn't)...
> def expmave(x, k):
>     def expmave_sub(a, b):
>         return a + k * (b - a)
>     return numpy.vectorize(expmave_sub).accumulate(x)
>     
> ---------------------------------------------
> Example 2 - forward fill a masked array:
>
> # naive brute force method...
> def forward_fill(x):
>     result = ma.array(x, copy=True)
>     for i in range(1, result.size):
>        if result[i] is ma.masked: result[i] = result[i-1]
>     return result
>
> # slicker method (if it worked, which it doesn't)...
> def forward_fill(x):
>     def forward_fill_sub(a, b):
>     	if b is ma.masked: return a
>     	else: return b
>     return numpy.vectorize(forward_fill_sub).accumulate(x)
> ---------------------------------------------
>
> Is their a good way to do these kinds of things without python looping? Or is that going to require writing a ufunc in C? Any help is greatly appreciated.
>
> Thanks,
>
> - Matt Knox
>   

Matt:  Here's a quick and dirty example of how to do this sort of thing 
in pyrex.  I do it all the time, and it works quite well.

# accumulator.pyx:
_doublesize = sizeof(double)
cdef extern from "Python.h":
    int PyObject_AsWriteBuffer(object, void **rbuf, Py_ssize_t *len)
    char *PyString_AsString(object)
def accumulator(object x, double k):
    cdef Py_ssize_t buflen
    cdef int ndim, i
    cdef double *xdata
    cdef void *xbuff
    # make a copy by casting to an array of doubles.
    x = x.astype('f8')
    # if buffer api is supported, get pointer to data buffers.
    if PyObject_AsWriteBuffer(x, &xbuff, &buflen) <> 0:
        raise RuntimeError('object does not support buffer API')
    ndim = buflen/_doublesize
    xdata = <double *>xbuff
    for i from 1 <= i < ndim:
        xdata[i] = xdata[i-1] + k * (xdata[i] - xdata[i-1])
    return x

# test.py
from accumulator import accumulator
from numpy import linspace
x = linspace(1.,10.,10)
k = 0.1
print x
x1 = accumulator(x,k)
print x1
def expmave(x, k):
    result = x.copy()
    for i in range(1, result.size):
        result[i] = result[i-1] + k * (result[i] - result[i-1])
    return result
x2 = expmave(x,k)
print x2 # should be the same as x1

# setup.py
import os
from distutils.core import setup, Extension
from Pyrex.Distutils import build_ext
setup(name = "accumulator",
  cmdclass          = {'build_ext': build_ext},
  keywords          = ["python","map projections","GIS","mapping","maps"],
  ext_modules = [Extension("accumulator",["accumulator.pyx"])])


to build, just do

python setup.py build_ext --inplace

then run test.py.

HTH,

-Jeff

-- 
Jeffrey S. Whitaker         Phone : (303)497-6313
NOAA/OAR/CDC  R/PSD1        FAX   : (303)497-6449
325 Broadway                Boulder, CO, USA 80305-3328