test_orthogonal_eval.py 9.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268
  1. import numpy as np
  2. from numpy.testing import assert_, assert_allclose
  3. import pytest
  4. from scipy.special import _ufuncs
  5. import scipy.special._orthogonal as orth
  6. from scipy.special._testutils import FuncData
  7. def test_eval_chebyt():
  8. n = np.arange(0, 10000, 7)
  9. x = 2*np.random.rand() - 1
  10. v1 = np.cos(n*np.arccos(x))
  11. v2 = _ufuncs.eval_chebyt(n, x)
  12. assert_(np.allclose(v1, v2, rtol=1e-15))
  13. def test_eval_genlaguerre_restriction():
  14. # check it returns nan for alpha <= -1
  15. assert_(np.isnan(_ufuncs.eval_genlaguerre(0, -1, 0)))
  16. assert_(np.isnan(_ufuncs.eval_genlaguerre(0.1, -1, 0)))
  17. def test_warnings():
  18. # ticket 1334
  19. with np.errstate(all='raise'):
  20. # these should raise no fp warnings
  21. _ufuncs.eval_legendre(1, 0)
  22. _ufuncs.eval_laguerre(1, 1)
  23. _ufuncs.eval_gegenbauer(1, 1, 0)
  24. class TestPolys:
  25. """
  26. Check that the eval_* functions agree with the constructed polynomials
  27. """
  28. def check_poly(self, func, cls, param_ranges=[], x_range=[], nn=10,
  29. nparam=10, nx=10, rtol=1e-8):
  30. np.random.seed(1234)
  31. dataset = []
  32. for n in np.arange(nn):
  33. params = [a + (b-a)*np.random.rand(nparam) for a,b in param_ranges]
  34. params = np.asarray(params).T
  35. if not param_ranges:
  36. params = [0]
  37. for p in params:
  38. if param_ranges:
  39. p = (n,) + tuple(p)
  40. else:
  41. p = (n,)
  42. x = x_range[0] + (x_range[1] - x_range[0])*np.random.rand(nx)
  43. x[0] = x_range[0] # always include domain start point
  44. x[1] = x_range[1] # always include domain end point
  45. poly = np.poly1d(cls(*p).coef)
  46. z = np.c_[np.tile(p, (nx,1)), x, poly(x)]
  47. dataset.append(z)
  48. dataset = np.concatenate(dataset, axis=0)
  49. def polyfunc(*p):
  50. p = (p[0].astype(int),) + p[1:]
  51. return func(*p)
  52. with np.errstate(all='raise'):
  53. ds = FuncData(polyfunc, dataset, list(range(len(param_ranges)+2)), -1,
  54. rtol=rtol)
  55. ds.check()
  56. def test_jacobi(self):
  57. self.check_poly(_ufuncs.eval_jacobi, orth.jacobi,
  58. param_ranges=[(-0.99, 10), (-0.99, 10)],
  59. x_range=[-1, 1], rtol=1e-5)
  60. def test_sh_jacobi(self):
  61. self.check_poly(_ufuncs.eval_sh_jacobi, orth.sh_jacobi,
  62. param_ranges=[(1, 10), (0, 1)], x_range=[0, 1],
  63. rtol=1e-5)
  64. def test_gegenbauer(self):
  65. self.check_poly(_ufuncs.eval_gegenbauer, orth.gegenbauer,
  66. param_ranges=[(-0.499, 10)], x_range=[-1, 1],
  67. rtol=1e-7)
  68. def test_chebyt(self):
  69. self.check_poly(_ufuncs.eval_chebyt, orth.chebyt,
  70. param_ranges=[], x_range=[-1, 1])
  71. def test_chebyu(self):
  72. self.check_poly(_ufuncs.eval_chebyu, orth.chebyu,
  73. param_ranges=[], x_range=[-1, 1])
  74. def test_chebys(self):
  75. self.check_poly(_ufuncs.eval_chebys, orth.chebys,
  76. param_ranges=[], x_range=[-2, 2])
  77. def test_chebyc(self):
  78. self.check_poly(_ufuncs.eval_chebyc, orth.chebyc,
  79. param_ranges=[], x_range=[-2, 2])
  80. def test_sh_chebyt(self):
  81. with np.errstate(all='ignore'):
  82. self.check_poly(_ufuncs.eval_sh_chebyt, orth.sh_chebyt,
  83. param_ranges=[], x_range=[0, 1])
  84. def test_sh_chebyu(self):
  85. self.check_poly(_ufuncs.eval_sh_chebyu, orth.sh_chebyu,
  86. param_ranges=[], x_range=[0, 1])
  87. def test_legendre(self):
  88. self.check_poly(_ufuncs.eval_legendre, orth.legendre,
  89. param_ranges=[], x_range=[-1, 1])
  90. def test_sh_legendre(self):
  91. with np.errstate(all='ignore'):
  92. self.check_poly(_ufuncs.eval_sh_legendre, orth.sh_legendre,
  93. param_ranges=[], x_range=[0, 1])
  94. def test_genlaguerre(self):
  95. self.check_poly(_ufuncs.eval_genlaguerre, orth.genlaguerre,
  96. param_ranges=[(-0.99, 10)], x_range=[0, 100])
  97. def test_laguerre(self):
  98. self.check_poly(_ufuncs.eval_laguerre, orth.laguerre,
  99. param_ranges=[], x_range=[0, 100])
  100. def test_hermite(self):
  101. self.check_poly(_ufuncs.eval_hermite, orth.hermite,
  102. param_ranges=[], x_range=[-100, 100])
  103. def test_hermitenorm(self):
  104. self.check_poly(_ufuncs.eval_hermitenorm, orth.hermitenorm,
  105. param_ranges=[], x_range=[-100, 100])
  106. class TestRecurrence:
  107. """
  108. Check that the eval_* functions sig='ld->d' and 'dd->d' agree.
  109. """
  110. def check_poly(self, func, param_ranges=[], x_range=[], nn=10,
  111. nparam=10, nx=10, rtol=1e-8):
  112. np.random.seed(1234)
  113. dataset = []
  114. for n in np.arange(nn):
  115. params = [a + (b-a)*np.random.rand(nparam) for a,b in param_ranges]
  116. params = np.asarray(params).T
  117. if not param_ranges:
  118. params = [0]
  119. for p in params:
  120. if param_ranges:
  121. p = (n,) + tuple(p)
  122. else:
  123. p = (n,)
  124. x = x_range[0] + (x_range[1] - x_range[0])*np.random.rand(nx)
  125. x[0] = x_range[0] # always include domain start point
  126. x[1] = x_range[1] # always include domain end point
  127. kw = dict(sig=(len(p)+1)*'d'+'->d')
  128. z = np.c_[np.tile(p, (nx,1)), x, func(*(p + (x,)), **kw)]
  129. dataset.append(z)
  130. dataset = np.concatenate(dataset, axis=0)
  131. def polyfunc(*p):
  132. p = (p[0].astype(int),) + p[1:]
  133. kw = dict(sig='l'+(len(p)-1)*'d'+'->d')
  134. return func(*p, **kw)
  135. with np.errstate(all='raise'):
  136. ds = FuncData(polyfunc, dataset, list(range(len(param_ranges)+2)), -1,
  137. rtol=rtol)
  138. ds.check()
  139. def test_jacobi(self):
  140. self.check_poly(_ufuncs.eval_jacobi,
  141. param_ranges=[(-0.99, 10), (-0.99, 10)],
  142. x_range=[-1, 1])
  143. def test_sh_jacobi(self):
  144. self.check_poly(_ufuncs.eval_sh_jacobi,
  145. param_ranges=[(1, 10), (0, 1)], x_range=[0, 1])
  146. def test_gegenbauer(self):
  147. self.check_poly(_ufuncs.eval_gegenbauer,
  148. param_ranges=[(-0.499, 10)], x_range=[-1, 1])
  149. def test_chebyt(self):
  150. self.check_poly(_ufuncs.eval_chebyt,
  151. param_ranges=[], x_range=[-1, 1])
  152. def test_chebyu(self):
  153. self.check_poly(_ufuncs.eval_chebyu,
  154. param_ranges=[], x_range=[-1, 1])
  155. def test_chebys(self):
  156. self.check_poly(_ufuncs.eval_chebys,
  157. param_ranges=[], x_range=[-2, 2])
  158. def test_chebyc(self):
  159. self.check_poly(_ufuncs.eval_chebyc,
  160. param_ranges=[], x_range=[-2, 2])
  161. def test_sh_chebyt(self):
  162. self.check_poly(_ufuncs.eval_sh_chebyt,
  163. param_ranges=[], x_range=[0, 1])
  164. def test_sh_chebyu(self):
  165. self.check_poly(_ufuncs.eval_sh_chebyu,
  166. param_ranges=[], x_range=[0, 1])
  167. def test_legendre(self):
  168. self.check_poly(_ufuncs.eval_legendre,
  169. param_ranges=[], x_range=[-1, 1])
  170. def test_sh_legendre(self):
  171. self.check_poly(_ufuncs.eval_sh_legendre,
  172. param_ranges=[], x_range=[0, 1])
  173. def test_genlaguerre(self):
  174. self.check_poly(_ufuncs.eval_genlaguerre,
  175. param_ranges=[(-0.99, 10)], x_range=[0, 100])
  176. def test_laguerre(self):
  177. self.check_poly(_ufuncs.eval_laguerre,
  178. param_ranges=[], x_range=[0, 100])
  179. def test_hermite(self):
  180. v = _ufuncs.eval_hermite(70, 1.0)
  181. a = -1.457076485701412e60
  182. assert_allclose(v, a)
  183. def test_hermite_domain():
  184. # Regression test for gh-11091.
  185. assert np.isnan(_ufuncs.eval_hermite(-1, 1.0))
  186. assert np.isnan(_ufuncs.eval_hermitenorm(-1, 1.0))
  187. @pytest.mark.parametrize("n", [0, 1, 2])
  188. @pytest.mark.parametrize("x", [0, 1, np.nan])
  189. def test_hermite_nan(n, x):
  190. # Regression test for gh-11369.
  191. assert np.isnan(_ufuncs.eval_hermite(n, x)) == np.any(np.isnan([n, x]))
  192. assert np.isnan(_ufuncs.eval_hermitenorm(n, x)) == np.any(np.isnan([n, x]))
  193. @pytest.mark.parametrize('n', [0, 1, 2, 3.2])
  194. @pytest.mark.parametrize('alpha', [1, np.nan])
  195. @pytest.mark.parametrize('x', [2, np.nan])
  196. def test_genlaguerre_nan(n, alpha, x):
  197. # Regression test for gh-11361.
  198. nan_laguerre = np.isnan(_ufuncs.eval_genlaguerre(n, alpha, x))
  199. nan_arg = np.any(np.isnan([n, alpha, x]))
  200. assert nan_laguerre == nan_arg
  201. @pytest.mark.parametrize('n', [0, 1, 2, 3.2])
  202. @pytest.mark.parametrize('alpha', [0.0, 1, np.nan])
  203. @pytest.mark.parametrize('x', [1e-6, 2, np.nan])
  204. def test_gegenbauer_nan(n, alpha, x):
  205. # Regression test for gh-11370.
  206. nan_gegenbauer = np.isnan(_ufuncs.eval_gegenbauer(n, alpha, x))
  207. nan_arg = np.any(np.isnan([n, alpha, x]))
  208. assert nan_gegenbauer == nan_arg