[Numpy-discussion] Generalized ufuncs?
Engel, Hans-Andreas
Hans-Andreas.Engel@deshaw....
Mon Aug 18 09:50:16 CDT 2008
"Robert Kern" <robert.kern@gmail.com> wrote in message
news:<3d375d730808172015t62739266hbee205b97bdc86c3@mail.gmail.com>...
> On Sun, Aug 17, 2008 at 21:55, Anne Archibald
<peridot.faceted@gmail.com> wrote:
> > 2008/8/17 Robert Kern <robert.kern@gmail.com>:
> >>
> >> I suggested that we move it to a branch for the time being so we
can
> >> play with it and come up with examples of its use. If you have
> >> examples that you have already written, I would love to see them.
I,
> >> for one, am amenable to seeing this in 1.2.0, but only if we push
back
> >> the release by at least a few weeks.
> >
> > This is not a worked example, but this is exactly what is needed to
> > make possible the "arrays of matrices" functions that were discussed
> > some time ago. For example, users frequently want to do something
like
> > multiply an array of matrices by an array of matrices; that is not
> > currently feasible without a very large temporary array (or a loop).
> >
> > But I think you were looking for examples of code using the
interface,
> > to see whether it worked well.
>
> I'll take what I can get. In order of increasing utility:
>
> 1. Descriptions of use cases (as above).
> 2. Mockups of the Python code demonstrating the use case (e.g.
> nonfunctional Python code showing a potential generalized ufunc
> with inputs and outputs).
> 3. The C code implementing the actual functionality for the use
case.
>
> --
> Robert Kern
>
> "I have come to believe that the whole world is an enigma, a harmless
> enigma that is made terrible by our own mad attempt to interpret it as
> though it had an underlying truth."
> -- Umberto Eco
Please find an example on "inner1d" in the following.
1. One use case for an inner function is described on
http://scipy.org/scipy/numpy/wiki/GeneralLoopingFunctions and
http://thread.gmane.org/gmane.comp.python.numeric.general/17694.
(Another one would be the "array of matrices" usage mentioned above; we
have called this "dot2d" with signature "(m,n),(n,p)->(m,p)" on
http://scipy.org/scipy/numpy/ticket/887: here the matrix multiplication
would occur with respect to the last two dimensions.)
2. The mockup python code would be:
>>> from numpy import *
>>> N = 10
>>> a = random.randn(3, 5, N)
>>> b = random.randn(5, N)
>>> # standard inner function
... inner(a, b).shape
(3, 5, 5)
>>> # new ufunc "inner1d" with signature "(i),(i)->()", satisfying
GeneralLoopingFunctions use case
... inner1d(a, b).shape
(3, 5)
3. Concrete implementation of inner1d in C:
/*
* This implements the function out = inner1d(in1, in2) with
* out[K] = sum_i { in1[K, i] * in2[K, i] }
* and multi-index K, as described on
* http://scipy.org/scipy/numpy/wiki/GeneralLoopingFunctions
* and on http://projects.scipy.org/scipy/numpy/ticket/887.
*/
static void
DOUBLE_inner1d(char **args, intp *dimensions, intp *steps, void *func)
{
/* Standard ufunc loop length and strides. */
intp dn = dimensions[0];
intp s0 = steps[0];
intp s1 = steps[1];
intp s2 = steps[2];
intp n;
/* Additional loop length and strides for dimension "i" in
elementary function. */
intp di = dimensions[1];
intp i_s1 = steps[3];
intp i_s2 = steps[4];
intp i;
/* Outer loop: equivalent to standard ufuncs */
for (n = 0; n < dn; n++, args[0] += s0, args[1] += s1, args[2] +=
s2) {
char *ip1 = args[0], *ip2 = args[1], *op = args[2];
/* Implement elementary function: out = sum_i { in1[i] * in2[i]
} */
npy_double sum = 0;
for (i = 0; i < di; i++) {
sum += (*(npy_double *)ip1) * (*(npy_double *)ip2);
ip1 += i_s1;
ip2 += i_s2;
}
*(npy_double *)op = sum;
}
}
/* Actually create the ufunc object */
static PyUFuncGenericFunction inner1d_functions[] = { DOUBLE_inner1d };
static void *inner1d_data[] = { (void *)NULL };
static char inner1d_sigs[] = { PyArray_DOUBLE, PyArray_DOUBLE,
PyArray_DOUBLE };
static void
addInner1d(PyObject *dictionary) {
PyObject *f = PyUFunc_FromFuncAndDataAndSignature(inner1d_functions,
inner1d_data, inner1d_sigs, 1,
2, 1, PyUFunc_None, "inner1d",
"inner on the last dimension and
broadcast on the other dimensions",
0,
"(i),(i)->()");
PyDict_SetItemString(dictionary, "inner1d", f);
}
More information about the Numpy-discussion
mailing list