test_reductions.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351
  1. from sympy.core.numbers import I
  2. from sympy.core.symbol import symbols
  3. from sympy.testing.pytest import raises
  4. from sympy.matrices import Matrix, zeros, eye
  5. from sympy.core.symbol import Symbol
  6. from sympy.core.numbers import Rational
  7. from sympy.functions.elementary.miscellaneous import sqrt
  8. from sympy.simplify.simplify import simplify
  9. from sympy.abc import x
  10. # Matrix tests
  11. def test_row_op():
  12. e = eye(3)
  13. raises(ValueError, lambda: e.elementary_row_op("abc"))
  14. raises(ValueError, lambda: e.elementary_row_op())
  15. raises(ValueError, lambda: e.elementary_row_op('n->kn', row=5, k=5))
  16. raises(ValueError, lambda: e.elementary_row_op('n->kn', row=-5, k=5))
  17. raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=1, row2=5))
  18. raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=5, row2=1))
  19. raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=-5, row2=1))
  20. raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=1, row2=-5))
  21. raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=5, k=5))
  22. raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=5, row2=1, k=5))
  23. raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=-5, row2=1, k=5))
  24. raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=-5, k=5))
  25. raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=1, k=5))
  26. # test various ways to set arguments
  27. assert e.elementary_row_op("n->kn", 0, 5) == Matrix([[5, 0, 0], [0, 1, 0], [0, 0, 1]])
  28. assert e.elementary_row_op("n->kn", 1, 5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  29. assert e.elementary_row_op("n->kn", row=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  30. assert e.elementary_row_op("n->kn", row1=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  31. assert e.elementary_row_op("n<->m", 0, 1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  32. assert e.elementary_row_op("n<->m", row1=0, row2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  33. assert e.elementary_row_op("n<->m", row=0, row2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  34. assert e.elementary_row_op("n->n+km", 0, 5, 1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
  35. assert e.elementary_row_op("n->n+km", row=0, k=5, row2=1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
  36. assert e.elementary_row_op("n->n+km", row1=0, k=5, row2=1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
  37. # make sure the matrix doesn't change size
  38. a = Matrix(2, 3, [0]*6)
  39. assert a.elementary_row_op("n->kn", 1, 5) == Matrix(2, 3, [0]*6)
  40. assert a.elementary_row_op("n<->m", 0, 1) == Matrix(2, 3, [0]*6)
  41. assert a.elementary_row_op("n->n+km", 0, 5, 1) == Matrix(2, 3, [0]*6)
  42. def test_col_op():
  43. e = eye(3)
  44. raises(ValueError, lambda: e.elementary_col_op("abc"))
  45. raises(ValueError, lambda: e.elementary_col_op())
  46. raises(ValueError, lambda: e.elementary_col_op('n->kn', col=5, k=5))
  47. raises(ValueError, lambda: e.elementary_col_op('n->kn', col=-5, k=5))
  48. raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=1, col2=5))
  49. raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=5, col2=1))
  50. raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=-5, col2=1))
  51. raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=1, col2=-5))
  52. raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=5, k=5))
  53. raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=5, col2=1, k=5))
  54. raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=-5, col2=1, k=5))
  55. raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=-5, k=5))
  56. raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=1, k=5))
  57. # test various ways to set arguments
  58. assert e.elementary_col_op("n->kn", 0, 5) == Matrix([[5, 0, 0], [0, 1, 0], [0, 0, 1]])
  59. assert e.elementary_col_op("n->kn", 1, 5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  60. assert e.elementary_col_op("n->kn", col=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  61. assert e.elementary_col_op("n->kn", col1=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
  62. assert e.elementary_col_op("n<->m", 0, 1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  63. assert e.elementary_col_op("n<->m", col1=0, col2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  64. assert e.elementary_col_op("n<->m", col=0, col2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
  65. assert e.elementary_col_op("n->n+km", 0, 5, 1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
  66. assert e.elementary_col_op("n->n+km", col=0, k=5, col2=1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
  67. assert e.elementary_col_op("n->n+km", col1=0, k=5, col2=1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
  68. # make sure the matrix doesn't change size
  69. a = Matrix(2, 3, [0]*6)
  70. assert a.elementary_col_op("n->kn", 1, 5) == Matrix(2, 3, [0]*6)
  71. assert a.elementary_col_op("n<->m", 0, 1) == Matrix(2, 3, [0]*6)
  72. assert a.elementary_col_op("n->n+km", 0, 5, 1) == Matrix(2, 3, [0]*6)
  73. def test_is_echelon():
  74. zro = zeros(3)
  75. ident = eye(3)
  76. assert zro.is_echelon
  77. assert ident.is_echelon
  78. a = Matrix(0, 0, [])
  79. assert a.is_echelon
  80. a = Matrix(2, 3, [3, 2, 1, 0, 0, 6])
  81. assert a.is_echelon
  82. a = Matrix(2, 3, [0, 0, 6, 3, 2, 1])
  83. assert not a.is_echelon
  84. x = Symbol('x')
  85. a = Matrix(3, 1, [x, 0, 0])
  86. assert a.is_echelon
  87. a = Matrix(3, 1, [x, x, 0])
  88. assert not a.is_echelon
  89. a = Matrix(3, 3, [0, 0, 0, 1, 2, 3, 0, 0, 0])
  90. assert not a.is_echelon
  91. def test_echelon_form():
  92. # echelon form is not unique, but the result
  93. # must be row-equivalent to the original matrix
  94. # and it must be in echelon form.
  95. a = zeros(3)
  96. e = eye(3)
  97. # we can assume the zero matrix and the identity matrix shouldn't change
  98. assert a.echelon_form() == a
  99. assert e.echelon_form() == e
  100. a = Matrix(0, 0, [])
  101. assert a.echelon_form() == a
  102. a = Matrix(1, 1, [5])
  103. assert a.echelon_form() == a
  104. # now we get to the real tests
  105. def verify_row_null_space(mat, rows, nulls):
  106. for v in nulls:
  107. assert all(t.is_zero for t in a_echelon*v)
  108. for v in rows:
  109. if not all(t.is_zero for t in v):
  110. assert not all(t.is_zero for t in a_echelon*v.transpose())
  111. a = Matrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
  112. nulls = [Matrix([
  113. [ 1],
  114. [-2],
  115. [ 1]])]
  116. rows = [a[i, :] for i in range(a.rows)]
  117. a_echelon = a.echelon_form()
  118. assert a_echelon.is_echelon
  119. verify_row_null_space(a, rows, nulls)
  120. a = Matrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 8])
  121. nulls = []
  122. rows = [a[i, :] for i in range(a.rows)]
  123. a_echelon = a.echelon_form()
  124. assert a_echelon.is_echelon
  125. verify_row_null_space(a, rows, nulls)
  126. a = Matrix(3, 3, [2, 1, 3, 0, 0, 0, 2, 1, 3])
  127. nulls = [Matrix([
  128. [Rational(-1, 2)],
  129. [ 1],
  130. [ 0]]),
  131. Matrix([
  132. [Rational(-3, 2)],
  133. [ 0],
  134. [ 1]])]
  135. rows = [a[i, :] for i in range(a.rows)]
  136. a_echelon = a.echelon_form()
  137. assert a_echelon.is_echelon
  138. verify_row_null_space(a, rows, nulls)
  139. # this one requires a row swap
  140. a = Matrix(3, 3, [2, 1, 3, 0, 0, 0, 1, 1, 3])
  141. nulls = [Matrix([
  142. [ 0],
  143. [ -3],
  144. [ 1]])]
  145. rows = [a[i, :] for i in range(a.rows)]
  146. a_echelon = a.echelon_form()
  147. assert a_echelon.is_echelon
  148. verify_row_null_space(a, rows, nulls)
  149. a = Matrix(3, 3, [0, 3, 3, 0, 2, 2, 0, 1, 1])
  150. nulls = [Matrix([
  151. [1],
  152. [0],
  153. [0]]),
  154. Matrix([
  155. [ 0],
  156. [-1],
  157. [ 1]])]
  158. rows = [a[i, :] for i in range(a.rows)]
  159. a_echelon = a.echelon_form()
  160. assert a_echelon.is_echelon
  161. verify_row_null_space(a, rows, nulls)
  162. a = Matrix(2, 3, [2, 2, 3, 3, 3, 0])
  163. nulls = [Matrix([
  164. [-1],
  165. [1],
  166. [0]])]
  167. rows = [a[i, :] for i in range(a.rows)]
  168. a_echelon = a.echelon_form()
  169. assert a_echelon.is_echelon
  170. verify_row_null_space(a, rows, nulls)
  171. def test_rref():
  172. e = Matrix(0, 0, [])
  173. assert e.rref(pivots=False) == e
  174. e = Matrix(1, 1, [1])
  175. a = Matrix(1, 1, [5])
  176. assert e.rref(pivots=False) == a.rref(pivots=False) == e
  177. a = Matrix(3, 1, [1, 2, 3])
  178. assert a.rref(pivots=False) == Matrix([[1], [0], [0]])
  179. a = Matrix(1, 3, [1, 2, 3])
  180. assert a.rref(pivots=False) == Matrix([[1, 2, 3]])
  181. a = Matrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
  182. assert a.rref(pivots=False) == Matrix([
  183. [1, 0, -1],
  184. [0, 1, 2],
  185. [0, 0, 0]])
  186. a = Matrix(3, 3, [1, 2, 3, 1, 2, 3, 1, 2, 3])
  187. b = Matrix(3, 3, [1, 2, 3, 0, 0, 0, 0, 0, 0])
  188. c = Matrix(3, 3, [0, 0, 0, 1, 2, 3, 0, 0, 0])
  189. d = Matrix(3, 3, [0, 0, 0, 0, 0, 0, 1, 2, 3])
  190. assert a.rref(pivots=False) == \
  191. b.rref(pivots=False) == \
  192. c.rref(pivots=False) == \
  193. d.rref(pivots=False) == b
  194. e = eye(3)
  195. z = zeros(3)
  196. assert e.rref(pivots=False) == e
  197. assert z.rref(pivots=False) == z
  198. a = Matrix([
  199. [ 0, 0, 1, 2, 2, -5, 3],
  200. [-1, 5, 2, 2, 1, -7, 5],
  201. [ 0, 0, -2, -3, -3, 8, -5],
  202. [-1, 5, 0, -1, -2, 1, 0]])
  203. mat, pivot_offsets = a.rref()
  204. assert mat == Matrix([
  205. [1, -5, 0, 0, 1, 1, -1],
  206. [0, 0, 1, 0, 0, -1, 1],
  207. [0, 0, 0, 1, 1, -2, 1],
  208. [0, 0, 0, 0, 0, 0, 0]])
  209. assert pivot_offsets == (0, 2, 3)
  210. a = Matrix([[Rational(1, 19), Rational(1, 5), 2, 3],
  211. [ 4, 5, 6, 7],
  212. [ 8, 9, 10, 11],
  213. [ 12, 13, 14, 15]])
  214. assert a.rref(pivots=False) == Matrix([
  215. [1, 0, 0, Rational(-76, 157)],
  216. [0, 1, 0, Rational(-5, 157)],
  217. [0, 0, 1, Rational(238, 157)],
  218. [0, 0, 0, 0]])
  219. x = Symbol('x')
  220. a = Matrix(2, 3, [x, 1, 1, sqrt(x), x, 1])
  221. for i, j in zip(a.rref(pivots=False),
  222. [1, 0, sqrt(x)*(-x + 1)/(-x**Rational(5, 2) + x),
  223. 0, 1, 1/(sqrt(x) + x + 1)]):
  224. assert simplify(i - j).is_zero
  225. def test_rref_rhs():
  226. a, b, c, d = symbols('a b c d')
  227. A = Matrix([[0, 0], [0, 0], [1, 2], [3, 4]])
  228. B = Matrix([a, b, c, d])
  229. assert A.rref_rhs(B) == (Matrix([
  230. [1, 0],
  231. [0, 1],
  232. [0, 0],
  233. [0, 0]]), Matrix([
  234. [ -2*c + d],
  235. [3*c/2 - d/2],
  236. [ a],
  237. [ b]]))
  238. def test_issue_17827():
  239. C = Matrix([
  240. [3, 4, -1, 1],
  241. [9, 12, -3, 3],
  242. [0, 2, 1, 3],
  243. [2, 3, 0, -2],
  244. [0, 3, 3, -5],
  245. [8, 15, 0, 6]
  246. ])
  247. # Tests for row/col within valid range
  248. D = C.elementary_row_op('n<->m', row1=2, row2=5)
  249. E = C.elementary_row_op('n->n+km', row1=5, row2=3, k=-4)
  250. F = C.elementary_row_op('n->kn', row=5, k=2)
  251. assert(D[5, :] == Matrix([[0, 2, 1, 3]]))
  252. assert(E[5, :] == Matrix([[0, 3, 0, 14]]))
  253. assert(F[5, :] == Matrix([[16, 30, 0, 12]]))
  254. # Tests for row/col out of range
  255. raises(ValueError, lambda: C.elementary_row_op('n<->m', row1=2, row2=6))
  256. raises(ValueError, lambda: C.elementary_row_op('n->kn', row=7, k=2))
  257. raises(ValueError, lambda: C.elementary_row_op('n->n+km', row1=-1, row2=5, k=2))
  258. def test_rank():
  259. m = Matrix([[1, 2], [x, 1 - 1/x]])
  260. assert m.rank() == 2
  261. n = Matrix(3, 3, range(1, 10))
  262. assert n.rank() == 2
  263. p = zeros(3)
  264. assert p.rank() == 0
  265. def test_issue_11434():
  266. ax, ay, bx, by, cx, cy, dx, dy, ex, ey, t0, t1 = \
  267. symbols('a_x a_y b_x b_y c_x c_y d_x d_y e_x e_y t_0 t_1')
  268. M = Matrix([[ax, ay, ax*t0, ay*t0, 0],
  269. [bx, by, bx*t0, by*t0, 0],
  270. [cx, cy, cx*t0, cy*t0, 1],
  271. [dx, dy, dx*t0, dy*t0, 1],
  272. [ex, ey, 2*ex*t1 - ex*t0, 2*ey*t1 - ey*t0, 0]])
  273. assert M.rank() == 4
  274. def test_rank_regression_from_so():
  275. # see:
  276. # https://stackoverflow.com/questions/19072700/why-does-sympy-give-me-the-wrong-answer-when-i-row-reduce-a-symbolic-matrix
  277. nu, lamb = symbols('nu, lambda')
  278. A = Matrix([[-3*nu, 1, 0, 0],
  279. [ 3*nu, -2*nu - 1, 2, 0],
  280. [ 0, 2*nu, (-1*nu) - lamb - 2, 3],
  281. [ 0, 0, nu + lamb, -3]])
  282. expected_reduced = Matrix([[1, 0, 0, 1/(nu**2*(-lamb - nu))],
  283. [0, 1, 0, 3/(nu*(-lamb - nu))],
  284. [0, 0, 1, 3/(-lamb - nu)],
  285. [0, 0, 0, 0]])
  286. expected_pivots = (0, 1, 2)
  287. reduced, pivots = A.rref()
  288. assert simplify(expected_reduced - reduced) == zeros(*A.shape)
  289. assert pivots == expected_pivots
  290. def test_issue_15872():
  291. A = Matrix([[1, 1, 1, 0], [-2, -1, 0, -1], [0, 0, -1, -1], [0, 0, 2, 1]])
  292. B = A - Matrix.eye(4) * I
  293. assert B.rank() == 3
  294. assert (B**2).rank() == 2
  295. assert (B**3).rank() == 2