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
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
