diff --git a/scipy/misc/__init__.py b/scipy/misc/__init__.py index 18083d2279af..cb2afff50cb1 100644 --- a/scipy/misc/__init__.py +++ b/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 diff --git a/scipy/misc/common.py b/scipy/misc/common.py index 652963f7c28c..cbd78df31894 100644 --- a/scipy/misc/common.py +++ b/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 diff --git a/scipy/misc/tests/test_common.py b/scipy/misc/tests/test_common.py index d29cc5e22774..0ec44b1805d6 100644 --- a/scipy/misc/tests/test_common.py +++ b/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))