test_numpy.py 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364
  1. import queue
  2. import threading
  3. import multiprocessing
  4. import numpy as np
  5. import pytest
  6. from numpy.random import random
  7. from numpy.testing import (
  8. assert_array_almost_equal, assert_array_equal, assert_allclose
  9. )
  10. from pytest import raises as assert_raises
  11. import scipy.fft as fft
  12. def fft1(x):
  13. L = len(x)
  14. phase = -2j*np.pi*(np.arange(L)/float(L))
  15. phase = np.arange(L).reshape(-1, 1) * phase
  16. return np.sum(x*np.exp(phase), axis=1)
  17. class TestFFTShift:
  18. def test_fft_n(self):
  19. assert_raises(ValueError, fft.fft, [1, 2, 3], 0)
  20. class TestFFT1D:
  21. def test_identity(self):
  22. maxlen = 512
  23. x = random(maxlen) + 1j*random(maxlen)
  24. xr = random(maxlen)
  25. for i in range(1,maxlen):
  26. assert_array_almost_equal(fft.ifft(fft.fft(x[0:i])), x[0:i],
  27. decimal=12)
  28. assert_array_almost_equal(fft.irfft(fft.rfft(xr[0:i]),i),
  29. xr[0:i], decimal=12)
  30. def test_fft(self):
  31. x = random(30) + 1j*random(30)
  32. expect = fft1(x)
  33. assert_array_almost_equal(expect, fft.fft(x))
  34. assert_array_almost_equal(expect, fft.fft(x, norm="backward"))
  35. assert_array_almost_equal(expect / np.sqrt(30),
  36. fft.fft(x, norm="ortho"))
  37. assert_array_almost_equal(expect / 30, fft.fft(x, norm="forward"))
  38. def test_ifft(self):
  39. x = random(30) + 1j*random(30)
  40. assert_array_almost_equal(x, fft.ifft(fft.fft(x)))
  41. for norm in ["backward", "ortho", "forward"]:
  42. assert_array_almost_equal(
  43. x, fft.ifft(fft.fft(x, norm=norm), norm=norm))
  44. def test_fft2(self):
  45. x = random((30, 20)) + 1j*random((30, 20))
  46. expect = fft.fft(fft.fft(x, axis=1), axis=0)
  47. assert_array_almost_equal(expect, fft.fft2(x))
  48. assert_array_almost_equal(expect, fft.fft2(x, norm="backward"))
  49. assert_array_almost_equal(expect / np.sqrt(30 * 20),
  50. fft.fft2(x, norm="ortho"))
  51. assert_array_almost_equal(expect / (30 * 20),
  52. fft.fft2(x, norm="forward"))
  53. def test_ifft2(self):
  54. x = random((30, 20)) + 1j*random((30, 20))
  55. expect = fft.ifft(fft.ifft(x, axis=1), axis=0)
  56. assert_array_almost_equal(expect, fft.ifft2(x))
  57. assert_array_almost_equal(expect, fft.ifft2(x, norm="backward"))
  58. assert_array_almost_equal(expect * np.sqrt(30 * 20),
  59. fft.ifft2(x, norm="ortho"))
  60. assert_array_almost_equal(expect * (30 * 20),
  61. fft.ifft2(x, norm="forward"))
  62. def test_fftn(self):
  63. x = random((30, 20, 10)) + 1j*random((30, 20, 10))
  64. expect = fft.fft(fft.fft(fft.fft(x, axis=2), axis=1), axis=0)
  65. assert_array_almost_equal(expect, fft.fftn(x))
  66. assert_array_almost_equal(expect, fft.fftn(x, norm="backward"))
  67. assert_array_almost_equal(expect / np.sqrt(30 * 20 * 10),
  68. fft.fftn(x, norm="ortho"))
  69. assert_array_almost_equal(expect / (30 * 20 * 10),
  70. fft.fftn(x, norm="forward"))
  71. def test_ifftn(self):
  72. x = random((30, 20, 10)) + 1j*random((30, 20, 10))
  73. expect = fft.ifft(fft.ifft(fft.ifft(x, axis=2), axis=1), axis=0)
  74. assert_array_almost_equal(expect, fft.ifftn(x))
  75. assert_array_almost_equal(expect, fft.ifftn(x, norm="backward"))
  76. assert_array_almost_equal(fft.ifftn(x) * np.sqrt(30 * 20 * 10),
  77. fft.ifftn(x, norm="ortho"))
  78. assert_array_almost_equal(expect * (30 * 20 * 10),
  79. fft.ifftn(x, norm="forward"))
  80. def test_rfft(self):
  81. x = random(29)
  82. for n in [x.size, 2*x.size]:
  83. for norm in [None, "backward", "ortho", "forward"]:
  84. assert_array_almost_equal(
  85. fft.fft(x, n=n, norm=norm)[:(n//2 + 1)],
  86. fft.rfft(x, n=n, norm=norm))
  87. assert_array_almost_equal(fft.rfft(x, n=n) / np.sqrt(n),
  88. fft.rfft(x, n=n, norm="ortho"))
  89. def test_irfft(self):
  90. x = random(30)
  91. assert_array_almost_equal(x, fft.irfft(fft.rfft(x)))
  92. for norm in ["backward", "ortho", "forward"]:
  93. assert_array_almost_equal(
  94. x, fft.irfft(fft.rfft(x, norm=norm), norm=norm))
  95. def test_rfft2(self):
  96. x = random((30, 20))
  97. expect = fft.fft2(x)[:, :11]
  98. assert_array_almost_equal(expect, fft.rfft2(x))
  99. assert_array_almost_equal(expect, fft.rfft2(x, norm="backward"))
  100. assert_array_almost_equal(expect / np.sqrt(30 * 20),
  101. fft.rfft2(x, norm="ortho"))
  102. assert_array_almost_equal(expect / (30 * 20),
  103. fft.rfft2(x, norm="forward"))
  104. def test_irfft2(self):
  105. x = random((30, 20))
  106. assert_array_almost_equal(x, fft.irfft2(fft.rfft2(x)))
  107. for norm in ["backward", "ortho", "forward"]:
  108. assert_array_almost_equal(
  109. x, fft.irfft2(fft.rfft2(x, norm=norm), norm=norm))
  110. def test_rfftn(self):
  111. x = random((30, 20, 10))
  112. expect = fft.fftn(x)[:, :, :6]
  113. assert_array_almost_equal(expect, fft.rfftn(x))
  114. assert_array_almost_equal(expect, fft.rfftn(x, norm="backward"))
  115. assert_array_almost_equal(expect / np.sqrt(30 * 20 * 10),
  116. fft.rfftn(x, norm="ortho"))
  117. assert_array_almost_equal(expect / (30 * 20 * 10),
  118. fft.rfftn(x, norm="forward"))
  119. def test_irfftn(self):
  120. x = random((30, 20, 10))
  121. assert_array_almost_equal(x, fft.irfftn(fft.rfftn(x)))
  122. for norm in ["backward", "ortho", "forward"]:
  123. assert_array_almost_equal(
  124. x, fft.irfftn(fft.rfftn(x, norm=norm), norm=norm))
  125. def test_hfft(self):
  126. x = random(14) + 1j*random(14)
  127. x_herm = np.concatenate((random(1), x, random(1)))
  128. x = np.concatenate((x_herm, x[::-1].conj()))
  129. expect = fft.fft(x)
  130. assert_array_almost_equal(expect, fft.hfft(x_herm))
  131. assert_array_almost_equal(expect, fft.hfft(x_herm, norm="backward"))
  132. assert_array_almost_equal(expect / np.sqrt(30),
  133. fft.hfft(x_herm, norm="ortho"))
  134. assert_array_almost_equal(expect / 30,
  135. fft.hfft(x_herm, norm="forward"))
  136. def test_ihfft(self):
  137. x = random(14) + 1j*random(14)
  138. x_herm = np.concatenate((random(1), x, random(1)))
  139. x = np.concatenate((x_herm, x[::-1].conj()))
  140. assert_array_almost_equal(x_herm, fft.ihfft(fft.hfft(x_herm)))
  141. for norm in ["backward", "ortho", "forward"]:
  142. assert_array_almost_equal(
  143. x_herm, fft.ihfft(fft.hfft(x_herm, norm=norm), norm=norm))
  144. def test_hfft2(self):
  145. x = random((30, 20))
  146. assert_array_almost_equal(x, fft.hfft2(fft.ihfft2(x)))
  147. for norm in ["backward", "ortho", "forward"]:
  148. assert_array_almost_equal(
  149. x, fft.hfft2(fft.ihfft2(x, norm=norm), norm=norm))
  150. def test_ihfft2(self):
  151. x = random((30, 20))
  152. expect = fft.ifft2(x)[:, :11]
  153. assert_array_almost_equal(expect, fft.ihfft2(x))
  154. assert_array_almost_equal(expect, fft.ihfft2(x, norm="backward"))
  155. assert_array_almost_equal(expect * np.sqrt(30 * 20),
  156. fft.ihfft2(x, norm="ortho"))
  157. assert_array_almost_equal(expect * (30 * 20),
  158. fft.ihfft2(x, norm="forward"))
  159. def test_hfftn(self):
  160. x = random((30, 20, 10))
  161. assert_array_almost_equal(x, fft.hfftn(fft.ihfftn(x)))
  162. for norm in ["backward", "ortho", "forward"]:
  163. assert_array_almost_equal(
  164. x, fft.hfftn(fft.ihfftn(x, norm=norm), norm=norm))
  165. def test_ihfftn(self):
  166. x = random((30, 20, 10))
  167. expect = fft.ifftn(x)[:, :, :6]
  168. assert_array_almost_equal(expect, fft.ihfftn(x))
  169. assert_array_almost_equal(expect, fft.ihfftn(x, norm="backward"))
  170. assert_array_almost_equal(expect * np.sqrt(30 * 20 * 10),
  171. fft.ihfftn(x, norm="ortho"))
  172. assert_array_almost_equal(expect * (30 * 20 * 10),
  173. fft.ihfftn(x, norm="forward"))
  174. @pytest.mark.parametrize("op", [fft.fftn, fft.ifftn,
  175. fft.rfftn, fft.irfftn,
  176. fft.hfftn, fft.ihfftn])
  177. def test_axes(self, op):
  178. x = random((30, 20, 10))
  179. axes = [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]
  180. for a in axes:
  181. op_tr = op(np.transpose(x, a))
  182. tr_op = np.transpose(op(x, axes=a), a)
  183. assert_array_almost_equal(op_tr, tr_op)
  184. @pytest.mark.parametrize("op", [fft.fft2, fft.ifft2,
  185. fft.rfft2, fft.irfft2,
  186. fft.hfft2, fft.ihfft2,
  187. fft.fftn, fft.ifftn,
  188. fft.rfftn, fft.irfftn,
  189. fft.hfftn, fft.ihfftn])
  190. def test_axes_subset_with_shape(self, op):
  191. x = random((16, 8, 4))
  192. axes = [(0, 1, 2), (0, 2, 1), (1, 2, 0)]
  193. for a in axes:
  194. # different shape on the first two axes
  195. shape = tuple([2*x.shape[ax] if ax in a[:2] else x.shape[ax]
  196. for ax in range(x.ndim)])
  197. # transform only the first two axes
  198. op_tr = op(np.transpose(x, a), s=shape[:2], axes=(0, 1))
  199. tr_op = np.transpose(op(x, s=shape[:2], axes=a[:2]), a)
  200. assert_array_almost_equal(op_tr, tr_op)
  201. def test_all_1d_norm_preserving(self):
  202. # verify that round-trip transforms are norm-preserving
  203. x = random(30)
  204. x_norm = np.linalg.norm(x)
  205. n = x.size * 2
  206. func_pairs = [(fft.fft, fft.ifft),
  207. (fft.rfft, fft.irfft),
  208. # hfft: order so the first function takes x.size samples
  209. # (necessary for comparison to x_norm above)
  210. (fft.ihfft, fft.hfft),
  211. ]
  212. for forw, back in func_pairs:
  213. for n in [x.size, 2*x.size]:
  214. for norm in ['backward', 'ortho', 'forward']:
  215. tmp = forw(x, n=n, norm=norm)
  216. tmp = back(tmp, n=n, norm=norm)
  217. assert_array_almost_equal(x_norm,
  218. np.linalg.norm(tmp))
  219. @pytest.mark.parametrize("dtype", [np.half, np.single, np.double,
  220. np.longdouble])
  221. def test_dtypes(self, dtype):
  222. # make sure that all input precisions are accepted
  223. x = random(30).astype(dtype)
  224. assert_array_almost_equal(fft.ifft(fft.fft(x)), x)
  225. assert_array_almost_equal(fft.irfft(fft.rfft(x)), x)
  226. assert_array_almost_equal(fft.hfft(fft.ihfft(x), len(x)), x)
  227. @pytest.mark.parametrize(
  228. "dtype",
  229. [np.float32, np.float64, np.longfloat,
  230. np.complex64, np.complex128, np.longcomplex])
  231. @pytest.mark.parametrize("order", ["F", 'non-contiguous'])
  232. @pytest.mark.parametrize(
  233. "fft",
  234. [fft.fft, fft.fft2, fft.fftn,
  235. fft.ifft, fft.ifft2, fft.ifftn])
  236. def test_fft_with_order(dtype, order, fft):
  237. # Check that FFT/IFFT produces identical results for C, Fortran and
  238. # non contiguous arrays
  239. rng = np.random.RandomState(42)
  240. X = rng.rand(8, 7, 13).astype(dtype, copy=False)
  241. if order == 'F':
  242. Y = np.asfortranarray(X)
  243. else:
  244. # Make a non contiguous array
  245. Y = X[::-1]
  246. X = np.ascontiguousarray(X[::-1])
  247. if fft.__name__.endswith('fft'):
  248. for axis in range(3):
  249. X_res = fft(X, axis=axis)
  250. Y_res = fft(Y, axis=axis)
  251. assert_array_almost_equal(X_res, Y_res)
  252. elif fft.__name__.endswith(('fft2', 'fftn')):
  253. axes = [(0, 1), (1, 2), (0, 2)]
  254. if fft.__name__.endswith('fftn'):
  255. axes.extend([(0,), (1,), (2,), None])
  256. for ax in axes:
  257. X_res = fft(X, axes=ax)
  258. Y_res = fft(Y, axes=ax)
  259. assert_array_almost_equal(X_res, Y_res)
  260. else:
  261. raise ValueError
  262. class TestFFTThreadSafe:
  263. threads = 16
  264. input_shape = (800, 200)
  265. def _test_mtsame(self, func, *args):
  266. def worker(args, q):
  267. q.put(func(*args))
  268. q = queue.Queue()
  269. expected = func(*args)
  270. # Spin off a bunch of threads to call the same function simultaneously
  271. t = [threading.Thread(target=worker, args=(args, q))
  272. for i in range(self.threads)]
  273. [x.start() for x in t]
  274. [x.join() for x in t]
  275. # Make sure all threads returned the correct value
  276. for i in range(self.threads):
  277. assert_array_equal(q.get(timeout=5), expected,
  278. 'Function returned wrong value in multithreaded context')
  279. def test_fft(self):
  280. a = np.ones(self.input_shape, dtype=np.complex128)
  281. self._test_mtsame(fft.fft, a)
  282. def test_ifft(self):
  283. a = np.full(self.input_shape, 1+0j)
  284. self._test_mtsame(fft.ifft, a)
  285. def test_rfft(self):
  286. a = np.ones(self.input_shape)
  287. self._test_mtsame(fft.rfft, a)
  288. def test_irfft(self):
  289. a = np.full(self.input_shape, 1+0j)
  290. self._test_mtsame(fft.irfft, a)
  291. def test_hfft(self):
  292. a = np.ones(self.input_shape, np.complex64)
  293. self._test_mtsame(fft.hfft, a)
  294. def test_ihfft(self):
  295. a = np.ones(self.input_shape)
  296. self._test_mtsame(fft.ihfft, a)
  297. @pytest.mark.parametrize("func", [fft.fft, fft.ifft, fft.rfft, fft.irfft])
  298. def test_multiprocess(func):
  299. # Test that fft still works after fork (gh-10422)
  300. with multiprocessing.Pool(2) as p:
  301. res = p.map(func, [np.ones(100) for _ in range(4)])
  302. expect = func(np.ones(100))
  303. for x in res:
  304. assert_allclose(x, expect)
  305. class TestIRFFTN:
  306. def test_not_last_axis_success(self):
  307. ar, ai = np.random.random((2, 16, 8, 32))
  308. a = ar + 1j*ai
  309. axes = (-2,)
  310. # Should not raise error
  311. fft.irfftn(a, axes=axes)