You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
347 lines
14 KiB
347 lines
14 KiB
5 months ago
|
from sympy.core.numbers import I
|
||
|
|
from sympy.core.symbol import symbols
|
||
|
|
from sympy.matrices.common import _MinimalMatrix, _CastableMatrix
|
||
|
|
from sympy.matrices.matrices import MatrixReductions
|
||
|
|
from sympy.testing.pytest import raises
|
||
|
|
from sympy.matrices import Matrix, zeros
|
||
|
|
from sympy.core.symbol import Symbol
|
||
|
|
from sympy.core.numbers import Rational
|
||
|
|
from sympy.functions.elementary.miscellaneous import sqrt
|
||
|
|
from sympy.simplify.simplify import simplify
|
||
|
|
from sympy.abc import x
|
||
|
|
|
||
|
|
class ReductionsOnlyMatrix(_MinimalMatrix, _CastableMatrix, MatrixReductions):
|
||
|
|
pass
|
||
|
|
|
||
|
|
def eye_Reductions(n):
|
||
|
|
return ReductionsOnlyMatrix(n, n, lambda i, j: int(i == j))
|
||
|
|
|
||
|
|
|
||
|
|
def zeros_Reductions(n):
|
||
|
|
return ReductionsOnlyMatrix(n, n, lambda i, j: 0)
|
||
|
|
|
||
|
|
# ReductionsOnlyMatrix tests
|
||
|
|
def test_row_op():
|
||
|
|
e = eye_Reductions(3)
|
||
|
|
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op("abc"))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op())
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->kn', row=5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->kn', row=-5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=1, row2=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=5, row2=1))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=-5, row2=1))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n<->m', row1=1, row2=-5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=5, row2=1, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=-5, row2=1, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=-5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_row_op('n->n+km', row1=1, row2=1, k=5))
|
||
|
|
|
||
|
|
# test various ways to set arguments
|
||
|
|
assert e.elementary_row_op("n->kn", 0, 5) == Matrix([[5, 0, 0], [0, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->kn", 1, 5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->kn", row=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->kn", row1=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n<->m", 0, 1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n<->m", row1=0, row2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n<->m", row=0, row2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->n+km", 0, 5, 1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->n+km", row=0, k=5, row2=1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_row_op("n->n+km", row1=0, k=5, row2=1) == Matrix([[1, 5, 0], [0, 1, 0], [0, 0, 1]])
|
||
|
|
|
||
|
|
# make sure the matrix doesn't change size
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_row_op("n->kn", 1, 5) == Matrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_row_op("n<->m", 0, 1) == Matrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_row_op("n->n+km", 0, 5, 1) == Matrix(2, 3, [0]*6)
|
||
|
|
|
||
|
|
|
||
|
|
def test_col_op():
|
||
|
|
e = eye_Reductions(3)
|
||
|
|
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op("abc"))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op())
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->kn', col=5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->kn', col=-5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=1, col2=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=5, col2=1))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=-5, col2=1))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n<->m', col1=1, col2=-5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=5, col2=1, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=-5, col2=1, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=-5, k=5))
|
||
|
|
raises(ValueError, lambda: e.elementary_col_op('n->n+km', col1=1, col2=1, k=5))
|
||
|
|
|
||
|
|
# test various ways to set arguments
|
||
|
|
assert e.elementary_col_op("n->kn", 0, 5) == Matrix([[5, 0, 0], [0, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->kn", 1, 5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->kn", col=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->kn", col1=1, k=5) == Matrix([[1, 0, 0], [0, 5, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n<->m", 0, 1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n<->m", col1=0, col2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n<->m", col=0, col2=1) == Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->n+km", 0, 5, 1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->n+km", col=0, k=5, col2=1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
|
||
|
|
assert e.elementary_col_op("n->n+km", col1=0, k=5, col2=1) == Matrix([[1, 0, 0], [5, 1, 0], [0, 0, 1]])
|
||
|
|
|
||
|
|
# make sure the matrix doesn't change size
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_col_op("n->kn", 1, 5) == Matrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_col_op("n<->m", 0, 1) == Matrix(2, 3, [0]*6)
|
||
|
|
assert a.elementary_col_op("n->n+km", 0, 5, 1) == Matrix(2, 3, [0]*6)
|
||
|
|
|
||
|
|
|
||
|
|
def test_is_echelon():
|
||
|
|
zro = zeros_Reductions(3)
|
||
|
|
ident = eye_Reductions(3)
|
||
|
|
|
||
|
|
assert zro.is_echelon
|
||
|
|
assert ident.is_echelon
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(0, 0, [])
|
||
|
|
assert a.is_echelon
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [3, 2, 1, 0, 0, 6])
|
||
|
|
assert a.is_echelon
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [0, 0, 6, 3, 2, 1])
|
||
|
|
assert not a.is_echelon
|
||
|
|
|
||
|
|
x = Symbol('x')
|
||
|
|
a = ReductionsOnlyMatrix(3, 1, [x, 0, 0])
|
||
|
|
assert a.is_echelon
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 1, [x, x, 0])
|
||
|
|
assert not a.is_echelon
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [0, 0, 0, 1, 2, 3, 0, 0, 0])
|
||
|
|
assert not a.is_echelon
|
||
|
|
|
||
|
|
|
||
|
|
def test_echelon_form():
|
||
|
|
# echelon form is not unique, but the result
|
||
|
|
# must be row-equivalent to the original matrix
|
||
|
|
# and it must be in echelon form.
|
||
|
|
|
||
|
|
a = zeros_Reductions(3)
|
||
|
|
e = eye_Reductions(3)
|
||
|
|
|
||
|
|
# we can assume the zero matrix and the identity matrix shouldn't change
|
||
|
|
assert a.echelon_form() == a
|
||
|
|
assert e.echelon_form() == e
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(0, 0, [])
|
||
|
|
assert a.echelon_form() == a
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(1, 1, [5])
|
||
|
|
assert a.echelon_form() == a
|
||
|
|
|
||
|
|
# now we get to the real tests
|
||
|
|
|
||
|
|
def verify_row_null_space(mat, rows, nulls):
|
||
|
|
for v in nulls:
|
||
|
|
assert all(t.is_zero for t in a_echelon*v)
|
||
|
|
for v in rows:
|
||
|
|
if not all(t.is_zero for t in v):
|
||
|
|
assert not all(t.is_zero for t in a_echelon*v.transpose())
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
|
||
|
|
nulls = [Matrix([
|
||
|
|
[ 1],
|
||
|
|
[-2],
|
||
|
|
[ 1]])]
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 8])
|
||
|
|
nulls = []
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [2, 1, 3, 0, 0, 0, 2, 1, 3])
|
||
|
|
nulls = [Matrix([
|
||
|
|
[Rational(-1, 2)],
|
||
|
|
[ 1],
|
||
|
|
[ 0]]),
|
||
|
|
Matrix([
|
||
|
|
[Rational(-3, 2)],
|
||
|
|
[ 0],
|
||
|
|
[ 1]])]
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
# this one requires a row swap
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [2, 1, 3, 0, 0, 0, 1, 1, 3])
|
||
|
|
nulls = [Matrix([
|
||
|
|
[ 0],
|
||
|
|
[ -3],
|
||
|
|
[ 1]])]
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [0, 3, 3, 0, 2, 2, 0, 1, 1])
|
||
|
|
nulls = [Matrix([
|
||
|
|
[1],
|
||
|
|
[0],
|
||
|
|
[0]]),
|
||
|
|
Matrix([
|
||
|
|
[ 0],
|
||
|
|
[-1],
|
||
|
|
[ 1]])]
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [2, 2, 3, 3, 3, 0])
|
||
|
|
nulls = [Matrix([
|
||
|
|
[-1],
|
||
|
|
[1],
|
||
|
|
[0]])]
|
||
|
|
rows = [a[i, :] for i in range(a.rows)]
|
||
|
|
a_echelon = a.echelon_form()
|
||
|
|
assert a_echelon.is_echelon
|
||
|
|
verify_row_null_space(a, rows, nulls)
|
||
|
|
|
||
|
|
|
||
|
|
def test_rref():
|
||
|
|
e = ReductionsOnlyMatrix(0, 0, [])
|
||
|
|
assert e.rref(pivots=False) == e
|
||
|
|
|
||
|
|
e = ReductionsOnlyMatrix(1, 1, [1])
|
||
|
|
a = ReductionsOnlyMatrix(1, 1, [5])
|
||
|
|
assert e.rref(pivots=False) == a.rref(pivots=False) == e
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 1, [1, 2, 3])
|
||
|
|
assert a.rref(pivots=False) == Matrix([[1], [0], [0]])
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(1, 3, [1, 2, 3])
|
||
|
|
assert a.rref(pivots=False) == Matrix([[1, 2, 3]])
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
|
||
|
|
assert a.rref(pivots=False) == Matrix([
|
||
|
|
[1, 0, -1],
|
||
|
|
[0, 1, 2],
|
||
|
|
[0, 0, 0]])
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix(3, 3, [1, 2, 3, 1, 2, 3, 1, 2, 3])
|
||
|
|
b = ReductionsOnlyMatrix(3, 3, [1, 2, 3, 0, 0, 0, 0, 0, 0])
|
||
|
|
c = ReductionsOnlyMatrix(3, 3, [0, 0, 0, 1, 2, 3, 0, 0, 0])
|
||
|
|
d = ReductionsOnlyMatrix(3, 3, [0, 0, 0, 0, 0, 0, 1, 2, 3])
|
||
|
|
assert a.rref(pivots=False) == \
|
||
|
|
b.rref(pivots=False) == \
|
||
|
|
c.rref(pivots=False) == \
|
||
|
|
d.rref(pivots=False) == b
|
||
|
|
|
||
|
|
e = eye_Reductions(3)
|
||
|
|
z = zeros_Reductions(3)
|
||
|
|
assert e.rref(pivots=False) == e
|
||
|
|
assert z.rref(pivots=False) == z
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix([
|
||
|
|
[ 0, 0, 1, 2, 2, -5, 3],
|
||
|
|
[-1, 5, 2, 2, 1, -7, 5],
|
||
|
|
[ 0, 0, -2, -3, -3, 8, -5],
|
||
|
|
[-1, 5, 0, -1, -2, 1, 0]])
|
||
|
|
mat, pivot_offsets = a.rref()
|
||
|
|
assert mat == Matrix([
|
||
|
|
[1, -5, 0, 0, 1, 1, -1],
|
||
|
|
[0, 0, 1, 0, 0, -1, 1],
|
||
|
|
[0, 0, 0, 1, 1, -2, 1],
|
||
|
|
[0, 0, 0, 0, 0, 0, 0]])
|
||
|
|
assert pivot_offsets == (0, 2, 3)
|
||
|
|
|
||
|
|
a = ReductionsOnlyMatrix([[Rational(1, 19), Rational(1, 5), 2, 3],
|
||
|
|
[ 4, 5, 6, 7],
|
||
|
|
[ 8, 9, 10, 11],
|
||
|
|
[ 12, 13, 14, 15]])
|
||
|
|
assert a.rref(pivots=False) == Matrix([
|
||
|
|
[1, 0, 0, Rational(-76, 157)],
|
||
|
|
[0, 1, 0, Rational(-5, 157)],
|
||
|
|
[0, 0, 1, Rational(238, 157)],
|
||
|
|
[0, 0, 0, 0]])
|
||
|
|
|
||
|
|
x = Symbol('x')
|
||
|
|
a = ReductionsOnlyMatrix(2, 3, [x, 1, 1, sqrt(x), x, 1])
|
||
|
|
for i, j in zip(a.rref(pivots=False),
|
||
|
|
[1, 0, sqrt(x)*(-x + 1)/(-x**Rational(5, 2) + x),
|
||
|
|
0, 1, 1/(sqrt(x) + x + 1)]):
|
||
|
|
assert simplify(i - j).is_zero
|
||
|
|
|
||
|
|
def test_issue_17827():
|
||
|
|
C = Matrix([
|
||
|
|
[3, 4, -1, 1],
|
||
|
|
[9, 12, -3, 3],
|
||
|
|
[0, 2, 1, 3],
|
||
|
|
[2, 3, 0, -2],
|
||
|
|
[0, 3, 3, -5],
|
||
|
|
[8, 15, 0, 6]
|
||
|
|
])
|
||
|
|
# Tests for row/col within valid range
|
||
|
|
D = C.elementary_row_op('n<->m', row1=2, row2=5)
|
||
|
|
E = C.elementary_row_op('n->n+km', row1=5, row2=3, k=-4)
|
||
|
|
F = C.elementary_row_op('n->kn', row=5, k=2)
|
||
|
|
assert(D[5, :] == Matrix([[0, 2, 1, 3]]))
|
||
|
|
assert(E[5, :] == Matrix([[0, 3, 0, 14]]))
|
||
|
|
assert(F[5, :] == Matrix([[16, 30, 0, 12]]))
|
||
|
|
# Tests for row/col out of range
|
||
|
|
raises(ValueError, lambda: C.elementary_row_op('n<->m', row1=2, row2=6))
|
||
|
|
raises(ValueError, lambda: C.elementary_row_op('n->kn', row=7, k=2))
|
||
|
|
raises(ValueError, lambda: C.elementary_row_op('n->n+km', row1=-1, row2=5, k=2))
|
||
|
|
|
||
|
|
def test_rank():
|
||
|
|
m = Matrix([[1, 2], [x, 1 - 1/x]])
|
||
|
|
assert m.rank() == 2
|
||
|
|
n = Matrix(3, 3, range(1, 10))
|
||
|
|
assert n.rank() == 2
|
||
|
|
p = zeros(3)
|
||
|
|
assert p.rank() == 0
|
||
|
|
|
||
|
|
def test_issue_11434():
|
||
|
|
ax, ay, bx, by, cx, cy, dx, dy, ex, ey, t0, t1 = \
|
||
|
|
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')
|
||
|
|
M = Matrix([[ax, ay, ax*t0, ay*t0, 0],
|
||
|
|
[bx, by, bx*t0, by*t0, 0],
|
||
|
|
[cx, cy, cx*t0, cy*t0, 1],
|
||
|
|
[dx, dy, dx*t0, dy*t0, 1],
|
||
|
|
[ex, ey, 2*ex*t1 - ex*t0, 2*ey*t1 - ey*t0, 0]])
|
||
|
|
assert M.rank() == 4
|
||
|
|
|
||
|
|
def test_rank_regression_from_so():
|
||
|
|
# see:
|
||
|
|
# https://stackoverflow.com/questions/19072700/why-does-sympy-give-me-the-wrong-answer-when-i-row-reduce-a-symbolic-matrix
|
||
|
|
|
||
|
|
nu, lamb = symbols('nu, lambda')
|
||
|
|
A = Matrix([[-3*nu, 1, 0, 0],
|
||
|
|
[ 3*nu, -2*nu - 1, 2, 0],
|
||
|
|
[ 0, 2*nu, (-1*nu) - lamb - 2, 3],
|
||
|
|
[ 0, 0, nu + lamb, -3]])
|
||
|
|
expected_reduced = Matrix([[1, 0, 0, 1/(nu**2*(-lamb - nu))],
|
||
|
|
[0, 1, 0, 3/(nu*(-lamb - nu))],
|
||
|
|
[0, 0, 1, 3/(-lamb - nu)],
|
||
|
|
[0, 0, 0, 0]])
|
||
|
|
expected_pivots = (0, 1, 2)
|
||
|
|
|
||
|
|
reduced, pivots = A.rref()
|
||
|
|
|
||
|
|
assert simplify(expected_reduced - reduced) == zeros(*A.shape)
|
||
|
|
assert pivots == expected_pivots
|
||
|
|
|
||
|
|
def test_issue_15872():
|
||
|
|
A = Matrix([[1, 1, 1, 0], [-2, -1, 0, -1], [0, 0, -1, -1], [0, 0, 2, 1]])
|
||
|
|
B = A - Matrix.eye(4) * I
|
||
|
|
assert B.rank() == 3
|
||
|
|
assert (B**2).rank() == 2
|
||
|
|
assert (B**3).rank() == 2
|