Added log-domain summation

appengine/findit/libs/math/logarithms.py

Index: appengine/findit/libs/math/logarithms.py
diff --git a/appengine/findit/libs/math/logarithms.py b/appengine/findit/libs/math/logarithms.py
new file mode 100644
index 0000000000000000000000000000000000000000..6e22fda13345c1f19f9530f50e67ca093baf7e13
--- /dev/null
+++ b/appengine/findit/libs/math/logarithms.py
@@ -0,0 +1,42 @@
+# Copyright 2016 The Chromium Authors. All rights reserved.

[Martin Barbella, 2016/12/02 20:07:02]

Please run gpylint on these files.

[wrengr, 2016/12/02 21:00:58]

Done.

+# Use of this source code is governed by a BSD-style license that can be
+# found in the LICENSE file.
+
+import math
+
+_LOG_ZERO = float('-inf')
+
+def logsumexp(xs):
+  """Efficiently and accurately compute a log-domain sum.
+
+  While ``math.log(math.fsum(math.exp(x) for x in xs)))`` accomplishes
+  the same end, it requires an additional ``len(xs) - 2`` logarithms
+  and introduces more opportunities for over/underflow and rounding
+  errors. Thus, this function is both more efficient and more accurate.
+
+  Args:
+    xs (collection of float): the collection of log-domain numbers to be
+      summed. The order of the elements doesn't matter, nor does the
+      concrete type of the collection (it could be ``list``, ``set``,
+      ``frozenset``, etc); the collection just needs to support ``len``,
+      ``max``, and ``iter``. Because we need the maximum before we
+      can start the iteration, this function requires (in principle)
+      two passes over the data; therefore we cannot take an iterator nor
+      generator. The length should be computed in constant time (else it
+      will introduce an unnecessary third pass over the data). Ideally
+      the maximum can also be computed in constant time (in which case
+      this function only needs to pass over the data once), though none
+      of Python's builtin types support that.
+
+  Returns:
+    The log-domain sum of ``xs``.
+  """
+  if not xs:
+    return _LOG_ZERO
+
+  maximum = max(xs)
+  if math.isinf(maximum):
+    return maximum
+
+  total = math.fsum(math.expm1(x - maximum) for x in xs)
+  return maximum + math.log1p(total + float(len(xs) - 1))