From d66b7bcf97f0018a104a3513dd88d19cca68adcd Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Mon, 11 Oct 2021 22:37:13 -0700 Subject: [PATCH 1/9] src/sage/matrix/matrix_double_dense.pyx: Update copyright according to 'git blame -w --date=format:%Y src/sage/matrix/matrix_double_dense.pyx | sort -k2' --- src/sage/matrix/matrix_double_dense.pyx | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_double_dense.pyx b/src/sage/matrix/matrix_double_dense.pyx index 8168239ff1c..fea7a93f88c 100644 --- a/src/sage/matrix/matrix_double_dense.pyx +++ b/src/sage/matrix/matrix_double_dense.pyx @@ -35,7 +35,27 @@ TESTS:: """ # **************************************************************************** -# Copyright (C) 2004,2005,2006 Joshua Kantor +# Copyright (C) 2004-2006 Joshua Kantor +# Copyright (C) 2008 Georg S. Weber +# Copyright (C) 2008-2011 Mike Hansen +# Copyright (C) 2008-2012 Jason Grout +# Copyright (C) 2009 Dag Sverre Seljebotn +# Copyright (C) 2009 Yann Laigle-Chapuy +# Copyright (C) 2009-2010 Florent Hivert +# Copyright (C) 2010-2012 Rob Beezer +# Copyright (C) 2011 Martin Raum +# Copyright (C) 2011-2012 J. H. Palmieri +# Copyright (C) 2011-2014 André Apitzsch +# Copyright (C) 2011-2018 Jeroen Demeyer +# Copyright (C) 2012 Kenneth Smith +# Copyright (C) 2016-2019 Frédéric Chapoton +# Copyright (C) 2017 Kiran Kedlaya +# Copyright (C) 2019 Chaman Agrawal +# Copyright (C) 2019-2021 Markus Wageringel +# Copyright (C) 2020 Michael Orlitzky +# Copyright (C) 2020 Victor Santos +# Copyright (C) 2021 Jonathan Kliem +# Copyright (C) 2021 Travis Scrimshaw # # Distributed under the terms of the GNU General Public License (GPL) # as published by the Free Software Foundation; either version 2 of From d33d173d96d8dc4c146d502a0a9c31b8db154899 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Mon, 11 Oct 2021 22:55:06 -0700 Subject: [PATCH 2/9] sage.matrix: Factor Matrix_double_dense through Matrix_numpy_dense, add Matrix_numpy_integer_dense --- src/sage/matrix/matrix_double_dense.pxd | 16 +- src/sage/matrix/matrix_double_dense.pyx | 403 +--------------- src/sage/matrix/matrix_numpy_dense.pxd | 12 + src/sage/matrix/matrix_numpy_dense.pyx | 448 ++++++++++++++++++ .../matrix/matrix_numpy_integer_dense.pxd | 6 + .../matrix/matrix_numpy_integer_dense.pyx | 58 +++ 6 files changed, 531 insertions(+), 412 deletions(-) create mode 100644 src/sage/matrix/matrix_numpy_dense.pxd create mode 100644 src/sage/matrix/matrix_numpy_dense.pyx create mode 100644 src/sage/matrix/matrix_numpy_integer_dense.pxd create mode 100644 src/sage/matrix/matrix_numpy_integer_dense.pyx diff --git a/src/sage/matrix/matrix_double_dense.pxd b/src/sage/matrix/matrix_double_dense.pxd index f30613822d2..db1feffa053 100644 --- a/src/sage/matrix/matrix_double_dense.pxd +++ b/src/sage/matrix/matrix_double_dense.pxd @@ -1,12 +1,6 @@ -from .matrix_dense cimport Matrix_dense -cimport numpy as cnumpy +from .matrix_numpy_dense cimport Matrix_numpy_dense -cdef class Matrix_double_dense(Matrix_dense): - cdef object _numpy_dtype - # cdef cnumpy.NPY_TYPES _numpy_dtypeint - cdef int _numpy_dtypeint - cdef object _python_dtype - cdef object _sage_dtype - cdef object _sage_vector_dtype - cdef Matrix_double_dense _new(self, int nrows=*, int ncols=*) - cdef cnumpy.ndarray _matrix_numpy + +cdef class Matrix_double_dense(Matrix_numpy_dense): + + pass diff --git a/src/sage/matrix/matrix_double_dense.pyx b/src/sage/matrix/matrix_double_dense.pyx index fea7a93f88c..84658065c44 100644 --- a/src/sage/matrix/matrix_double_dense.pyx +++ b/src/sage/matrix/matrix_double_dense.pyx @@ -68,15 +68,9 @@ import math import sage.rings.real_double import sage.rings.complex_double -from .matrix cimport Matrix -from .args cimport MatrixArgs_init -from sage.structure.coerce cimport coercion_model -from sage.structure.element import is_Matrix -from sage.structure.element cimport ModuleElement,Vector +from sage.structure.element cimport Vector from .constructor import matrix cimport sage.structure.element -from .matrix_space import MatrixSpace -from sage.misc.decorators import rename_keyword cimport numpy as cnumpy @@ -87,7 +81,7 @@ scipy = None cnumpy.import_array() -cdef class Matrix_double_dense(Matrix_dense): +cdef class Matrix_double_dense(Matrix_numpy_dense): """ Base class for matrices over the Real Double Field and the Complex Double Field. These are supposed to be fast matrix operations @@ -131,26 +125,6 @@ cdef class Matrix_double_dense(Matrix_dense): 6694819972852100501 # 64-bit 1829383573 # 32-bit """ - def __create_matrix__(self): - """ - Create a new uninitialized numpy matrix to hold the data for the class. - - This function assumes that self._numpy_dtypeint and - self._nrows and self._ncols have already been initialized. - - EXAMPLES: - - In this example, we throw away the current matrix and make a - new uninitialized matrix representing the data for the class:: - - sage: a = matrix(RDF, 3, range(9)) - sage: a.__create_matrix__() - """ - cdef cnumpy.npy_intp dims[2] - dims[0] = self._nrows - dims[1] = self._ncols - self._matrix_numpy = cnumpy.PyArray_SimpleNew(2, dims, self._numpy_dtypeint) - return def LU_valid(self): r""" @@ -167,133 +141,6 @@ cdef class Matrix_double_dense(Matrix_dense): """ return self.fetch('PLU_factors') is not None - def __init__(self, parent, entries=None, copy=None, bint coerce=True): - r""" - Fill the matrix with entries. - - The numpy matrix must have already been allocated. - - INPUT: - - - ``parent`` -- a matrix space over ``RDF`` - - - ``entries`` -- see :func:`matrix` - - - ``copy`` -- ignored (for backwards compatibility) - - - ``coerce`` -- if ``True`` (the default), convert elements to the - base ring before passing them to NumPy. If ``False``, pass the - elements to NumPy as given. - - EXAMPLES:: - - sage: matrix(RDF,3,range(9)) - [0.0 1.0 2.0] - [3.0 4.0 5.0] - [6.0 7.0 8.0] - sage: matrix(CDF,3,3,2) - [2.0 0.0 0.0] - [0.0 2.0 0.0] - [0.0 0.0 2.0] - - TESTS:: - - sage: matrix(RDF,3,0) - [] - sage: matrix(RDF,3,3,0) - [0.0 0.0 0.0] - [0.0 0.0 0.0] - [0.0 0.0 0.0] - sage: matrix(RDF,3,3,1) - [1.0 0.0 0.0] - [0.0 1.0 0.0] - [0.0 0.0 1.0] - sage: matrix(RDF,3,3,2) - [2.0 0.0 0.0] - [0.0 2.0 0.0] - [0.0 0.0 2.0] - sage: matrix(CDF,3,0) - [] - sage: matrix(CDF,3,3,0) - [0.0 0.0 0.0] - [0.0 0.0 0.0] - [0.0 0.0 0.0] - sage: matrix(CDF,3,3,1) - [1.0 0.0 0.0] - [0.0 1.0 0.0] - [0.0 0.0 1.0] - sage: matrix(CDF,3,3,range(9)) - [0.0 1.0 2.0] - [3.0 4.0 5.0] - [6.0 7.0 8.0] - sage: matrix(CDF,2,2,[CDF(1+I)*j for j in range(4)]) - [ 0.0 1.0 + 1.0*I] - [2.0 + 2.0*I 3.0 + 3.0*I] - """ - ma = MatrixArgs_init(parent, entries) - cdef long i, j - it = ma.iter(coerce) - for i in range(ma.nrows): - for j in range(ma.ncols): - self.set_unsafe(i, j, next(it)) - - cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, object value): - """ - Set the (i,j) entry to value without any bounds checking, - mutability checking, etc. - """ - # We assume that Py_ssize_t is the same as cnumpy.npy_intp - - # We must patch the ndarrayobject.h file so that the SETITEM - # macro does not have a semicolon at the end for this to work. - # Cython wraps the macro in a function that converts the - # returned int to a python object, which leads to compilation - # errors because after preprocessing you get something that - # looks like "););". This is bug - # http://scipy.org/scipy/numpy/ticket/918 - - # We call the self._python_dtype function on the value since - # numpy does not know how to deal with complex numbers other - # than the built-in complex number type. - cdef int status - status = cnumpy.PyArray_SETITEM(self._matrix_numpy, - cnumpy.PyArray_GETPTR2(self._matrix_numpy, i, j), - self._python_dtype(value)) - #TODO: Throw an error if status == -1 - - cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): - """ - Get the (i,j) entry without any bounds checking, etc. - """ - # We assume that Py_ssize_t is the same as cnumpy.npy_intp - return self._sage_dtype(cnumpy.PyArray_GETITEM(self._matrix_numpy, - cnumpy.PyArray_GETPTR2(self._matrix_numpy, i, j))) - - cdef Matrix_double_dense _new(self, int nrows=-1, int ncols=-1): - """ - Return a new uninitialized matrix with same parent as ``self``. - - INPUT: - - - nrows -- (default self._nrows) number of rows in returned matrix - - ncols -- (default self._ncols) number of columns in returned matrix - - """ - cdef Matrix_double_dense m - if nrows == -1 and ncols == -1: - nrows = self._nrows - ncols = self._ncols - parent = self._parent - else: - if nrows == -1: - nrows = self._nrows - if ncols == -1: - ncols = self._ncols - parent = self.matrix_space(nrows, ncols) - m = self.__class__.__new__(self.__class__, parent, None, None, None) - return m - - ######################################################################## # LEVEL 2 functionality # * def _pickle @@ -500,43 +347,6 @@ cdef class Matrix_double_dense(Matrix_dense): raise ZeroDivisionError("input matrix must be nonsingular") return M - def __copy__(self): - r""" - Return a new copy of this matrix. - - EXAMPLES:: - - sage: a = matrix(RDF,1,3, [1,2,-3]) - sage: a - [ 1.0 2.0 -3.0] - sage: b = a.__copy__() - sage: b - [ 1.0 2.0 -3.0] - sage: b is a - False - sage: b == a - True - sage: b[0,0] = 3 - sage: a[0,0] # note that a hasn't changed - 1.0 - - :: - - sage: copy(MatrixSpace(RDF,0,0,sparse=False).zero_matrix()) - [] - """ - if self._nrows == 0 or self._ncols == 0: - # Create a brand new empty matrix. This is needed to prevent a - # recursive loop: a copy of zero_matrix is asked otherwise. - return self.__class__(self.parent(), [], self._nrows, self._ncols) - - cdef Matrix_double_dense A - A = self._new(self._nrows, self._ncols) - A._matrix_numpy = self._matrix_numpy.copy() - if self._subdivisions is not None: - A.subdivide(*self.subdivisions()) - return A - # def _list(self): # def _dict(self): @@ -2014,48 +1824,6 @@ cdef class Matrix_double_dense(Matrix_dense): return sage.rings.real_double.RDF(sum(numpy.log(abs(numpy.diag(U._matrix_numpy))))) - def transpose(self): - """ - Return the transpose of this matrix, without changing self. - - EXAMPLES:: - - sage: m = matrix(RDF,2,3,range(6)); m - [0.0 1.0 2.0] - [3.0 4.0 5.0] - sage: m2 = m.transpose() - sage: m[0,0] = 2 - sage: m2 #note that m2 hasn't changed - [0.0 3.0] - [1.0 4.0] - [2.0 5.0] - - ``.T`` is a convenient shortcut for the transpose:: - - sage: m.T - [2.0 3.0] - [1.0 4.0] - [2.0 5.0] - - sage: m = matrix(RDF,0,3); m - [] - sage: m.transpose() - [] - sage: m.transpose().parent() - Full MatrixSpace of 3 by 0 dense matrices over Real Double Field - - """ - if self._nrows == 0 or self._ncols == 0: - return self.new_matrix(self._ncols, self._nrows) - - cdef Matrix_double_dense trans - trans = self._new(self._ncols, self._nrows) - trans._matrix_numpy = self._matrix_numpy.transpose().copy() - if self._subdivisions is not None: - row_divs, col_divs = self.subdivisions() - trans.subdivide(col_divs, row_divs) - return trans - def conjugate(self): r""" Return the conjugate of this matrix, i.e. the matrix whose entries are @@ -2429,55 +2197,6 @@ cdef class Matrix_double_dense(Matrix_dense): self.cache('QR_factors', QR) return QR - def is_symmetric(self, tol = 1e-12): - """ - Return whether this matrix is symmetric, to the given tolerance. - - EXAMPLES:: - - sage: m = matrix(RDF,2,2,range(4)); m - [0.0 1.0] - [2.0 3.0] - sage: m.is_symmetric() - False - sage: m[1,0] = 1.1; m - [0.0 1.0] - [1.1 3.0] - sage: m.is_symmetric() - False - - The tolerance inequality is strict: - sage: m.is_symmetric(tol=0.1) - False - sage: m.is_symmetric(tol=0.11) - True - - TESTS: - - Complex entries are supported (:trac:`27831`). :: - - sage: a = matrix(CDF, [(21, 0.6 + 18.5*i), (0.6 - 18.5*i, 21)]) - sage: a.is_symmetric() - False - """ - cdef Py_ssize_t i, j - tol = float(tol) - key = 'symmetric_%s'%tol - b = self.fetch(key) - if not b is None: - return b - if self._nrows != self._ncols: - self.cache(key, False) - return False - b = True - for i from 0 < i < self._nrows: - for j from 0 <= j < i: - if abs(self.get_unsafe(i,j) - self.get_unsafe(j,i)) > tol: - b = False - break - self.cache(key, b) - return b - def is_unitary(self, tol=1e-12, algorithm='orthonormal'): r""" Return ``True`` if the columns of the matrix are an orthonormal basis. @@ -2668,42 +2387,6 @@ cdef class Matrix_double_dense(Matrix_dense): self.cache(key, unitary) return unitary - def _is_lower_triangular(self, tol): - r""" - Return ``True`` if the entries above the diagonal are all zero. - - INPUT: - - - ``tol`` - the largest value of the absolute value of the - difference between two matrix entries for which they will - still be considered equal. - - OUTPUT: - - Returns ``True`` if each entry above the diagonal (entries - with a row index less than the column index) is zero. - - EXAMPLES:: - - sage: A = matrix(RDF, [[ 2.0, 0.0, 0.0], - ....: [ 1.0, 3.0, 0.0], - ....: [-4.0, 2.0, -1.0]]) - sage: A._is_lower_triangular(1.0e-17) - True - sage: A[1,2] = 10^-13 - sage: A._is_lower_triangular(1.0e-14) - False - """ - global numpy - if numpy is None: - import numpy - cdef Py_ssize_t i, j - for i in range(self._nrows): - for j in range(i+1, self._ncols): - if abs(self.get_unsafe(i,j)) > tol: - return False - return True - def _is_hermitian(self, tol = 1e-12, algorithm='orthonormal', skew=False): r""" Return ``True`` if the matrix is (skew-)Hermitian. @@ -3960,88 +3643,6 @@ cdef class Matrix_double_dense(Matrix_dense): ans = numpy.dot(self._matrix_numpy, v_numpy) return M(ans) - def numpy(self, dtype=None): - """ - This method returns a copy of the matrix as a numpy array. It - uses the numpy C/api so is very fast. - - INPUT: - - - ``dtype`` - The desired data-type for the array. If not given, - then the type will be determined as the minimum type required - to hold the objects in the sequence. - - EXAMPLES:: - - sage: m = matrix(RDF,[[1,2],[3,4]]) - sage: n = m.numpy() - sage: import numpy - sage: numpy.linalg.eig(n) - (array([-0.37228132, 5.37228132]), array([[-0.82456484, -0.41597356], - [ 0.56576746, -0.90937671]])) - sage: m = matrix(RDF, 2, range(6)); m - [0.0 1.0 2.0] - [3.0 4.0 5.0] - sage: m.numpy() - array([[0., 1., 2.], - [3., 4., 5.]]) - - Alternatively, numpy automatically calls this function (via - the magic :meth:`__array__` method) to convert Sage matrices - to numpy arrays:: - - sage: import numpy - sage: m = matrix(RDF, 2, range(6)); m - [0.0 1.0 2.0] - [3.0 4.0 5.0] - sage: numpy.array(m) - array([[0., 1., 2.], - [3., 4., 5.]]) - sage: numpy.array(m).dtype - dtype('float64') - sage: m = matrix(CDF, 2, range(6)); m - [0.0 1.0 2.0] - [3.0 4.0 5.0] - sage: numpy.array(m) - array([[0.+0.j, 1.+0.j, 2.+0.j], - [3.+0.j, 4.+0.j, 5.+0.j]]) - sage: numpy.array(m).dtype - dtype('complex128') - - TESTS:: - - sage: m = matrix(RDF,0,5,[]); m - [] - sage: m.numpy() - array([], shape=(0, 5), dtype=float64) - sage: m = matrix(RDF,5,0,[]); m - [] - sage: m.numpy() - array([], shape=(5, 0), dtype=float64) - """ - import numpy as np - if dtype is None or self._numpy_dtype == np.dtype(dtype): - return self._matrix_numpy.copy() - else: - return Matrix_dense.numpy(self, dtype=dtype) - - def _replace_self_with_numpy(self,numpy_matrix): - """ - - EXAMPLES:: - - sage: import numpy - sage: a = numpy.array([[1,2],[3,4]], 'float64') - sage: m = matrix(RDF,2,2,0) - sage: m._replace_self_with_numpy(a) - sage: m - [1.0 2.0] - [3.0 4.0] - """ - if (self._matrix_numpy).shape != (numpy_matrix).shape: - raise ValueError("matrix shapes are not the same") - self._matrix_numpy = numpy_matrix.astype(self._numpy_dtype) - def _replace_self_with_numpy32(self,numpy_matrix): """ diff --git a/src/sage/matrix/matrix_numpy_dense.pxd b/src/sage/matrix/matrix_numpy_dense.pxd new file mode 100644 index 00000000000..a0ec36c9228 --- /dev/null +++ b/src/sage/matrix/matrix_numpy_dense.pxd @@ -0,0 +1,12 @@ +from .matrix_dense cimport Matrix_dense +cimport numpy as cnumpy + + +cdef class Matrix_numpy_dense(Matrix_dense): + cdef object _numpy_dtype + cdef int _numpy_dtypeint + cdef object _python_dtype + cdef object _sage_dtype + cdef object _sage_vector_dtype + cdef Matrix_numpy_dense _new(self, int nrows=*, int ncols=*) + cdef cnumpy.ndarray _matrix_numpy diff --git a/src/sage/matrix/matrix_numpy_dense.pyx b/src/sage/matrix/matrix_numpy_dense.pyx new file mode 100644 index 00000000000..5047d423c94 --- /dev/null +++ b/src/sage/matrix/matrix_numpy_dense.pyx @@ -0,0 +1,448 @@ +""" +Dense matrices using a NumPy backend + +AUTHORS: + +- Jason Grout, Sep 2008: switch to NumPy backend, factored out the Matrix_double_dense class + +- Josh Kantor + +- William Stein: many bug fixes and touch ups. +""" + +# **************************************************************************** +# Copyright (C) 2004-2006 Joshua Kantor +# Copyright (C) 2008 Georg S. Weber +# Copyright (C) 2008-2011 Mike Hansen +# Copyright (C) 2008-2012 Jason Grout +# Copyright (C) 2009 Dag Sverre Seljebotn +# Copyright (C) 2009 Yann Laigle-Chapuy +# Copyright (C) 2009-2010 Florent Hivert +# Copyright (C) 2010-2012 Rob Beezer +# Copyright (C) 2011 Martin Raum +# Copyright (C) 2011-2012 J. H. Palmieri +# Copyright (C) 2011-2014 André Apitzsch +# Copyright (C) 2011-2018 Jeroen Demeyer +# Copyright (C) 2012 Kenneth Smith +# Copyright (C) 2016-2019 Frédéric Chapoton +# Copyright (C) 2017 Kiran Kedlaya +# Copyright (C) 2019 Chaman Agrawal +# Copyright (C) 2019-2021 Markus Wageringel +# Copyright (C) 2020 Michael Orlitzky +# Copyright (C) 2020 Victor Santos +# Copyright (C) 2021 Jonathan Kliem +# Copyright (C) 2021 Travis Scrimshaw +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +from .matrix cimport Matrix +from .args cimport MatrixArgs_init +cimport sage.structure.element + +cimport numpy as cnumpy + +numpy = None +scipy = None + +# This is for the Numpy C API to work +cnumpy.import_array() + + +cdef class Matrix_numpy_dense(Matrix_dense): + + def __create_matrix__(self): + """ + Create a new uninitialized numpy matrix to hold the data for the class. + + This function assumes that self._numpy_dtypeint and + self._nrows and self._ncols have already been initialized. + + EXAMPLES: + + In this example, we throw away the current matrix and make a + new uninitialized matrix representing the data for the class:: + + sage: a = matrix(RDF, 3, range(9)) + sage: a.__create_matrix__() + """ + cdef cnumpy.npy_intp dims[2] + dims[0] = self._nrows + dims[1] = self._ncols + self._matrix_numpy = cnumpy.PyArray_SimpleNew(2, dims, self._numpy_dtypeint) + return + + def __init__(self, parent, entries=None, copy=None, bint coerce=True): + r""" + Fill the matrix with entries. + + The numpy matrix must have already been allocated. + + INPUT: + + - ``parent`` -- a matrix space over ``RDF`` + + - ``entries`` -- see :func:`matrix` + + - ``copy`` -- ignored (for backwards compatibility) + + - ``coerce`` -- if ``True`` (the default), convert elements to the + base ring before passing them to NumPy. If ``False``, pass the + elements to NumPy as given. + + EXAMPLES:: + + sage: matrix(RDF,3,range(9)) + [0.0 1.0 2.0] + [3.0 4.0 5.0] + [6.0 7.0 8.0] + sage: matrix(CDF,3,3,2) + [2.0 0.0 0.0] + [0.0 2.0 0.0] + [0.0 0.0 2.0] + + TESTS:: + + sage: matrix(RDF,3,0) + [] + sage: matrix(RDF,3,3,0) + [0.0 0.0 0.0] + [0.0 0.0 0.0] + [0.0 0.0 0.0] + sage: matrix(RDF,3,3,1) + [1.0 0.0 0.0] + [0.0 1.0 0.0] + [0.0 0.0 1.0] + sage: matrix(RDF,3,3,2) + [2.0 0.0 0.0] + [0.0 2.0 0.0] + [0.0 0.0 2.0] + sage: matrix(CDF,3,0) + [] + sage: matrix(CDF,3,3,0) + [0.0 0.0 0.0] + [0.0 0.0 0.0] + [0.0 0.0 0.0] + sage: matrix(CDF,3,3,1) + [1.0 0.0 0.0] + [0.0 1.0 0.0] + [0.0 0.0 1.0] + sage: matrix(CDF,3,3,range(9)) + [0.0 1.0 2.0] + [3.0 4.0 5.0] + [6.0 7.0 8.0] + sage: matrix(CDF,2,2,[CDF(1+I)*j for j in range(4)]) + [ 0.0 1.0 + 1.0*I] + [2.0 + 2.0*I 3.0 + 3.0*I] + """ + ma = MatrixArgs_init(parent, entries) + cdef long i, j + it = ma.iter(coerce) + for i in range(ma.nrows): + for j in range(ma.ncols): + self.set_unsafe(i, j, next(it)) + + cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, object value): + """ + Set the (i,j) entry to value without any bounds checking, + mutability checking, etc. + """ + # We assume that Py_ssize_t is the same as cnumpy.npy_intp + + # We must patch the ndarrayobject.h file so that the SETITEM + # macro does not have a semicolon at the end for this to work. + # Cython wraps the macro in a function that converts the + # returned int to a python object, which leads to compilation + # errors because after preprocessing you get something that + # looks like "););". This is bug + # http://scipy.org/scipy/numpy/ticket/918 + + # We call the self._python_dtype function on the value since + # numpy does not know how to deal with complex numbers other + # than the built-in complex number type. + cdef int status + status = cnumpy.PyArray_SETITEM(self._matrix_numpy, + cnumpy.PyArray_GETPTR2(self._matrix_numpy, i, j), + self._python_dtype(value)) + #TODO: Throw an error if status == -1 + + cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): + """ + Get the (i,j) entry without any bounds checking, etc. + """ + # We assume that Py_ssize_t is the same as cnumpy.npy_intp + return self._sage_dtype(cnumpy.PyArray_GETITEM(self._matrix_numpy, + cnumpy.PyArray_GETPTR2(self._matrix_numpy, i, j))) + + cdef Matrix_numpy_dense _new(self, int nrows=-1, int ncols=-1): + """ + Return a new uninitialized matrix with same parent as ``self``. + + INPUT: + + - nrows -- (default self._nrows) number of rows in returned matrix + - ncols -- (default self._ncols) number of columns in returned matrix + + """ + cdef Matrix_numpy_dense m + if nrows == -1 and ncols == -1: + nrows = self._nrows + ncols = self._ncols + parent = self._parent + else: + if nrows == -1: + nrows = self._nrows + if ncols == -1: + ncols = self._ncols + parent = self.matrix_space(nrows, ncols) + m = self.__class__.__new__(self.__class__, parent, None, None, None) + return m + + def __copy__(self): + r""" + Return a new copy of this matrix. + + EXAMPLES:: + + sage: a = matrix(RDF,1,3, [1,2,-3]) + sage: a + [ 1.0 2.0 -3.0] + sage: b = a.__copy__() + sage: b + [ 1.0 2.0 -3.0] + sage: b is a + False + sage: b == a + True + sage: b[0,0] = 3 + sage: a[0,0] # note that a hasn't changed + 1.0 + + :: + + sage: copy(MatrixSpace(RDF,0,0,sparse=False).zero_matrix()) + [] + """ + if self._nrows == 0 or self._ncols == 0: + # Create a brand new empty matrix. This is needed to prevent a + # recursive loop: a copy of zero_matrix is asked otherwise. + return self.__class__(self.parent(), [], self._nrows, self._ncols) + + cdef Matrix_numpy_dense A + A = self._new(self._nrows, self._ncols) + A._matrix_numpy = self._matrix_numpy.copy() + if self._subdivisions is not None: + A.subdivide(*self.subdivisions()) + return A + + def transpose(self): + """ + Return the transpose of this matrix, without changing self. + + EXAMPLES:: + + sage: m = matrix(RDF,2,3,range(6)); m + [0.0 1.0 2.0] + [3.0 4.0 5.0] + sage: m2 = m.transpose() + sage: m[0,0] = 2 + sage: m2 #note that m2 hasn't changed + [0.0 3.0] + [1.0 4.0] + [2.0 5.0] + + ``.T`` is a convenient shortcut for the transpose:: + + sage: m.T + [2.0 3.0] + [1.0 4.0] + [2.0 5.0] + + sage: m = matrix(RDF,0,3); m + [] + sage: m.transpose() + [] + sage: m.transpose().parent() + Full MatrixSpace of 3 by 0 dense matrices over Real Double Field + + """ + if self._nrows == 0 or self._ncols == 0: + return self.new_matrix(self._ncols, self._nrows) + + cdef Matrix_numpy_dense trans + trans = self._new(self._ncols, self._nrows) + trans._matrix_numpy = self._matrix_numpy.transpose().copy() + if self._subdivisions is not None: + row_divs, col_divs = self.subdivisions() + trans.subdivide(col_divs, row_divs) + return trans + + def is_symmetric(self, tol = 1e-12): + """ + Return whether this matrix is symmetric, to the given tolerance. + + EXAMPLES:: + + sage: m = matrix(RDF,2,2,range(4)); m + [0.0 1.0] + [2.0 3.0] + sage: m.is_symmetric() + False + sage: m[1,0] = 1.1; m + [0.0 1.0] + [1.1 3.0] + sage: m.is_symmetric() + False + + The tolerance inequality is strict: + sage: m.is_symmetric(tol=0.1) + False + sage: m.is_symmetric(tol=0.11) + True + + TESTS: + + Complex entries are supported (:trac:`27831`). :: + + sage: a = matrix(CDF, [(21, 0.6 + 18.5*i), (0.6 - 18.5*i, 21)]) + sage: a.is_symmetric() + False + """ + cdef Py_ssize_t i, j + tol = float(tol) + key = 'symmetric_%s'%tol + b = self.fetch(key) + if not b is None: + return b + if self._nrows != self._ncols: + self.cache(key, False) + return False + b = True + for i from 0 < i < self._nrows: + for j from 0 <= j < i: + if abs(self.get_unsafe(i,j) - self.get_unsafe(j,i)) > tol: + b = False + break + self.cache(key, b) + return b + + def _is_lower_triangular(self, tol): + r""" + Return ``True`` if the entries above the diagonal are all zero. + + INPUT: + + - ``tol`` - the largest value of the absolute value of the + difference between two matrix entries for which they will + still be considered equal. + + OUTPUT: + + Returns ``True`` if each entry above the diagonal (entries + with a row index less than the column index) is zero. + + EXAMPLES:: + + sage: A = matrix(RDF, [[ 2.0, 0.0, 0.0], + ....: [ 1.0, 3.0, 0.0], + ....: [-4.0, 2.0, -1.0]]) + sage: A._is_lower_triangular(1.0e-17) + True + sage: A[1,2] = 10^-13 + sage: A._is_lower_triangular(1.0e-14) + False + """ + global numpy + if numpy is None: + import numpy + cdef Py_ssize_t i, j + for i in range(self._nrows): + for j in range(i+1, self._ncols): + if abs(self.get_unsafe(i,j)) > tol: + return False + return True + + def numpy(self, dtype=None): + """ + This method returns a copy of the matrix as a numpy array. It + uses the numpy C/api so is very fast. + + INPUT: + + - ``dtype`` - The desired data-type for the array. If not given, + then the type will be determined as the minimum type required + to hold the objects in the sequence. + + EXAMPLES:: + + sage: m = matrix(RDF,[[1,2],[3,4]]) + sage: n = m.numpy() + sage: import numpy + sage: numpy.linalg.eig(n) + (array([-0.37228132, 5.37228132]), array([[-0.82456484, -0.41597356], + [ 0.56576746, -0.90937671]])) + sage: m = matrix(RDF, 2, range(6)); m + [0.0 1.0 2.0] + [3.0 4.0 5.0] + sage: m.numpy() + array([[0., 1., 2.], + [3., 4., 5.]]) + + Alternatively, numpy automatically calls this function (via + the magic :meth:`__array__` method) to convert Sage matrices + to numpy arrays:: + + sage: import numpy + sage: m = matrix(RDF, 2, range(6)); m + [0.0 1.0 2.0] + [3.0 4.0 5.0] + sage: numpy.array(m) + array([[0., 1., 2.], + [3., 4., 5.]]) + sage: numpy.array(m).dtype + dtype('float64') + sage: m = matrix(CDF, 2, range(6)); m + [0.0 1.0 2.0] + [3.0 4.0 5.0] + sage: numpy.array(m) + array([[0.+0.j, 1.+0.j, 2.+0.j], + [3.+0.j, 4.+0.j, 5.+0.j]]) + sage: numpy.array(m).dtype + dtype('complex128') + + TESTS:: + + sage: m = matrix(RDF,0,5,[]); m + [] + sage: m.numpy() + array([], shape=(0, 5), dtype=float64) + sage: m = matrix(RDF,5,0,[]); m + [] + sage: m.numpy() + array([], shape=(5, 0), dtype=float64) + """ + import numpy as np + if dtype is None or self._numpy_dtype == np.dtype(dtype): + return self._matrix_numpy.copy() + else: + return Matrix_dense.numpy(self, dtype=dtype) + + def _replace_self_with_numpy(self,numpy_matrix): + """ + + EXAMPLES:: + + sage: import numpy + sage: a = numpy.array([[1,2],[3,4]], 'float64') + sage: m = matrix(RDF,2,2,0) + sage: m._replace_self_with_numpy(a) + sage: m + [1.0 2.0] + [3.0 4.0] + """ + if (self._matrix_numpy).shape != (numpy_matrix).shape: + raise ValueError("matrix shapes are not the same") + self._matrix_numpy = numpy_matrix.astype(self._numpy_dtype) diff --git a/src/sage/matrix/matrix_numpy_integer_dense.pxd b/src/sage/matrix/matrix_numpy_integer_dense.pxd new file mode 100644 index 00000000000..c74a2c7511c --- /dev/null +++ b/src/sage/matrix/matrix_numpy_integer_dense.pxd @@ -0,0 +1,6 @@ +from .matrix_numpy_dense cimport Matrix_numpy_dense + + +cdef class Matrix_numpy_integer_dense(Matrix_numpy_dense): + + pass diff --git a/src/sage/matrix/matrix_numpy_integer_dense.pyx b/src/sage/matrix/matrix_numpy_integer_dense.pyx new file mode 100644 index 00000000000..94919b5b428 --- /dev/null +++ b/src/sage/matrix/matrix_numpy_integer_dense.pyx @@ -0,0 +1,58 @@ +r""" +Dense integer matrices using a NumPy backend + + +EXAMPLES:: + + sage: from sage.matrix.matrix_numpy_integer_dense import Matrix_numpy_integer_dense + sage: M = Matrix_numpy_integer_dense(MatrixSpace(ZZ, 2, 3)); M + [0 0 0] + [0 0 0] + sage: M[1,2] = 47 + sage: M + [ 0 0 0] + [ 0 0 47] + sage: M.numpy() + array([[ 0, 0, 0], + [ 0, 0, 47]]) +""" + +# **************************************************************************** +# Copyright (C) 2021 Matthias Koeppe +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +from sage.rings.integer_ring import ZZ + +cimport numpy as cnumpy + +numpy = None +scipy = None + +# This is for the Numpy C API to work +cnumpy.import_array() + + +cdef class Matrix_numpy_integer_dense(Matrix_numpy_dense): + r""" + TESTS:: + + sage: from sage.matrix.matrix_numpy_integer_dense import Matrix_numpy_integer_dense + sage: M = Matrix_numpy_integer_dense(MatrixSpace(ZZ, 2, 3)); M + sage: TestSuite(M).run() + """ + + def __cinit__(self): + global numpy + if numpy is None: + import numpy + self._numpy_dtype = numpy.dtype('int64') + self._numpy_dtypeint = cnumpy.NPY_INT64 + self._python_dtype = int + self._sage_dtype = ZZ + self.__create_matrix__() + return From 98ced3cbac17fc6ef5ffc7831693831f647912b9 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 12 Oct 2021 12:39:34 -0700 Subject: [PATCH 3/9] src/sage/matrix/matrix_numpy_integer_dense.pyx: Fix doctest --- src/sage/matrix/matrix_numpy_integer_dense.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_numpy_integer_dense.pyx b/src/sage/matrix/matrix_numpy_integer_dense.pyx index 94919b5b428..e832386f618 100644 --- a/src/sage/matrix/matrix_numpy_integer_dense.pyx +++ b/src/sage/matrix/matrix_numpy_integer_dense.pyx @@ -42,7 +42,7 @@ cdef class Matrix_numpy_integer_dense(Matrix_numpy_dense): TESTS:: sage: from sage.matrix.matrix_numpy_integer_dense import Matrix_numpy_integer_dense - sage: M = Matrix_numpy_integer_dense(MatrixSpace(ZZ, 2, 3)); M + sage: M = Matrix_numpy_integer_dense(MatrixSpace(ZZ, 2, 3)) sage: TestSuite(M).run() """ From cde4e5463a5b0bbc08e9cbcdd181d6f9488683b5 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 12 Oct 2021 14:35:26 -0700 Subject: [PATCH 4/9] src/sage/modules/vector_double_dense.pyx: Update copyright according to 'git blame -w --date=format:%Y src/sage/modules/vector_double_dense.pyx | sort -k2' --- src/sage/modules/vector_double_dense.pyx | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/sage/modules/vector_double_dense.pyx b/src/sage/modules/vector_double_dense.pyx index 83334163bea..668929037cf 100644 --- a/src/sage/modules/vector_double_dense.pyx +++ b/src/sage/modules/vector_double_dense.pyx @@ -32,13 +32,20 @@ AUTHORS: """ #***************************************************************************** -# Copyright (C) 2006 William Stein +# Copyright (C) 2006-2010 William Stein +# Copyright (C) 2009 Alexandru Ghitza +# Copyright (C) 2020 Antonio Rojas +# Copyright (C) 2017 Frédéric Chapoton +# Copyright (C) 2008-2009 Jason Grout +# Copyright (C) 2014-2016 Jeroen Demeyer +# Copyright (C) 2011 Mike Hansen +# Copyright (C) 2011 Rob Beezer # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. -# http://www.gnu.org/licenses/ +# https://www.gnu.org/licenses/ #***************************************************************************** cimport numpy From 2ffb118b1aa1e5b001f128af0783f163dd06cb9c Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 12 Oct 2021 15:23:45 -0700 Subject: [PATCH 5/9] Vector_double_dense: Factor through new class Vector_numpy_dense --- src/sage/modules/vector_double_dense.pxd | 14 +- src/sage/modules/vector_double_dense.pyx | 258 +------------------ src/sage/modules/vector_numpy_dense.pxd | 12 + src/sage/modules/vector_numpy_dense.pyx | 302 +++++++++++++++++++++++ 4 files changed, 319 insertions(+), 267 deletions(-) create mode 100644 src/sage/modules/vector_numpy_dense.pxd create mode 100644 src/sage/modules/vector_numpy_dense.pyx diff --git a/src/sage/modules/vector_double_dense.pxd b/src/sage/modules/vector_double_dense.pxd index fbe3e04f25c..f29b4597d15 100644 --- a/src/sage/modules/vector_double_dense.pxd +++ b/src/sage/modules/vector_double_dense.pxd @@ -1,12 +1,4 @@ -from .free_module_element cimport FreeModuleElement -cimport numpy +from .vector_numpy_dense cimport Vector_numpy_dense -cdef class Vector_double_dense(FreeModuleElement): - cdef object _numpy_dtype - cdef int _numpy_dtypeint - cdef object _python_dtype - cdef object _sage_dtype - cdef object _sage_vector_dtype - cdef numpy.ndarray _vector_numpy - cdef Vector_double_dense _new(self, numpy.ndarray vector_numpy) - cdef _replace_self_with_numpy(self, numpy.ndarray numpy_array) +cdef class Vector_double_dense(Vector_numpy_dense): + pass diff --git a/src/sage/modules/vector_double_dense.pyx b/src/sage/modules/vector_double_dense.pyx index 668929037cf..9c4efc48b4c 100644 --- a/src/sage/modules/vector_double_dense.pyx +++ b/src/sage/modules/vector_double_dense.pyx @@ -50,9 +50,8 @@ AUTHORS: cimport numpy import numpy -from .free_module_element import FreeModuleElement -from sage.structure.element cimport Element, ModuleElement, RingElement, Vector +from sage.structure.element cimport Element, Vector from sage.rings.real_double import RDF from sage.rings.complex_double import CDF @@ -61,7 +60,7 @@ from sage.rings.complex_double import CDF numpy.import_array() -cdef class Vector_double_dense(FreeModuleElement): +cdef class Vector_double_dense(Vector_numpy_dense): """ Base class for vectors over the Real Double Field and the Complex Double Field. These are supposed to be fast vector operations @@ -79,204 +78,6 @@ cdef class Vector_double_dense(FreeModuleElement): sage: v*v 30.0 """ - def __cinit__(self, parent, entries, coerce=True, copy=True): - """ - Set up a new vector - - EXAMPLES:: - - sage: v = vector(RDF, range(3)) - sage: v.__new__(v.__class__, v.parent(), [1,1,1]) # random output - (3.00713073107e-261, 3.06320700422e-322, 2.86558074588e-322) - sage: v = vector(CDF, range(3)) - sage: v.__new__(v.__class__, v.parent(), [1,1,1]) # random output - (-2.26770549592e-39, 5.1698223615e-252*I, -5.9147262602e-62 + 4.63145528786e-258*I) - """ - self._is_immutable = 0 - self._degree = parent.degree() - self._parent = parent - - cdef Vector_double_dense _new(self, numpy.ndarray vector_numpy): - """ - Return a new vector with same parent as self. - """ - cdef Vector_double_dense v - v = self.__class__.__new__(self.__class__,self._parent,None,None,None) - v._is_immutable = 0 - v._parent = self._parent - v._degree = self._parent.degree() - - v._vector_numpy = vector_numpy - return v - - def __create_vector__(self): - """ - Create a new uninitialized numpy array to hold the data for the class. - - This function assumes that self._numpy_dtypeint and - self._nrows and self._ncols have already been initialized. - - EXAMPLES: - - In this example, we throw away the current array and make a - new uninitialized array representing the data for the class. :: - - sage: a=vector(RDF, range(9)) - sage: a.__create_vector__() - """ - cdef numpy.npy_intp dims[1] - dims[0] = self._degree - self._vector_numpy = numpy.PyArray_SimpleNew(1, dims, self._numpy_dtypeint) - return - - cdef bint is_dense_c(self): - """ - Return True (i.e., 1) if self is dense. - """ - return 1 - - cdef bint is_sparse_c(self): - """ - Return True (i.e., 1) if self is sparse. - """ - return 0 - - - def __copy__(self, copy=True): - """ - Return a copy of the vector - - EXAMPLES:: - - sage: a = vector(RDF, range(9)) - sage: a == copy(a) - True - """ - if self._degree == 0: - return self - from copy import copy - return self._new(copy(self._vector_numpy)) - - def __init__(self, parent, entries, coerce = True, copy = True): - """ - Fill the vector with entries. - - The numpy array must have already been allocated. - - EXAMPLES:: - - sage: vector(RDF, range(9)) - (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) - sage: vector(CDF, 5) - (0.0, 0.0, 0.0, 0.0, 0.0) - - TESTS:: - - sage: vector(CDF, 0) - () - sage: vector(RDF, 0) - () - sage: vector(CDF, 4) - (0.0, 0.0, 0.0, 0.0) - sage: vector(RDF, 4) - (0.0, 0.0, 0.0, 0.0) - sage: vector(CDF, [CDF(1+I)*j for j in range(4)]) - (0.0, 1.0 + 1.0*I, 2.0 + 2.0*I, 3.0 + 3.0*I) - sage: vector(RDF, 4, range(4)) - (0.0, 1.0, 2.0, 3.0) - - sage: V = RDF^2 - sage: V.element_class(V, 5) - Traceback (most recent call last): - ... - TypeError: entries must be a list or 0 - sage: V.element_class(V, 0) - (0.0, 0.0) - """ - cdef Py_ssize_t i,j - if isinstance(entries,(tuple, list)): - if len(entries)!=self._degree: - raise TypeError("entries has wrong length") - - if coerce: - for i from 0<=i`_ - of the returned array. - - EXAMPLES:: - - sage: v = vector(CDF,4,range(4)) - sage: v.numpy() - array([0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j]) - sage: v = vector(CDF,0) - sage: v.numpy() - array([], dtype=complex128) - sage: v = vector(RDF,4,range(4)) - sage: v.numpy() - array([0., 1., 2., 3.]) - sage: v = vector(RDF,0) - sage: v.numpy() - array([], dtype=float64) - - A numpy dtype may be requested manually:: - - sage: import numpy - sage: v = vector(CDF, 3, range(3)) - sage: v.numpy() - array([0.+0.j, 1.+0.j, 2.+0.j]) - sage: v.numpy(dtype=numpy.float64) - array([0., 1., 2.]) - sage: v.numpy(dtype=numpy.float32) - array([0., 1., 2.], dtype=float32) - """ - if dtype is None or dtype is self._vector_numpy.dtype: - from copy import copy - return copy(self._vector_numpy) - else: - return super(Vector_double_dense, self).numpy(dtype) diff --git a/src/sage/modules/vector_numpy_dense.pxd b/src/sage/modules/vector_numpy_dense.pxd new file mode 100644 index 00000000000..b019bc8ebac --- /dev/null +++ b/src/sage/modules/vector_numpy_dense.pxd @@ -0,0 +1,12 @@ +from .free_module_element cimport FreeModuleElement +cimport numpy + +cdef class Vector_numpy_dense(FreeModuleElement): + cdef object _numpy_dtype + cdef int _numpy_dtypeint + cdef object _python_dtype + cdef object _sage_dtype + cdef object _sage_vector_dtype + cdef numpy.ndarray _vector_numpy + cdef Vector_numpy_dense _new(self, numpy.ndarray vector_numpy) + cdef _replace_self_with_numpy(self, numpy.ndarray numpy_array) diff --git a/src/sage/modules/vector_numpy_dense.pyx b/src/sage/modules/vector_numpy_dense.pyx new file mode 100644 index 00000000000..1594f817c9c --- /dev/null +++ b/src/sage/modules/vector_numpy_dense.pyx @@ -0,0 +1,302 @@ +r""" +Dense vectors using a NumPy backend. + +AUTHORS: + +- Jason Grout, Oct 2008: switch to numpy backend, factored out + ``Vector_double_dense`` class +- Josh Kantor +- William Stein +""" + +#***************************************************************************** +# Copyright (C) 2006-2010 William Stein +# Copyright (C) 2009 Alexandru Ghitza +# Copyright (C) 2020 Antonio Rojas +# Copyright (C) 2017 Frédéric Chapoton +# Copyright (C) 2008-2009 Jason Grout +# Copyright (C) 2014-2016 Jeroen Demeyer +# Copyright (C) 2011 Mike Hansen +# Copyright (C) 2011 Rob Beezer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + +cimport numpy +import numpy +from .free_module_element import FreeModuleElement + +# This is for the NumPy C API (the PyArray... functions) to work +numpy.import_array() + + +cdef class Vector_numpy_dense(FreeModuleElement): + """ + Base class for vectors implemented using numpy arrays. + + This class cannot be instantiated on its own. The numpy vector + creation depends on several variables that are set in the + subclasses. + + EXAMPLES:: + + sage: v = vector(RDF, [1,2,3,4]); v + (1.0, 2.0, 3.0, 4.0) + sage: v*v + 30.0 + """ + + def __cinit__(self, parent, entries, coerce=True, copy=True): + """ + Set up a new vector + + EXAMPLES:: + + sage: v = vector(RDF, range(3)) + sage: v.__new__(v.__class__, v.parent(), [1,1,1]) # random output + (3.00713073107e-261, 3.06320700422e-322, 2.86558074588e-322) + sage: v = vector(CDF, range(3)) + sage: v.__new__(v.__class__, v.parent(), [1,1,1]) # random output + (-2.26770549592e-39, 5.1698223615e-252*I, -5.9147262602e-62 + 4.63145528786e-258*I) + """ + self._is_immutable = 0 + self._degree = parent.degree() + self._parent = parent + + cdef Vector_numpy_dense _new(self, numpy.ndarray vector_numpy): + """ + Return a new vector with same parent as self. + """ + cdef Vector_numpy_dense v + v = self.__class__.__new__(self.__class__,self._parent,None,None,None) + v._is_immutable = 0 + v._parent = self._parent + v._degree = self._parent.degree() + + v._vector_numpy = vector_numpy + return v + + def __create_vector__(self): + """ + Create a new uninitialized numpy array to hold the data for the class. + + This function assumes that self._numpy_dtypeint and + self._nrows and self._ncols have already been initialized. + + EXAMPLES: + + In this example, we throw away the current array and make a + new uninitialized array representing the data for the class. :: + + sage: a=vector(RDF, range(9)) + sage: a.__create_vector__() + """ + cdef numpy.npy_intp dims[1] + dims[0] = self._degree + self._vector_numpy = numpy.PyArray_SimpleNew(1, dims, self._numpy_dtypeint) + return + + cdef bint is_dense_c(self): + """ + Return True (i.e., 1) if self is dense. + """ + return 1 + + cdef bint is_sparse_c(self): + """ + Return True (i.e., 1) if self is sparse. + """ + return 0 + + + def __copy__(self, copy=True): + """ + Return a copy of the vector + + EXAMPLES:: + + sage: a = vector(RDF, range(9)) + sage: a == copy(a) + True + """ + if self._degree == 0: + return self + from copy import copy + return self._new(copy(self._vector_numpy)) + + def __init__(self, parent, entries, coerce = True, copy = True): + """ + Fill the vector with entries. + + The numpy array must have already been allocated. + + EXAMPLES:: + + sage: vector(RDF, range(9)) + (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) + sage: vector(CDF, 5) + (0.0, 0.0, 0.0, 0.0, 0.0) + + TESTS:: + + sage: vector(CDF, 0) + () + sage: vector(RDF, 0) + () + sage: vector(CDF, 4) + (0.0, 0.0, 0.0, 0.0) + sage: vector(RDF, 4) + (0.0, 0.0, 0.0, 0.0) + sage: vector(CDF, [CDF(1+I)*j for j in range(4)]) + (0.0, 1.0 + 1.0*I, 2.0 + 2.0*I, 3.0 + 3.0*I) + sage: vector(RDF, 4, range(4)) + (0.0, 1.0, 2.0, 3.0) + + sage: V = RDF^2 + sage: V.element_class(V, 5) + Traceback (most recent call last): + ... + TypeError: entries must be a list or 0 + sage: V.element_class(V, 0) + (0.0, 0.0) + """ + cdef Py_ssize_t i,j + if isinstance(entries,(tuple, list)): + if len(entries)!=self._degree: + raise TypeError("entries has wrong length") + + if coerce: + for i from 0<=i`_ + of the returned array. + + EXAMPLES:: + + sage: v = vector(CDF,4,range(4)) + sage: v.numpy() + array([0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j]) + sage: v = vector(CDF,0) + sage: v.numpy() + array([], dtype=complex128) + sage: v = vector(RDF,4,range(4)) + sage: v.numpy() + array([0., 1., 2., 3.]) + sage: v = vector(RDF,0) + sage: v.numpy() + array([], dtype=float64) + + A numpy dtype may be requested manually:: + + sage: import numpy + sage: v = vector(CDF, 3, range(3)) + sage: v.numpy() + array([0.+0.j, 1.+0.j, 2.+0.j]) + sage: v.numpy(dtype=numpy.float64) + array([0., 1., 2.]) + sage: v.numpy(dtype=numpy.float32) + array([0., 1., 2.], dtype=float32) + """ + if dtype is None or dtype is self._vector_numpy.dtype: + from copy import copy + return copy(self._vector_numpy) + else: + return super(Vector_numpy_dense, self).numpy(dtype) From 5b288c35d91cbbf1aaca136cfa6cb5a35b5bbe1f Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 12 Oct 2021 15:33:23 -0700 Subject: [PATCH 6/9] sage.modules.vector_numpy_integer_dense: New --- .../modules/vector_numpy_integer_dense.pxd | 5 +++ .../modules/vector_numpy_integer_dense.pyx | 44 +++++++++++++++++++ 2 files changed, 49 insertions(+) create mode 100644 src/sage/modules/vector_numpy_integer_dense.pxd create mode 100644 src/sage/modules/vector_numpy_integer_dense.pyx diff --git a/src/sage/modules/vector_numpy_integer_dense.pxd b/src/sage/modules/vector_numpy_integer_dense.pxd new file mode 100644 index 00000000000..e39b90bbcb9 --- /dev/null +++ b/src/sage/modules/vector_numpy_integer_dense.pxd @@ -0,0 +1,5 @@ +from .vector_numpy_dense cimport Vector_numpy_dense + +cdef class Vector_numpy_integer_dense(Vector_numpy_dense): + + pass diff --git a/src/sage/modules/vector_numpy_integer_dense.pyx b/src/sage/modules/vector_numpy_integer_dense.pyx new file mode 100644 index 00000000000..b73a21eb849 --- /dev/null +++ b/src/sage/modules/vector_numpy_integer_dense.pyx @@ -0,0 +1,44 @@ +r""" +Dense integer vectors using a NumPy backend. + +EXAMPLES:: + + sage: from sage.modules.vector_numpy_integer_dense import Vector_numpy_integer_dense + sage: v = Vector_numpy_integer_dense(FreeModule(ZZ, 3), [0, 0, 0]); v + (0, 0, 0) + sage: v[1] = 42 + sage: v + (0, 42, 0) + sage: v.numpy() + array([ 0, 42, 0]) +""" + +# **************************************************************************** +# Copyright (C) 2021 Matthias Koeppe +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +cimport numpy +import numpy + +from sage.structure.element cimport Element, Vector + +from sage.rings.integer_ring import ZZ + +# This is for the NumPy C API (the PyArray... functions) to work +numpy.import_array() + + +cdef class Vector_numpy_integer_dense(Vector_numpy_dense): + + def __cinit__(self, parent, entries, coerce=True, copy=True): + self._numpy_dtype = numpy.dtype('int64') + self._numpy_dtypeint = numpy.NPY_INT64 + self._python_dtype = int + self._sage_dtype = ZZ + self.__create_vector__() + return From 5323b255cfc7f6e461b02a2fb37aec6dd32e7aad Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 16 Nov 2021 12:27:22 -0800 Subject: [PATCH 7/9] src/sage/matrix/matrix_numpy_dense.pyx: Returns -> Return --- src/sage/matrix/matrix_numpy_dense.pyx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix_numpy_dense.pyx b/src/sage/matrix/matrix_numpy_dense.pyx index 5047d423c94..76f1b82ed01 100644 --- a/src/sage/matrix/matrix_numpy_dense.pyx +++ b/src/sage/matrix/matrix_numpy_dense.pyx @@ -341,7 +341,7 @@ cdef class Matrix_numpy_dense(Matrix_dense): OUTPUT: - Returns ``True`` if each entry above the diagonal (entries + Return ``True`` if each entry above the diagonal (entries with a row index less than the column index) is zero. EXAMPLES:: @@ -367,8 +367,9 @@ cdef class Matrix_numpy_dense(Matrix_dense): def numpy(self, dtype=None): """ - This method returns a copy of the matrix as a numpy array. It - uses the numpy C/api so is very fast. + Return a copy of the matrix as a numpy array. + + It uses the numpy C/api so is very fast. INPUT: From 93f3925a4b0beaef33d48608dfb3a61c26419d49 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Tue, 15 Feb 2022 20:45:44 -0800 Subject: [PATCH 8/9] src/sage/matrix/matrix_numpy_dense.pyx: Fix doc markup, whitespace --- src/sage/matrix/matrix_numpy_dense.pyx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix_numpy_dense.pyx b/src/sage/matrix/matrix_numpy_dense.pyx index 76f1b82ed01..215e8f46d76 100644 --- a/src/sage/matrix/matrix_numpy_dense.pyx +++ b/src/sage/matrix/matrix_numpy_dense.pyx @@ -280,7 +280,7 @@ cdef class Matrix_numpy_dense(Matrix_dense): trans.subdivide(col_divs, row_divs) return trans - def is_symmetric(self, tol = 1e-12): + def is_symmetric(self, tol=1e-12): """ Return whether this matrix is symmetric, to the given tolerance. @@ -297,7 +297,8 @@ cdef class Matrix_numpy_dense(Matrix_dense): sage: m.is_symmetric() False - The tolerance inequality is strict: + The tolerance inequality is strict:: + sage: m.is_symmetric(tol=0.1) False sage: m.is_symmetric(tol=0.11) @@ -431,7 +432,7 @@ cdef class Matrix_numpy_dense(Matrix_dense): else: return Matrix_dense.numpy(self, dtype=dtype) - def _replace_self_with_numpy(self,numpy_matrix): + def _replace_self_with_numpy(self, numpy_matrix): """ EXAMPLES:: From b0a1a04f20e05e50869e73895eb12692f98362a2 Mon Sep 17 00:00:00 2001 From: Matthias Koeppe Date: Sun, 20 Feb 2022 09:51:31 -0800 Subject: [PATCH 9/9] src/sage/{matrix/matrix,modules/vector}_numpy_integer_dense.pyx: Fix doctest output in 32bit --- src/sage/matrix/matrix_numpy_integer_dense.pyx | 4 ++-- src/sage/modules/vector_numpy_integer_dense.pyx | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix_numpy_integer_dense.pyx b/src/sage/matrix/matrix_numpy_integer_dense.pyx index e832386f618..bc605df9b92 100644 --- a/src/sage/matrix/matrix_numpy_integer_dense.pyx +++ b/src/sage/matrix/matrix_numpy_integer_dense.pyx @@ -13,8 +13,8 @@ EXAMPLES:: [ 0 0 0] [ 0 0 47] sage: M.numpy() - array([[ 0, 0, 0], - [ 0, 0, 47]]) + array([[ 0, 0, 0], [ 0, 0, 47]]) # 64-bit + array([[ 0, 0, 0], [ 0, 0, 47]], dtype=int64) # 32-bit """ # **************************************************************************** diff --git a/src/sage/modules/vector_numpy_integer_dense.pyx b/src/sage/modules/vector_numpy_integer_dense.pyx index b73a21eb849..faeb8531771 100644 --- a/src/sage/modules/vector_numpy_integer_dense.pyx +++ b/src/sage/modules/vector_numpy_integer_dense.pyx @@ -10,7 +10,8 @@ EXAMPLES:: sage: v (0, 42, 0) sage: v.numpy() - array([ 0, 42, 0]) + array([ 0, 42, 0]) # 64-bit + array([ 0, 42, 0], dtype=int64) # 32-bit """ # ****************************************************************************