"Robert Kern" <[EMAIL PROTECTED]> wrote in message news:<[EMAIL PROTECTED]>... > On Sun, Aug 17, 2008 at 21:55, Anne Archibald <[EMAIL PROTECTED]> wrote: > > 2008/8/17 Robert Kern <[EMAIL PROTECTED]>: > >> > >> 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); } _______________________________________________ Numpy-discussion mailing list [email protected] http://projects.scipy.org/mailman/listinfo/numpy-discussion
