[Numpy-discussion] generalized ufunc problem

[David Warde-Farley]

Hi Warren,

Thanks for the reply. I actually have a very similar version in my own  
code; I was just hoping to figure out the generalized ufunc  
architecture. There aren't many examples of actual uses of this  
capability in NumPy, so I wanted to try and exercise it a bit.  
logsumexp is kind of a perfect scenario for generalized ufuncs as it  
requires operations over an entire dimension (the max operation).

Cheers,

David



On 21-Jan-10, at 7:30 PM, Warren Weckesser wrote:

> 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 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);
>>
>> ....
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
>
>
> from numpy import *
> from scipy.maxentropy import logsumexp
>
> from my_logsumexp import my_logsumexp
>
>
> if __name__ == "__main__":
>
>    #
>    #--------------------------------------------
>    # 1D tests.
>    #--------------------------------------------
>    x = array([1.0,2.0,3.0])
>    p = logsumexp(x)
>    q = my_logsumexp(x)
>    assert p == q
>
>    x = array([1.0,2.0,3.0])
>    p = logsumexp(x)
>    q = my_logsumexp(x, axis=0)
>    assert p == q
>
>    #--------------------------------------------
>    # A 2D test.
>    #--------------------------------------------
>
>    a = array([[1.0,2.0,3.0],[0.1,0.2,0.3]])
>    q = my_logsumexp(a, axis=0)
>    assert allclose(q[0], logsumexp(a[:,0]))
>    assert allclose(q[1], logsumexp(a[:,1]))
>    assert allclose(q[2], logsumexp(a[:,2]))
>    q = my_logsumexp(a, axis=1)
>    assert allclose(q[0], logsumexp(a[0]))
>    assert allclose(q[1], logsumexp(a[1]))
>
>    #--------------------------------------------
>    # A 3D test.
>    #--------------------------------------------
>    L = 3
>    M = 4
>    N = 5
>    w = random.random((L, M, N))
>    q0 = empty((M,N))
>    for j in range(M):
>        for k in range(N):
>            q0[j,k] = logsumexp(w[:,j,k])
>    q1 = empty((L,N))
>    for i in range(L):
>        for k in range(N):
>            q1[i,k] = logsumexp(w[i,:,k])
>    q2 = empty((L,M))
>    for i in range(L):
>        for j in range(M):
>            q2[i,j] = logsumexp(w[i,j,:])
>
>    assert allclose(q0, my_logsumexp(w, axis=0))
>    assert allclose(q1, my_logsumexp(w, axis=1))
>    assert allclose(q2, my_logsumexp(w, axis=2))
>    #-------------------------------------------- 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion