[Numpy-discussion] Generalised ufuncs branch

[Engel, Hans-Andreas]

"Stéfan van der Walt" <stefan at sun.ac.za> wrote in message news:<9457e7c80808151645v457e155apa192bdbcf3f91d47 at mail.gmail.com>...
Hi all,

I have moved the generalised ufuncs functionality off to a branch:

http://svn.scipy.org/svn/numpy/branches/gen_ufuncs

Please try it out and give us your feedback.  We shall also pound on
it at the sprint during SciPy'08, and thereafter decide how and when
to best integrate it into NumPy.

Thanks to Wenjie Fu and Hans-Andreas Engel for taking the time to
think this issue through and to submit such a high-quality patch.

Regards
Stéfan
----------


Motivated by Stefan's and Chuck's off-line comments, I have attached an example implementation of a generalized ufunc:


/*
 *  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);
}