test_special_matrices.py 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690
  1. import pytest
  2. import numpy as np
  3. from numpy import arange, add, array, eye, copy, sqrt
  4. from numpy.testing import (assert_equal, assert_array_equal,
  5. assert_array_almost_equal, assert_allclose)
  6. from pytest import raises as assert_raises
  7. from scipy.fft import fft
  8. from scipy.special import comb
  9. from scipy.linalg import (toeplitz, hankel, circulant, hadamard, leslie, dft,
  10. companion, tri, triu, tril, kron, block_diag,
  11. helmert, hilbert, invhilbert, pascal, invpascal,
  12. fiedler, fiedler_companion, eigvals,
  13. convolution_matrix)
  14. from numpy.linalg import cond
  15. def get_mat(n):
  16. data = arange(n)
  17. data = add.outer(data, data)
  18. return data
  19. class TestTri:
  20. def test_basic(self):
  21. assert_equal(tri(4), array([[1, 0, 0, 0],
  22. [1, 1, 0, 0],
  23. [1, 1, 1, 0],
  24. [1, 1, 1, 1]]))
  25. assert_equal(tri(4, dtype='f'), array([[1, 0, 0, 0],
  26. [1, 1, 0, 0],
  27. [1, 1, 1, 0],
  28. [1, 1, 1, 1]], 'f'))
  29. def test_diag(self):
  30. assert_equal(tri(4, k=1), array([[1, 1, 0, 0],
  31. [1, 1, 1, 0],
  32. [1, 1, 1, 1],
  33. [1, 1, 1, 1]]))
  34. assert_equal(tri(4, k=-1), array([[0, 0, 0, 0],
  35. [1, 0, 0, 0],
  36. [1, 1, 0, 0],
  37. [1, 1, 1, 0]]))
  38. def test_2d(self):
  39. assert_equal(tri(4, 3), array([[1, 0, 0],
  40. [1, 1, 0],
  41. [1, 1, 1],
  42. [1, 1, 1]]))
  43. assert_equal(tri(3, 4), array([[1, 0, 0, 0],
  44. [1, 1, 0, 0],
  45. [1, 1, 1, 0]]))
  46. def test_diag2d(self):
  47. assert_equal(tri(3, 4, k=2), array([[1, 1, 1, 0],
  48. [1, 1, 1, 1],
  49. [1, 1, 1, 1]]))
  50. assert_equal(tri(4, 3, k=-2), array([[0, 0, 0],
  51. [0, 0, 0],
  52. [1, 0, 0],
  53. [1, 1, 0]]))
  54. class TestTril:
  55. def test_basic(self):
  56. a = (100*get_mat(5)).astype('l')
  57. b = a.copy()
  58. for k in range(5):
  59. for l in range(k+1, 5):
  60. b[k, l] = 0
  61. assert_equal(tril(a), b)
  62. def test_diag(self):
  63. a = (100*get_mat(5)).astype('f')
  64. b = a.copy()
  65. for k in range(5):
  66. for l in range(k+3, 5):
  67. b[k, l] = 0
  68. assert_equal(tril(a, k=2), b)
  69. b = a.copy()
  70. for k in range(5):
  71. for l in range(max((k-1, 0)), 5):
  72. b[k, l] = 0
  73. assert_equal(tril(a, k=-2), b)
  74. class TestTriu:
  75. def test_basic(self):
  76. a = (100*get_mat(5)).astype('l')
  77. b = a.copy()
  78. for k in range(5):
  79. for l in range(k+1, 5):
  80. b[l, k] = 0
  81. assert_equal(triu(a), b)
  82. def test_diag(self):
  83. a = (100*get_mat(5)).astype('f')
  84. b = a.copy()
  85. for k in range(5):
  86. for l in range(max((k-1, 0)), 5):
  87. b[l, k] = 0
  88. assert_equal(triu(a, k=2), b)
  89. b = a.copy()
  90. for k in range(5):
  91. for l in range(k+3, 5):
  92. b[l, k] = 0
  93. assert_equal(triu(a, k=-2), b)
  94. class TestToeplitz:
  95. def test_basic(self):
  96. y = toeplitz([1, 2, 3])
  97. assert_array_equal(y, [[1, 2, 3], [2, 1, 2], [3, 2, 1]])
  98. y = toeplitz([1, 2, 3], [1, 4, 5])
  99. assert_array_equal(y, [[1, 4, 5], [2, 1, 4], [3, 2, 1]])
  100. def test_complex_01(self):
  101. data = (1.0 + arange(3.0)) * (1.0 + 1.0j)
  102. x = copy(data)
  103. t = toeplitz(x)
  104. # Calling toeplitz should not change x.
  105. assert_array_equal(x, data)
  106. # According to the docstring, x should be the first column of t.
  107. col0 = t[:, 0]
  108. assert_array_equal(col0, data)
  109. assert_array_equal(t[0, 1:], data[1:].conj())
  110. def test_scalar_00(self):
  111. """Scalar arguments still produce a 2D array."""
  112. t = toeplitz(10)
  113. assert_array_equal(t, [[10]])
  114. t = toeplitz(10, 20)
  115. assert_array_equal(t, [[10]])
  116. def test_scalar_01(self):
  117. c = array([1, 2, 3])
  118. t = toeplitz(c, 1)
  119. assert_array_equal(t, [[1], [2], [3]])
  120. def test_scalar_02(self):
  121. c = array([1, 2, 3])
  122. t = toeplitz(c, array(1))
  123. assert_array_equal(t, [[1], [2], [3]])
  124. def test_scalar_03(self):
  125. c = array([1, 2, 3])
  126. t = toeplitz(c, array([1]))
  127. assert_array_equal(t, [[1], [2], [3]])
  128. def test_scalar_04(self):
  129. r = array([10, 2, 3])
  130. t = toeplitz(1, r)
  131. assert_array_equal(t, [[1, 2, 3]])
  132. class TestHankel:
  133. def test_basic(self):
  134. y = hankel([1, 2, 3])
  135. assert_array_equal(y, [[1, 2, 3], [2, 3, 0], [3, 0, 0]])
  136. y = hankel([1, 2, 3], [3, 4, 5])
  137. assert_array_equal(y, [[1, 2, 3], [2, 3, 4], [3, 4, 5]])
  138. class TestCirculant:
  139. def test_basic(self):
  140. y = circulant([1, 2, 3])
  141. assert_array_equal(y, [[1, 3, 2], [2, 1, 3], [3, 2, 1]])
  142. class TestHadamard:
  143. def test_basic(self):
  144. y = hadamard(1)
  145. assert_array_equal(y, [[1]])
  146. y = hadamard(2, dtype=float)
  147. assert_array_equal(y, [[1.0, 1.0], [1.0, -1.0]])
  148. y = hadamard(4)
  149. assert_array_equal(y, [[1, 1, 1, 1],
  150. [1, -1, 1, -1],
  151. [1, 1, -1, -1],
  152. [1, -1, -1, 1]])
  153. assert_raises(ValueError, hadamard, 0)
  154. assert_raises(ValueError, hadamard, 5)
  155. class TestLeslie:
  156. def test_bad_shapes(self):
  157. assert_raises(ValueError, leslie, [[1, 1], [2, 2]], [3, 4, 5])
  158. assert_raises(ValueError, leslie, [3, 4, 5], [[1, 1], [2, 2]])
  159. assert_raises(ValueError, leslie, [1, 2], [1, 2])
  160. assert_raises(ValueError, leslie, [1], [])
  161. def test_basic(self):
  162. a = leslie([1, 2, 3], [0.25, 0.5])
  163. expected = array([[1.0, 2.0, 3.0],
  164. [0.25, 0.0, 0.0],
  165. [0.0, 0.5, 0.0]])
  166. assert_array_equal(a, expected)
  167. class TestCompanion:
  168. def test_bad_shapes(self):
  169. assert_raises(ValueError, companion, [[1, 1], [2, 2]])
  170. assert_raises(ValueError, companion, [0, 4, 5])
  171. assert_raises(ValueError, companion, [1])
  172. assert_raises(ValueError, companion, [])
  173. def test_basic(self):
  174. c = companion([1, 2, 3])
  175. expected = array([
  176. [-2.0, -3.0],
  177. [1.0, 0.0]])
  178. assert_array_equal(c, expected)
  179. c = companion([2.0, 5.0, -10.0])
  180. expected = array([
  181. [-2.5, 5.0],
  182. [1.0, 0.0]])
  183. assert_array_equal(c, expected)
  184. class TestBlockDiag:
  185. def test_basic(self):
  186. x = block_diag(eye(2), [[1, 2], [3, 4], [5, 6]], [[1, 2, 3]])
  187. assert_array_equal(x, [[1, 0, 0, 0, 0, 0, 0],
  188. [0, 1, 0, 0, 0, 0, 0],
  189. [0, 0, 1, 2, 0, 0, 0],
  190. [0, 0, 3, 4, 0, 0, 0],
  191. [0, 0, 5, 6, 0, 0, 0],
  192. [0, 0, 0, 0, 1, 2, 3]])
  193. def test_dtype(self):
  194. x = block_diag([[1.5]])
  195. assert_equal(x.dtype, float)
  196. x = block_diag([[True]])
  197. assert_equal(x.dtype, bool)
  198. def test_mixed_dtypes(self):
  199. actual = block_diag([[1]], [[1j]])
  200. desired = np.array([[1, 0], [0, 1j]])
  201. assert_array_equal(actual, desired)
  202. def test_scalar_and_1d_args(self):
  203. a = block_diag(1)
  204. assert_equal(a.shape, (1, 1))
  205. assert_array_equal(a, [[1]])
  206. a = block_diag([2, 3], 4)
  207. assert_array_equal(a, [[2, 3, 0], [0, 0, 4]])
  208. def test_bad_arg(self):
  209. assert_raises(ValueError, block_diag, [[[1]]])
  210. def test_no_args(self):
  211. a = block_diag()
  212. assert_equal(a.ndim, 2)
  213. assert_equal(a.nbytes, 0)
  214. def test_empty_matrix_arg(self):
  215. # regression test for gh-4596: check the shape of the result
  216. # for empty matrix inputs. Empty matrices are no longer ignored
  217. # (gh-4908) it is viewed as a shape (1, 0) matrix.
  218. a = block_diag([[1, 0], [0, 1]],
  219. [],
  220. [[2, 3], [4, 5], [6, 7]])
  221. assert_array_equal(a, [[1, 0, 0, 0],
  222. [0, 1, 0, 0],
  223. [0, 0, 0, 0],
  224. [0, 0, 2, 3],
  225. [0, 0, 4, 5],
  226. [0, 0, 6, 7]])
  227. def test_zerosized_matrix_arg(self):
  228. # test for gh-4908: check the shape of the result for
  229. # zero-sized matrix inputs, i.e. matrices with shape (0,n) or (n,0).
  230. # note that [[]] takes shape (1,0)
  231. a = block_diag([[1, 0], [0, 1]],
  232. [[]],
  233. [[2, 3], [4, 5], [6, 7]],
  234. np.zeros([0, 2], dtype='int32'))
  235. assert_array_equal(a, [[1, 0, 0, 0, 0, 0],
  236. [0, 1, 0, 0, 0, 0],
  237. [0, 0, 0, 0, 0, 0],
  238. [0, 0, 2, 3, 0, 0],
  239. [0, 0, 4, 5, 0, 0],
  240. [0, 0, 6, 7, 0, 0]])
  241. class TestKron:
  242. def test_basic(self):
  243. a = kron(array([[1, 2], [3, 4]]), array([[1, 1, 1]]))
  244. assert_array_equal(a, array([[1, 1, 1, 2, 2, 2],
  245. [3, 3, 3, 4, 4, 4]]))
  246. m1 = array([[1, 2], [3, 4]])
  247. m2 = array([[10], [11]])
  248. a = kron(m1, m2)
  249. expected = array([[10, 20],
  250. [11, 22],
  251. [30, 40],
  252. [33, 44]])
  253. assert_array_equal(a, expected)
  254. class TestHelmert:
  255. def test_orthogonality(self):
  256. for n in range(1, 7):
  257. H = helmert(n, full=True)
  258. Id = np.eye(n)
  259. assert_allclose(H.dot(H.T), Id, atol=1e-12)
  260. assert_allclose(H.T.dot(H), Id, atol=1e-12)
  261. def test_subspace(self):
  262. for n in range(2, 7):
  263. H_full = helmert(n, full=True)
  264. H_partial = helmert(n)
  265. for U in H_full[1:, :].T, H_partial.T:
  266. C = np.eye(n) - np.full((n, n), 1 / n)
  267. assert_allclose(U.dot(U.T), C)
  268. assert_allclose(U.T.dot(U), np.eye(n-1), atol=1e-12)
  269. class TestHilbert:
  270. def test_basic(self):
  271. h3 = array([[1.0, 1/2., 1/3.],
  272. [1/2., 1/3., 1/4.],
  273. [1/3., 1/4., 1/5.]])
  274. assert_array_almost_equal(hilbert(3), h3)
  275. assert_array_equal(hilbert(1), [[1.0]])
  276. h0 = hilbert(0)
  277. assert_equal(h0.shape, (0, 0))
  278. class TestInvHilbert:
  279. def test_basic(self):
  280. invh1 = array([[1]])
  281. assert_array_equal(invhilbert(1, exact=True), invh1)
  282. assert_array_equal(invhilbert(1), invh1)
  283. invh2 = array([[4, -6],
  284. [-6, 12]])
  285. assert_array_equal(invhilbert(2, exact=True), invh2)
  286. assert_array_almost_equal(invhilbert(2), invh2)
  287. invh3 = array([[9, -36, 30],
  288. [-36, 192, -180],
  289. [30, -180, 180]])
  290. assert_array_equal(invhilbert(3, exact=True), invh3)
  291. assert_array_almost_equal(invhilbert(3), invh3)
  292. invh4 = array([[16, -120, 240, -140],
  293. [-120, 1200, -2700, 1680],
  294. [240, -2700, 6480, -4200],
  295. [-140, 1680, -4200, 2800]])
  296. assert_array_equal(invhilbert(4, exact=True), invh4)
  297. assert_array_almost_equal(invhilbert(4), invh4)
  298. invh5 = array([[25, -300, 1050, -1400, 630],
  299. [-300, 4800, -18900, 26880, -12600],
  300. [1050, -18900, 79380, -117600, 56700],
  301. [-1400, 26880, -117600, 179200, -88200],
  302. [630, -12600, 56700, -88200, 44100]])
  303. assert_array_equal(invhilbert(5, exact=True), invh5)
  304. assert_array_almost_equal(invhilbert(5), invh5)
  305. invh17 = array([
  306. [289, -41616, 1976760, -46124400, 629598060, -5540462928,
  307. 33374693352, -143034400080, 446982500250, -1033026222800,
  308. 1774926873720, -2258997839280, 2099709530100, -1384423866000,
  309. 613101997800, -163493866080, 19835652870],
  310. [-41616, 7990272, -426980160, 10627061760, -151103534400,
  311. 1367702848512, -8410422724704, 36616806420480, -115857864064800,
  312. 270465047424000, -468580694662080, 600545887119360,
  313. -561522320049600, 372133135180800, -165537539406000,
  314. 44316454993920, -5395297580640],
  315. [1976760, -426980160, 24337869120, -630981792000, 9228108708000,
  316. -85267724461920, 532660105897920, -2348052711713280,
  317. 7504429831470000, -17664748409880000, 30818191841236800,
  318. -39732544853164800, 37341234283298400, -24857330514030000,
  319. 11100752642520000, -2982128117299200, 364182586693200],
  320. [-46124400, 10627061760, -630981792000, 16826181120000,
  321. -251209625940000, 2358021022156800, -14914482965141760,
  322. 66409571644416000, -214015221119700000, 507295338950400000,
  323. -890303319857952000, 1153715376477081600, -1089119333262870000,
  324. 727848632044800000, -326170262829600000, 87894302404608000,
  325. -10763618673376800],
  326. [629598060, -151103534400, 9228108708000,
  327. -251209625940000, 3810012660090000, -36210360321495360,
  328. 231343968720664800, -1038687206500944000, 3370739732635275000,
  329. -8037460526495400000, 14178080368737885600, -18454939322943942000,
  330. 17489975175339030000, -11728977435138600000, 5272370630081100000,
  331. -1424711708039692800, 174908803442373000],
  332. [-5540462928, 1367702848512, -85267724461920, 2358021022156800,
  333. -36210360321495360, 347619459086355456, -2239409617216035264,
  334. 10124803292907663360, -33052510749726468000,
  335. 79217210949138662400, -140362995650505067440,
  336. 183420385176741672960, -174433352415381259200,
  337. 117339159519533952000, -52892422160973595200,
  338. 14328529177999196160, -1763080738699119840],
  339. [33374693352, -8410422724704, 532660105897920,
  340. -14914482965141760, 231343968720664800, -2239409617216035264,
  341. 14527452132196331328, -66072377044391477760,
  342. 216799987176909536400, -521925895055522958000,
  343. 928414062734059661760, -1217424500995626443520,
  344. 1161358898976091015200, -783401860847777371200,
  345. 354015418167362952000, -96120549902411274240,
  346. 11851820521255194480],
  347. [-143034400080, 36616806420480, -2348052711713280,
  348. 66409571644416000, -1038687206500944000, 10124803292907663360,
  349. -66072377044391477760, 302045152202932469760,
  350. -995510145200094810000, 2405996923185123840000,
  351. -4294704507885446054400, 5649058909023744614400,
  352. -5403874060541811254400, 3654352703663101440000,
  353. -1655137020003255360000, 450325202737117593600,
  354. -55630994283442749600],
  355. [446982500250, -115857864064800, 7504429831470000,
  356. -214015221119700000, 3370739732635275000, -33052510749726468000,
  357. 216799987176909536400, -995510145200094810000,
  358. 3293967392206196062500, -7988661659013106500000,
  359. 14303908928401362270000, -18866974090684772052000,
  360. 18093328327706957325000, -12263364009096700500000,
  361. 5565847995255512250000, -1517208935002984080000,
  362. 187754605706619279900],
  363. [-1033026222800, 270465047424000, -17664748409880000,
  364. 507295338950400000, -8037460526495400000, 79217210949138662400,
  365. -521925895055522958000, 2405996923185123840000,
  366. -7988661659013106500000, 19434404971634224000000,
  367. -34894474126569249192000, 46141453390504792320000,
  368. -44349976506971935800000, 30121928988527376000000,
  369. -13697025107665828500000, 3740200989399948902400,
  370. -463591619028689580000],
  371. [1774926873720, -468580694662080,
  372. 30818191841236800, -890303319857952000, 14178080368737885600,
  373. -140362995650505067440, 928414062734059661760,
  374. -4294704507885446054400, 14303908928401362270000,
  375. -34894474126569249192000, 62810053427824648545600,
  376. -83243376594051600326400, 80177044485212743068000,
  377. -54558343880470209780000, 24851882355348879230400,
  378. -6797096028813368678400, 843736746632215035600],
  379. [-2258997839280, 600545887119360, -39732544853164800,
  380. 1153715376477081600, -18454939322943942000, 183420385176741672960,
  381. -1217424500995626443520, 5649058909023744614400,
  382. -18866974090684772052000, 46141453390504792320000,
  383. -83243376594051600326400, 110552468520163390156800,
  384. -106681852579497947388000, 72720410752415168870400,
  385. -33177973900974346080000, 9087761081682520473600,
  386. -1129631016152221783200],
  387. [2099709530100, -561522320049600, 37341234283298400,
  388. -1089119333262870000, 17489975175339030000,
  389. -174433352415381259200, 1161358898976091015200,
  390. -5403874060541811254400, 18093328327706957325000,
  391. -44349976506971935800000, 80177044485212743068000,
  392. -106681852579497947388000, 103125790826848015808400,
  393. -70409051543137015800000, 32171029219823375700000,
  394. -8824053728865840192000, 1098252376814660067000],
  395. [-1384423866000, 372133135180800,
  396. -24857330514030000, 727848632044800000, -11728977435138600000,
  397. 117339159519533952000, -783401860847777371200,
  398. 3654352703663101440000, -12263364009096700500000,
  399. 30121928988527376000000, -54558343880470209780000,
  400. 72720410752415168870400, -70409051543137015800000,
  401. 48142941226076592000000, -22027500987368499000000,
  402. 6049545098753157120000, -753830033789944188000],
  403. [613101997800, -165537539406000,
  404. 11100752642520000, -326170262829600000, 5272370630081100000,
  405. -52892422160973595200, 354015418167362952000,
  406. -1655137020003255360000, 5565847995255512250000,
  407. -13697025107665828500000, 24851882355348879230400,
  408. -33177973900974346080000, 32171029219823375700000,
  409. -22027500987368499000000, 10091416708498869000000,
  410. -2774765838662800128000, 346146444087219270000],
  411. [-163493866080, 44316454993920, -2982128117299200,
  412. 87894302404608000, -1424711708039692800,
  413. 14328529177999196160, -96120549902411274240,
  414. 450325202737117593600, -1517208935002984080000,
  415. 3740200989399948902400, -6797096028813368678400,
  416. 9087761081682520473600, -8824053728865840192000,
  417. 6049545098753157120000, -2774765838662800128000,
  418. 763806510427609497600, -95382575704033754400],
  419. [19835652870, -5395297580640, 364182586693200, -10763618673376800,
  420. 174908803442373000, -1763080738699119840, 11851820521255194480,
  421. -55630994283442749600, 187754605706619279900,
  422. -463591619028689580000, 843736746632215035600,
  423. -1129631016152221783200, 1098252376814660067000,
  424. -753830033789944188000, 346146444087219270000,
  425. -95382575704033754400, 11922821963004219300]
  426. ])
  427. assert_array_equal(invhilbert(17, exact=True), invh17)
  428. assert_allclose(invhilbert(17), invh17.astype(float), rtol=1e-12)
  429. def test_inverse(self):
  430. for n in range(1, 10):
  431. a = hilbert(n)
  432. b = invhilbert(n)
  433. # The Hilbert matrix is increasingly badly conditioned,
  434. # so take that into account in the test
  435. c = cond(a)
  436. assert_allclose(a.dot(b), eye(n), atol=1e-15*c, rtol=1e-15*c)
  437. class TestPascal:
  438. cases = [
  439. (1, array([[1]]), array([[1]])),
  440. (2, array([[1, 1],
  441. [1, 2]]),
  442. array([[1, 0],
  443. [1, 1]])),
  444. (3, array([[1, 1, 1],
  445. [1, 2, 3],
  446. [1, 3, 6]]),
  447. array([[1, 0, 0],
  448. [1, 1, 0],
  449. [1, 2, 1]])),
  450. (4, array([[1, 1, 1, 1],
  451. [1, 2, 3, 4],
  452. [1, 3, 6, 10],
  453. [1, 4, 10, 20]]),
  454. array([[1, 0, 0, 0],
  455. [1, 1, 0, 0],
  456. [1, 2, 1, 0],
  457. [1, 3, 3, 1]])),
  458. ]
  459. def check_case(self, n, sym, low):
  460. assert_array_equal(pascal(n), sym)
  461. assert_array_equal(pascal(n, kind='lower'), low)
  462. assert_array_equal(pascal(n, kind='upper'), low.T)
  463. assert_array_almost_equal(pascal(n, exact=False), sym)
  464. assert_array_almost_equal(pascal(n, exact=False, kind='lower'), low)
  465. assert_array_almost_equal(pascal(n, exact=False, kind='upper'), low.T)
  466. def test_cases(self):
  467. for n, sym, low in self.cases:
  468. self.check_case(n, sym, low)
  469. def test_big(self):
  470. p = pascal(50)
  471. assert p[-1, -1] == comb(98, 49, exact=True)
  472. def test_threshold(self):
  473. # Regression test. An early version of `pascal` returned an
  474. # array of type np.uint64 for n=35, but that data type is too small
  475. # to hold p[-1, -1]. The second assert_equal below would fail
  476. # because p[-1, -1] overflowed.
  477. p = pascal(34)
  478. assert_equal(2*p.item(-1, -2), p.item(-1, -1), err_msg="n = 34")
  479. p = pascal(35)
  480. assert_equal(2.*p.item(-1, -2), 1.*p.item(-1, -1), err_msg="n = 35")
  481. def test_invpascal():
  482. def check_invpascal(n, kind, exact):
  483. ip = invpascal(n, kind=kind, exact=exact)
  484. p = pascal(n, kind=kind, exact=exact)
  485. # Matrix-multiply ip and p, and check that we get the identity matrix.
  486. # We can't use the simple expression e = ip.dot(p), because when
  487. # n < 35 and exact is True, p.dtype is np.uint64 and ip.dtype is
  488. # np.int64. The product of those dtypes is np.float64, which loses
  489. # precision when n is greater than 18. Instead we'll cast both to
  490. # object arrays, and then multiply.
  491. e = ip.astype(object).dot(p.astype(object))
  492. assert_array_equal(e, eye(n), err_msg="n=%d kind=%r exact=%r" %
  493. (n, kind, exact))
  494. kinds = ['symmetric', 'lower', 'upper']
  495. ns = [1, 2, 5, 18]
  496. for n in ns:
  497. for kind in kinds:
  498. for exact in [True, False]:
  499. check_invpascal(n, kind, exact)
  500. ns = [19, 34, 35, 50]
  501. for n in ns:
  502. for kind in kinds:
  503. check_invpascal(n, kind, True)
  504. def test_dft():
  505. m = dft(2)
  506. expected = array([[1.0, 1.0], [1.0, -1.0]])
  507. assert_array_almost_equal(m, expected)
  508. m = dft(2, scale='n')
  509. assert_array_almost_equal(m, expected/2.0)
  510. m = dft(2, scale='sqrtn')
  511. assert_array_almost_equal(m, expected/sqrt(2.0))
  512. x = array([0, 1, 2, 3, 4, 5, 0, 1])
  513. m = dft(8)
  514. mx = m.dot(x)
  515. fx = fft(x)
  516. assert_array_almost_equal(mx, fx)
  517. def test_fiedler():
  518. f = fiedler([])
  519. assert_equal(f.size, 0)
  520. f = fiedler([123.])
  521. assert_array_equal(f, np.array([[0.]]))
  522. f = fiedler(np.arange(1, 7))
  523. des = np.array([[0, 1, 2, 3, 4, 5],
  524. [1, 0, 1, 2, 3, 4],
  525. [2, 1, 0, 1, 2, 3],
  526. [3, 2, 1, 0, 1, 2],
  527. [4, 3, 2, 1, 0, 1],
  528. [5, 4, 3, 2, 1, 0]])
  529. assert_array_equal(f, des)
  530. def test_fiedler_companion():
  531. fc = fiedler_companion([])
  532. assert_equal(fc.size, 0)
  533. fc = fiedler_companion([1.])
  534. assert_equal(fc.size, 0)
  535. fc = fiedler_companion([1., 2.])
  536. assert_array_equal(fc, np.array([[-2.]]))
  537. fc = fiedler_companion([1e-12, 2., 3.])
  538. assert_array_almost_equal(fc, companion([1e-12, 2., 3.]))
  539. with assert_raises(ValueError):
  540. fiedler_companion([0, 1, 2])
  541. fc = fiedler_companion([1., -16., 86., -176., 105.])
  542. assert_array_almost_equal(eigvals(fc),
  543. np.array([7., 5., 3., 1.]))
  544. class TestConvolutionMatrix:
  545. """
  546. Test convolution_matrix vs. numpy.convolve for various parameters.
  547. """
  548. def create_vector(self, n, cpx):
  549. """Make a complex or real test vector of length n."""
  550. x = np.linspace(-2.5, 2.2, n)
  551. if cpx:
  552. x = x + 1j*np.linspace(-1.5, 3.1, n)
  553. return x
  554. def test_bad_n(self):
  555. # n must be a positive integer
  556. with pytest.raises(ValueError, match='n must be a positive integer'):
  557. convolution_matrix([1, 2, 3], 0)
  558. def test_bad_first_arg(self):
  559. # first arg must be a 1d array, otherwise ValueError
  560. with pytest.raises(ValueError, match='one-dimensional'):
  561. convolution_matrix(1, 4)
  562. def test_empty_first_arg(self):
  563. # first arg must have at least one value
  564. with pytest.raises(ValueError, match=r'len\(a\)'):
  565. convolution_matrix([], 4)
  566. def test_bad_mode(self):
  567. # mode must be in ('full', 'valid', 'same')
  568. with pytest.raises(ValueError, match='mode.*must be one of'):
  569. convolution_matrix((1, 1), 4, mode='invalid argument')
  570. @pytest.mark.parametrize('cpx', [False, True])
  571. @pytest.mark.parametrize('na', [1, 2, 9])
  572. @pytest.mark.parametrize('nv', [1, 2, 9])
  573. @pytest.mark.parametrize('mode', [None, 'full', 'valid', 'same'])
  574. def test_against_numpy_convolve(self, cpx, na, nv, mode):
  575. a = self.create_vector(na, cpx)
  576. v = self.create_vector(nv, cpx)
  577. if mode is None:
  578. y1 = np.convolve(v, a)
  579. A = convolution_matrix(a, nv)
  580. else:
  581. y1 = np.convolve(v, a, mode)
  582. A = convolution_matrix(a, nv, mode)
  583. y2 = A @ v
  584. assert_array_almost_equal(y1, y2)