Ticket #10974: trac_10974_stack-augment-general.patch

File trac_10974_stack-augment-general.patch, 17.2 KB (added by rbeezer, 8 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  
    984984    ############################################################################################
    985985    # Building matrices out of other matrices, rows, or columns
    986986    ############################################################################################
    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``.
    990991
    991             [ self  ]
    992             [ other ]
     992        INPUT:
    993993
    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 1-row 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 11|12 13|14]
     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.  ::
    9951171
    9961172            sage: M = Matrix(QQ, 2, 3, range(6))
    9971173            sage: N = Matrix(QQ, 1, 3, [10,11,12])
     
    9991175            [ 0  1  2]
    10001176            [ 3  4  5]
    10011177            [10 11 12]
     1178
     1179        AUTHOR:
     1180
     1181        - Rob Beezer (2011-03-19) - rewritten to mirror code for :meth:`augment`
    10021182        """
    1003         cdef Py_ssize_t r, c
     1183        from sage.matrix.constructor import matrix
    10041184
    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')
    10071189
    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))
    10101195        if not (self._base_ring is other.base_ring()):
    10111196            other = other.change_ring(self._base_ring)
    10121197
    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
    10151214        return Z
    10161215
     1216
    10171217    def matrix_from_columns(self, columns):
    10181218        """
    10191219        Return the matrix constructed from self using columns with indices
     
    15141714
    15151715        Augmenting with a matrix. ::
    15161716
    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))
    15191719            sage: A.augment(B)
    15201720            [ 0  1  2  3  0  1  2]
    15211721            [ 4  5  6  7  3  4  5]
     
    15361736            sage: A.augment(B)
    15371737            Traceback (most recent call last):
    15381738            ...
    1539             TypeError: number of rows must be the same
     1739            TypeError: number of rows must be the same, 2 != 3
    15401740
    15411741            sage: v = vector(RR, [100, 200, 300])
    15421742            sage: A.augment(v)
    15431743            Traceback (most recent call last):
    15441744            ...
    1545             TypeError: number of rows must be the same
     1745            TypeError: number of rows must be the same, 2 != 3
    15461746
    15471747        Setting ``subdivide`` to ``True`` will, in its simplest form,
    15481748        add a subdivision between ``self`` and ``right``. ::
    15491749
    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))
    15521752            sage: A.augment(B, subdivide=True)
    15531753            [ 0  1  2  3| 0  1  2  3  4]
    15541754            [ 4  5  6  7| 5  6  7  8  9]
     
    15581758        if subdivisions are requested.  (So multiple augmentations can
    15591759        be recorded.) ::
    15601760
    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])
    15651765            sage: A.augment(B, subdivide=True)
    15661766            [0|1|0 1|2]
    15671767            [2|3|3 4|5]
     
    15711771        Otherwise, this information is discarded and must be managed
    15721772        separately. ::
    15731773
    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)
    15781778            sage: A.augment(B, subdivide=True)
    15791779            [0 1|0 1 2]
    15801780            [---+-----]
     
    15821782            [4 5|6 7 8]
    15831783            [---+-----]
    15841784
    1585             sage: A.subdivide([1,2],None)
     1785            sage: A.subdivide([1,2], None)
    15861786            sage: A.augment(B, subdivide=True)
    15871787            [0 1|0 1 2]
    15881788            [2 3|3 4 5]
     
    16381838        from sage.matrix.constructor import matrix
    16391839
    16401840        if hasattr(right, '_vector_'):
    1641             right = matrix(right.degree(), right.list())
     1841            right = right.column()
    16421842        if not isinstance(right, sage.matrix.matrix1.Matrix):
    16431843            raise TypeError("a matrix must be augmented with another matrix, or a vector")
    16441844
     
    16461846        other = right
    16471847
    16481848        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))
    16501850        if not (self._base_ring is other.base_ring()):
    16511851            other = other.change_ring(self._base_ring)
    16521852
     
    16641864                Z.set_unsafe(r, c+nc, other.get_unsafe(r,c))
    16651865
    16661866        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)
    16771868
    16781869        return Z
    16791870
  • sage/matrix/matrix2.pyx

    diff -r 6a679959b54b -r cc0a175bb6a0 sage/matrix/matrix2.pyx
    a b  
    52025202            raise IndexError, "Submatrix %s,%s has no entry %s,%s"%(i,j, x, y)
    52035203        return self[self.subdivisions[0][i] + x , self.subdivisions[1][j] + y]
    52045204
     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 11|12 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
    52055300    def get_subdivisions(self):
    52065301        """
    52075302        Returns the current subdivision of self.