[Numpy-discussion] generalized ufunc problem

[Warren Weckesser]

David,

I haven't tried creating a ufunc before, so I can't help you with that, 
but since you are working on logsumexp, you might be interested in the 
version I posted here in October:

    http://mail.scipy.org/pipermail/scipy-user/2009-October/022931.html

and the attached tests.


Warren


David Warde-Farley wrote:
> 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@_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);
>
> ....
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>   

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: my_logsumexp_test.py
Url: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20100121/6c64ab74/attachment.pl