Ticket #10974: trac_10974_stackaugmentgeneral.patch
File trac_10974_stackaugmentgeneral.patch, 17.2 KB (added by , 9 years ago) 


sage/matrix/matrix1.pyx
# HG changeset patch # User Rob Beezer <beezer@ups.edu> # Date 1300742627 25200 # Node ID cc0a175bb6a0b06d7a2d3f114cc433585aa93d2f # Parent 6a679959b54b7975af1841df4dd77d28b0faef28 10974: Overhaul matrix stack, augment for general classes diff r 6a679959b54b r cc0a175bb6a0 sage/matrix/matrix1.pyx
a b 984 984 ############################################################################################ 985 985 # Building matrices out of other matrices, rows, or columns 986 986 ############################################################################################ 987 def stack(self, other): 988 """ 989 Return the augmented matrix self on top of other:: 987 def stack(self, bottom, subdivide=False): 988 r""" 989 Returns a new matrix formed by appending the matrix 990 (or vector) ``bottom`` beneath ``self``. 990 991 991 [ self ] 992 [ other ] 992 INPUT: 993 993 994 EXAMPLES:: 994  ``bottom``  a matrix, vector or free module element, whose 995 dimensions are compatible with ``self``. 996 997  ``subdivide``  default: ``False``  request the resulting 998 matrix to have a new subdivision, separating ``self`` from ``bottom``. 999 1000 OUTPUT: 1001 1002 A new matrix formed by appending ``bottom`` beneath ``self``. 1003 If ``bottom`` is a vector (or free module element) then in this context 1004 it is appropriate to consider it as a row vector. (The code first 1005 converts a vector to a 1row matrix.) 1006 1007 If ``subdivide`` is ``True`` then any row subdivisions for 1008 the two matrices are preserved, and a new subdivision is added 1009 between ``self`` and ``bottom``. If the column divisions are 1010 identical, then they are preserved, otherwise they are discarded. 1011 When ``subdivide`` is ``False`` there is no subdivision information 1012 in the result. 1013 1014 .. warning:: 1015 If ``subdivide`` is ``True`` then unequal column subdivisions 1016 will be discarded, since it would be ambiguous how to interpret 1017 them. If the subdivision behavior is not what you need, 1018 you can manage subdivisions yourself with methods like 1019 :meth:`~sage.matrix.matrix2.Matrix.get_subdivisions` 1020 and 1021 :meth:`~sage.matrix.matrix2.Matrix.subdivide`. 1022 You might also find :func:`~sage.matrix.constructor.block_matrix` 1023 or 1024 :func:`~sage.matrix.constructor.block_diagonal_matrix` 1025 useful and simpler in some instances. 1026 1027 1028 EXAMPLES: 1029 1030 Stacking with a matrix. :: 1031 1032 sage: A = matrix(QQ, 4, 3, range(12)) 1033 sage: B = matrix(QQ, 3, 3, range(9)) 1034 sage: A.stack(B) 1035 [ 0 1 2] 1036 [ 3 4 5] 1037 [ 6 7 8] 1038 [ 9 10 11] 1039 [ 0 1 2] 1040 [ 3 4 5] 1041 [ 6 7 8] 1042 1043 Stacking with a vector. :: 1044 1045 sage: A = matrix(QQ, 3, 2, [0, 2, 4, 6, 8, 10]) 1046 sage: v = vector(QQ, 2, [100, 200]) 1047 sage: A.stack(v) 1048 [ 0 2] 1049 [ 4 6] 1050 [ 8 10] 1051 [100 200] 1052 1053 Errors are raised if the sizes are incompatible. :: 1054 1055 sage: A = matrix(RR, [[1, 2],[3, 4]]) 1056 sage: B = matrix(RR, [[10, 20, 30], [40, 50, 60]]) 1057 sage: A.stack(B) 1058 Traceback (most recent call last): 1059 ... 1060 TypeError: number of columns must be the same, 2 != 3 1061 1062 sage: v = vector(RR, [100, 200, 300]) 1063 sage: A.stack(v) 1064 Traceback (most recent call last): 1065 ... 1066 TypeError: number of columns must be the same, 2 != 3 1067 1068 Setting ``subdivide`` to ``True`` will, in its simplest form, 1069 add a subdivision between ``self`` and ``bottom``. :: 1070 1071 sage: A = matrix(QQ, 2, 5, range(10)) 1072 sage: B = matrix(QQ, 3, 5, range(15)) 1073 sage: A.stack(B, subdivide=True) 1074 [ 0 1 2 3 4] 1075 [ 5 6 7 8 9] 1076 [] 1077 [ 0 1 2 3 4] 1078 [ 5 6 7 8 9] 1079 [10 11 12 13 14] 1080 1081 Row subdivisions are preserved by stacking, and enriched, 1082 if subdivisions are requested. (So multiple stackings can 1083 be recorded.) :: 1084 1085 sage: A = matrix(QQ, 2, 4, range(8)) 1086 sage: A.subdivide([1], None) 1087 sage: B = matrix(QQ, 3, 4, range(12)) 1088 sage: B.subdivide([2], None) 1089 sage: A.stack(B, subdivide=True) 1090 [ 0 1 2 3] 1091 [] 1092 [ 4 5 6 7] 1093 [] 1094 [ 0 1 2 3] 1095 [ 4 5 6 7] 1096 [] 1097 [ 8 9 10 11] 1098 1099 Column subdivisions can be preserved, but only if they are identical. 1100 Otherwise, this information is discarded and must be managed 1101 separately. :: 1102 1103 sage: A = matrix(QQ, 2, 5, range(10)) 1104 sage: A.subdivide(None, [2,4]) 1105 sage: B = matrix(QQ, 3, 5, range(15)) 1106 sage: B.subdivide(None, [2,4]) 1107 sage: A.stack(B, subdivide=True) 1108 [ 0 1 2 3 4] 1109 [ 5 6 7 8 9] 1110 [++] 1111 [ 0 1 2 3 4] 1112 [ 5 6 7 8 9] 1113 [10 1112 1314] 1114 1115 sage: A.subdivide(None, [1,2]) 1116 sage: A.stack(B, subdivide=True) 1117 [ 0 1 2 3 4] 1118 [ 5 6 7 8 9] 1119 [] 1120 [ 0 1 2 3 4] 1121 [ 5 6 7 8 9] 1122 [10 11 12 13 14] 1123 1124 The result retains the base ring of ``self`` by coercing 1125 the elements of ``bottom`` into the base ring of ``self``. :: 1126 1127 sage: A = matrix(QQ, 1, 2, [1,2]) 1128 sage: B = matrix(RR, 1, 2, [sin(1.1), sin(2.2)]) 1129 sage: C = A.stack(B); C 1130 [ 1 2] 1131 [183017397/205358938 106580492/131825561] 1132 sage: C.parent() 1133 Full MatrixSpace of 2 by 2 dense matrices over Rational Field 1134 1135 sage: D = B.stack(A); D 1136 [0.891207360061435 0.808496403819590] 1137 [ 1.00000000000000 2.00000000000000] 1138 sage: D.parent() 1139 Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision 1140 1141 Sometimes it is not possible to coerce into the base ring 1142 of ``self``. A solution is to change the base ring of ``self`` 1143 to a more expansive ring. Here we mix the rationals with 1144 a ring of polynomials with rational coefficients. :: 1145 1146 sage: R = PolynomialRing(QQ, 'y') 1147 sage: A = matrix(QQ, 1, 2, [1,2]) 1148 sage: B = matrix(R, 1, 2, ['y', 'y^2']) 1149 1150 sage: C = B.stack(A); C 1151 [ y y^2] 1152 [ 1 2] 1153 sage: C.parent() 1154 Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field 1155 1156 sage: D = A.stack(B) 1157 Traceback (most recent call last): 1158 ... 1159 TypeError: not a constant polynomial 1160 1161 sage: E = A.change_ring(R) 1162 sage: F = E.stack(B); F 1163 [ 1 2] 1164 [ y y^2] 1165 sage: F.parent() 1166 Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field 1167 1168 TESTS: 1169 1170 A legacy test from the original implementation. :: 995 1171 996 1172 sage: M = Matrix(QQ, 2, 3, range(6)) 997 1173 sage: N = Matrix(QQ, 1, 3, [10,11,12]) … … 999 1175 [ 0 1 2] 1000 1176 [ 3 4 5] 1001 1177 [10 11 12] 1178 1179 AUTHOR: 1180 1181  Rob Beezer (20110319)  rewritten to mirror code for :meth:`augment` 1002 1182 """ 1003 cdef Py_ssize_t r, c1183 from sage.matrix.constructor import matrix 1004 1184 1005 if not isinstance(other, Matrix): 1006 raise TypeError, "other must be a matrix" 1185 if hasattr(bottom, '_vector_'): 1186 bottom = bottom.row() 1187 if not isinstance(bottom, sage.matrix.matrix1.Matrix): 1188 raise TypeError('a matrix must be stacked with another matrix, or a vector') 1007 1189 1008 if self._ncols != other.ncols(): 1009 raise TypeError, "number of columns must be the same" 1190 cdef Matrix other 1191 other = bottom 1192 1193 if self._ncols != other._ncols: 1194 raise TypeError('number of columns must be the same, {0} != {1}'.format(self._ncols, other._ncols)) 1010 1195 if not (self._base_ring is other.base_ring()): 1011 1196 other = other.change_ring(self._base_ring) 1012 1197 1013 v = self.list() + other.list() 1014 Z = self.new_matrix(nrows = self._nrows + other.nrows(), entries=v, coerce=False, copy=False) 1198 cdef Matrix Z 1199 Z = self.new_matrix(nrows = self._nrows + other._nrows) 1200 1201 cdef Py_ssize_t r, c 1202 for r from 0 <= r < self._nrows: 1203 for c from 0 <= c < self._ncols: 1204 Z.set_unsafe(r,c, self.get_unsafe(r,c)) 1205 nr = self.nrows() 1206 1207 for r from 0 <= r < other._nrows: 1208 for c from 0 <= c < other._ncols: 1209 Z.set_unsafe(r+nr, c, other.get_unsafe(r,c)) 1210 1211 if subdivide: 1212 Z._subdivide_on_stack(self, other) 1213 1015 1214 return Z 1016 1215 1216 1017 1217 def matrix_from_columns(self, columns): 1018 1218 """ 1019 1219 Return the matrix constructed from self using columns with indices … … 1514 1714 1515 1715 Augmenting with a matrix. :: 1516 1716 1517 sage: A = matrix( ZZ, 3, range(12))1518 sage: B = matrix( ZZ, 3, range(9))1717 sage: A = matrix(QQ, 3, range(12)) 1718 sage: B = matrix(QQ, 3, range(9)) 1519 1719 sage: A.augment(B) 1520 1720 [ 0 1 2 3 0 1 2] 1521 1721 [ 4 5 6 7 3 4 5] … … 1536 1736 sage: A.augment(B) 1537 1737 Traceback (most recent call last): 1538 1738 ... 1539 TypeError: number of rows must be the same 1739 TypeError: number of rows must be the same, 2 != 3 1540 1740 1541 1741 sage: v = vector(RR, [100, 200, 300]) 1542 1742 sage: A.augment(v) 1543 1743 Traceback (most recent call last): 1544 1744 ... 1545 TypeError: number of rows must be the same 1745 TypeError: number of rows must be the same, 2 != 3 1546 1746 1547 1747 Setting ``subdivide`` to ``True`` will, in its simplest form, 1548 1748 add a subdivision between ``self`` and ``right``. :: 1549 1749 1550 sage: A = matrix( 3, range(12))1551 sage: B = matrix( 3, range(15))1750 sage: A = matrix(QQ, 3, range(12)) 1751 sage: B = matrix(QQ, 3, range(15)) 1552 1752 sage: A.augment(B, subdivide=True) 1553 1753 [ 0 1 2 3 0 1 2 3 4] 1554 1754 [ 4 5 6 7 5 6 7 8 9] … … 1558 1758 if subdivisions are requested. (So multiple augmentations can 1559 1759 be recorded.) :: 1560 1760 1561 sage: A = matrix( 3, range(6))1562 sage: A.subdivide(None, [1])1563 sage: B = matrix( 3, range(9))1564 sage: B.subdivide(None, [2])1761 sage: A = matrix(QQ, 3, range(6)) 1762 sage: A.subdivide(None, [1]) 1763 sage: B = matrix(QQ, 3, range(9)) 1764 sage: B.subdivide(None, [2]) 1565 1765 sage: A.augment(B, subdivide=True) 1566 1766 [010 12] 1567 1767 [233 45] … … 1571 1771 Otherwise, this information is discarded and must be managed 1572 1772 separately. :: 1573 1773 1574 sage: A = matrix( 3, range(6))1575 sage: A.subdivide([1,3], None)1576 sage: B = matrix( 3, range(9))1577 sage: B.subdivide([1,3], None)1774 sage: A = matrix(QQ, 3, range(6)) 1775 sage: A.subdivide([1,3], None) 1776 sage: B = matrix(QQ, 3, range(9)) 1777 sage: B.subdivide([1,3], None) 1578 1778 sage: A.augment(B, subdivide=True) 1579 1779 [0 10 1 2] 1580 1780 [+] … … 1582 1782 [4 56 7 8] 1583 1783 [+] 1584 1784 1585 sage: A.subdivide([1,2], None)1785 sage: A.subdivide([1,2], None) 1586 1786 sage: A.augment(B, subdivide=True) 1587 1787 [0 10 1 2] 1588 1788 [2 33 4 5] … … 1638 1838 from sage.matrix.constructor import matrix 1639 1839 1640 1840 if hasattr(right, '_vector_'): 1641 right = matrix(right.degree(), right.list())1841 right = right.column() 1642 1842 if not isinstance(right, sage.matrix.matrix1.Matrix): 1643 1843 raise TypeError("a matrix must be augmented with another matrix, or a vector") 1644 1844 … … 1646 1846 other = right 1647 1847 1648 1848 if self._nrows != other._nrows: 1649 raise TypeError , "number of rows must be the same"1849 raise TypeError('number of rows must be the same, {0} != {1}'.format(self._nrows, other._nrows)) 1650 1850 if not (self._base_ring is other.base_ring()): 1651 1851 other = other.change_ring(self._base_ring) 1652 1852 … … 1664 1864 Z.set_unsafe(r, c+nc, other.get_unsafe(r,c)) 1665 1865 1666 1866 if subdivide: 1667 self_subs_rows, self_subs_cols = self.get_subdivisions() 1668 other_subs_rows, other_subs_cols = other.get_subdivisions() 1669 if self_subs_rows == other_subs_rows: 1670 z_subs_rows = self_subs_rows 1671 else: 1672 z_subs_rows = None 1673 z_subs_cols = self_subs_cols + [nc] 1674 for col in other_subs_cols: 1675 z_subs_cols.append(col+nc) 1676 Z.subdivide(z_subs_rows, z_subs_cols) 1867 Z._subdivide_on_augment(self, other) 1677 1868 1678 1869 return Z 1679 1870 
sage/matrix/matrix2.pyx
diff r 6a679959b54b r cc0a175bb6a0 sage/matrix/matrix2.pyx
a b 5202 5202 raise IndexError, "Submatrix %s,%s has no entry %s,%s"%(i,j, x, y) 5203 5203 return self[self.subdivisions[0][i] + x , self.subdivisions[1][j] + y] 5204 5204 5205 def _subdivide_on_augment(self, left, right): 5206 r""" 5207 Helper method to manage subdivisions when augmenting a matrix. 5208 5209 INPUT: 5210 5211  ``left``, ``right``  two matrices, such that if ``left`` is 5212 augmented by placing ``right`` on the right side of ``left``, 5213 then the result is ``self``. It is the responsibility of the 5214 calling routine to ensure this condition holds. 5215 5216 OUTPUT: 5217 5218 ``None``. A new subdivision is created between ``left`` and 5219 ``right`` for ``self``. If possible, row subdivisions are 5220 preserved in ``self``, but if the two sets of row subdivisions 5221 are incompatible, they are removed. 5222 5223 EXAMPLE:: 5224 5225 sage: A = matrix(QQ, 3, 4, range(12)) 5226 sage: B = matrix(QQ, 3, 6, range(18)) 5227 sage: C = A.augment(B) 5228 sage: C._subdivide_on_augment(A, B) 5229 sage: C 5230 [ 0 1 2 3 0 1 2 3 4 5] 5231 [ 4 5 6 7 6 7 8 9 10 11] 5232 [ 8 9 10 1112 13 14 15 16 17] 5233 5234 More descriptive, but indirect, doctests are at 5235 :meth:`sage.matrix.matrix1.Matrix.augment`. 5236 """ 5237 left_rows, left_cols = left.get_subdivisions() 5238 right_rows, right_cols = right.get_subdivisions() 5239 if left_rows == right_rows: 5240 self_rows = left_rows 5241 else: 5242 self_rows = None 5243 nc = left.ncols() 5244 self_cols = left_cols + [nc] 5245 for col in right_cols: 5246 self_cols.append(col+nc) 5247 self.subdivide(self_rows, self_cols) 5248 return None 5249 5250 def _subdivide_on_stack(self, top, bottom): 5251 r""" 5252 Helper method to manage subdivisions when stacking a matrix. 5253 5254 INPUT: 5255 5256  ``top``, ``bottom``  two matrices, such that if ``top`` is 5257 stacked by placing ``top`` above ``bottom``, then the result 5258 is ``self``. It is the responsibility of the calling routine 5259 to ensure this condition holds. 5260 5261 OUTPUT: 5262 5263 ``None``. A new subdivision is created between ``top`` and 5264 ``bottom`` for ``self``. If possible, column subdivisions are 5265 preserved in ``self``, but if the two sets of solumn subdivisions 5266 are incompatible, they are removed. 5267 5268 EXAMPLE:: 5269 5270 sage: A = matrix(QQ, 3, 2, range(6)) 5271 sage: B = matrix(QQ, 4, 2, range(8)) 5272 sage: C = A.stack(B) 5273 sage: C._subdivide_on_stack(A, B) 5274 sage: C 5275 [0 1] 5276 [2 3] 5277 [4 5] 5278 [] 5279 [0 1] 5280 [2 3] 5281 [4 5] 5282 [6 7] 5283 5284 More descriptive, but indirect, doctests are at 5285 :meth:`sage.matrix.matrix1.Matrix.augment`. 5286 """ 5287 top_rows, top_cols = top.get_subdivisions() 5288 bottom_rows, bottom_cols = bottom.get_subdivisions() 5289 if top_cols == bottom_cols: 5290 self_cols = top_cols 5291 else: 5292 self_cols = None 5293 nr = top.nrows() 5294 self_rows = top_rows + [nr] 5295 for row in bottom_rows: 5296 self_rows.append(row+nr) 5297 self.subdivide(self_rows, self_cols) 5298 return None 5299 5205 5300 def get_subdivisions(self): 5206 5301 """ 5207 5302 Returns the current subdivision of self.