Skip to content
Browse files

Merge pull request #363 from pv/sparse-eye

MAINT: sparse: refactor identity() into eye()

Makes also identity() return a dia_matrix by default. Closes Trac #1767
  • Loading branch information...
2 parents 2a649c6 + 0b4b881 commit eb17505e73a9e6721ddd3165e2c038f565a89930 @pv pv committed
View
4 scipy/sparse/base.py
@@ -361,8 +361,8 @@ def __pow__(self, other):
raise ValueError('exponent must be >= 0')
if other == 0:
- from construct import identity
- return identity( self.shape[0], dtype=self.dtype )
+ from construct import eye
+ return eye( self.shape[0], dtype=self.dtype )
elif other == 1:
return self.copy()
else:
View
8 scipy/sparse/benchmarks/bench_sparse.py
@@ -47,7 +47,7 @@ class BenchmarkSparse(TestCase):
def bench_arithmetic(self):
matrices = []
- #matrices.append( ('A','Identity', sparse.identity(500**2,format='csr')) )
+ #matrices.append( ('A','Identity', sparse.eye(500**2,format='csr')) )
matrices.append( ('A','Poisson5pt', poisson2d(250,format='csr')) )
matrices.append( ('B','Poisson5pt^2', poisson2d(250,format='csr')**2) )
@@ -128,8 +128,8 @@ def bench_sort(self):
def bench_matvec(self):
matrices = []
- matrices.append(('Identity', sparse.identity(10**4,format='dia')))
- matrices.append(('Identity', sparse.identity(10**4,format='csr')))
+ matrices.append(('Identity', sparse.eye(10**4,format='dia')))
+ matrices.append(('Identity', sparse.eye(10**4,format='csr')))
matrices.append(('Poisson5pt', poisson2d(300,format='lil')))
matrices.append(('Poisson5pt', poisson2d(300,format='dok')))
matrices.append(('Poisson5pt', poisson2d(300,format='dia')))
@@ -217,7 +217,7 @@ def bench_construction(self):
"""build matrices by inserting single values"""
matrices = []
matrices.append( ('Empty',csr_matrix((10000,10000))) )
- matrices.append( ('Identity',sparse.identity(10000)) )
+ matrices.append( ('Identity',sparse.eye(10000)) )
matrices.append( ('Poisson5pt', poisson2d(100)) )
print
View
71 scipy/sparse/construct.py
@@ -209,35 +209,60 @@ def identity(n, dtype='d', format=None):
with 3 stored elements (1 diagonals) in DIAgonal format>
"""
+ return eye(n, n, dtype=dtype, format=format)
- if format in ['csr','csc']:
- indptr = np.arange(n+1, dtype=np.intc)
- indices = np.arange(n, dtype=np.intc)
- data = np.ones(n, dtype=dtype)
- cls = eval('%s_matrix' % format)
- return cls((data,indices,indptr),(n,n))
- elif format == 'coo':
- row = np.arange(n, dtype=np.intc)
- col = np.arange(n, dtype=np.intc)
- data = np.ones(n, dtype=dtype)
- return coo_matrix((data,(row,col)),(n,n))
- elif format == 'dia':
- data = np.ones(n, dtype=dtype)
- diags = [0]
- return dia_matrix((data,diags), shape=(n,n))
- else:
- return identity(n, dtype=dtype, format='csr').asformat(format)
-
+def eye(m, n=None, k=0, dtype=float, format=None):
+ """Sparse matrix with ones on diagonal
-def eye(m, n, k=0, dtype='d', format=None):
- """eye(m, n) returns a sparse (m x n) matrix where the k-th diagonal
+ Returns a sparse (m x n) matrix where the k-th diagonal
is all ones and everything else is zeros.
+
+ Parameters
+ ----------
+ n : integer
+ Number of rows in the matrix.
+ m : integer, optional
+ Number of columns. Default: n
+ k : integer, optional
+ Diagonal to place ones on. Default: 0 (main diagonal)
+ dtype :
+ Data type of the matrix
+ format : string
+ Sparse format of the result, e.g. format="csr", etc.
+
+ Examples
+ --------
+ >>> from scipy import sparse
+ >>> sparse.eye(3).todense()
+ matrix([[ 1., 0., 0.],
+ [ 0., 1., 0.],
+ [ 0., 0., 1.]])
+ >>> sparse.eye(3, dtype=np.int8)
+ <3x3 sparse matrix of type '<type 'numpy.int8'>'
+ with 3 stored elements (1 diagonals) in DIAgonal format>
+
"""
+ if n is None:
+ n = m
m,n = int(m),int(n)
+
+ if m == n and k == 0:
+ # fast branch for special formats
+ if format in ['csr', 'csc']:
+ indptr = np.arange(n+1, dtype=np.intc)
+ indices = np.arange(n, dtype=np.intc)
+ data = np.ones(n, dtype=dtype)
+ cls = {'csr': csr_matrix, 'csc': csc_matrix}[format]
+ return cls((data,indices,indptr),(n,n))
+ elif format == 'coo':
+ row = np.arange(n, dtype=np.intc)
+ col = np.arange(n, dtype=np.intc)
+ data = np.ones(n, dtype=dtype)
+ return coo_matrix((data,(row,col)),(n,n))
+
diags = np.ones((1, max(0, min(m + k, n))), dtype=dtype)
return spdiags(diags, k, m, n).asformat(format)
-
def kron(A, B, format=None):
"""kronecker product of sparse matrices A and B
@@ -355,8 +380,8 @@ def kronsum(A, B, format=None):
dtype = upcast(A.dtype, B.dtype)
- L = kron(identity(B.shape[0],dtype=dtype), A, format=format)
- R = kron(B, identity(A.shape[0],dtype=dtype), format=format)
+ L = kron(eye(B.shape[0],dtype=dtype), A, format=format)
+ R = kron(B, eye(A.shape[0],dtype=dtype), format=format)
return (L+R).asformat(format) #since L + R is not always same format
View
8 scipy/sparse/linalg/eigen/arpack/arpack.py
@@ -47,7 +47,7 @@
import _arpack
import numpy as np
from scipy.sparse.linalg.interface import aslinearoperator, LinearOperator
-from scipy.sparse import identity, csc_matrix, csr_matrix, \
+from scipy.sparse import eye, csc_matrix, csr_matrix, \
isspmatrix, isspmatrix_csr
from scipy.linalg import lu_factor, lu_solve
from scipy.sparse.sputils import isdense
@@ -1035,7 +1035,7 @@ def get_OPinv_matvec(A, M, sigma, symmetric=False, tol=0):
A.flat[::A.shape[1] + 1] -= sigma
return LuInv(A).matvec
elif isspmatrix(A):
- A = A - sigma * identity(A.shape[0])
+ A = A - sigma * eye(A.shape[0])
if symmetric and isspmatrix_csr(A):
A = A.T
return SpLuInv(A.tocsc()).matvec
@@ -1181,7 +1181,7 @@ def eigs(A, k=6, M=None, sigma=None, which='LM', v0=None,
--------
Find 6 eigenvectors of the identity matrix:
- >>> id = np.identity(13)
+ >>> id = np.eye(13)
>>> vals, vecs = sp.sparse.linalg.eigs(id, k=6)
>>> vals
array([ 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])
@@ -1425,7 +1425,7 @@ def eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None,
Examples
--------
- >>> id = np.identity(13)
+ >>> id = np.eye(13)
>>> vals, vecs = sp.sparse.linalg.eigsh(id, k=6)
>>> vals
array([ 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])
View
8 scipy/sparse/linalg/isolve/tests/test_lsmr.py
@@ -15,7 +15,7 @@
Dept of MS&E, Stanford University.
"""
-from numpy import arange, concatenate, identity, zeros, ones, sqrt, \
+from numpy import arange, concatenate, eye, zeros, ones, sqrt, \
transpose, hstack
from numpy.linalg import norm
from numpy.testing import run_module_suite, assert_almost_equal
@@ -36,17 +36,17 @@ def assertCompatibleSystem(self, A, xtrue):
assert_almost_equal(norm(x - xtrue), 0, 6)
def testIdentityACase1(self):
- A = identity(self.n)
+ A = eye(self.n)
xtrue = zeros((self.n, 1))
self.assertCompatibleSystem(A, xtrue)
def testIdentityACase2(self):
- A = identity(self.n)
+ A = eye(self.n)
xtrue = ones((self.n,1))
self.assertCompatibleSystem(A, xtrue)
def testIdentityACase3(self):
- A = identity(self.n)
+ A = eye(self.n)
xtrue = transpose(arange(self.n,0,-1))
self.assertCompatibleSystem(A, xtrue)
View
10 scipy/sparse/linalg/tests/test_matfuncs.py
@@ -6,7 +6,7 @@
"""
import numpy as np
-from numpy import array, identity, dot, sqrt, double, exp, random
+from numpy import array, eye, dot, sqrt, double, exp, random
from numpy.testing import TestCase, run_module_suite, assert_array_almost_equal, \
assert_array_almost_equal_nulp
@@ -28,8 +28,8 @@ def test_padecases_dtype(self):
for dtype in [np.float32, np.float64, np.complex64, np.complex128]:
# test double-precision cases
for scale in [1e-2, 1e-1, 5e-1, 1, 10]:
- a = scale * identity(3, dtype=dtype)
- e = exp(scale) * identity(3, dtype=dtype)
+ a = scale * eye(3, dtype=dtype)
+ e = exp(scale) * eye(3, dtype=dtype)
assert_array_almost_equal_nulp(expm(a), e, nulp=100)
def test_padecases_dtype_sparse(self):
@@ -38,7 +38,7 @@ def test_padecases_dtype_sparse(self):
# test double-precision cases
for scale in [1e-2, 1e-1, 5e-1, 1, 10]:
a = scale * speye(3, 3, dtype=dtype, format='csc')
- e = exp(scale) * identity(3, dtype=dtype)
+ e = exp(scale) * eye(3, dtype=dtype)
assert_array_almost_equal_nulp(expm(a).toarray(), e, nulp=100)
def test_logm_consistency(self):
@@ -47,7 +47,7 @@ def test_logm_consistency(self):
for n in range(1, 10):
for scale in [1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1, 1e2]:
# make logm(a) be of a given scale
- a = (identity(n) + random.rand(n, n) * scale).astype(dtype)
+ a = (eye(n) + random.rand(n, n) * scale).astype(dtype)
if np.iscomplexobj(a):
a = a + 1j * random.rand(n, n) * scale
assert_array_almost_equal(expm(logm(a)), a)
View
15 scipy/sparse/tests/test_construct.py
@@ -207,6 +207,21 @@ def test_eye(self):
for n in [3, 5]:
for k in range(-5,6):
assert_equal(construct.eye(m, n, k=k).toarray(), np.eye(m, n, k=k))
+ if m == n:
+ assert_equal(construct.eye(m, k=k).toarray(), np.eye(m, n, k=k))
+
+ def test_eye_one(self):
+ assert_equal(construct.eye(1).toarray(), [[1]])
+ assert_equal(construct.eye(2).toarray(), [[1,0],[0,1]])
+
+ I = construct.eye(3, dtype='int8', format='dia')
+ assert_equal( I.dtype, np.dtype('int8') )
+ assert_equal( I.format, 'dia' )
+
+ for fmt in sparse_formats:
+ I = construct.eye( 3, format=fmt )
+ assert_equal( I.format, fmt )
+ assert_equal( I.toarray(), [[1,0,0],[0,1,0],[0,0,1]])
def test_kron(self):
cases = []

0 comments on commit eb17505

Please sign in to comment.
Something went wrong with that request. Please try again.