test_numpy.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351
  1. from sympy.concrete.summations import Sum
  2. from sympy.core.mod import Mod
  3. from sympy.core.relational import (Equality, Unequality)
  4. from sympy.core.symbol import Symbol
  5. from sympy.functions.elementary.miscellaneous import sqrt
  6. from sympy.functions.elementary.piecewise import Piecewise
  7. from sympy.functions.special.gamma_functions import polygamma
  8. from sympy.functions.special.error_functions import (Si, Ci)
  9. from sympy.matrices.expressions.blockmatrix import BlockMatrix
  10. from sympy.matrices.expressions.matexpr import MatrixSymbol
  11. from sympy.matrices.expressions.special import Identity
  12. from sympy.utilities.lambdify import lambdify
  13. from sympy.abc import x, i, j, a, b, c, d
  14. from sympy.core import Pow
  15. from sympy.codegen.matrix_nodes import MatrixSolve
  16. from sympy.codegen.numpy_nodes import logaddexp, logaddexp2
  17. from sympy.codegen.cfunctions import log1p, expm1, hypot, log10, exp2, log2, Sqrt
  18. from sympy.tensor.array import Array
  19. from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct, ArrayAdd, \
  20. PermuteDims, ArrayDiagonal
  21. from sympy.printing.numpy import NumPyPrinter, SciPyPrinter, _numpy_known_constants, \
  22. _numpy_known_functions, _scipy_known_constants, _scipy_known_functions
  23. from sympy.tensor.array.expressions.from_matrix_to_array import convert_matrix_to_array
  24. from sympy.testing.pytest import skip, raises
  25. from sympy.external import import_module
  26. np = import_module('numpy')
  27. if np:
  28. deafult_float_info = np.finfo(np.array([]).dtype)
  29. NUMPY_DEFAULT_EPSILON = deafult_float_info.eps
  30. def test_numpy_piecewise_regression():
  31. """
  32. NumPyPrinter needs to print Piecewise()'s choicelist as a list to avoid
  33. breaking compatibility with numpy 1.8. This is not necessary in numpy 1.9+.
  34. See gh-9747 and gh-9749 for details.
  35. """
  36. printer = NumPyPrinter()
  37. p = Piecewise((1, x < 0), (0, True))
  38. assert printer.doprint(p) == \
  39. 'numpy.select([numpy.less(x, 0),True], [1,0], default=numpy.nan)'
  40. assert printer.module_imports == {'numpy': {'select', 'less', 'nan'}}
  41. def test_numpy_logaddexp():
  42. lae = logaddexp(a, b)
  43. assert NumPyPrinter().doprint(lae) == 'numpy.logaddexp(a, b)'
  44. lae2 = logaddexp2(a, b)
  45. assert NumPyPrinter().doprint(lae2) == 'numpy.logaddexp2(a, b)'
  46. def test_sum():
  47. if not np:
  48. skip("NumPy not installed")
  49. s = Sum(x ** i, (i, a, b))
  50. f = lambdify((a, b, x), s, 'numpy')
  51. a_, b_ = 0, 10
  52. x_ = np.linspace(-1, +1, 10)
  53. assert np.allclose(f(a_, b_, x_), sum(x_ ** i_ for i_ in range(a_, b_ + 1)))
  54. s = Sum(i * x, (i, a, b))
  55. f = lambdify((a, b, x), s, 'numpy')
  56. a_, b_ = 0, 10
  57. x_ = np.linspace(-1, +1, 10)
  58. assert np.allclose(f(a_, b_, x_), sum(i_ * x_ for i_ in range(a_, b_ + 1)))
  59. def test_multiple_sums():
  60. if not np:
  61. skip("NumPy not installed")
  62. s = Sum((x + j) * i, (i, a, b), (j, c, d))
  63. f = lambdify((a, b, c, d, x), s, 'numpy')
  64. a_, b_ = 0, 10
  65. c_, d_ = 11, 21
  66. x_ = np.linspace(-1, +1, 10)
  67. assert np.allclose(f(a_, b_, c_, d_, x_),
  68. sum((x_ + j_) * i_ for i_ in range(a_, b_ + 1) for j_ in range(c_, d_ + 1)))
  69. def test_codegen_einsum():
  70. if not np:
  71. skip("NumPy not installed")
  72. M = MatrixSymbol("M", 2, 2)
  73. N = MatrixSymbol("N", 2, 2)
  74. cg = convert_matrix_to_array(M * N)
  75. f = lambdify((M, N), cg, 'numpy')
  76. ma = np.array([[1, 2], [3, 4]])
  77. mb = np.array([[1,-2], [-1, 3]])
  78. assert (f(ma, mb) == np.matmul(ma, mb)).all()
  79. def test_codegen_extra():
  80. if not np:
  81. skip("NumPy not installed")
  82. M = MatrixSymbol("M", 2, 2)
  83. N = MatrixSymbol("N", 2, 2)
  84. P = MatrixSymbol("P", 2, 2)
  85. Q = MatrixSymbol("Q", 2, 2)
  86. ma = np.array([[1, 2], [3, 4]])
  87. mb = np.array([[1,-2], [-1, 3]])
  88. mc = np.array([[2, 0], [1, 2]])
  89. md = np.array([[1,-1], [4, 7]])
  90. cg = ArrayTensorProduct(M, N)
  91. f = lambdify((M, N), cg, 'numpy')
  92. assert (f(ma, mb) == np.einsum(ma, [0, 1], mb, [2, 3])).all()
  93. cg = ArrayAdd(M, N)
  94. f = lambdify((M, N), cg, 'numpy')
  95. assert (f(ma, mb) == ma+mb).all()
  96. cg = ArrayAdd(M, N, P)
  97. f = lambdify((M, N, P), cg, 'numpy')
  98. assert (f(ma, mb, mc) == ma+mb+mc).all()
  99. cg = ArrayAdd(M, N, P, Q)
  100. f = lambdify((M, N, P, Q), cg, 'numpy')
  101. assert (f(ma, mb, mc, md) == ma+mb+mc+md).all()
  102. cg = PermuteDims(M, [1, 0])
  103. f = lambdify((M,), cg, 'numpy')
  104. assert (f(ma) == ma.T).all()
  105. cg = PermuteDims(ArrayTensorProduct(M, N), [1, 2, 3, 0])
  106. f = lambdify((M, N), cg, 'numpy')
  107. assert (f(ma, mb) == np.transpose(np.einsum(ma, [0, 1], mb, [2, 3]), (1, 2, 3, 0))).all()
  108. cg = ArrayDiagonal(ArrayTensorProduct(M, N), (1, 2))
  109. f = lambdify((M, N), cg, 'numpy')
  110. assert (f(ma, mb) == np.diagonal(np.einsum(ma, [0, 1], mb, [2, 3]), axis1=1, axis2=2)).all()
  111. def test_relational():
  112. if not np:
  113. skip("NumPy not installed")
  114. e = Equality(x, 1)
  115. f = lambdify((x,), e)
  116. x_ = np.array([0, 1, 2])
  117. assert np.array_equal(f(x_), [False, True, False])
  118. e = Unequality(x, 1)
  119. f = lambdify((x,), e)
  120. x_ = np.array([0, 1, 2])
  121. assert np.array_equal(f(x_), [True, False, True])
  122. e = (x < 1)
  123. f = lambdify((x,), e)
  124. x_ = np.array([0, 1, 2])
  125. assert np.array_equal(f(x_), [True, False, False])
  126. e = (x <= 1)
  127. f = lambdify((x,), e)
  128. x_ = np.array([0, 1, 2])
  129. assert np.array_equal(f(x_), [True, True, False])
  130. e = (x > 1)
  131. f = lambdify((x,), e)
  132. x_ = np.array([0, 1, 2])
  133. assert np.array_equal(f(x_), [False, False, True])
  134. e = (x >= 1)
  135. f = lambdify((x,), e)
  136. x_ = np.array([0, 1, 2])
  137. assert np.array_equal(f(x_), [False, True, True])
  138. def test_mod():
  139. if not np:
  140. skip("NumPy not installed")
  141. e = Mod(a, b)
  142. f = lambdify((a, b), e)
  143. a_ = np.array([0, 1, 2, 3])
  144. b_ = 2
  145. assert np.array_equal(f(a_, b_), [0, 1, 0, 1])
  146. a_ = np.array([0, 1, 2, 3])
  147. b_ = np.array([2, 2, 2, 2])
  148. assert np.array_equal(f(a_, b_), [0, 1, 0, 1])
  149. a_ = np.array([2, 3, 4, 5])
  150. b_ = np.array([2, 3, 4, 5])
  151. assert np.array_equal(f(a_, b_), [0, 0, 0, 0])
  152. def test_pow():
  153. if not np:
  154. skip('NumPy not installed')
  155. expr = Pow(2, -1, evaluate=False)
  156. f = lambdify([], expr, 'numpy')
  157. assert f() == 0.5
  158. def test_expm1():
  159. if not np:
  160. skip("NumPy not installed")
  161. f = lambdify((a,), expm1(a), 'numpy')
  162. assert abs(f(1e-10) - 1e-10 - 5e-21) <= 1e-10 * NUMPY_DEFAULT_EPSILON
  163. def test_log1p():
  164. if not np:
  165. skip("NumPy not installed")
  166. f = lambdify((a,), log1p(a), 'numpy')
  167. assert abs(f(1e-99) - 1e-99) <= 1e-99 * NUMPY_DEFAULT_EPSILON
  168. def test_hypot():
  169. if not np:
  170. skip("NumPy not installed")
  171. assert abs(lambdify((a, b), hypot(a, b), 'numpy')(3, 4) - 5) <= NUMPY_DEFAULT_EPSILON
  172. def test_log10():
  173. if not np:
  174. skip("NumPy not installed")
  175. assert abs(lambdify((a,), log10(a), 'numpy')(100) - 2) <= NUMPY_DEFAULT_EPSILON
  176. def test_exp2():
  177. if not np:
  178. skip("NumPy not installed")
  179. assert abs(lambdify((a,), exp2(a), 'numpy')(5) - 32) <= NUMPY_DEFAULT_EPSILON
  180. def test_log2():
  181. if not np:
  182. skip("NumPy not installed")
  183. assert abs(lambdify((a,), log2(a), 'numpy')(256) - 8) <= NUMPY_DEFAULT_EPSILON
  184. def test_Sqrt():
  185. if not np:
  186. skip("NumPy not installed")
  187. assert abs(lambdify((a,), Sqrt(a), 'numpy')(4) - 2) <= NUMPY_DEFAULT_EPSILON
  188. def test_sqrt():
  189. if not np:
  190. skip("NumPy not installed")
  191. assert abs(lambdify((a,), sqrt(a), 'numpy')(4) - 2) <= NUMPY_DEFAULT_EPSILON
  192. def test_matsolve():
  193. if not np:
  194. skip("NumPy not installed")
  195. M = MatrixSymbol("M", 3, 3)
  196. x = MatrixSymbol("x", 3, 1)
  197. expr = M**(-1) * x + x
  198. matsolve_expr = MatrixSolve(M, x) + x
  199. f = lambdify((M, x), expr)
  200. f_matsolve = lambdify((M, x), matsolve_expr)
  201. m0 = np.array([[1, 2, 3], [3, 2, 5], [5, 6, 7]])
  202. assert np.linalg.matrix_rank(m0) == 3
  203. x0 = np.array([3, 4, 5])
  204. assert np.allclose(f_matsolve(m0, x0), f(m0, x0))
  205. def test_16857():
  206. if not np:
  207. skip("NumPy not installed")
  208. a_1 = MatrixSymbol('a_1', 10, 3)
  209. a_2 = MatrixSymbol('a_2', 10, 3)
  210. a_3 = MatrixSymbol('a_3', 10, 3)
  211. a_4 = MatrixSymbol('a_4', 10, 3)
  212. A = BlockMatrix([[a_1, a_2], [a_3, a_4]])
  213. assert A.shape == (20, 6)
  214. printer = NumPyPrinter()
  215. assert printer.doprint(A) == 'numpy.block([[a_1, a_2], [a_3, a_4]])'
  216. def test_issue_17006():
  217. if not np:
  218. skip("NumPy not installed")
  219. M = MatrixSymbol("M", 2, 2)
  220. f = lambdify(M, M + Identity(2))
  221. ma = np.array([[1, 2], [3, 4]])
  222. mr = np.array([[2, 2], [3, 5]])
  223. assert (f(ma) == mr).all()
  224. from sympy.core.symbol import symbols
  225. n = symbols('n', integer=True)
  226. N = MatrixSymbol("M", n, n)
  227. raises(NotImplementedError, lambda: lambdify(N, N + Identity(n)))
  228. def test_numpy_array():
  229. assert NumPyPrinter().doprint(Array(((1, 2), (3, 5)))) == 'numpy.array([[1, 2], [3, 5]])'
  230. assert NumPyPrinter().doprint(Array((1, 2))) == 'numpy.array((1, 2))'
  231. def test_numpy_known_funcs_consts():
  232. assert _numpy_known_constants['NaN'] == 'numpy.nan'
  233. assert _numpy_known_constants['EulerGamma'] == 'numpy.euler_gamma'
  234. assert _numpy_known_functions['acos'] == 'numpy.arccos'
  235. assert _numpy_known_functions['log'] == 'numpy.log'
  236. def test_scipy_known_funcs_consts():
  237. assert _scipy_known_constants['GoldenRatio'] == 'scipy.constants.golden_ratio'
  238. assert _scipy_known_constants['Pi'] == 'scipy.constants.pi'
  239. assert _scipy_known_functions['erf'] == 'scipy.special.erf'
  240. assert _scipy_known_functions['factorial'] == 'scipy.special.factorial'
  241. def test_numpy_print_methods():
  242. prntr = NumPyPrinter()
  243. assert hasattr(prntr, '_print_acos')
  244. assert hasattr(prntr, '_print_log')
  245. def test_scipy_print_methods():
  246. prntr = SciPyPrinter()
  247. assert hasattr(prntr, '_print_acos')
  248. assert hasattr(prntr, '_print_log')
  249. assert hasattr(prntr, '_print_erf')
  250. assert hasattr(prntr, '_print_factorial')
  251. assert hasattr(prntr, '_print_chebyshevt')
  252. k = Symbol('k', integer=True, nonnegative=True)
  253. x = Symbol('x', real=True)
  254. assert prntr.doprint(polygamma(k, x)) == "scipy.special.polygamma(k, x)"
  255. assert prntr.doprint(Si(x)) == "scipy.special.sici(x)[0]"
  256. assert prntr.doprint(Ci(x)) == "scipy.special.sici(x)[1]"