test_logsumexp.py 6.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  1. import numpy as np
  2. from numpy.testing import (assert_almost_equal, assert_equal, assert_allclose,
  3. assert_array_almost_equal, assert_)
  4. from scipy.special import logsumexp, softmax
  5. def test_logsumexp():
  6. # Test whether logsumexp() function correctly handles large inputs.
  7. a = np.arange(200)
  8. desired = np.log(np.sum(np.exp(a)))
  9. assert_almost_equal(logsumexp(a), desired)
  10. # Now test with large numbers
  11. b = [1000, 1000]
  12. desired = 1000.0 + np.log(2.0)
  13. assert_almost_equal(logsumexp(b), desired)
  14. n = 1000
  15. b = np.full(n, 10000, dtype='float64')
  16. desired = 10000.0 + np.log(n)
  17. assert_almost_equal(logsumexp(b), desired)
  18. x = np.array([1e-40] * 1000000)
  19. logx = np.log(x)
  20. X = np.vstack([x, x])
  21. logX = np.vstack([logx, logx])
  22. assert_array_almost_equal(np.exp(logsumexp(logX)), X.sum())
  23. assert_array_almost_equal(np.exp(logsumexp(logX, axis=0)), X.sum(axis=0))
  24. assert_array_almost_equal(np.exp(logsumexp(logX, axis=1)), X.sum(axis=1))
  25. # Handling special values properly
  26. assert_equal(logsumexp(np.inf), np.inf)
  27. assert_equal(logsumexp(-np.inf), -np.inf)
  28. assert_equal(logsumexp(np.nan), np.nan)
  29. assert_equal(logsumexp([-np.inf, -np.inf]), -np.inf)
  30. # Handling an array with different magnitudes on the axes
  31. assert_array_almost_equal(logsumexp([[1e10, 1e-10],
  32. [-1e10, -np.inf]], axis=-1),
  33. [1e10, -1e10])
  34. # Test keeping dimensions
  35. assert_array_almost_equal(logsumexp([[1e10, 1e-10],
  36. [-1e10, -np.inf]],
  37. axis=-1,
  38. keepdims=True),
  39. [[1e10], [-1e10]])
  40. # Test multiple axes
  41. assert_array_almost_equal(logsumexp([[1e10, 1e-10],
  42. [-1e10, -np.inf]],
  43. axis=(-1,-2)),
  44. 1e10)
  45. def test_logsumexp_b():
  46. a = np.arange(200)
  47. b = np.arange(200, 0, -1)
  48. desired = np.log(np.sum(b*np.exp(a)))
  49. assert_almost_equal(logsumexp(a, b=b), desired)
  50. a = [1000, 1000]
  51. b = [1.2, 1.2]
  52. desired = 1000 + np.log(2 * 1.2)
  53. assert_almost_equal(logsumexp(a, b=b), desired)
  54. x = np.array([1e-40] * 100000)
  55. b = np.linspace(1, 1000, 100000)
  56. logx = np.log(x)
  57. X = np.vstack((x, x))
  58. logX = np.vstack((logx, logx))
  59. B = np.vstack((b, b))
  60. assert_array_almost_equal(np.exp(logsumexp(logX, b=B)), (B * X).sum())
  61. assert_array_almost_equal(np.exp(logsumexp(logX, b=B, axis=0)),
  62. (B * X).sum(axis=0))
  63. assert_array_almost_equal(np.exp(logsumexp(logX, b=B, axis=1)),
  64. (B * X).sum(axis=1))
  65. def test_logsumexp_sign():
  66. a = [1,1,1]
  67. b = [1,-1,-1]
  68. r, s = logsumexp(a, b=b, return_sign=True)
  69. assert_almost_equal(r,1)
  70. assert_equal(s,-1)
  71. def test_logsumexp_sign_zero():
  72. a = [1,1]
  73. b = [1,-1]
  74. r, s = logsumexp(a, b=b, return_sign=True)
  75. assert_(not np.isfinite(r))
  76. assert_(not np.isnan(r))
  77. assert_(r < 0)
  78. assert_equal(s,0)
  79. def test_logsumexp_sign_shape():
  80. a = np.ones((1,2,3,4))
  81. b = np.ones_like(a)
  82. r, s = logsumexp(a, axis=2, b=b, return_sign=True)
  83. assert_equal(r.shape, s.shape)
  84. assert_equal(r.shape, (1,2,4))
  85. r, s = logsumexp(a, axis=(1,3), b=b, return_sign=True)
  86. assert_equal(r.shape, s.shape)
  87. assert_equal(r.shape, (1,3))
  88. def test_logsumexp_shape():
  89. a = np.ones((1, 2, 3, 4))
  90. b = np.ones_like(a)
  91. r = logsumexp(a, axis=2, b=b)
  92. assert_equal(r.shape, (1, 2, 4))
  93. r = logsumexp(a, axis=(1, 3), b=b)
  94. assert_equal(r.shape, (1, 3))
  95. def test_logsumexp_b_zero():
  96. a = [1,10000]
  97. b = [1,0]
  98. assert_almost_equal(logsumexp(a, b=b), 1)
  99. def test_logsumexp_b_shape():
  100. a = np.zeros((4,1,2,1))
  101. b = np.ones((3,1,5))
  102. logsumexp(a, b=b)
  103. def test_softmax_fixtures():
  104. assert_allclose(softmax([1000, 0, 0, 0]), np.array([1, 0, 0, 0]),
  105. rtol=1e-13)
  106. assert_allclose(softmax([1, 1]), np.array([.5, .5]), rtol=1e-13)
  107. assert_allclose(softmax([0, 1]), np.array([1, np.e])/(1 + np.e),
  108. rtol=1e-13)
  109. # Expected value computed using mpmath (with mpmath.mp.dps = 200) and then
  110. # converted to float.
  111. x = np.arange(4)
  112. expected = np.array([0.03205860328008499,
  113. 0.08714431874203256,
  114. 0.23688281808991013,
  115. 0.6439142598879722])
  116. assert_allclose(softmax(x), expected, rtol=1e-13)
  117. # Translation property. If all the values are changed by the same amount,
  118. # the softmax result does not change.
  119. assert_allclose(softmax(x + 100), expected, rtol=1e-13)
  120. # When axis=None, softmax operates on the entire array, and preserves
  121. # the shape.
  122. assert_allclose(softmax(x.reshape(2, 2)), expected.reshape(2, 2),
  123. rtol=1e-13)
  124. def test_softmax_multi_axes():
  125. assert_allclose(softmax([[1000, 0], [1000, 0]], axis=0),
  126. np.array([[.5, .5], [.5, .5]]), rtol=1e-13)
  127. assert_allclose(softmax([[1000, 0], [1000, 0]], axis=1),
  128. np.array([[1, 0], [1, 0]]), rtol=1e-13)
  129. # Expected value computed using mpmath (with mpmath.mp.dps = 200) and then
  130. # converted to float.
  131. x = np.array([[-25, 0, 25, 50],
  132. [1, 325, 749, 750]])
  133. expected = np.array([[2.678636961770877e-33,
  134. 1.9287498479371314e-22,
  135. 1.3887943864771144e-11,
  136. 0.999999999986112],
  137. [0.0,
  138. 1.9444526359919372e-185,
  139. 0.2689414213699951,
  140. 0.7310585786300048]])
  141. assert_allclose(softmax(x, axis=1), expected, rtol=1e-13)
  142. assert_allclose(softmax(x.T, axis=0), expected.T, rtol=1e-13)
  143. # 3-d input, with a tuple for the axis.
  144. x3d = x.reshape(2, 2, 2)
  145. assert_allclose(softmax(x3d, axis=(1, 2)), expected.reshape(2, 2, 2),
  146. rtol=1e-13)