Skip to content
This repository has been archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
Refactor stacking, use coercion
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdemeyer
jdemeyer committed Jan 7, 2015
1 parent 0f1cbae commit e3d7bbd
Show file tree
Hide file tree
Showing 7 changed files with 159 additions and 183 deletions.
2 changes: 1 addition & 1 deletion src/sage/matrix/matrix1.pxd
View file
@@ -1,4 +1,4 @@
cimport matrix0

cdef class Matrix(matrix0.Matrix):
pass
cdef _stack_impl(self, bottom)
168 changes: 115 additions & 53 deletions src/sage/matrix/matrix1.pyx
View file
Expand Up @@ -23,6 +23,7 @@ include "sage/ext/python.pxi"

import sage.modules.free_module


cdef class Matrix(matrix0.Matrix):
###################################################
# Coercion to Various Systems
Expand Down Expand Up @@ -1116,8 +1117,10 @@ cdef class Matrix(matrix0.Matrix):
###########################################################################
def stack(self, bottom, subdivide=False):
r"""
Returns a new matrix formed by appending the matrix
(or vector) ``bottom`` beneath ``self``.
Return the matrix ``self`` on top of ``bottom``::
[ self ]
[ bottom ]
INPUT:
Expand All @@ -1142,6 +1145,7 @@ cdef class Matrix(matrix0.Matrix):
in the result.
.. warning::
If ``subdivide`` is ``True`` then unequal column subdivisions
will be discarded, since it would be ambiguous how to interpret
them. If the subdivision behavior is not what you need,
Expand Down Expand Up @@ -1187,13 +1191,13 @@ cdef class Matrix(matrix0.Matrix):
sage: A.stack(B)
Traceback (most recent call last):
...
TypeError: number of columns must be the same, 2 != 3
TypeError: number of columns must be the same, not 2 and 3
sage: v = vector(RR, [100, 200, 300])
sage: A.stack(v)
Traceback (most recent call last):
...
TypeError: number of columns must be the same, 2 != 3
TypeError: number of columns must be the same, not 2 and 3
Setting ``subdivide`` to ``True`` will, in its simplest form,
add a subdivision between ``self`` and ``bottom``. ::
Expand Down Expand Up @@ -1251,49 +1255,70 @@ cdef class Matrix(matrix0.Matrix):
[ 5 6 7 8 9]
[10 11 12 13 14]
The result retains the base ring of ``self`` by coercing
the elements of ``bottom`` into the base ring of ``self``. ::
The base ring of the result is the common parent for the base
rings of ``self`` and ``bottom``. In particular, the parent for
``A.stack(B)`` and ``B.stack(A)`` should be equal::
sage: A = matrix(QQ, 1, 2, [1,2])
sage: B = matrix(RR, 1, 2, [sin(1.1), sin(2.2)])
sage: C = A.stack(B); C
[ 1 2]
[183017397/205358938 106580492/131825561]
[ 1.00000000000000 2.00000000000000]
[0.891207360061435 0.808496403819590]
sage: C.parent()
Full MatrixSpace of 2 by 2 dense matrices over Rational Field
Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision
sage: D = B.stack(A); D
[0.891207360061435 0.808496403819590]
[ 1.00000000000000 2.00000000000000]
sage: D.parent()
Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision
Sometimes it is not possible to coerce into the base ring
of ``self``. A solution is to change the base ring of ``self``
to a more expansive ring. Here we mix the rationals with
a ring of polynomials with rational coefficients. ::
::
sage: R = PolynomialRing(QQ, 'y')
sage: A = matrix(QQ, 1, 2, [1,2])
sage: B = matrix(R, 1, 2, ['y', 'y^2'])
sage: R.<y> = PolynomialRing(ZZ)
sage: A = matrix(QQ, 1, 2, [1, 2/3])
sage: B = matrix(R, 1, 2, [y, y^2])
sage: C = B.stack(A); C
sage: C = A.stack(B); C
[ 1 2/3]
[ y y^2]
[ 1 2]
sage: C.parent()
Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field
sage: D = A.stack(B)
Traceback (most recent call last):
...
TypeError: not a constant polynomial
Stacking a dense matrix atop a sparse one returns a sparse
matrix::
sage: E = A.change_ring(R)
sage: F = E.stack(B); F
[ 1 2]
[ y y^2]
sage: F.parent()
Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field
sage: M = Matrix(ZZ, 2, 3, range(6), sparse=False)
sage: N = diagonal_matrix([10,11,12], sparse=True)
sage: P = M.stack(N); P
[ 0 1 2]
[ 3 4 5]
[10 0 0]
[ 0 11 0]
[ 0 0 12]
sage: P.is_sparse()
True
sage: P = N.stack(M); P
[10 0 0]
[ 0 11 0]
[ 0 0 12]
[ 0 1 2]
[ 3 4 5]
sage: P.is_sparse()
True
One can stack matrices over different rings (:trac:`16399`). ::
sage: M = Matrix(ZZ, 2, 3, range(6))
sage: N = Matrix(QQ, 1, 3, [10,11,12])
sage: M.stack(N)
[ 0 1 2]
[ 3 4 5]
[10 11 12]
sage: N.stack(M)
[10 11 12]
[ 0 1 2]
[ 3 4 5]
TESTS:
Expand All @@ -1306,43 +1331,80 @@ cdef class Matrix(matrix0.Matrix):
[ 3 4 5]
[10 11 12]
AUTHOR:
Non-matrices fail gracefully::
- Rob Beezer (2011-03-19) - rewritten to mirror code for :meth:`augment`
"""
from sage.matrix.constructor import matrix
sage: M.stack(polygen(QQ))
Traceback (most recent call last):
...
TypeError: a matrix must be stacked with another matrix or a vector
AUTHORS:
if hasattr(bottom, '_vector_'):
bottom = bottom.row()
if not isinstance(bottom, sage.matrix.matrix1.Matrix):
raise TypeError('a matrix must be stacked with another matrix, or '
'a vector')
- Rob Beezer (2011-03-19): rewritten to mirror code for :meth:`augment`
- Jeroen Demeyer (2015-01-06): refactor, see :trac:`16399`.
Put all boilerplate in one place (here) and put the actual
type-dependent implementation in ``_stack_impl``.
"""
cdef Matrix other
other = bottom
if isinstance(bottom, Matrix):
other = <Matrix>bottom
else:
if hasattr(bottom, '_vector_'):
bottom = bottom.row()
else:
raise TypeError('a matrix must be stacked with '
'another matrix or a vector')
other = <Matrix?>bottom

if self._ncols != other._ncols:
raise TypeError('number of columns must be the same, '
'{0} != {1}'.format(self._ncols, other._ncols))
if not (self._base_ring is other.base_ring()):
other = other.change_ring(self._base_ring)
raise TypeError("number of columns must be the same, not %s and %s" %
(self.ncols(), bottom.ncols()) )

top_ring = self._base_ring
bottom_ring = other._base_ring
if top_ring is not bottom_ring:
from sage.structure.element import get_coercion_model
coercion_model = get_coercion_model()
R = coercion_model.common_parent(top_ring, bottom_ring)
if top_ring is not R:
self = self.change_ring(R)
if bottom_ring is not R:
other = other.change_ring(R)

if type(self) is not type(other):
# If one of the matrices is sparse, return a sparse matrix
if self.is_sparse_c() and not other.is_sparse_c():
other = other.sparse_matrix()
elif other.is_sparse_c() and not self.is_sparse_c():
self = self.sparse_matrix()

Z = self._stack_impl(other)
if subdivide:
Z._subdivide_on_stack(self, other)
return Z

cdef _stack_impl(self, bottom):
"""
Implementation of :meth:`stack`.
Assume that ``self`` and ``other`` are compatible in the sense
that they have the same base ring and that both are either
dense or sparse.
"""
cdef Matrix other = <Matrix>bottom
cdef Matrix Z
Z = self.new_matrix(nrows = self._nrows + other._nrows)
Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols)

cdef Py_ssize_t r, c
for r from 0 <= r < self._nrows:
for c from 0 <= c < self._ncols:
Z.set_unsafe(r,c, self.get_unsafe(r,c))
nr = self.nrows()

for r from 0 <= r < other._nrows:
for c from 0 <= c < other._ncols:
cdef Py_ssize_t nr = self._nrows
for r in range(self._nrows):
for c in range(self._ncols):
Z.set_unsafe(r, c, self.get_unsafe(r,c))
for r in range(other._nrows):
for c in range(other._ncols):
Z.set_unsafe(r+nr, c, other.get_unsafe(r,c))

if subdivide:
Z._subdivide_on_stack(self, other)

return Z

def augment(self, right, subdivide=False):
Expand Down
84 changes: 16 additions & 68 deletions src/sage/matrix/matrix_integer_dense.pyx
View file
Expand Up @@ -95,7 +95,7 @@ from sage.rings.integer_ring cimport IntegerRing_class
from sage.rings.finite_rings.integer_mod_ring import IntegerModRing
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint
from sage.structure.element cimport ModuleElement, RingElement, Element, Vector, CoercionModel
from sage.structure.element cimport ModuleElement, RingElement, Element, Vector
from sage.structure.element import is_Vector
from sage.structure.sequence import Sequence

Expand Down Expand Up @@ -4679,9 +4679,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse
#################################################################
# operations with matrices
#################################################################
def stack(self, bottom, subdivide=False):
cdef _stack_impl(self, bottom):
r"""
Return the matrix ``self`` on top of ``bottom``:
Return the matrix ``self`` on top of ``bottom``::
[ self ]
[ bottom ]
Expand Down Expand Up @@ -4718,73 +4718,21 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse
[-----------]
[ 0 1 2 3]
[ 4 5 6 7]
"""
cdef Matrix_integer_dense other = <Matrix_integer_dense>bottom
cdef Matrix_integer_dense Z
Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols)

TESTS:
Stacking a dense matrix atop a sparse one should work::
sage: M = Matrix(ZZ, 2, 3, range(6))
sage: M.is_sparse()
False
sage: N = diagonal_matrix([10,11,12], sparse=True)
sage: N.is_sparse()
True
sage: P = M.stack(N); P
[ 0 1 2]
[ 3 4 5]
[10 0 0]
[ 0 11 0]
[ 0 0 12]
sage: P.is_sparse()
False
One can stack matrices over different rings (:trac:`16399`). ::
sage: M = Matrix(ZZ, 2, 3, range(6))
sage: N = Matrix(QQ, 1, 3, [10,11,12])
sage: M.stack(N)
[0 1 2 3 4]
[5 6 7 8 9]
[0 1 2 3 4]
sage: N.stack(M)
?
sage: M2 = Matrix(ZZ['x'], 2, 3, range(6))
sage: N.stack(M2)
?
"""
if hasattr(bottom, '_vector_'):
bottom = bottom.row()
if self._ncols != bottom.ncols():
raise TypeError("number of columns must be the same")
cdef CoercionModel coercion_model

top_ring = self._base_ring
bottom_ring = bottom.base_ring()
if not (top_ring is bottom_ring):
if top_ring.has_coerce_map_from(bottom_ring):
bottom = bottom.change_ring(top_ring)
elif bottom_ring.has_coerce_map_from(top_ring):
new_top = self.change_ring(bottom_ring)
return new_top.stack(bottom, subdivide=subdivide):
else:
from sage.structure.element import get_coercion_model
coercion_model = get_coercion_model()
com_ring = coercion_model.common_parent(top_ring, bottom_ring)
self = self.change_ring(com_ring) # allowed ?
bottom = other.change_ring(com_ring)
cdef Py_ssize_t r, c
cdef Py_ssize_t nr = self._nrows
for r in range(self._nrows):
for c in range(self._ncols):
fmpz_set(fmpz_mat_entry(Z._matrix, r, c),fmpz_mat_entry(self._matrix, r, c))
for r in range(other._nrows):
for c in range(other._ncols):
fmpz_set(fmpz_mat_entry(Z._matrix, r+nr, c),fmpz_mat_entry(other._matrix, r, c))

cdef Matrix_integer_dense other = bottom.dense_matrix()
cdef Matrix_integer_dense M
M = self.new_matrix(nrows = self._nrows + other._nrows, ncols = self.ncols())
cdef Py_ssize_t i, j, k
for j from 0 <= j < self._ncols:
for i from 0 <= i < self._nrows:
fmpz_set(fmpz_mat_entry(M._matrix,i,j),fmpz_mat_entry(self._matrix,i,j))
for i from 0 <= i < other._nrows:
fmpz_set(fmpz_mat_entry(M._matrix,self._nrows + i,j),fmpz_mat_entry(other._matrix,i,j))
if subdivide:
M._subdivide_on_stack(self, other)
return M
return Z

def augment(self, right, subdivide=False):
r"""
Expand Down

0 comments on commit e3d7bbd

Please sign in to comment.