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 >> @t...@_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 >> > > > > 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@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion