[Numpy-discussion] Generalized ufuncs?

[Engel, Hans-Andreas]

"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);
}