2008/5/11 Robert Kern <[EMAIL PROTECTED]>:
> Perhaps, but ufuncs only allow 0 or 1 for this value, currently.
That's a shame, minus infinity is the identity for maximum too.
> Also, I was wrong about using PyUFunc_ff_f. Instead, use PyUFunc_ff_f_As_dd_d.
Hmm. Well, I tried implementing logsum(), and it turns out neither
PyUFunc_dd_d nor any of the other functions - mostly flagged as being
part of the ufunc API and described in the automatically-generated
ufunc_api.txt - are exported. All are static and not described in any
header I could find. That said, I know scipy defines its own ufuncs,
so it must either reimplement these or have some way I hadn't thought
of to get at them.
I've attached a patch to add the ufunc logsum to numpy. It's a bit
nasty for several reasons:
* I wasn't sure where to put it, so I put it in _compiled_base.c in
numpy/lib/src. I was very hesitant to touch umathmodule.c.src with its
sui generis macro language.
* I'm not sure what to do about log1p - it seems to be available in
spite of HAVE_LOG1P not being defined. In any case if it doesn't exist
it seems crazy to implement it again here.
* since PyUFunc_dd_d does not seem to be exported, I just cut-n-pasted
its code here. Obviously not a solution, but what are extension
writers supposed to do?
Do we want a ufunc version of logsum() in numpy? I got the impression
that in spite of some people's doubts as to its utility there was a
strong contingent of users who do want it. This plus a logdot() would
cover most of the operations one might want to do on numbers in log
representation. (Maybe subtraction?). It could be written in pure
python, for the most part, and reduce() would save a few exps and logs
at the cost of a temporary, but accumulate() remains a challenge. (I
don't expect many people feel the need for reduceat()).
Anne
P.S. it was surprisingly difficult to persuade numpy to find my tests.
I suppose that's what the proposed transition to nose is for? -A
Index: numpy/lib/src/_compiled_base.c
===================================================================
--- numpy/lib/src/_compiled_base.c (revision 5155)
+++ numpy/lib/src/_compiled_base.c (working copy)
@@ -1,6 +1,7 @@
#include "Python.h"
#include "structmember.h"
#include "numpy/noprefix.h"
+#include "numpy/ufuncobject.h"
static PyObject *ErrorObject;
#define Py_Try(BOOLEAN) {if (!(BOOLEAN)) goto fail;}
@@ -835,16 +836,68 @@
return;
}
+#if !defined(HAVE_LOG1P) && 0
+/* logsum ufunc */
+static double log1p(double x) /* is this really necessary? */
+{
+ double u = 1. + x;
+ if (u == 1.0) {
+ return x;
+ } else {
+ return log(u) * x / (u-1.);
+ }
+}
+#endif HAVE_LOG1P
+
+
+static double logsum(double x, double y) {
+ if (x>y) {
+ return x+log1p(exp(y-x));
+ } else {
+ return y+log1p(exp(x-y));
+ }
+};
+static void
+logsum_dd_d(char **args, intp *dimensions, intp *steps)
+{
+ intp n = dimensions[0];
+ intp is1 = steps[0];
+ intp is2 = steps[1];
+ intp os = steps[2];
+ char *ip1 = args[0];
+ char *ip2 = args[1];
+ char *op = args[2];
+ intp i;
+
+ for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op += os) {
+ double *in1 = (double *)ip1;
+ double *in2 = (double *)ip2;
+ double *out = (double *)op;
+
+ *out = logsum(*in1, *in2);
+ }
+}
+
+
+static char logsum_sigs[] = {
+ NPY_FLOAT64, NPY_FLOAT64, NPY_FLOAT64
+};
+static PyUFuncGenericFunction logsum_functions[] = {logsum_dd_d};
+static void* logsum_data[] = {logsum};
+static char logsum_doc[] = "Compute log(exp(x)+exp(y)) without numeric overflows and underflows";
+
+
/* Initialization function for the module (*must* be called init<name>) */
PyMODINIT_FUNC init_compiled_base(void) {
- PyObject *m, *d, *s;
+ PyObject *m, *d, *s, *f;
/* Create the module and add the functions */
m = Py_InitModule("_compiled_base", methods);
/* Import the array objects */
import_array();
+ import_umath();
/* Add some symbolic constants to the module */
d = PyModule_GetDict(m);
@@ -857,6 +910,23 @@
PyDict_SetItemString(d, "error", ErrorObject);
Py_DECREF(ErrorObject);
+ /* set up logsum ufunc */
+ f = PyUFunc_FromFuncAndData(
+ logsum_functions,
+ logsum_data,
+ logsum_sigs,
+ 1, // The number of type signatures.
+ 2, // The number of inputs.
+ 1, // The number of outputs.
+ PyUFunc_None, // The identity element for reduction.
+ // No good one to use for this function,
+ // unfortunately.
+ "logsum", // The name of the ufunc.
+ logsum_doc,
+ 0 // Dummy for API backwards compatibility.
+ );
+ PyDict_SetItemString(d, "logsum", f);
+ Py_DECREF(f);
/* define PyGetSetDescr_Type and PyMemberDescr_Type */
define_types();
Index: numpy/lib/tests/test_ufunclike.py
===================================================================
--- numpy/lib/tests/test_ufunclike.py (revision 5155)
+++ numpy/lib/tests/test_ufunclike.py (working copy)
@@ -58,9 +58,45 @@
"""
from numpy.testing import *
+#set_package_path() # is this needed?
+import numpy as np
+#restore_path()
class TestDocs(NumpyTestCase):
def check_doctests(self): return self.rundocs()
+class TestLogsum(NumpyTestCase):
+ def check_normal_values(self):
+ for x,y in [(1,1),
+ (1,2),
+ (1,1e-20),
+ (1,1+1e-12)]:
+ assert_almost_equal(x+y,np.exp(np.logsum(np.log(x),np.log(y))))
+ assert_equal(np.exp(np.logsum(np.log(y),np.log(x))),
+ np.exp(np.logsum(np.log(x),np.log(y))))
+
+ def check_tiny_values(self):
+ const = -1000
+ for x,y in [(1,1),
+ (1,2),
+ (1e20,1e-20),
+ (1,1+1e-12)]:
+ assert_almost_equal(const+np.log(x+y),
+ np.logsum(np.log(x)+const,np.log(y)+const))
+
+ def check_reduce_accumulate(self):
+ a = [1,10,2,0.1,0.01,100]
+ assert_almost_equal(np.add.reduce(a),
+ np.exp(np.logsum.reduce(np.log(a))))
+ assert_almost_equal(np.add.accumulate(a),
+ np.exp(np.logsum.accumulate(np.log(a))))
+
+ def check_zeros(self):
+ x = 0
+ y = 1
+ assert_almost_equal(x+y,np.exp(np.logsum(np.log(x),np.log(y))))
+ assert_equal(np.exp(np.logsum(np.log(y),np.log(x))),
+ np.exp(np.logsum(np.log(x),np.log(y))))
+
if __name__ == "__main__":
NumpyTest().run()
Index: numpy/lib/function_base.py
===================================================================
--- numpy/lib/function_base.py (revision 5155)
+++ numpy/lib/function_base.py (working copy)
@@ -9,7 +9,7 @@
'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning',
'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc',
'add_docstring', 'meshgrid', 'delete', 'insert', 'append',
- 'interp'
+ 'interp', 'logsum'
]
import warnings
@@ -25,7 +25,7 @@
from numpy.core.numerictypes import typecodes, number
from numpy.lib.shape_base import atleast_1d, atleast_2d
from numpy.lib.twodim_base import diag
-from _compiled_base import _insert, add_docstring
+from _compiled_base import _insert, add_docstring, logsum
from _compiled_base import digitize, bincount, interp as compiled_interp
from arraysetops import setdiff1d
import numpy as np
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
