Neal Becker wrote:
> I need to write functions that apply element-wise numpy arrays, and should
> be independent of the number of dimensions.
>
> numpy has a mechanism for generic (dimension-independent) iteration over
> elements of the array, but is there any way to use this (or implement
> something similar) in cython?
There's two methods: ufuncs and multi-iterators. ufuncs is the most
common one -- search for ufunc on the Cython wiki, there should be an
example somewhere.
I just did this using multi-iterators (which I had to do because I
needed some context -- a temporary buffer -- during the calculation).
Attaching my code for calculating a Wigner 3j symbol elementwise with
full broadcasting.
(It is not as fast as it could be -- it is possible to extract out the
inner loop and have it as a C loop rather than rely on the iterator,
which is faster. I didn't bother yet though, mainly because each Wigner
symbol calculation is rather heavy.)
The NumPy docs is rather good in this area, check there next! (search
for the multiiterator functions I use + ufuncs)
Also you need a newer numpy.pxd than the one shipped with Cython for
both ufuncs and multiiterators. Just replace the one in Cython/Includes
with the attachment. I will include it in the next Cython release.
[Edit: The attachments were too large, but they're up at
sage.math.washington.edu/home/dagss as numpy.pxd and wigner.pyx. The
main idea is:
cdef np.broadcast it = np.broadcast(l1, l2, l3, m1, m2, m3)
shape = (<object>it).shape
if out is None:
out = np.empty(shape, dtype)
else:
if out.shape != shape:
raise ValueError(u"out argument has incorrect shape")
if out.dtype != dtype:
raise ValueError(u"out argument has incorrect dtype")
it = np.broadcast(l1, l2, l3, m1, m2, m3, out)
with nogil:
while np.PyArray_MultiIter_NOTDONE(it):
l1val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 0))[0]
l2val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 1))[0]
l3val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 2))[0]
m1val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 3))[0]
m2val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 4))[0]
m3val = (<dtype_t*>np.PyArray_MultiIter_DATA(it, 4))[0]
...make outval...
(<dtype_t*>np.PyArray_MultiIter_DATA(it, 6))[0] = outval
np.PyArray_MultiIter_NEXT(it)
]
--
Dag Sverre
_______________________________________________
Cython-dev mailing list
[email protected]
http://codespeak.net/mailman/listinfo/cython-dev
