[Numpy-discussion] generalized ufunc problem

[David Warde-Farley]

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,  
adding, and then adding back E. Functionally the same  
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 };


/* and inside addUfuncs() */
...

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

....