1""" 2Basic methods common to all matrices to be used 3when creating more advanced matrices (e.g., matrices over rings, 4etc.). 5""" 6 7from collections import defaultdict 8from collections.abc import Iterable 9from inspect import isfunction 10from functools import reduce 11 12from sympy.core.logic import FuzzyBool 13from sympy.assumptions.refine import refine 14from sympy.core import SympifyError, Add 15from sympy.core.basic import Atom 16from sympy.core.compatibility import as_int, is_sequence 17from sympy.core.decorators import call_highest_priority 18from sympy.core.kind import Kind, NumberKind 19from sympy.core.logic import fuzzy_and 20from sympy.core.singleton import S 21from sympy.core.symbol import Symbol 22from sympy.core.sympify import sympify 23from sympy.functions import Abs 24from sympy.polys.polytools import Poly 25from sympy.simplify import simplify as _simplify 26from sympy.simplify.simplify import dotprodsimp as _dotprodsimp 27from sympy.utilities.exceptions import SymPyDeprecationWarning 28from sympy.utilities.iterables import flatten 29from sympy.utilities.misc import filldedent 30from sympy.tensor.array import NDimArray 31 32from .utilities import _get_intermediate_simp_bool 33 34 35class MatrixError(Exception): 36 pass 37 38 39class ShapeError(ValueError, MatrixError): 40 """Wrong matrix shape""" 41 pass 42 43 44class NonSquareMatrixError(ShapeError): 45 pass 46 47 48class NonInvertibleMatrixError(ValueError, MatrixError): 49 """The matrix in not invertible (division by multidimensional zero error).""" 50 pass 51 52 53class NonPositiveDefiniteMatrixError(ValueError, MatrixError): 54 """The matrix is not a positive-definite matrix.""" 55 pass 56 57 58class MatrixRequired: 59 """All subclasses of matrix objects must implement the 60 required matrix properties listed here.""" 61 rows = None # type: int 62 cols = None # type: int 63 _simplify = None 64 65 @classmethod 66 def _new(cls, *args, **kwargs): 67 """`_new` must, at minimum, be callable as 68 `_new(rows, cols, mat) where mat is a flat list of the 69 elements of the matrix.""" 70 raise NotImplementedError("Subclasses must implement this.") 71 72 def __eq__(self, other): 73 raise NotImplementedError("Subclasses must implement this.") 74 75 def __getitem__(self, key): 76 """Implementations of __getitem__ should accept ints, in which 77 case the matrix is indexed as a flat list, tuples (i,j) in which 78 case the (i,j) entry is returned, slices, or mixed tuples (a,b) 79 where a and b are any combintion of slices and integers.""" 80 raise NotImplementedError("Subclasses must implement this.") 81 82 def __len__(self): 83 """The total number of entries in the matrix.""" 84 raise NotImplementedError("Subclasses must implement this.") 85 86 @property 87 def shape(self): 88 raise NotImplementedError("Subclasses must implement this.") 89 90 91class MatrixShaping(MatrixRequired): 92 """Provides basic matrix shaping and extracting of submatrices""" 93 94 def _eval_col_del(self, col): 95 def entry(i, j): 96 return self[i, j] if j < col else self[i, j + 1] 97 return self._new(self.rows, self.cols - 1, entry) 98 99 def _eval_col_insert(self, pos, other): 100 101 def entry(i, j): 102 if j < pos: 103 return self[i, j] 104 elif pos <= j < pos + other.cols: 105 return other[i, j - pos] 106 return self[i, j - other.cols] 107 108 return self._new(self.rows, self.cols + other.cols, 109 lambda i, j: entry(i, j)) 110 111 def _eval_col_join(self, other): 112 rows = self.rows 113 114 def entry(i, j): 115 if i < rows: 116 return self[i, j] 117 return other[i - rows, j] 118 119 return classof(self, other)._new(self.rows + other.rows, self.cols, 120 lambda i, j: entry(i, j)) 121 122 def _eval_extract(self, rowsList, colsList): 123 mat = list(self) 124 cols = self.cols 125 indices = (i * cols + j for i in rowsList for j in colsList) 126 return self._new(len(rowsList), len(colsList), 127 list(mat[i] for i in indices)) 128 129 def _eval_get_diag_blocks(self): 130 sub_blocks = [] 131 132 def recurse_sub_blocks(M): 133 i = 1 134 while i <= M.shape[0]: 135 if i == 1: 136 to_the_right = M[0, i:] 137 to_the_bottom = M[i:, 0] 138 else: 139 to_the_right = M[:i, i:] 140 to_the_bottom = M[i:, :i] 141 if any(to_the_right) or any(to_the_bottom): 142 i += 1 143 continue 144 else: 145 sub_blocks.append(M[:i, :i]) 146 if M.shape == M[:i, :i].shape: 147 return 148 else: 149 recurse_sub_blocks(M[i:, i:]) 150 return 151 152 recurse_sub_blocks(self) 153 return sub_blocks 154 155 def _eval_row_del(self, row): 156 def entry(i, j): 157 return self[i, j] if i < row else self[i + 1, j] 158 return self._new(self.rows - 1, self.cols, entry) 159 160 def _eval_row_insert(self, pos, other): 161 entries = list(self) 162 insert_pos = pos * self.cols 163 entries[insert_pos:insert_pos] = list(other) 164 return self._new(self.rows + other.rows, self.cols, entries) 165 166 def _eval_row_join(self, other): 167 cols = self.cols 168 169 def entry(i, j): 170 if j < cols: 171 return self[i, j] 172 return other[i, j - cols] 173 174 return classof(self, other)._new(self.rows, self.cols + other.cols, 175 lambda i, j: entry(i, j)) 176 177 def _eval_tolist(self): 178 return [list(self[i,:]) for i in range(self.rows)] 179 180 def _eval_todok(self): 181 dok = {} 182 rows, cols = self.shape 183 for i in range(rows): 184 for j in range(cols): 185 val = self[i, j] 186 if val != self.zero: 187 dok[i, j] = val 188 return dok 189 190 def _eval_vec(self): 191 rows = self.rows 192 193 def entry(n, _): 194 # we want to read off the columns first 195 j = n // rows 196 i = n - j * rows 197 return self[i, j] 198 199 return self._new(len(self), 1, entry) 200 201 def _eval_vech(self, diagonal): 202 c = self.cols 203 v = [] 204 if diagonal: 205 for j in range(c): 206 for i in range(j, c): 207 v.append(self[i, j]) 208 else: 209 for j in range(c): 210 for i in range(j + 1, c): 211 v.append(self[i, j]) 212 return self._new(len(v), 1, v) 213 214 def col_del(self, col): 215 """Delete the specified column.""" 216 if col < 0: 217 col += self.cols 218 if not 0 <= col < self.cols: 219 raise IndexError("Column {} is out of range.".format(col)) 220 return self._eval_col_del(col) 221 222 def col_insert(self, pos, other): 223 """Insert one or more columns at the given column position. 224 225 Examples 226 ======== 227 228 >>> from sympy import zeros, ones 229 >>> M = zeros(3) 230 >>> V = ones(3, 1) 231 >>> M.col_insert(1, V) 232 Matrix([ 233 [0, 1, 0, 0], 234 [0, 1, 0, 0], 235 [0, 1, 0, 0]]) 236 237 See Also 238 ======== 239 240 col 241 row_insert 242 """ 243 # Allows you to build a matrix even if it is null matrix 244 if not self: 245 return type(self)(other) 246 247 pos = as_int(pos) 248 249 if pos < 0: 250 pos = self.cols + pos 251 if pos < 0: 252 pos = 0 253 elif pos > self.cols: 254 pos = self.cols 255 256 if self.rows != other.rows: 257 raise ShapeError( 258 "`self` and `other` must have the same number of rows.") 259 260 return self._eval_col_insert(pos, other) 261 262 def col_join(self, other): 263 """Concatenates two matrices along self's last and other's first row. 264 265 Examples 266 ======== 267 268 >>> from sympy import zeros, ones 269 >>> M = zeros(3) 270 >>> V = ones(1, 3) 271 >>> M.col_join(V) 272 Matrix([ 273 [0, 0, 0], 274 [0, 0, 0], 275 [0, 0, 0], 276 [1, 1, 1]]) 277 278 See Also 279 ======== 280 281 col 282 row_join 283 """ 284 # A null matrix can always be stacked (see #10770) 285 if self.rows == 0 and self.cols != other.cols: 286 return self._new(0, other.cols, []).col_join(other) 287 288 if self.cols != other.cols: 289 raise ShapeError( 290 "`self` and `other` must have the same number of columns.") 291 return self._eval_col_join(other) 292 293 def col(self, j): 294 """Elementary column selector. 295 296 Examples 297 ======== 298 299 >>> from sympy import eye 300 >>> eye(2).col(0) 301 Matrix([ 302 [1], 303 [0]]) 304 305 See Also 306 ======== 307 308 row 309 col_del 310 col_join 311 col_insert 312 """ 313 return self[:, j] 314 315 def extract(self, rowsList, colsList): 316 """Return a submatrix by specifying a list of rows and columns. 317 Negative indices can be given. All indices must be in the range 318 -n <= i < n where n is the number of rows or columns. 319 320 Examples 321 ======== 322 323 >>> from sympy import Matrix 324 >>> m = Matrix(4, 3, range(12)) 325 >>> m 326 Matrix([ 327 [0, 1, 2], 328 [3, 4, 5], 329 [6, 7, 8], 330 [9, 10, 11]]) 331 >>> m.extract([0, 1, 3], [0, 1]) 332 Matrix([ 333 [0, 1], 334 [3, 4], 335 [9, 10]]) 336 337 Rows or columns can be repeated: 338 339 >>> m.extract([0, 0, 1], [-1]) 340 Matrix([ 341 [2], 342 [2], 343 [5]]) 344 345 Every other row can be taken by using range to provide the indices: 346 347 >>> m.extract(range(0, m.rows, 2), [-1]) 348 Matrix([ 349 [2], 350 [8]]) 351 352 RowsList or colsList can also be a list of booleans, in which case 353 the rows or columns corresponding to the True values will be selected: 354 355 >>> m.extract([0, 1, 2, 3], [True, False, True]) 356 Matrix([ 357 [0, 2], 358 [3, 5], 359 [6, 8], 360 [9, 11]]) 361 """ 362 363 if not is_sequence(rowsList) or not is_sequence(colsList): 364 raise TypeError("rowsList and colsList must be iterable") 365 # ensure rowsList and colsList are lists of integers 366 if rowsList and all(isinstance(i, bool) for i in rowsList): 367 rowsList = [index for index, item in enumerate(rowsList) if item] 368 if colsList and all(isinstance(i, bool) for i in colsList): 369 colsList = [index for index, item in enumerate(colsList) if item] 370 371 # ensure everything is in range 372 rowsList = [a2idx(k, self.rows) for k in rowsList] 373 colsList = [a2idx(k, self.cols) for k in colsList] 374 375 return self._eval_extract(rowsList, colsList) 376 377 def get_diag_blocks(self): 378 """Obtains the square sub-matrices on the main diagonal of a square matrix. 379 380 Useful for inverting symbolic matrices or solving systems of 381 linear equations which may be decoupled by having a block diagonal 382 structure. 383 384 Examples 385 ======== 386 387 >>> from sympy import Matrix 388 >>> from sympy.abc import x, y, z 389 >>> A = Matrix([[1, 3, 0, 0], [y, z*z, 0, 0], [0, 0, x, 0], [0, 0, 0, 0]]) 390 >>> a1, a2, a3 = A.get_diag_blocks() 391 >>> a1 392 Matrix([ 393 [1, 3], 394 [y, z**2]]) 395 >>> a2 396 Matrix([[x]]) 397 >>> a3 398 Matrix([[0]]) 399 400 """ 401 return self._eval_get_diag_blocks() 402 403 @classmethod 404 def hstack(cls, *args): 405 """Return a matrix formed by joining args horizontally (i.e. 406 by repeated application of row_join). 407 408 Examples 409 ======== 410 411 >>> from sympy.matrices import Matrix, eye 412 >>> Matrix.hstack(eye(2), 2*eye(2)) 413 Matrix([ 414 [1, 0, 2, 0], 415 [0, 1, 0, 2]]) 416 """ 417 if len(args) == 0: 418 return cls._new() 419 420 kls = type(args[0]) 421 return reduce(kls.row_join, args) 422 423 def reshape(self, rows, cols): 424 """Reshape the matrix. Total number of elements must remain the same. 425 426 Examples 427 ======== 428 429 >>> from sympy import Matrix 430 >>> m = Matrix(2, 3, lambda i, j: 1) 431 >>> m 432 Matrix([ 433 [1, 1, 1], 434 [1, 1, 1]]) 435 >>> m.reshape(1, 6) 436 Matrix([[1, 1, 1, 1, 1, 1]]) 437 >>> m.reshape(3, 2) 438 Matrix([ 439 [1, 1], 440 [1, 1], 441 [1, 1]]) 442 443 """ 444 if self.rows * self.cols != rows * cols: 445 raise ValueError("Invalid reshape parameters %d %d" % (rows, cols)) 446 return self._new(rows, cols, lambda i, j: self[i * cols + j]) 447 448 def row_del(self, row): 449 """Delete the specified row.""" 450 if row < 0: 451 row += self.rows 452 if not 0 <= row < self.rows: 453 raise IndexError("Row {} is out of range.".format(row)) 454 455 return self._eval_row_del(row) 456 457 def row_insert(self, pos, other): 458 """Insert one or more rows at the given row position. 459 460 Examples 461 ======== 462 463 >>> from sympy import zeros, ones 464 >>> M = zeros(3) 465 >>> V = ones(1, 3) 466 >>> M.row_insert(1, V) 467 Matrix([ 468 [0, 0, 0], 469 [1, 1, 1], 470 [0, 0, 0], 471 [0, 0, 0]]) 472 473 See Also 474 ======== 475 476 row 477 col_insert 478 """ 479 # Allows you to build a matrix even if it is null matrix 480 if not self: 481 return self._new(other) 482 483 pos = as_int(pos) 484 485 if pos < 0: 486 pos = self.rows + pos 487 if pos < 0: 488 pos = 0 489 elif pos > self.rows: 490 pos = self.rows 491 492 if self.cols != other.cols: 493 raise ShapeError( 494 "`self` and `other` must have the same number of columns.") 495 496 return self._eval_row_insert(pos, other) 497 498 def row_join(self, other): 499 """Concatenates two matrices along self's last and rhs's first column 500 501 Examples 502 ======== 503 504 >>> from sympy import zeros, ones 505 >>> M = zeros(3) 506 >>> V = ones(3, 1) 507 >>> M.row_join(V) 508 Matrix([ 509 [0, 0, 0, 1], 510 [0, 0, 0, 1], 511 [0, 0, 0, 1]]) 512 513 See Also 514 ======== 515 516 row 517 col_join 518 """ 519 # A null matrix can always be stacked (see #10770) 520 if self.cols == 0 and self.rows != other.rows: 521 return self._new(other.rows, 0, []).row_join(other) 522 523 if self.rows != other.rows: 524 raise ShapeError( 525 "`self` and `rhs` must have the same number of rows.") 526 return self._eval_row_join(other) 527 528 def diagonal(self, k=0): 529 """Returns the kth diagonal of self. The main diagonal 530 corresponds to `k=0`; diagonals above and below correspond to 531 `k > 0` and `k < 0`, respectively. The values of `self[i, j]` 532 for which `j - i = k`, are returned in order of increasing 533 `i + j`, starting with `i + j = |k|`. 534 535 Examples 536 ======== 537 538 >>> from sympy import Matrix 539 >>> m = Matrix(3, 3, lambda i, j: j - i); m 540 Matrix([ 541 [ 0, 1, 2], 542 [-1, 0, 1], 543 [-2, -1, 0]]) 544 >>> _.diagonal() 545 Matrix([[0, 0, 0]]) 546 >>> m.diagonal(1) 547 Matrix([[1, 1]]) 548 >>> m.diagonal(-2) 549 Matrix([[-2]]) 550 551 Even though the diagonal is returned as a Matrix, the element 552 retrieval can be done with a single index: 553 554 >>> Matrix.diag(1, 2, 3).diagonal()[1] # instead of [0, 1] 555 2 556 557 See Also 558 ======== 559 diag - to create a diagonal matrix 560 """ 561 rv = [] 562 k = as_int(k) 563 r = 0 if k > 0 else -k 564 c = 0 if r else k 565 while True: 566 if r == self.rows or c == self.cols: 567 break 568 rv.append(self[r, c]) 569 r += 1 570 c += 1 571 if not rv: 572 raise ValueError(filldedent(''' 573 The %s diagonal is out of range [%s, %s]''' % ( 574 k, 1 - self.rows, self.cols - 1))) 575 return self._new(1, len(rv), rv) 576 577 def row(self, i): 578 """Elementary row selector. 579 580 Examples 581 ======== 582 583 >>> from sympy import eye 584 >>> eye(2).row(0) 585 Matrix([[1, 0]]) 586 587 See Also 588 ======== 589 590 col 591 row_del 592 row_join 593 row_insert 594 """ 595 return self[i, :] 596 597 @property 598 def shape(self): 599 """The shape (dimensions) of the matrix as the 2-tuple (rows, cols). 600 601 Examples 602 ======== 603 604 >>> from sympy.matrices import zeros 605 >>> M = zeros(2, 3) 606 >>> M.shape 607 (2, 3) 608 >>> M.rows 609 2 610 >>> M.cols 611 3 612 """ 613 return (self.rows, self.cols) 614 615 def todok(self): 616 """Return the matrix as dictionary of keys. 617 618 Examples 619 ======== 620 621 >>> from sympy import Matrix 622 >>> M = Matrix.eye(3) 623 >>> M.todok() 624 {(0, 0): 1, (1, 1): 1, (2, 2): 1} 625 """ 626 return self._eval_todok() 627 628 def tolist(self): 629 """Return the Matrix as a nested Python list. 630 631 Examples 632 ======== 633 634 >>> from sympy import Matrix, ones 635 >>> m = Matrix(3, 3, range(9)) 636 >>> m 637 Matrix([ 638 [0, 1, 2], 639 [3, 4, 5], 640 [6, 7, 8]]) 641 >>> m.tolist() 642 [[0, 1, 2], [3, 4, 5], [6, 7, 8]] 643 >>> ones(3, 0).tolist() 644 [[], [], []] 645 646 When there are no rows then it will not be possible to tell how 647 many columns were in the original matrix: 648 649 >>> ones(0, 3).tolist() 650 [] 651 652 """ 653 if not self.rows: 654 return [] 655 if not self.cols: 656 return [[] for i in range(self.rows)] 657 return self._eval_tolist() 658 659 def todod(M): 660 """Returns matrix as dict of dicts containing non-zero elements of the Matrix 661 662 Examples 663 ======== 664 665 >>> from sympy import Matrix 666 >>> A = Matrix([[0, 1],[0, 3]]) 667 >>> A 668 Matrix([ 669 [0, 1], 670 [0, 3]]) 671 >>> A.todod() 672 {0: {1: 1}, 1: {1: 3}} 673 674 675 """ 676 rowsdict = {} 677 Mlol = M.tolist() 678 for i, Mi in enumerate(Mlol): 679 row = {j: Mij for j, Mij in enumerate(Mi) if Mij} 680 if row: 681 rowsdict[i] = row 682 return rowsdict 683 684 def vec(self): 685 """Return the Matrix converted into a one column matrix by stacking columns 686 687 Examples 688 ======== 689 690 >>> from sympy import Matrix 691 >>> m=Matrix([[1, 3], [2, 4]]) 692 >>> m 693 Matrix([ 694 [1, 3], 695 [2, 4]]) 696 >>> m.vec() 697 Matrix([ 698 [1], 699 [2], 700 [3], 701 [4]]) 702 703 See Also 704 ======== 705 706 vech 707 """ 708 return self._eval_vec() 709 710 def vech(self, diagonal=True, check_symmetry=True): 711 """Reshapes the matrix into a column vector by stacking the 712 elements in the lower triangle. 713 714 Parameters 715 ========== 716 717 diagonal : bool, optional 718 If ``True``, it includes the diagonal elements. 719 720 check_symmetry : bool, optional 721 If ``True``, it checks whether the matrix is symmetric. 722 723 Examples 724 ======== 725 726 >>> from sympy import Matrix 727 >>> m=Matrix([[1, 2], [2, 3]]) 728 >>> m 729 Matrix([ 730 [1, 2], 731 [2, 3]]) 732 >>> m.vech() 733 Matrix([ 734 [1], 735 [2], 736 [3]]) 737 >>> m.vech(diagonal=False) 738 Matrix([[2]]) 739 740 Notes 741 ===== 742 743 This should work for symmetric matrices and ``vech`` can 744 represent symmetric matrices in vector form with less size than 745 ``vec``. 746 747 See Also 748 ======== 749 750 vec 751 """ 752 if not self.is_square: 753 raise NonSquareMatrixError 754 755 if check_symmetry and not self.is_symmetric(): 756 raise ValueError("The matrix is not symmetric.") 757 758 return self._eval_vech(diagonal) 759 760 @classmethod 761 def vstack(cls, *args): 762 """Return a matrix formed by joining args vertically (i.e. 763 by repeated application of col_join). 764 765 Examples 766 ======== 767 768 >>> from sympy.matrices import Matrix, eye 769 >>> Matrix.vstack(eye(2), 2*eye(2)) 770 Matrix([ 771 [1, 0], 772 [0, 1], 773 [2, 0], 774 [0, 2]]) 775 """ 776 if len(args) == 0: 777 return cls._new() 778 779 kls = type(args[0]) 780 return reduce(kls.col_join, args) 781 782 783class MatrixSpecial(MatrixRequired): 784 """Construction of special matrices""" 785 786 @classmethod 787 def _eval_diag(cls, rows, cols, diag_dict): 788 """diag_dict is a defaultdict containing 789 all the entries of the diagonal matrix.""" 790 def entry(i, j): 791 return diag_dict[(i, j)] 792 return cls._new(rows, cols, entry) 793 794 @classmethod 795 def _eval_eye(cls, rows, cols): 796 vals = [cls.zero]*(rows*cols) 797 vals[::cols+1] = [cls.one]*min(rows, cols) 798 return cls._new(rows, cols, vals, copy=False) 799 800 @classmethod 801 def _eval_jordan_block(cls, rows, cols, eigenvalue, band='upper'): 802 if band == 'lower': 803 def entry(i, j): 804 if i == j: 805 return eigenvalue 806 elif j + 1 == i: 807 return cls.one 808 return cls.zero 809 else: 810 def entry(i, j): 811 if i == j: 812 return eigenvalue 813 elif i + 1 == j: 814 return cls.one 815 return cls.zero 816 return cls._new(rows, cols, entry) 817 818 @classmethod 819 def _eval_ones(cls, rows, cols): 820 def entry(i, j): 821 return cls.one 822 return cls._new(rows, cols, entry) 823 824 @classmethod 825 def _eval_zeros(cls, rows, cols): 826 return cls._new(rows, cols, [cls.zero]*(rows*cols), copy=False) 827 828 @classmethod 829 def _eval_wilkinson(cls, n): 830 def entry(i, j): 831 return cls.one if i + 1 == j else cls.zero 832 833 D = cls._new(2*n + 1, 2*n + 1, entry) 834 835 wminus = cls.diag([i for i in range(-n, n + 1)], unpack=True) + D + D.T 836 wplus = abs(cls.diag([i for i in range(-n, n + 1)], unpack=True)) + D + D.T 837 838 return wminus, wplus 839 840 @classmethod 841 def diag(kls, *args, strict=False, unpack=True, rows=None, cols=None, **kwargs): 842 """Returns a matrix with the specified diagonal. 843 If matrices are passed, a block-diagonal matrix 844 is created (i.e. the "direct sum" of the matrices). 845 846 kwargs 847 ====== 848 849 rows : rows of the resulting matrix; computed if 850 not given. 851 852 cols : columns of the resulting matrix; computed if 853 not given. 854 855 cls : class for the resulting matrix 856 857 unpack : bool which, when True (default), unpacks a single 858 sequence rather than interpreting it as a Matrix. 859 860 strict : bool which, when False (default), allows Matrices to 861 have variable-length rows. 862 863 Examples 864 ======== 865 866 >>> from sympy.matrices import Matrix 867 >>> Matrix.diag(1, 2, 3) 868 Matrix([ 869 [1, 0, 0], 870 [0, 2, 0], 871 [0, 0, 3]]) 872 873 The current default is to unpack a single sequence. If this is 874 not desired, set `unpack=False` and it will be interpreted as 875 a matrix. 876 877 >>> Matrix.diag([1, 2, 3]) == Matrix.diag(1, 2, 3) 878 True 879 880 When more than one element is passed, each is interpreted as 881 something to put on the diagonal. Lists are converted to 882 matrices. Filling of the diagonal always continues from 883 the bottom right hand corner of the previous item: this 884 will create a block-diagonal matrix whether the matrices 885 are square or not. 886 887 >>> col = [1, 2, 3] 888 >>> row = [[4, 5]] 889 >>> Matrix.diag(col, row) 890 Matrix([ 891 [1, 0, 0], 892 [2, 0, 0], 893 [3, 0, 0], 894 [0, 4, 5]]) 895 896 When `unpack` is False, elements within a list need not all be 897 of the same length. Setting `strict` to True would raise a 898 ValueError for the following: 899 900 >>> Matrix.diag([[1, 2, 3], [4, 5], [6]], unpack=False) 901 Matrix([ 902 [1, 2, 3], 903 [4, 5, 0], 904 [6, 0, 0]]) 905 906 The type of the returned matrix can be set with the ``cls`` 907 keyword. 908 909 >>> from sympy.matrices import ImmutableMatrix 910 >>> from sympy.utilities.misc import func_name 911 >>> func_name(Matrix.diag(1, cls=ImmutableMatrix)) 912 'ImmutableDenseMatrix' 913 914 A zero dimension matrix can be used to position the start of 915 the filling at the start of an arbitrary row or column: 916 917 >>> from sympy import ones 918 >>> r2 = ones(0, 2) 919 >>> Matrix.diag(r2, 1, 2) 920 Matrix([ 921 [0, 0, 1, 0], 922 [0, 0, 0, 2]]) 923 924 See Also 925 ======== 926 eye 927 diagonal - to extract a diagonal 928 .dense.diag 929 .expressions.blockmatrix.BlockMatrix 930 .sparsetools.banded - to create multi-diagonal matrices 931 """ 932 from sympy.matrices.matrices import MatrixBase 933 from sympy.matrices.dense import Matrix 934 from sympy.matrices.sparse import SparseMatrix 935 klass = kwargs.get('cls', kls) 936 if unpack and len(args) == 1 and is_sequence(args[0]) and \ 937 not isinstance(args[0], MatrixBase): 938 args = args[0] 939 940 # fill a default dict with the diagonal entries 941 diag_entries = defaultdict(int) 942 rmax = cmax = 0 # keep track of the biggest index seen 943 for m in args: 944 if isinstance(m, list): 945 if strict: 946 # if malformed, Matrix will raise an error 947 _ = Matrix(m) 948 r, c = _.shape 949 m = _.tolist() 950 else: 951 r, c, smat = SparseMatrix._handle_creation_inputs(m) 952 for (i, j), _ in smat.items(): 953 diag_entries[(i + rmax, j + cmax)] = _ 954 m = [] # to skip process below 955 elif hasattr(m, 'shape'): # a Matrix 956 # convert to list of lists 957 r, c = m.shape 958 m = m.tolist() 959 else: # in this case, we're a single value 960 diag_entries[(rmax, cmax)] = m 961 rmax += 1 962 cmax += 1 963 continue 964 # process list of lists 965 for i in range(len(m)): 966 for j, _ in enumerate(m[i]): 967 diag_entries[(i + rmax, j + cmax)] = _ 968 rmax += r 969 cmax += c 970 if rows is None: 971 rows, cols = cols, rows 972 if rows is None: 973 rows, cols = rmax, cmax 974 else: 975 cols = rows if cols is None else cols 976 if rows < rmax or cols < cmax: 977 raise ValueError(filldedent(''' 978 The constructed matrix is {} x {} but a size of {} x {} 979 was specified.'''.format(rmax, cmax, rows, cols))) 980 return klass._eval_diag(rows, cols, diag_entries) 981 982 @classmethod 983 def eye(kls, rows, cols=None, **kwargs): 984 """Returns an identity matrix. 985 986 Args 987 ==== 988 989 rows : rows of the matrix 990 cols : cols of the matrix (if None, cols=rows) 991 992 kwargs 993 ====== 994 cls : class of the returned matrix 995 """ 996 if cols is None: 997 cols = rows 998 if rows < 0 or cols < 0: 999 raise ValueError("Cannot create a {} x {} matrix. " 1000 "Both dimensions must be positive".format(rows, cols)) 1001 klass = kwargs.get('cls', kls) 1002 rows, cols = as_int(rows), as_int(cols) 1003 1004 return klass._eval_eye(rows, cols) 1005 1006 @classmethod 1007 def jordan_block(kls, size=None, eigenvalue=None, *, band='upper', **kwargs): 1008 """Returns a Jordan block 1009 1010 Parameters 1011 ========== 1012 1013 size : Integer, optional 1014 Specifies the shape of the Jordan block matrix. 1015 1016 eigenvalue : Number or Symbol 1017 Specifies the value for the main diagonal of the matrix. 1018 1019 .. note:: 1020 The keyword ``eigenval`` is also specified as an alias 1021 of this keyword, but it is not recommended to use. 1022 1023 We may deprecate the alias in later release. 1024 1025 band : 'upper' or 'lower', optional 1026 Specifies the position of the off-diagonal to put `1` s on. 1027 1028 cls : Matrix, optional 1029 Specifies the matrix class of the output form. 1030 1031 If it is not specified, the class type where the method is 1032 being executed on will be returned. 1033 1034 rows, cols : Integer, optional 1035 Specifies the shape of the Jordan block matrix. See Notes 1036 section for the details of how these key works. 1037 1038 .. note:: 1039 This feature will be deprecated in the future. 1040 1041 1042 Returns 1043 ======= 1044 1045 Matrix 1046 A Jordan block matrix. 1047 1048 Raises 1049 ====== 1050 1051 ValueError 1052 If insufficient arguments are given for matrix size 1053 specification, or no eigenvalue is given. 1054 1055 Examples 1056 ======== 1057 1058 Creating a default Jordan block: 1059 1060 >>> from sympy import Matrix 1061 >>> from sympy.abc import x 1062 >>> Matrix.jordan_block(4, x) 1063 Matrix([ 1064 [x, 1, 0, 0], 1065 [0, x, 1, 0], 1066 [0, 0, x, 1], 1067 [0, 0, 0, x]]) 1068 1069 Creating an alternative Jordan block matrix where `1` is on 1070 lower off-diagonal: 1071 1072 >>> Matrix.jordan_block(4, x, band='lower') 1073 Matrix([ 1074 [x, 0, 0, 0], 1075 [1, x, 0, 0], 1076 [0, 1, x, 0], 1077 [0, 0, 1, x]]) 1078 1079 Creating a Jordan block with keyword arguments 1080 1081 >>> Matrix.jordan_block(size=4, eigenvalue=x) 1082 Matrix([ 1083 [x, 1, 0, 0], 1084 [0, x, 1, 0], 1085 [0, 0, x, 1], 1086 [0, 0, 0, x]]) 1087 1088 Notes 1089 ===== 1090 1091 .. note:: 1092 This feature will be deprecated in the future. 1093 1094 The keyword arguments ``size``, ``rows``, ``cols`` relates to 1095 the Jordan block size specifications. 1096 1097 If you want to create a square Jordan block, specify either 1098 one of the three arguments. 1099 1100 If you want to create a rectangular Jordan block, specify 1101 ``rows`` and ``cols`` individually. 1102 1103 +--------------------------------+---------------------+ 1104 | Arguments Given | Matrix Shape | 1105 +----------+----------+----------+----------+----------+ 1106 | size | rows | cols | rows | cols | 1107 +==========+==========+==========+==========+==========+ 1108 | size | Any | size | size | 1109 +----------+----------+----------+----------+----------+ 1110 | | None | ValueError | 1111 | +----------+----------+----------+----------+ 1112 | None | rows | None | rows | rows | 1113 | +----------+----------+----------+----------+ 1114 | | None | cols | cols | cols | 1115 + +----------+----------+----------+----------+ 1116 | | rows | cols | rows | cols | 1117 +----------+----------+----------+----------+----------+ 1118 1119 References 1120 ========== 1121 1122 .. [1] https://en.wikipedia.org/wiki/Jordan_matrix 1123 """ 1124 if 'rows' in kwargs or 'cols' in kwargs: 1125 SymPyDeprecationWarning( 1126 feature="Keyword arguments 'rows' or 'cols'", 1127 issue=16102, 1128 useinstead="a more generic banded matrix constructor", 1129 deprecated_since_version="1.4" 1130 ).warn() 1131 1132 klass = kwargs.pop('cls', kls) 1133 rows = kwargs.pop('rows', None) 1134 cols = kwargs.pop('cols', None) 1135 1136 eigenval = kwargs.get('eigenval', None) 1137 if eigenvalue is None and eigenval is None: 1138 raise ValueError("Must supply an eigenvalue") 1139 elif eigenvalue != eigenval and None not in (eigenval, eigenvalue): 1140 raise ValueError( 1141 "Inconsistent values are given: 'eigenval'={}, " 1142 "'eigenvalue'={}".format(eigenval, eigenvalue)) 1143 else: 1144 if eigenval is not None: 1145 eigenvalue = eigenval 1146 1147 if (size, rows, cols) == (None, None, None): 1148 raise ValueError("Must supply a matrix size") 1149 1150 if size is not None: 1151 rows, cols = size, size 1152 elif rows is not None and cols is None: 1153 cols = rows 1154 elif cols is not None and rows is None: 1155 rows = cols 1156 1157 rows, cols = as_int(rows), as_int(cols) 1158 1159 return klass._eval_jordan_block(rows, cols, eigenvalue, band) 1160 1161 @classmethod 1162 def ones(kls, rows, cols=None, **kwargs): 1163 """Returns a matrix of ones. 1164 1165 Args 1166 ==== 1167 1168 rows : rows of the matrix 1169 cols : cols of the matrix (if None, cols=rows) 1170 1171 kwargs 1172 ====== 1173 cls : class of the returned matrix 1174 """ 1175 if cols is None: 1176 cols = rows 1177 klass = kwargs.get('cls', kls) 1178 rows, cols = as_int(rows), as_int(cols) 1179 1180 return klass._eval_ones(rows, cols) 1181 1182 @classmethod 1183 def zeros(kls, rows, cols=None, **kwargs): 1184 """Returns a matrix of zeros. 1185 1186 Args 1187 ==== 1188 1189 rows : rows of the matrix 1190 cols : cols of the matrix (if None, cols=rows) 1191 1192 kwargs 1193 ====== 1194 cls : class of the returned matrix 1195 """ 1196 if cols is None: 1197 cols = rows 1198 if rows < 0 or cols < 0: 1199 raise ValueError("Cannot create a {} x {} matrix. " 1200 "Both dimensions must be positive".format(rows, cols)) 1201 klass = kwargs.get('cls', kls) 1202 rows, cols = as_int(rows), as_int(cols) 1203 1204 return klass._eval_zeros(rows, cols) 1205 1206 @classmethod 1207 def companion(kls, poly): 1208 """Returns a companion matrix of a polynomial. 1209 1210 Examples 1211 ======== 1212 1213 >>> from sympy import Matrix, Poly, Symbol, symbols 1214 >>> x = Symbol('x') 1215 >>> c0, c1, c2, c3, c4 = symbols('c0:5') 1216 >>> p = Poly(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + x**5, x) 1217 >>> Matrix.companion(p) 1218 Matrix([ 1219 [0, 0, 0, 0, -c0], 1220 [1, 0, 0, 0, -c1], 1221 [0, 1, 0, 0, -c2], 1222 [0, 0, 1, 0, -c3], 1223 [0, 0, 0, 1, -c4]]) 1224 """ 1225 poly = kls._sympify(poly) 1226 if not isinstance(poly, Poly): 1227 raise ValueError("{} must be a Poly instance.".format(poly)) 1228 if not poly.is_monic: 1229 raise ValueError("{} must be a monic polynomial.".format(poly)) 1230 if not poly.is_univariate: 1231 raise ValueError( 1232 "{} must be a univariate polynomial.".format(poly)) 1233 1234 size = poly.degree() 1235 if not size >= 1: 1236 raise ValueError( 1237 "{} must have degree not less than 1.".format(poly)) 1238 1239 coeffs = poly.all_coeffs() 1240 def entry(i, j): 1241 if j == size - 1: 1242 return -coeffs[-1 - i] 1243 elif i == j + 1: 1244 return kls.one 1245 return kls.zero 1246 return kls._new(size, size, entry) 1247 1248 1249 @classmethod 1250 def wilkinson(kls, n, **kwargs): 1251 """Returns two square Wilkinson Matrix of size 2*n + 1 1252 $W_{2n + 1}^-, W_{2n + 1}^+ =$ Wilkinson(n) 1253 1254 Examples 1255 ======== 1256 1257 >>> from sympy.matrices import Matrix 1258 >>> wminus, wplus = Matrix.wilkinson(3) 1259 >>> wminus 1260 Matrix([ 1261 [-3, 1, 0, 0, 0, 0, 0], 1262 [ 1, -2, 1, 0, 0, 0, 0], 1263 [ 0, 1, -1, 1, 0, 0, 0], 1264 [ 0, 0, 1, 0, 1, 0, 0], 1265 [ 0, 0, 0, 1, 1, 1, 0], 1266 [ 0, 0, 0, 0, 1, 2, 1], 1267 [ 0, 0, 0, 0, 0, 1, 3]]) 1268 >>> wplus 1269 Matrix([ 1270 [3, 1, 0, 0, 0, 0, 0], 1271 [1, 2, 1, 0, 0, 0, 0], 1272 [0, 1, 1, 1, 0, 0, 0], 1273 [0, 0, 1, 0, 1, 0, 0], 1274 [0, 0, 0, 1, 1, 1, 0], 1275 [0, 0, 0, 0, 1, 2, 1], 1276 [0, 0, 0, 0, 0, 1, 3]]) 1277 1278 References 1279 ========== 1280 1281 .. [1] https://blogs.mathworks.com/cleve/2013/04/15/wilkinsons-matrices-2/ 1282 .. [2] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press, Oxford, 1965, 662 pp. 1283 1284 """ 1285 klass = kwargs.get('cls', kls) 1286 n = as_int(n) 1287 return klass._eval_wilkinson(n) 1288 1289class MatrixProperties(MatrixRequired): 1290 """Provides basic properties of a matrix.""" 1291 1292 def _eval_atoms(self, *types): 1293 result = set() 1294 for i in self: 1295 result.update(i.atoms(*types)) 1296 return result 1297 1298 def _eval_free_symbols(self): 1299 return set().union(*(i.free_symbols for i in self if i)) 1300 1301 def _eval_has(self, *patterns): 1302 return any(a.has(*patterns) for a in self) 1303 1304 def _eval_is_anti_symmetric(self, simpfunc): 1305 if not all(simpfunc(self[i, j] + self[j, i]).is_zero for i in range(self.rows) for j in range(self.cols)): 1306 return False 1307 return True 1308 1309 def _eval_is_diagonal(self): 1310 for i in range(self.rows): 1311 for j in range(self.cols): 1312 if i != j and self[i, j]: 1313 return False 1314 return True 1315 1316 # _eval_is_hermitian is called by some general sympy 1317 # routines and has a different *args signature. Make 1318 # sure the names don't clash by adding `_matrix_` in name. 1319 def _eval_is_matrix_hermitian(self, simpfunc): 1320 mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i].conjugate())) 1321 return mat.is_zero_matrix 1322 1323 def _eval_is_Identity(self) -> FuzzyBool: 1324 def dirac(i, j): 1325 if i == j: 1326 return 1 1327 return 0 1328 1329 return all(self[i, j] == dirac(i, j) 1330 for i in range(self.rows) 1331 for j in range(self.cols)) 1332 1333 def _eval_is_lower_hessenberg(self): 1334 return all(self[i, j].is_zero 1335 for i in range(self.rows) 1336 for j in range(i + 2, self.cols)) 1337 1338 def _eval_is_lower(self): 1339 return all(self[i, j].is_zero 1340 for i in range(self.rows) 1341 for j in range(i + 1, self.cols)) 1342 1343 def _eval_is_symbolic(self): 1344 return self.has(Symbol) 1345 1346 def _eval_is_symmetric(self, simpfunc): 1347 mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i])) 1348 return mat.is_zero_matrix 1349 1350 def _eval_is_zero_matrix(self): 1351 if any(i.is_zero == False for i in self): 1352 return False 1353 if any(i.is_zero is None for i in self): 1354 return None 1355 return True 1356 1357 def _eval_is_upper_hessenberg(self): 1358 return all(self[i, j].is_zero 1359 for i in range(2, self.rows) 1360 for j in range(min(self.cols, (i - 1)))) 1361 1362 def _eval_values(self): 1363 return [i for i in self if not i.is_zero] 1364 1365 def _has_positive_diagonals(self): 1366 diagonal_entries = (self[i, i] for i in range(self.rows)) 1367 return fuzzy_and(x.is_positive for x in diagonal_entries) 1368 1369 def _has_nonnegative_diagonals(self): 1370 diagonal_entries = (self[i, i] for i in range(self.rows)) 1371 return fuzzy_and(x.is_nonnegative for x in diagonal_entries) 1372 1373 def atoms(self, *types): 1374 """Returns the atoms that form the current object. 1375 1376 Examples 1377 ======== 1378 1379 >>> from sympy.abc import x, y 1380 >>> from sympy.matrices import Matrix 1381 >>> Matrix([[x]]) 1382 Matrix([[x]]) 1383 >>> _.atoms() 1384 {x} 1385 >>> Matrix([[x, y], [y, x]]) 1386 Matrix([ 1387 [x, y], 1388 [y, x]]) 1389 >>> _.atoms() 1390 {x, y} 1391 """ 1392 1393 types = tuple(t if isinstance(t, type) else type(t) for t in types) 1394 if not types: 1395 types = (Atom,) 1396 return self._eval_atoms(*types) 1397 1398 @property 1399 def free_symbols(self): 1400 """Returns the free symbols within the matrix. 1401 1402 Examples 1403 ======== 1404 1405 >>> from sympy.abc import x 1406 >>> from sympy.matrices import Matrix 1407 >>> Matrix([[x], [1]]).free_symbols 1408 {x} 1409 """ 1410 return self._eval_free_symbols() 1411 1412 def has(self, *patterns): 1413 """Test whether any subexpression matches any of the patterns. 1414 1415 Examples 1416 ======== 1417 1418 >>> from sympy import Matrix, SparseMatrix, Float 1419 >>> from sympy.abc import x, y 1420 >>> A = Matrix(((1, x), (0.2, 3))) 1421 >>> B = SparseMatrix(((1, x), (0.2, 3))) 1422 >>> A.has(x) 1423 True 1424 >>> A.has(y) 1425 False 1426 >>> A.has(Float) 1427 True 1428 >>> B.has(x) 1429 True 1430 >>> B.has(y) 1431 False 1432 >>> B.has(Float) 1433 True 1434 """ 1435 return self._eval_has(*patterns) 1436 1437 def is_anti_symmetric(self, simplify=True): 1438 """Check if matrix M is an antisymmetric matrix, 1439 that is, M is a square matrix with all M[i, j] == -M[j, i]. 1440 1441 When ``simplify=True`` (default), the sum M[i, j] + M[j, i] is 1442 simplified before testing to see if it is zero. By default, 1443 the SymPy simplify function is used. To use a custom function 1444 set simplify to a function that accepts a single argument which 1445 returns a simplified expression. To skip simplification, set 1446 simplify to False but note that although this will be faster, 1447 it may induce false negatives. 1448 1449 Examples 1450 ======== 1451 1452 >>> from sympy import Matrix, symbols 1453 >>> m = Matrix(2, 2, [0, 1, -1, 0]) 1454 >>> m 1455 Matrix([ 1456 [ 0, 1], 1457 [-1, 0]]) 1458 >>> m.is_anti_symmetric() 1459 True 1460 >>> x, y = symbols('x y') 1461 >>> m = Matrix(2, 3, [0, 0, x, -y, 0, 0]) 1462 >>> m 1463 Matrix([ 1464 [ 0, 0, x], 1465 [-y, 0, 0]]) 1466 >>> m.is_anti_symmetric() 1467 False 1468 1469 >>> from sympy.abc import x, y 1470 >>> m = Matrix(3, 3, [0, x**2 + 2*x + 1, y, 1471 ... -(x + 1)**2 , 0, x*y, 1472 ... -y, -x*y, 0]) 1473 1474 Simplification of matrix elements is done by default so even 1475 though two elements which should be equal and opposite wouldn't 1476 pass an equality test, the matrix is still reported as 1477 anti-symmetric: 1478 1479 >>> m[0, 1] == -m[1, 0] 1480 False 1481 >>> m.is_anti_symmetric() 1482 True 1483 1484 If 'simplify=False' is used for the case when a Matrix is already 1485 simplified, this will speed things up. Here, we see that without 1486 simplification the matrix does not appear anti-symmetric: 1487 1488 >>> m.is_anti_symmetric(simplify=False) 1489 False 1490 1491 But if the matrix were already expanded, then it would appear 1492 anti-symmetric and simplification in the is_anti_symmetric routine 1493 is not needed: 1494 1495 >>> m = m.expand() 1496 >>> m.is_anti_symmetric(simplify=False) 1497 True 1498 """ 1499 # accept custom simplification 1500 simpfunc = simplify 1501 if not isfunction(simplify): 1502 simpfunc = _simplify if simplify else lambda x: x 1503 1504 if not self.is_square: 1505 return False 1506 return self._eval_is_anti_symmetric(simpfunc) 1507 1508 def is_diagonal(self): 1509 """Check if matrix is diagonal, 1510 that is matrix in which the entries outside the main diagonal are all zero. 1511 1512 Examples 1513 ======== 1514 1515 >>> from sympy import Matrix, diag 1516 >>> m = Matrix(2, 2, [1, 0, 0, 2]) 1517 >>> m 1518 Matrix([ 1519 [1, 0], 1520 [0, 2]]) 1521 >>> m.is_diagonal() 1522 True 1523 1524 >>> m = Matrix(2, 2, [1, 1, 0, 2]) 1525 >>> m 1526 Matrix([ 1527 [1, 1], 1528 [0, 2]]) 1529 >>> m.is_diagonal() 1530 False 1531 1532 >>> m = diag(1, 2, 3) 1533 >>> m 1534 Matrix([ 1535 [1, 0, 0], 1536 [0, 2, 0], 1537 [0, 0, 3]]) 1538 >>> m.is_diagonal() 1539 True 1540 1541 See Also 1542 ======== 1543 1544 is_lower 1545 is_upper 1546 sympy.matrices.matrices.MatrixEigen.is_diagonalizable 1547 diagonalize 1548 """ 1549 return self._eval_is_diagonal() 1550 1551 @property 1552 def is_weakly_diagonally_dominant(self): 1553 r"""Tests if the matrix is row weakly diagonally dominant. 1554 1555 Explanation 1556 =========== 1557 1558 A $n, n$ matrix $A$ is row weakly diagonally dominant if 1559 1560 .. math:: 1561 \left|A_{i, i}\right| \ge \sum_{j = 0, j \neq i}^{n-1} 1562 \left|A_{i, j}\right| \quad {\text{for all }} 1563 i \in \{ 0, ..., n-1 \} 1564 1565 Examples 1566 ======== 1567 1568 >>> from sympy.matrices import Matrix 1569 >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]]) 1570 >>> A.is_weakly_diagonally_dominant 1571 True 1572 1573 >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]]) 1574 >>> A.is_weakly_diagonally_dominant 1575 False 1576 1577 >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]]) 1578 >>> A.is_weakly_diagonally_dominant 1579 True 1580 1581 Notes 1582 ===== 1583 1584 If you want to test whether a matrix is column diagonally 1585 dominant, you can apply the test after transposing the matrix. 1586 """ 1587 if not self.is_square: 1588 return False 1589 1590 rows, cols = self.shape 1591 1592 def test_row(i): 1593 summation = self.zero 1594 for j in range(cols): 1595 if i != j: 1596 summation += Abs(self[i, j]) 1597 return (Abs(self[i, i]) - summation).is_nonnegative 1598 1599 return fuzzy_and(test_row(i) for i in range(rows)) 1600 1601 @property 1602 def is_strongly_diagonally_dominant(self): 1603 r"""Tests if the matrix is row strongly diagonally dominant. 1604 1605 Explanation 1606 =========== 1607 1608 A $n, n$ matrix $A$ is row strongly diagonally dominant if 1609 1610 .. math:: 1611 \left|A_{i, i}\right| > \sum_{j = 0, j \neq i}^{n-1} 1612 \left|A_{i, j}\right| \quad {\text{for all }} 1613 i \in \{ 0, ..., n-1 \} 1614 1615 Examples 1616 ======== 1617 1618 >>> from sympy.matrices import Matrix 1619 >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]]) 1620 >>> A.is_strongly_diagonally_dominant 1621 False 1622 1623 >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]]) 1624 >>> A.is_strongly_diagonally_dominant 1625 False 1626 1627 >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]]) 1628 >>> A.is_strongly_diagonally_dominant 1629 True 1630 1631 Notes 1632 ===== 1633 1634 If you want to test whether a matrix is column diagonally 1635 dominant, you can apply the test after transposing the matrix. 1636 """ 1637 if not self.is_square: 1638 return False 1639 1640 rows, cols = self.shape 1641 1642 def test_row(i): 1643 summation = self.zero 1644 for j in range(cols): 1645 if i != j: 1646 summation += Abs(self[i, j]) 1647 return (Abs(self[i, i]) - summation).is_positive 1648 1649 return fuzzy_and(test_row(i) for i in range(rows)) 1650 1651 @property 1652 def is_hermitian(self): 1653 """Checks if the matrix is Hermitian. 1654 1655 In a Hermitian matrix element i,j is the complex conjugate of 1656 element j,i. 1657 1658 Examples 1659 ======== 1660 1661 >>> from sympy.matrices import Matrix 1662 >>> from sympy import I 1663 >>> from sympy.abc import x 1664 >>> a = Matrix([[1, I], [-I, 1]]) 1665 >>> a 1666 Matrix([ 1667 [ 1, I], 1668 [-I, 1]]) 1669 >>> a.is_hermitian 1670 True 1671 >>> a[0, 0] = 2*I 1672 >>> a.is_hermitian 1673 False 1674 >>> a[0, 0] = x 1675 >>> a.is_hermitian 1676 >>> a[0, 1] = a[1, 0]*I 1677 >>> a.is_hermitian 1678 False 1679 """ 1680 if not self.is_square: 1681 return False 1682 1683 return self._eval_is_matrix_hermitian(_simplify) 1684 1685 @property 1686 def is_Identity(self) -> FuzzyBool: 1687 if not self.is_square: 1688 return False 1689 return self._eval_is_Identity() 1690 1691 @property 1692 def is_lower_hessenberg(self): 1693 r"""Checks if the matrix is in the lower-Hessenberg form. 1694 1695 The lower hessenberg matrix has zero entries 1696 above the first superdiagonal. 1697 1698 Examples 1699 ======== 1700 1701 >>> from sympy.matrices import Matrix 1702 >>> a = Matrix([[1, 2, 0, 0], [5, 2, 3, 0], [3, 4, 3, 7], [5, 6, 1, 1]]) 1703 >>> a 1704 Matrix([ 1705 [1, 2, 0, 0], 1706 [5, 2, 3, 0], 1707 [3, 4, 3, 7], 1708 [5, 6, 1, 1]]) 1709 >>> a.is_lower_hessenberg 1710 True 1711 1712 See Also 1713 ======== 1714 1715 is_upper_hessenberg 1716 is_lower 1717 """ 1718 return self._eval_is_lower_hessenberg() 1719 1720 @property 1721 def is_lower(self): 1722 """Check if matrix is a lower triangular matrix. True can be returned 1723 even if the matrix is not square. 1724 1725 Examples 1726 ======== 1727 1728 >>> from sympy import Matrix 1729 >>> m = Matrix(2, 2, [1, 0, 0, 1]) 1730 >>> m 1731 Matrix([ 1732 [1, 0], 1733 [0, 1]]) 1734 >>> m.is_lower 1735 True 1736 1737 >>> m = Matrix(4, 3, [0, 0, 0, 2, 0, 0, 1, 4 , 0, 6, 6, 5]) 1738 >>> m 1739 Matrix([ 1740 [0, 0, 0], 1741 [2, 0, 0], 1742 [1, 4, 0], 1743 [6, 6, 5]]) 1744 >>> m.is_lower 1745 True 1746 1747 >>> from sympy.abc import x, y 1748 >>> m = Matrix(2, 2, [x**2 + y, y**2 + x, 0, x + y]) 1749 >>> m 1750 Matrix([ 1751 [x**2 + y, x + y**2], 1752 [ 0, x + y]]) 1753 >>> m.is_lower 1754 False 1755 1756 See Also 1757 ======== 1758 1759 is_upper 1760 is_diagonal 1761 is_lower_hessenberg 1762 """ 1763 return self._eval_is_lower() 1764 1765 @property 1766 def is_square(self): 1767 """Checks if a matrix is square. 1768 1769 A matrix is square if the number of rows equals the number of columns. 1770 The empty matrix is square by definition, since the number of rows and 1771 the number of columns are both zero. 1772 1773 Examples 1774 ======== 1775 1776 >>> from sympy import Matrix 1777 >>> a = Matrix([[1, 2, 3], [4, 5, 6]]) 1778 >>> b = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) 1779 >>> c = Matrix([]) 1780 >>> a.is_square 1781 False 1782 >>> b.is_square 1783 True 1784 >>> c.is_square 1785 True 1786 """ 1787 return self.rows == self.cols 1788 1789 def is_symbolic(self): 1790 """Checks if any elements contain Symbols. 1791 1792 Examples 1793 ======== 1794 1795 >>> from sympy.matrices import Matrix 1796 >>> from sympy.abc import x, y 1797 >>> M = Matrix([[x, y], [1, 0]]) 1798 >>> M.is_symbolic() 1799 True 1800 1801 """ 1802 return self._eval_is_symbolic() 1803 1804 def is_symmetric(self, simplify=True): 1805 """Check if matrix is symmetric matrix, 1806 that is square matrix and is equal to its transpose. 1807 1808 By default, simplifications occur before testing symmetry. 1809 They can be skipped using 'simplify=False'; while speeding things a bit, 1810 this may however induce false negatives. 1811 1812 Examples 1813 ======== 1814 1815 >>> from sympy import Matrix 1816 >>> m = Matrix(2, 2, [0, 1, 1, 2]) 1817 >>> m 1818 Matrix([ 1819 [0, 1], 1820 [1, 2]]) 1821 >>> m.is_symmetric() 1822 True 1823 1824 >>> m = Matrix(2, 2, [0, 1, 2, 0]) 1825 >>> m 1826 Matrix([ 1827 [0, 1], 1828 [2, 0]]) 1829 >>> m.is_symmetric() 1830 False 1831 1832 >>> m = Matrix(2, 3, [0, 0, 0, 0, 0, 0]) 1833 >>> m 1834 Matrix([ 1835 [0, 0, 0], 1836 [0, 0, 0]]) 1837 >>> m.is_symmetric() 1838 False 1839 1840 >>> from sympy.abc import x, y 1841 >>> m = Matrix(3, 3, [1, x**2 + 2*x + 1, y, (x + 1)**2 , 2, 0, y, 0, 3]) 1842 >>> m 1843 Matrix([ 1844 [ 1, x**2 + 2*x + 1, y], 1845 [(x + 1)**2, 2, 0], 1846 [ y, 0, 3]]) 1847 >>> m.is_symmetric() 1848 True 1849 1850 If the matrix is already simplified, you may speed-up is_symmetric() 1851 test by using 'simplify=False'. 1852 1853 >>> bool(m.is_symmetric(simplify=False)) 1854 False 1855 >>> m1 = m.expand() 1856 >>> m1.is_symmetric(simplify=False) 1857 True 1858 """ 1859 simpfunc = simplify 1860 if not isfunction(simplify): 1861 simpfunc = _simplify if simplify else lambda x: x 1862 1863 if not self.is_square: 1864 return False 1865 1866 return self._eval_is_symmetric(simpfunc) 1867 1868 @property 1869 def is_upper_hessenberg(self): 1870 """Checks if the matrix is the upper-Hessenberg form. 1871 1872 The upper hessenberg matrix has zero entries 1873 below the first subdiagonal. 1874 1875 Examples 1876 ======== 1877 1878 >>> from sympy.matrices import Matrix 1879 >>> a = Matrix([[1, 4, 2, 3], [3, 4, 1, 7], [0, 2, 3, 4], [0, 0, 1, 3]]) 1880 >>> a 1881 Matrix([ 1882 [1, 4, 2, 3], 1883 [3, 4, 1, 7], 1884 [0, 2, 3, 4], 1885 [0, 0, 1, 3]]) 1886 >>> a.is_upper_hessenberg 1887 True 1888 1889 See Also 1890 ======== 1891 1892 is_lower_hessenberg 1893 is_upper 1894 """ 1895 return self._eval_is_upper_hessenberg() 1896 1897 @property 1898 def is_upper(self): 1899 """Check if matrix is an upper triangular matrix. True can be returned 1900 even if the matrix is not square. 1901 1902 Examples 1903 ======== 1904 1905 >>> from sympy import Matrix 1906 >>> m = Matrix(2, 2, [1, 0, 0, 1]) 1907 >>> m 1908 Matrix([ 1909 [1, 0], 1910 [0, 1]]) 1911 >>> m.is_upper 1912 True 1913 1914 >>> m = Matrix(4, 3, [5, 1, 9, 0, 4 , 6, 0, 0, 5, 0, 0, 0]) 1915 >>> m 1916 Matrix([ 1917 [5, 1, 9], 1918 [0, 4, 6], 1919 [0, 0, 5], 1920 [0, 0, 0]]) 1921 >>> m.is_upper 1922 True 1923 1924 >>> m = Matrix(2, 3, [4, 2, 5, 6, 1, 1]) 1925 >>> m 1926 Matrix([ 1927 [4, 2, 5], 1928 [6, 1, 1]]) 1929 >>> m.is_upper 1930 False 1931 1932 See Also 1933 ======== 1934 1935 is_lower 1936 is_diagonal 1937 is_upper_hessenberg 1938 """ 1939 return all(self[i, j].is_zero 1940 for i in range(1, self.rows) 1941 for j in range(min(i, self.cols))) 1942 1943 @property 1944 def is_zero_matrix(self): 1945 """Checks if a matrix is a zero matrix. 1946 1947 A matrix is zero if every element is zero. A matrix need not be square 1948 to be considered zero. The empty matrix is zero by the principle of 1949 vacuous truth. For a matrix that may or may not be zero (e.g. 1950 contains a symbol), this will be None 1951 1952 Examples 1953 ======== 1954 1955 >>> from sympy import Matrix, zeros 1956 >>> from sympy.abc import x 1957 >>> a = Matrix([[0, 0], [0, 0]]) 1958 >>> b = zeros(3, 4) 1959 >>> c = Matrix([[0, 1], [0, 0]]) 1960 >>> d = Matrix([]) 1961 >>> e = Matrix([[x, 0], [0, 0]]) 1962 >>> a.is_zero_matrix 1963 True 1964 >>> b.is_zero_matrix 1965 True 1966 >>> c.is_zero_matrix 1967 False 1968 >>> d.is_zero_matrix 1969 True 1970 >>> e.is_zero_matrix 1971 """ 1972 return self._eval_is_zero_matrix() 1973 1974 def values(self): 1975 """Return non-zero values of self.""" 1976 return self._eval_values() 1977 1978 1979class MatrixOperations(MatrixRequired): 1980 """Provides basic matrix shape and elementwise 1981 operations. Should not be instantiated directly.""" 1982 1983 def _eval_adjoint(self): 1984 return self.transpose().conjugate() 1985 1986 def _eval_applyfunc(self, f): 1987 out = self._new(self.rows, self.cols, [f(x) for x in self]) 1988 return out 1989 1990 def _eval_as_real_imag(self): # type: ignore 1991 from sympy.functions.elementary.complexes import re, im 1992 1993 return (self.applyfunc(re), self.applyfunc(im)) 1994 1995 def _eval_conjugate(self): 1996 return self.applyfunc(lambda x: x.conjugate()) 1997 1998 def _eval_permute_cols(self, perm): 1999 # apply the permutation to a list 2000 mapping = list(perm) 2001 2002 def entry(i, j): 2003 return self[i, mapping[j]] 2004 2005 return self._new(self.rows, self.cols, entry) 2006 2007 def _eval_permute_rows(self, perm): 2008 # apply the permutation to a list 2009 mapping = list(perm) 2010 2011 def entry(i, j): 2012 return self[mapping[i], j] 2013 2014 return self._new(self.rows, self.cols, entry) 2015 2016 def _eval_trace(self): 2017 return sum(self[i, i] for i in range(self.rows)) 2018 2019 def _eval_transpose(self): 2020 return self._new(self.cols, self.rows, lambda i, j: self[j, i]) 2021 2022 def adjoint(self): 2023 """Conjugate transpose or Hermitian conjugation.""" 2024 return self._eval_adjoint() 2025 2026 def applyfunc(self, f): 2027 """Apply a function to each element of the matrix. 2028 2029 Examples 2030 ======== 2031 2032 >>> from sympy import Matrix 2033 >>> m = Matrix(2, 2, lambda i, j: i*2+j) 2034 >>> m 2035 Matrix([ 2036 [0, 1], 2037 [2, 3]]) 2038 >>> m.applyfunc(lambda i: 2*i) 2039 Matrix([ 2040 [0, 2], 2041 [4, 6]]) 2042 2043 """ 2044 if not callable(f): 2045 raise TypeError("`f` must be callable.") 2046 2047 return self._eval_applyfunc(f) 2048 2049 def as_real_imag(self, deep=True, **hints): 2050 """Returns a tuple containing the (real, imaginary) part of matrix.""" 2051 # XXX: Ignoring deep and hints... 2052 return self._eval_as_real_imag() 2053 2054 def conjugate(self): 2055 """Return the by-element conjugation. 2056 2057 Examples 2058 ======== 2059 2060 >>> from sympy.matrices import SparseMatrix 2061 >>> from sympy import I 2062 >>> a = SparseMatrix(((1, 2 + I), (3, 4), (I, -I))) 2063 >>> a 2064 Matrix([ 2065 [1, 2 + I], 2066 [3, 4], 2067 [I, -I]]) 2068 >>> a.C 2069 Matrix([ 2070 [ 1, 2 - I], 2071 [ 3, 4], 2072 [-I, I]]) 2073 2074 See Also 2075 ======== 2076 2077 transpose: Matrix transposition 2078 H: Hermite conjugation 2079 sympy.matrices.matrices.MatrixBase.D: Dirac conjugation 2080 """ 2081 return self._eval_conjugate() 2082 2083 def doit(self, **kwargs): 2084 return self.applyfunc(lambda x: x.doit()) 2085 2086 def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None, verbose=False): 2087 """Apply evalf() to each element of self.""" 2088 options = {'subs':subs, 'maxn':maxn, 'chop':chop, 'strict':strict, 2089 'quad':quad, 'verbose':verbose} 2090 return self.applyfunc(lambda i: i.evalf(n, **options)) 2091 2092 def expand(self, deep=True, modulus=None, power_base=True, power_exp=True, 2093 mul=True, log=True, multinomial=True, basic=True, **hints): 2094 """Apply core.function.expand to each entry of the matrix. 2095 2096 Examples 2097 ======== 2098 2099 >>> from sympy.abc import x 2100 >>> from sympy.matrices import Matrix 2101 >>> Matrix(1, 1, [x*(x+1)]) 2102 Matrix([[x*(x + 1)]]) 2103 >>> _.expand() 2104 Matrix([[x**2 + x]]) 2105 2106 """ 2107 return self.applyfunc(lambda x: x.expand( 2108 deep, modulus, power_base, power_exp, mul, log, multinomial, basic, 2109 **hints)) 2110 2111 @property 2112 def H(self): 2113 """Return Hermite conjugate. 2114 2115 Examples 2116 ======== 2117 2118 >>> from sympy import Matrix, I 2119 >>> m = Matrix((0, 1 + I, 2, 3)) 2120 >>> m 2121 Matrix([ 2122 [ 0], 2123 [1 + I], 2124 [ 2], 2125 [ 3]]) 2126 >>> m.H 2127 Matrix([[0, 1 - I, 2, 3]]) 2128 2129 See Also 2130 ======== 2131 2132 conjugate: By-element conjugation 2133 sympy.matrices.matrices.MatrixBase.D: Dirac conjugation 2134 """ 2135 return self.T.C 2136 2137 def permute(self, perm, orientation='rows', direction='forward'): 2138 r"""Permute the rows or columns of a matrix by the given list of 2139 swaps. 2140 2141 Parameters 2142 ========== 2143 2144 perm : Permutation, list, or list of lists 2145 A representation for the permutation. 2146 2147 If it is ``Permutation``, it is used directly with some 2148 resizing with respect to the matrix size. 2149 2150 If it is specified as list of lists, 2151 (e.g., ``[[0, 1], [0, 2]]``), then the permutation is formed 2152 from applying the product of cycles. The direction how the 2153 cyclic product is applied is described in below. 2154 2155 If it is specified as a list, the list should represent 2156 an array form of a permutation. (e.g., ``[1, 2, 0]``) which 2157 would would form the swapping function 2158 `0 \mapsto 1, 1 \mapsto 2, 2\mapsto 0`. 2159 2160 orientation : 'rows', 'cols' 2161 A flag to control whether to permute the rows or the columns 2162 2163 direction : 'forward', 'backward' 2164 A flag to control whether to apply the permutations from 2165 the start of the list first, or from the back of the list 2166 first. 2167 2168 For example, if the permutation specification is 2169 ``[[0, 1], [0, 2]]``, 2170 2171 If the flag is set to ``'forward'``, the cycle would be 2172 formed as `0 \mapsto 2, 2 \mapsto 1, 1 \mapsto 0`. 2173 2174 If the flag is set to ``'backward'``, the cycle would be 2175 formed as `0 \mapsto 1, 1 \mapsto 2, 2 \mapsto 0`. 2176 2177 If the argument ``perm`` is not in a form of list of lists, 2178 this flag takes no effect. 2179 2180 Examples 2181 ======== 2182 2183 >>> from sympy.matrices import eye 2184 >>> M = eye(3) 2185 >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='forward') 2186 Matrix([ 2187 [0, 0, 1], 2188 [1, 0, 0], 2189 [0, 1, 0]]) 2190 2191 >>> from sympy.matrices import eye 2192 >>> M = eye(3) 2193 >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='backward') 2194 Matrix([ 2195 [0, 1, 0], 2196 [0, 0, 1], 2197 [1, 0, 0]]) 2198 2199 Notes 2200 ===== 2201 2202 If a bijective function 2203 `\sigma : \mathbb{N}_0 \rightarrow \mathbb{N}_0` denotes the 2204 permutation. 2205 2206 If the matrix `A` is the matrix to permute, represented as 2207 a horizontal or a vertical stack of vectors: 2208 2209 .. math:: 2210 A = 2211 \begin{bmatrix} 2212 a_0 \\ a_1 \\ \vdots \\ a_{n-1} 2213 \end{bmatrix} = 2214 \begin{bmatrix} 2215 \alpha_0 & \alpha_1 & \cdots & \alpha_{n-1} 2216 \end{bmatrix} 2217 2218 If the matrix `B` is the result, the permutation of matrix rows 2219 is defined as: 2220 2221 .. math:: 2222 B := \begin{bmatrix} 2223 a_{\sigma(0)} \\ a_{\sigma(1)} \\ \vdots \\ a_{\sigma(n-1)} 2224 \end{bmatrix} 2225 2226 And the permutation of matrix columns is defined as: 2227 2228 .. math:: 2229 B := \begin{bmatrix} 2230 \alpha_{\sigma(0)} & \alpha_{\sigma(1)} & 2231 \cdots & \alpha_{\sigma(n-1)} 2232 \end{bmatrix} 2233 """ 2234 from sympy.combinatorics import Permutation 2235 2236 # allow british variants and `columns` 2237 if direction == 'forwards': 2238 direction = 'forward' 2239 if direction == 'backwards': 2240 direction = 'backward' 2241 if orientation == 'columns': 2242 orientation = 'cols' 2243 2244 if direction not in ('forward', 'backward'): 2245 raise TypeError("direction='{}' is an invalid kwarg. " 2246 "Try 'forward' or 'backward'".format(direction)) 2247 if orientation not in ('rows', 'cols'): 2248 raise TypeError("orientation='{}' is an invalid kwarg. " 2249 "Try 'rows' or 'cols'".format(orientation)) 2250 2251 if not isinstance(perm, (Permutation, Iterable)): 2252 raise ValueError( 2253 "{} must be a list, a list of lists, " 2254 "or a SymPy permutation object.".format(perm)) 2255 2256 # ensure all swaps are in range 2257 max_index = self.rows if orientation == 'rows' else self.cols 2258 if not all(0 <= t <= max_index for t in flatten(list(perm))): 2259 raise IndexError("`swap` indices out of range.") 2260 2261 if perm and not isinstance(perm, Permutation) and \ 2262 isinstance(perm[0], Iterable): 2263 if direction == 'forward': 2264 perm = list(reversed(perm)) 2265 perm = Permutation(perm, size=max_index+1) 2266 else: 2267 perm = Permutation(perm, size=max_index+1) 2268 2269 if orientation == 'rows': 2270 return self._eval_permute_rows(perm) 2271 if orientation == 'cols': 2272 return self._eval_permute_cols(perm) 2273 2274 def permute_cols(self, swaps, direction='forward'): 2275 """Alias for 2276 ``self.permute(swaps, orientation='cols', direction=direction)`` 2277 2278 See Also 2279 ======== 2280 2281 permute 2282 """ 2283 return self.permute(swaps, orientation='cols', direction=direction) 2284 2285 def permute_rows(self, swaps, direction='forward'): 2286 """Alias for 2287 ``self.permute(swaps, orientation='rows', direction=direction)`` 2288 2289 See Also 2290 ======== 2291 2292 permute 2293 """ 2294 return self.permute(swaps, orientation='rows', direction=direction) 2295 2296 def refine(self, assumptions=True): 2297 """Apply refine to each element of the matrix. 2298 2299 Examples 2300 ======== 2301 2302 >>> from sympy import Symbol, Matrix, Abs, sqrt, Q 2303 >>> x = Symbol('x') 2304 >>> Matrix([[Abs(x)**2, sqrt(x**2)],[sqrt(x**2), Abs(x)**2]]) 2305 Matrix([ 2306 [ Abs(x)**2, sqrt(x**2)], 2307 [sqrt(x**2), Abs(x)**2]]) 2308 >>> _.refine(Q.real(x)) 2309 Matrix([ 2310 [ x**2, Abs(x)], 2311 [Abs(x), x**2]]) 2312 2313 """ 2314 return self.applyfunc(lambda x: refine(x, assumptions)) 2315 2316 def replace(self, F, G, map=False, simultaneous=True, exact=None): 2317 """Replaces Function F in Matrix entries with Function G. 2318 2319 Examples 2320 ======== 2321 2322 >>> from sympy import symbols, Function, Matrix 2323 >>> F, G = symbols('F, G', cls=Function) 2324 >>> M = Matrix(2, 2, lambda i, j: F(i+j)) ; M 2325 Matrix([ 2326 [F(0), F(1)], 2327 [F(1), F(2)]]) 2328 >>> N = M.replace(F,G) 2329 >>> N 2330 Matrix([ 2331 [G(0), G(1)], 2332 [G(1), G(2)]]) 2333 """ 2334 return self.applyfunc( 2335 lambda x: x.replace(F, G, map=map, simultaneous=simultaneous, exact=exact)) 2336 2337 def rot90(self, k=1): 2338 """Rotates Matrix by 90 degrees 2339 2340 Parameters 2341 ========== 2342 2343 k : int 2344 Specifies how many times the matrix is rotated by 90 degrees 2345 (clockwise when positive, counter-clockwise when negative). 2346 2347 Examples 2348 ======== 2349 2350 >>> from sympy import Matrix, symbols 2351 >>> A = Matrix(2, 2, symbols('a:d')) 2352 >>> A 2353 Matrix([ 2354 [a, b], 2355 [c, d]]) 2356 2357 Rotating the matrix clockwise one time: 2358 2359 >>> A.rot90(1) 2360 Matrix([ 2361 [c, a], 2362 [d, b]]) 2363 2364 Rotating the matrix anticlockwise two times: 2365 2366 >>> A.rot90(-2) 2367 Matrix([ 2368 [d, c], 2369 [b, a]]) 2370 """ 2371 2372 mod = k%4 2373 if mod == 0: 2374 return self 2375 if mod == 1: 2376 return self[::-1, ::].T 2377 if mod == 2: 2378 return self[::-1, ::-1] 2379 if mod == 3: 2380 return self[::, ::-1].T 2381 2382 def simplify(self, **kwargs): 2383 """Apply simplify to each element of the matrix. 2384 2385 Examples 2386 ======== 2387 2388 >>> from sympy.abc import x, y 2389 >>> from sympy import sin, cos 2390 >>> from sympy.matrices import SparseMatrix 2391 >>> SparseMatrix(1, 1, [x*sin(y)**2 + x*cos(y)**2]) 2392 Matrix([[x*sin(y)**2 + x*cos(y)**2]]) 2393 >>> _.simplify() 2394 Matrix([[x]]) 2395 """ 2396 return self.applyfunc(lambda x: x.simplify(**kwargs)) 2397 2398 def subs(self, *args, **kwargs): # should mirror core.basic.subs 2399 """Return a new matrix with subs applied to each entry. 2400 2401 Examples 2402 ======== 2403 2404 >>> from sympy.abc import x, y 2405 >>> from sympy.matrices import SparseMatrix, Matrix 2406 >>> SparseMatrix(1, 1, [x]) 2407 Matrix([[x]]) 2408 >>> _.subs(x, y) 2409 Matrix([[y]]) 2410 >>> Matrix(_).subs(y, x) 2411 Matrix([[x]]) 2412 """ 2413 2414 if len(args) == 1 and not isinstance(args[0], (dict, set)) and iter(args[0]) and not is_sequence(args[0]): 2415 args = (list(args[0]),) 2416 2417 return self.applyfunc(lambda x: x.subs(*args, **kwargs)) 2418 2419 def trace(self): 2420 """ 2421 Returns the trace of a square matrix i.e. the sum of the 2422 diagonal elements. 2423 2424 Examples 2425 ======== 2426 2427 >>> from sympy import Matrix 2428 >>> A = Matrix(2, 2, [1, 2, 3, 4]) 2429 >>> A.trace() 2430 5 2431 2432 """ 2433 if self.rows != self.cols: 2434 raise NonSquareMatrixError() 2435 return self._eval_trace() 2436 2437 def transpose(self): 2438 """ 2439 Returns the transpose of the matrix. 2440 2441 Examples 2442 ======== 2443 2444 >>> from sympy import Matrix 2445 >>> A = Matrix(2, 2, [1, 2, 3, 4]) 2446 >>> A.transpose() 2447 Matrix([ 2448 [1, 3], 2449 [2, 4]]) 2450 2451 >>> from sympy import Matrix, I 2452 >>> m=Matrix(((1, 2+I), (3, 4))) 2453 >>> m 2454 Matrix([ 2455 [1, 2 + I], 2456 [3, 4]]) 2457 >>> m.transpose() 2458 Matrix([ 2459 [ 1, 3], 2460 [2 + I, 4]]) 2461 >>> m.T == m.transpose() 2462 True 2463 2464 See Also 2465 ======== 2466 2467 conjugate: By-element conjugation 2468 2469 """ 2470 return self._eval_transpose() 2471 2472 @property 2473 def T(self): 2474 '''Matrix transposition''' 2475 return self.transpose() 2476 2477 @property 2478 def C(self): 2479 '''By-element conjugation''' 2480 return self.conjugate() 2481 2482 def n(self, *args, **kwargs): 2483 """Apply evalf() to each element of self.""" 2484 return self.evalf(*args, **kwargs) 2485 2486 def xreplace(self, rule): # should mirror core.basic.xreplace 2487 """Return a new matrix with xreplace applied to each entry. 2488 2489 Examples 2490 ======== 2491 2492 >>> from sympy.abc import x, y 2493 >>> from sympy.matrices import SparseMatrix, Matrix 2494 >>> SparseMatrix(1, 1, [x]) 2495 Matrix([[x]]) 2496 >>> _.xreplace({x: y}) 2497 Matrix([[y]]) 2498 >>> Matrix(_).xreplace({y: x}) 2499 Matrix([[x]]) 2500 """ 2501 return self.applyfunc(lambda x: x.xreplace(rule)) 2502 2503 def _eval_simplify(self, **kwargs): 2504 # XXX: We can't use self.simplify here as mutable subclasses will 2505 # override simplify and have it return None 2506 return MatrixOperations.simplify(self, **kwargs) 2507 2508 def _eval_trigsimp(self, **opts): 2509 from sympy.simplify import trigsimp 2510 return self.applyfunc(lambda x: trigsimp(x, **opts)) 2511 2512 def upper_triangular(self, k=0): 2513 """returns the elements on and above the kth diagonal of a matrix. 2514 If k is not specified then simply returns upper-triangular portion 2515 of a matrix 2516 2517 Examples 2518 ======== 2519 2520 >>> from sympy import ones 2521 >>> A = ones(4) 2522 >>> A.upper_triangular() 2523 Matrix([ 2524 [1, 1, 1, 1], 2525 [0, 1, 1, 1], 2526 [0, 0, 1, 1], 2527 [0, 0, 0, 1]]) 2528 2529 >>> A.upper_triangular(2) 2530 Matrix([ 2531 [0, 0, 1, 1], 2532 [0, 0, 0, 1], 2533 [0, 0, 0, 0], 2534 [0, 0, 0, 0]]) 2535 2536 >>> A.upper_triangular(-1) 2537 Matrix([ 2538 [1, 1, 1, 1], 2539 [1, 1, 1, 1], 2540 [0, 1, 1, 1], 2541 [0, 0, 1, 1]]) 2542 2543 """ 2544 2545 def entry(i, j): 2546 return self[i, j] if i + k <= j else self.zero 2547 2548 return self._new(self.rows, self.cols, entry) 2549 2550 2551 def lower_triangular(self, k=0): 2552 """returns the elements on and below the kth diagonal of a matrix. 2553 If k is not specified then simply returns lower-triangular portion 2554 of a matrix 2555 2556 Examples 2557 ======== 2558 2559 >>> from sympy import ones 2560 >>> A = ones(4) 2561 >>> A.lower_triangular() 2562 Matrix([ 2563 [1, 0, 0, 0], 2564 [1, 1, 0, 0], 2565 [1, 1, 1, 0], 2566 [1, 1, 1, 1]]) 2567 2568 >>> A.lower_triangular(-2) 2569 Matrix([ 2570 [0, 0, 0, 0], 2571 [0, 0, 0, 0], 2572 [1, 0, 0, 0], 2573 [1, 1, 0, 0]]) 2574 2575 >>> A.lower_triangular(1) 2576 Matrix([ 2577 [1, 1, 0, 0], 2578 [1, 1, 1, 0], 2579 [1, 1, 1, 1], 2580 [1, 1, 1, 1]]) 2581 2582 """ 2583 2584 def entry(i, j): 2585 return self[i, j] if i + k >= j else self.zero 2586 2587 return self._new(self.rows, self.cols, entry) 2588 2589 2590 2591class MatrixArithmetic(MatrixRequired): 2592 """Provides basic matrix arithmetic operations. 2593 Should not be instantiated directly.""" 2594 2595 _op_priority = 10.01 2596 2597 def _eval_Abs(self): 2598 return self._new(self.rows, self.cols, lambda i, j: Abs(self[i, j])) 2599 2600 def _eval_add(self, other): 2601 return self._new(self.rows, self.cols, 2602 lambda i, j: self[i, j] + other[i, j]) 2603 2604 def _eval_matrix_mul(self, other): 2605 def entry(i, j): 2606 vec = [self[i,k]*other[k,j] for k in range(self.cols)] 2607 try: 2608 return Add(*vec) 2609 except (TypeError, SympifyError): 2610 # Some matrices don't work with `sum` or `Add` 2611 # They don't work with `sum` because `sum` tries to add `0` 2612 # Fall back to a safe way to multiply if the `Add` fails. 2613 return reduce(lambda a, b: a + b, vec) 2614 2615 return self._new(self.rows, other.cols, entry) 2616 2617 def _eval_matrix_mul_elementwise(self, other): 2618 return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other[i,j]) 2619 2620 def _eval_matrix_rmul(self, other): 2621 def entry(i, j): 2622 return sum(other[i,k]*self[k,j] for k in range(other.cols)) 2623 return self._new(other.rows, self.cols, entry) 2624 2625 def _eval_pow_by_recursion(self, num): 2626 if num == 1: 2627 return self 2628 2629 if num % 2 == 1: 2630 a, b = self, self._eval_pow_by_recursion(num - 1) 2631 else: 2632 a = b = self._eval_pow_by_recursion(num // 2) 2633 2634 return a.multiply(b) 2635 2636 def _eval_pow_by_cayley(self, exp): 2637 from sympy.discrete.recurrences import linrec_coeffs 2638 row = self.shape[0] 2639 p = self.charpoly() 2640 2641 coeffs = (-p).all_coeffs()[1:] 2642 coeffs = linrec_coeffs(coeffs, exp) 2643 new_mat = self.eye(row) 2644 ans = self.zeros(row) 2645 2646 for i in range(row): 2647 ans += coeffs[i]*new_mat 2648 new_mat *= self 2649 2650 return ans 2651 2652 def _eval_pow_by_recursion_dotprodsimp(self, num, prevsimp=None): 2653 if prevsimp is None: 2654 prevsimp = [True]*len(self) 2655 2656 if num == 1: 2657 return self 2658 2659 if num % 2 == 1: 2660 a, b = self, self._eval_pow_by_recursion_dotprodsimp(num - 1, 2661 prevsimp=prevsimp) 2662 else: 2663 a = b = self._eval_pow_by_recursion_dotprodsimp(num // 2, 2664 prevsimp=prevsimp) 2665 2666 m = a.multiply(b, dotprodsimp=False) 2667 lenm = len(m) 2668 elems = [None]*lenm 2669 2670 for i in range(lenm): 2671 if prevsimp[i]: 2672 elems[i], prevsimp[i] = _dotprodsimp(m[i], withsimp=True) 2673 else: 2674 elems[i] = m[i] 2675 2676 return m._new(m.rows, m.cols, elems) 2677 2678 def _eval_scalar_mul(self, other): 2679 return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other) 2680 2681 def _eval_scalar_rmul(self, other): 2682 return self._new(self.rows, self.cols, lambda i, j: other*self[i,j]) 2683 2684 def _eval_Mod(self, other): 2685 from sympy import Mod 2686 return self._new(self.rows, self.cols, lambda i, j: Mod(self[i, j], other)) 2687 2688 # python arithmetic functions 2689 def __abs__(self): 2690 """Returns a new matrix with entry-wise absolute values.""" 2691 return self._eval_Abs() 2692 2693 @call_highest_priority('__radd__') 2694 def __add__(self, other): 2695 """Return self + other, raising ShapeError if shapes don't match.""" 2696 if isinstance(other, NDimArray): # Matrix and array addition is currently not implemented 2697 return NotImplemented 2698 other = _matrixify(other) 2699 # matrix-like objects can have shapes. This is 2700 # our first sanity check. 2701 if hasattr(other, 'shape'): 2702 if self.shape != other.shape: 2703 raise ShapeError("Matrix size mismatch: %s + %s" % ( 2704 self.shape, other.shape)) 2705 2706 # honest sympy matrices defer to their class's routine 2707 if getattr(other, 'is_Matrix', False): 2708 # call the highest-priority class's _eval_add 2709 a, b = self, other 2710 if a.__class__ != classof(a, b): 2711 b, a = a, b 2712 return a._eval_add(b) 2713 # Matrix-like objects can be passed to CommonMatrix routines directly. 2714 if getattr(other, 'is_MatrixLike', False): 2715 return MatrixArithmetic._eval_add(self, other) 2716 2717 raise TypeError('cannot add %s and %s' % (type(self), type(other))) 2718 2719 @call_highest_priority('__rtruediv__') 2720 def __truediv__(self, other): 2721 return self * (self.one / other) 2722 2723 @call_highest_priority('__rmatmul__') 2724 def __matmul__(self, other): 2725 other = _matrixify(other) 2726 if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False): 2727 return NotImplemented 2728 2729 return self.__mul__(other) 2730 2731 def __mod__(self, other): 2732 return self.applyfunc(lambda x: x % other) 2733 2734 @call_highest_priority('__rmul__') 2735 def __mul__(self, other): 2736 """Return self*other where other is either a scalar or a matrix 2737 of compatible dimensions. 2738 2739 Examples 2740 ======== 2741 2742 >>> from sympy.matrices import Matrix 2743 >>> A = Matrix([[1, 2, 3], [4, 5, 6]]) 2744 >>> 2*A == A*2 == Matrix([[2, 4, 6], [8, 10, 12]]) 2745 True 2746 >>> B = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) 2747 >>> A*B 2748 Matrix([ 2749 [30, 36, 42], 2750 [66, 81, 96]]) 2751 >>> B*A 2752 Traceback (most recent call last): 2753 ... 2754 ShapeError: Matrices size mismatch. 2755 >>> 2756 2757 See Also 2758 ======== 2759 2760 matrix_multiply_elementwise 2761 """ 2762 2763 return self.multiply(other) 2764 2765 def multiply(self, other, dotprodsimp=None): 2766 """Same as __mul__() but with optional simplification. 2767 2768 Parameters 2769 ========== 2770 2771 dotprodsimp : bool, optional 2772 Specifies whether intermediate term algebraic simplification is used 2773 during matrix multiplications to control expression blowup and thus 2774 speed up calculation. Default is off. 2775 """ 2776 2777 isimpbool = _get_intermediate_simp_bool(False, dotprodsimp) 2778 other = _matrixify(other) 2779 # matrix-like objects can have shapes. This is 2780 # our first sanity check. Double check other is not explicitly not a Matrix. 2781 if (hasattr(other, 'shape') and len(other.shape) == 2 and 2782 (getattr(other, 'is_Matrix', True) or 2783 getattr(other, 'is_MatrixLike', True))): 2784 if self.shape[1] != other.shape[0]: 2785 raise ShapeError("Matrix size mismatch: %s * %s." % ( 2786 self.shape, other.shape)) 2787 2788 # honest sympy matrices defer to their class's routine 2789 if getattr(other, 'is_Matrix', False): 2790 m = self._eval_matrix_mul(other) 2791 if isimpbool: 2792 return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m]) 2793 return m 2794 2795 # Matrix-like objects can be passed to CommonMatrix routines directly. 2796 if getattr(other, 'is_MatrixLike', False): 2797 return MatrixArithmetic._eval_matrix_mul(self, other) 2798 2799 # if 'other' is not iterable then scalar multiplication. 2800 if not isinstance(other, Iterable): 2801 try: 2802 return self._eval_scalar_mul(other) 2803 except TypeError: 2804 pass 2805 2806 return NotImplemented 2807 2808 def multiply_elementwise(self, other): 2809 """Return the Hadamard product (elementwise product) of A and B 2810 2811 Examples 2812 ======== 2813 2814 >>> from sympy.matrices import Matrix 2815 >>> A = Matrix([[0, 1, 2], [3, 4, 5]]) 2816 >>> B = Matrix([[1, 10, 100], [100, 10, 1]]) 2817 >>> A.multiply_elementwise(B) 2818 Matrix([ 2819 [ 0, 10, 200], 2820 [300, 40, 5]]) 2821 2822 See Also 2823 ======== 2824 2825 sympy.matrices.matrices.MatrixBase.cross 2826 sympy.matrices.matrices.MatrixBase.dot 2827 multiply 2828 """ 2829 if self.shape != other.shape: 2830 raise ShapeError("Matrix shapes must agree {} != {}".format(self.shape, other.shape)) 2831 2832 return self._eval_matrix_mul_elementwise(other) 2833 2834 def __neg__(self): 2835 return self._eval_scalar_mul(-1) 2836 2837 @call_highest_priority('__rpow__') 2838 def __pow__(self, exp): 2839 """Return self**exp a scalar or symbol.""" 2840 2841 return self.pow(exp) 2842 2843 2844 def pow(self, exp, method=None): 2845 r"""Return self**exp a scalar or symbol. 2846 2847 Parameters 2848 ========== 2849 2850 method : multiply, mulsimp, jordan, cayley 2851 If multiply then it returns exponentiation using recursion. 2852 If jordan then Jordan form exponentiation will be used. 2853 If cayley then the exponentiation is done using Cayley-Hamilton 2854 theorem. 2855 If mulsimp then the exponentiation is done using recursion 2856 with dotprodsimp. This specifies whether intermediate term 2857 algebraic simplification is used during naive matrix power to 2858 control expression blowup and thus speed up calculation. 2859 If None, then it heuristically decides which method to use. 2860 2861 """ 2862 2863 if method is not None and method not in ['multiply', 'mulsimp', 'jordan', 'cayley']: 2864 raise TypeError('No such method') 2865 if self.rows != self.cols: 2866 raise NonSquareMatrixError() 2867 a = self 2868 jordan_pow = getattr(a, '_matrix_pow_by_jordan_blocks', None) 2869 exp = sympify(exp) 2870 2871 if exp.is_zero: 2872 return a._new(a.rows, a.cols, lambda i, j: int(i == j)) 2873 if exp == 1: 2874 return a 2875 2876 diagonal = getattr(a, 'is_diagonal', None) 2877 if diagonal is not None and diagonal(): 2878 return a._new(a.rows, a.cols, lambda i, j: a[i,j]**exp if i == j else 0) 2879 2880 if exp.is_Number and exp % 1 == 0: 2881 if a.rows == 1: 2882 return a._new([[a[0]**exp]]) 2883 if exp < 0: 2884 exp = -exp 2885 a = a.inv() 2886 # When certain conditions are met, 2887 # Jordan block algorithm is faster than 2888 # computation by recursion. 2889 if method == 'jordan': 2890 try: 2891 return jordan_pow(exp) 2892 except MatrixError: 2893 if method == 'jordan': 2894 raise 2895 2896 elif method == 'cayley': 2897 if not exp.is_Number or exp % 1 != 0: 2898 raise ValueError("cayley method is only valid for integer powers") 2899 return a._eval_pow_by_cayley(exp) 2900 2901 elif method == "mulsimp": 2902 if not exp.is_Number or exp % 1 != 0: 2903 raise ValueError("mulsimp method is only valid for integer powers") 2904 return a._eval_pow_by_recursion_dotprodsimp(exp) 2905 2906 elif method == "multiply": 2907 if not exp.is_Number or exp % 1 != 0: 2908 raise ValueError("multiply method is only valid for integer powers") 2909 return a._eval_pow_by_recursion(exp) 2910 2911 elif method is None and exp.is_Number and exp % 1 == 0: 2912 # Decide heuristically which method to apply 2913 if a.rows == 2 and exp > 100000: 2914 return jordan_pow(exp) 2915 elif _get_intermediate_simp_bool(True, None): 2916 return a._eval_pow_by_recursion_dotprodsimp(exp) 2917 elif exp > 10000: 2918 return a._eval_pow_by_cayley(exp) 2919 else: 2920 return a._eval_pow_by_recursion(exp) 2921 2922 if jordan_pow: 2923 try: 2924 return jordan_pow(exp) 2925 except NonInvertibleMatrixError: 2926 # Raised by jordan_pow on zero determinant matrix unless exp is 2927 # definitely known to be a non-negative integer. 2928 # Here we raise if n is definitely not a non-negative integer 2929 # but otherwise we can leave this as an unevaluated MatPow. 2930 if exp.is_integer is False or exp.is_nonnegative is False: 2931 raise 2932 2933 from sympy.matrices.expressions import MatPow 2934 return MatPow(a, exp) 2935 2936 @call_highest_priority('__add__') 2937 def __radd__(self, other): 2938 return self + other 2939 2940 @call_highest_priority('__matmul__') 2941 def __rmatmul__(self, other): 2942 other = _matrixify(other) 2943 if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False): 2944 return NotImplemented 2945 2946 return self.__rmul__(other) 2947 2948 @call_highest_priority('__mul__') 2949 def __rmul__(self, other): 2950 return self.rmultiply(other) 2951 2952 def rmultiply(self, other, dotprodsimp=None): 2953 """Same as __rmul__() but with optional simplification. 2954 2955 Parameters 2956 ========== 2957 2958 dotprodsimp : bool, optional 2959 Specifies whether intermediate term algebraic simplification is used 2960 during matrix multiplications to control expression blowup and thus 2961 speed up calculation. Default is off. 2962 """ 2963 isimpbool = _get_intermediate_simp_bool(False, dotprodsimp) 2964 other = _matrixify(other) 2965 # matrix-like objects can have shapes. This is 2966 # our first sanity check. Double check other is not explicitly not a Matrix. 2967 if (hasattr(other, 'shape') and len(other.shape) == 2 and 2968 (getattr(other, 'is_Matrix', True) or 2969 getattr(other, 'is_MatrixLike', True))): 2970 if self.shape[0] != other.shape[1]: 2971 raise ShapeError("Matrix size mismatch.") 2972 2973 # honest sympy matrices defer to their class's routine 2974 if getattr(other, 'is_Matrix', False): 2975 m = self._eval_matrix_rmul(other) 2976 if isimpbool: 2977 return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m]) 2978 return m 2979 # Matrix-like objects can be passed to CommonMatrix routines directly. 2980 if getattr(other, 'is_MatrixLike', False): 2981 return MatrixArithmetic._eval_matrix_rmul(self, other) 2982 2983 # if 'other' is not iterable then scalar multiplication. 2984 if not isinstance(other, Iterable): 2985 try: 2986 return self._eval_scalar_rmul(other) 2987 except TypeError: 2988 pass 2989 2990 return NotImplemented 2991 2992 @call_highest_priority('__sub__') 2993 def __rsub__(self, a): 2994 return (-self) + a 2995 2996 @call_highest_priority('__rsub__') 2997 def __sub__(self, a): 2998 return self + (-a) 2999 3000class MatrixCommon(MatrixArithmetic, MatrixOperations, MatrixProperties, 3001 MatrixSpecial, MatrixShaping): 3002 """All common matrix operations including basic arithmetic, shaping, 3003 and special matrices like `zeros`, and `eye`.""" 3004 _diff_wrt = True # type: bool 3005 3006 3007class _MinimalMatrix: 3008 """Class providing the minimum functionality 3009 for a matrix-like object and implementing every method 3010 required for a `MatrixRequired`. This class does not have everything 3011 needed to become a full-fledged SymPy object, but it will satisfy the 3012 requirements of anything inheriting from `MatrixRequired`. If you wish 3013 to make a specialized matrix type, make sure to implement these 3014 methods and properties with the exception of `__init__` and `__repr__` 3015 which are included for convenience.""" 3016 3017 is_MatrixLike = True 3018 _sympify = staticmethod(sympify) 3019 _class_priority = 3 3020 zero = S.Zero 3021 one = S.One 3022 3023 is_Matrix = True 3024 is_MatrixExpr = False 3025 3026 @classmethod 3027 def _new(cls, *args, **kwargs): 3028 return cls(*args, **kwargs) 3029 3030 def __init__(self, rows, cols=None, mat=None, copy=False): 3031 if isfunction(mat): 3032 # if we passed in a function, use that to populate the indices 3033 mat = list(mat(i, j) for i in range(rows) for j in range(cols)) 3034 if cols is None and mat is None: 3035 mat = rows 3036 rows, cols = getattr(mat, 'shape', (rows, cols)) 3037 try: 3038 # if we passed in a list of lists, flatten it and set the size 3039 if cols is None and mat is None: 3040 mat = rows 3041 cols = len(mat[0]) 3042 rows = len(mat) 3043 mat = [x for l in mat for x in l] 3044 except (IndexError, TypeError): 3045 pass 3046 self.mat = tuple(self._sympify(x) for x in mat) 3047 self.rows, self.cols = rows, cols 3048 if self.rows is None or self.cols is None: 3049 raise NotImplementedError("Cannot initialize matrix with given parameters") 3050 3051 def __getitem__(self, key): 3052 def _normalize_slices(row_slice, col_slice): 3053 """Ensure that row_slice and col_slice don't have 3054 `None` in their arguments. Any integers are converted 3055 to slices of length 1""" 3056 if not isinstance(row_slice, slice): 3057 row_slice = slice(row_slice, row_slice + 1, None) 3058 row_slice = slice(*row_slice.indices(self.rows)) 3059 3060 if not isinstance(col_slice, slice): 3061 col_slice = slice(col_slice, col_slice + 1, None) 3062 col_slice = slice(*col_slice.indices(self.cols)) 3063 3064 return (row_slice, col_slice) 3065 3066 def _coord_to_index(i, j): 3067 """Return the index in _mat corresponding 3068 to the (i,j) position in the matrix. """ 3069 return i * self.cols + j 3070 3071 if isinstance(key, tuple): 3072 i, j = key 3073 if isinstance(i, slice) or isinstance(j, slice): 3074 # if the coordinates are not slices, make them so 3075 # and expand the slices so they don't contain `None` 3076 i, j = _normalize_slices(i, j) 3077 3078 rowsList, colsList = list(range(self.rows))[i], \ 3079 list(range(self.cols))[j] 3080 indices = (i * self.cols + j for i in rowsList for j in 3081 colsList) 3082 return self._new(len(rowsList), len(colsList), 3083 list(self.mat[i] for i in indices)) 3084 3085 # if the key is a tuple of ints, change 3086 # it to an array index 3087 key = _coord_to_index(i, j) 3088 return self.mat[key] 3089 3090 def __eq__(self, other): 3091 try: 3092 classof(self, other) 3093 except TypeError: 3094 return False 3095 return ( 3096 self.shape == other.shape and list(self) == list(other)) 3097 3098 def __len__(self): 3099 return self.rows*self.cols 3100 3101 def __repr__(self): 3102 return "_MinimalMatrix({}, {}, {})".format(self.rows, self.cols, 3103 self.mat) 3104 3105 @property 3106 def shape(self): 3107 return (self.rows, self.cols) 3108 3109 3110class _CastableMatrix: # this is needed here ONLY FOR TESTS. 3111 def as_mutable(self): 3112 return self 3113 3114 def as_immutable(self): 3115 return self 3116 3117 3118class _MatrixWrapper: 3119 """Wrapper class providing the minimum functionality for a matrix-like 3120 object: .rows, .cols, .shape, indexability, and iterability. CommonMatrix 3121 math operations should work on matrix-like objects. This one is intended for 3122 matrix-like objects which use the same indexing format as SymPy with respect 3123 to returning matrix elements instead of rows for non-tuple indexes. 3124 """ 3125 3126 is_Matrix = False # needs to be here because of __getattr__ 3127 is_MatrixLike = True 3128 3129 def __init__(self, mat, shape): 3130 self.mat = mat 3131 self.shape = shape 3132 self.rows, self.cols = shape 3133 3134 def __getitem__(self, key): 3135 if isinstance(key, tuple): 3136 return sympify(self.mat.__getitem__(key)) 3137 3138 return sympify(self.mat.__getitem__((key // self.rows, key % self.cols))) 3139 3140 def __iter__(self): # supports numpy.matrix and numpy.array 3141 mat = self.mat 3142 cols = self.cols 3143 3144 return iter(sympify(mat[r, c]) for r in range(self.rows) for c in range(cols)) 3145 3146 3147class MatrixKind(Kind): 3148 """ 3149 Kind for all matrices in SymPy. 3150 3151 Basic class for this kind is ``MatrixBase`` and ``MatrixExpr``, 3152 but any expression representing the matrix can have this. 3153 3154 Parameters 3155 ========== 3156 3157 element_kind : Kind 3158 Kind of the element. Default is :obj:NumberKind `<sympy.core.kind.NumberKind>`, 3159 which means that the matrix contains only numbers. 3160 3161 Examples 3162 ======== 3163 3164 Any instance of matrix class has ``MatrixKind``. 3165 3166 >>> from sympy import MatrixSymbol 3167 >>> A = MatrixSymbol('A', 2,2) 3168 >>> A.kind 3169 MatrixKind(NumberKind) 3170 3171 Although expression representing a matrix may be not instance of 3172 matrix class, it will have ``MatrixKind`` as well. 3173 3174 >>> from sympy import Integral 3175 >>> from sympy.matrices.expressions import MatrixExpr 3176 >>> from sympy.abc import x 3177 >>> intM = Integral(A, x) 3178 >>> isinstance(intM, MatrixExpr) 3179 False 3180 >>> intM.kind 3181 MatrixKind(NumberKind) 3182 3183 Use ``isinstance()`` to check for ``MatrixKind` without specifying 3184 the element kind. Use ``is`` with specifying the element kind. 3185 3186 >>> from sympy import Matrix 3187 >>> from sympy.matrices import MatrixKind 3188 >>> from sympy.core.kind import NumberKind 3189 >>> M = Matrix([1, 2]) 3190 >>> isinstance(M.kind, MatrixKind) 3191 True 3192 >>> M.kind is MatrixKind(NumberKind) 3193 True 3194 3195 See Also 3196 ======== 3197 3198 shape : Function to return the shape of objects with ``MatrixKind``. 3199 3200 """ 3201 def __new__(cls, element_kind=NumberKind): 3202 obj = super().__new__(cls, element_kind) 3203 obj.element_kind = element_kind 3204 return obj 3205 3206 def __repr__(self): 3207 return "MatrixKind(%s)" % self.element_kind 3208 3209 3210def _matrixify(mat): 3211 """If `mat` is a Matrix or is matrix-like, 3212 return a Matrix or MatrixWrapper object. Otherwise 3213 `mat` is passed through without modification.""" 3214 3215 if getattr(mat, 'is_Matrix', False) or getattr(mat, 'is_MatrixLike', False): 3216 return mat 3217 3218 if not(getattr(mat, 'is_Matrix', True) or getattr(mat, 'is_MatrixLike', True)): 3219 return mat 3220 3221 shape = None 3222 3223 if hasattr(mat, 'shape'): # numpy, scipy.sparse 3224 if len(mat.shape) == 2: 3225 shape = mat.shape 3226 elif hasattr(mat, 'rows') and hasattr(mat, 'cols'): # mpmath 3227 shape = (mat.rows, mat.cols) 3228 3229 if shape: 3230 return _MatrixWrapper(mat, shape) 3231 3232 return mat 3233 3234 3235def a2idx(j, n=None): 3236 """Return integer after making positive and validating against n.""" 3237 if type(j) is not int: 3238 jindex = getattr(j, '__index__', None) 3239 if jindex is not None: 3240 j = jindex() 3241 else: 3242 raise IndexError("Invalid index a[%r]" % (j,)) 3243 if n is not None: 3244 if j < 0: 3245 j += n 3246 if not (j >= 0 and j < n): 3247 raise IndexError("Index out of range: a[%s]" % (j,)) 3248 return int(j) 3249 3250 3251def classof(A, B): 3252 """ 3253 Get the type of the result when combining matrices of different types. 3254 3255 Currently the strategy is that immutability is contagious. 3256 3257 Examples 3258 ======== 3259 3260 >>> from sympy import Matrix, ImmutableMatrix 3261 >>> from sympy.matrices.common import classof 3262 >>> M = Matrix([[1, 2], [3, 4]]) # a Mutable Matrix 3263 >>> IM = ImmutableMatrix([[1, 2], [3, 4]]) 3264 >>> classof(M, IM) 3265 <class 'sympy.matrices.immutable.ImmutableDenseMatrix'> 3266 """ 3267 priority_A = getattr(A, '_class_priority', None) 3268 priority_B = getattr(B, '_class_priority', None) 3269 if None not in (priority_A, priority_B): 3270 if A._class_priority > B._class_priority: 3271 return A.__class__ 3272 else: 3273 return B.__class__ 3274 3275 try: 3276 import numpy 3277 except ImportError: 3278 pass 3279 else: 3280 if isinstance(A, numpy.ndarray): 3281 return B.__class__ 3282 if isinstance(B, numpy.ndarray): 3283 return A.__class__ 3284 3285 raise TypeError("Incompatible classes %s, %s" % (A.__class__, B.__class__)) 3286