REF: Remove logsumexp_b and add argument to logsumexp

scipy/misc/__init__.py

@@ -31,7 +31,6 @@
    info - Get help information for a function, class, or module
    lena - Get classic image processing example image Lena
    logsumexp - Compute the log of the sum of exponentials of input elements
-   logsumexp_b - Compute the log of the sum of scaled exponentials of inputs
    pade - Pade approximation to function as the ratio of two polynomials
    radon -
    toimage - Takes a numpy array and returns a PIL image

scipy/misc/common.py

@@ -8,13 +8,12 @@
                   r_, rollaxis, sum
 
 __all__ = ['logsumexp', 'factorial','factorial2','factorialk','comb',
-           'central_diff_weights', 'derivative', 'pade', 'lena',
-           'logsumexp_b']
+           'central_diff_weights', 'derivative', 'pade', 'lena']
 
 # XXX: the factorial functions could move to scipy.special, and the others
 # to numpy perhaps?
 
-def logsumexp(a, axis=None):
+def logsumexp(a, axis=None, b=None):
     """Compute the log of the sum of exponentials of input elements.
 
     Parameters
@@ -24,12 +23,17 @@ def logsumexp(a, axis=None):
     axis : int, optional
         Axis over which the sum is taken. By default `axis` is None,
         and all elements are summed.
+    b : array-like, optional
+        Scaling factor for exp(`a`) must be of the same shape as `a` or
+        broadcastable to `a`.
+
 
     Returns
     -------
     res : ndarray
         The result, ``np.log(np.sum(np.exp(a)))`` calculated in a numerically
-        more stable way.
+        more stable way. If `b` is given then ``np.log(np.sum(b*np.exp(a)))``
+        is returned.
 
     See Also
     --------
@@ -50,59 +54,30 @@ def logsumexp(a, axis=None):
     >>> logsumexp(a)
     9.4586297444267107
 
-    """
-    a = asarray(a)
-    if axis is None:
-        a = a.ravel()
-    else:
-        a = rollaxis(a, axis)
-    a_max = a.max(axis=0)
-    out = log(sum(exp(a - a_max), axis=0))
-    out += a_max
-    return out
-
-def logsumexp_b(a, b, axis=None):
-    """
-    Compute the log of the sum of scaled exponentials of inputs.
-
-    Parameters
-    ----------
-    a : array-like
-        Input array
-    b : array-like
-        Scaling factor for exp(`a`) must be of the same shape as `a` or
-        broadcastable to `a`.
-    axis : int, optional
-        Axis over which the sum is taken. By default `axis` is None, and
-        all elements are summed.
-
-    Returns
-    -------
-    res : ndarray
-        The result, ``np.log(np.sum(b*np.exp(a)))`` calculated in a numerically
-        more stable way.
+    With weights
 
-    Examples
-    --------
-    >>> import numpy as np
-    >>> from scipy.misc import logsumexp_b
     >>> a = np.arange(10)
     >>> b = np.arange(10, 0, -1)
-    >>> logsumexp_b(a, b)
+    >>> logsumexp(a, b=b)
     9.9170178533034665
     >>> np.log(np.sum(b*np.exp(a)))
     9.9170178533034647
     """
     a = asarray(a)
-    b = asarray(b)
     if axis is None:
         a = a.ravel()
-        b = b.ravel()
     else:
         a = rollaxis(a, axis)
-        b = rollaxis(b, axis)
     a_max = a.max(axis=0)
-    out = log(sum(b * exp(a - a_max), axis=0))
+    if b is not None:
+        b = asarray(b)
+        if axis is None:
+            b = b.ravel()
+        else:
+            b = rollaxis(b, axis)
+        out = log(sum(b * exp(a - a_max), axis=0))
+    else:
+        out = log(sum(exp(a - a_max), axis=0))
     out += a_max
     return out
 

scipy/misc/tests/test_common.py

@@ -2,7 +2,7 @@
 from numpy.testing import assert_array_equal, assert_almost_equal, \
                           assert_array_almost_equal
 
-from scipy.misc import pade, logsumexp, logsumexp_b
+from scipy.misc import pade, logsumexp
 
 
 def test_pade_trivial():
@@ -60,12 +60,12 @@ def test_logsumexp_b():
     a = np.arange(200)
     b = np.arange(200, 0, -1)
     desired = np.log(np.sum(b*np.exp(a)))
-    assert_almost_equal(logsumexp_b(a, b), desired)
+    assert_almost_equal(logsumexp(a, b=b), desired)
 
     a = [1000, 1000]
     b = [1.2, 1.2]
     desired = 1000 + np.log(2 * 1.2)
-    assert_almost_equal(logsumexp_b(a, b), desired)
+    assert_almost_equal(logsumexp(a, b=b), desired)
 
     x = np.array([1e-40] * 100000)
     b = np.linspace(1, 1000, 1e5)
@@ -74,9 +74,9 @@ def test_logsumexp_b():
     X = np.vstack((x, x))
     logX = np.vstack((logx, logx))
     B = np.vstack((b, b))
-    assert_array_almost_equal(np.exp(logsumexp_b(logX, B)), (B * X).sum())
-    assert_array_almost_equal(np.exp(logsumexp_b(logX, B, axis=0)),
+    assert_array_almost_equal(np.exp(logsumexp(logX, b=B)), (B * X).sum())
+    assert_array_almost_equal(np.exp(logsumexp(logX, b=B, axis=0)),
                                 (B * X).sum(axis=0))
-    assert_array_almost_equal(np.exp(logsumexp_b(logX, B, axis=1)),
+    assert_array_almost_equal(np.exp(logsumexp(logX, b=B, axis=1)),
                                 (B * X).sum(axis=1))