# [Numpy-discussion] generalized ufunc problem

David Warde-Farley dwf at cs.toronto.edu
Thu Jan 21 19:13:17 EST 2010

```I decided to take a crack at adding a generalized ufunc for logsumexp,
i.e. collapsed an array along the last dimension by subtracting the
maximum element E along that dimension, taking the exponential,
logaddexp.reduce() but presumably faster and less prone to error
accumulation.

I added the following to umath_tests.c.src and got everything to
compile, but for some reason it doesn't give me the behaviour I'm
looking for. I was expecting a (500, 50) array to be collapsed down to
a (500,) array. Is that not what the signature calls for?

Thanks,

David

char *logsumexp_signature = "(i)->()";

/**begin repeat

#TYPE=LONG,DOUBLE#
#typ=npy_long, npy_double#
#EXPFUN=expl, exp#
#LOGFUN=logl, log#
*/

/*
*  This implements the function
*        out[n] = sum_i { in1[n, i] * in2[n, i] }.
*/

static void
@TYPE at _logsumexp(char **args, intp *dimensions, intp *steps, void
*NPY_UNUSED(func))
{
INIT_OUTER_LOOP_3
intp di = dimensions[0];
intp i;
intp is1=steps[0];
BEGIN_OUTER_LOOP_3
char *ip1=args[0], *op=args[1];
@typ@ max = (*(@typ@ *)ip1);
@typ@ sum = 0;

for (i = 0; i < di; i++) {
max = max < (*(@typ@ *)ip1) ? (*(@typ@ *)ip1) : max;
ip1 += is1;
}
ip1 = args[0];
for (i = 0; i < di; i++) {
sum += @EXPFUN@((*(@typ@ *)ip1) - max);
ip1 += is1;
}
*(@typ@ *)op = @LOGFUN@(sum + max);
END_OUTER_LOOP
}

/**end repeat**/

static PyUFuncGenericFunction logsumexp_functions[] =
{ LONG_logsumexp, DOUBLE_logsumexp };
static void * logsumexp_data[] = { (void *)NULL, (void *)NULL };
static char logsumexp_signatures[] = { PyArray_LONG, PyArray_LONG,
PyArray_DOUBLE, PyArray_DOUBLE };

...

f = PyUFunc_FromFuncAndData(logsumexp_functions, logsumexp_data,
logsumexp_signatures, 2,                                    1, 1,
PyUFunc_None, "logsumexp",                                    "inner1d
with a weight argument \\n""\"          \"\n""                  \\
\"(i)->()\\\" \\n""", 0);    PyDict_SetItemString(dictionary,
"logsumexp", f);    Py_DECREF(f);

....

```