[Scipy-svn] r3641 - trunk/scipy/sparse
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Dec 13 15:25:45 EST 2007
Author: wnbell
Date: 2007-12-13 14:25:41 -0600 (Thu, 13 Dec 2007)
New Revision: 3641
Added:
trunk/scipy/sparse/construct.py
trunk/scipy/sparse/spfuncs.py
Modified:
trunk/scipy/sparse/__init__.py
trunk/scipy/sparse/info.py
trunk/scipy/sparse/sparse.py
Log:
intial split of sparse module into different files
Modified: trunk/scipy/sparse/__init__.py
===================================================================
--- trunk/scipy/sparse/__init__.py 2007-12-13 16:45:21 UTC (rev 3640)
+++ trunk/scipy/sparse/__init__.py 2007-12-13 20:25:41 UTC (rev 3641)
@@ -3,6 +3,8 @@
from info import __doc__
from sparse import *
+from construct import *
+from spfuncs import *
__all__ = filter(lambda s:not s.startswith('_'),dir())
from numpy.testing import NumpyTest
Added: trunk/scipy/sparse/construct.py
===================================================================
--- trunk/scipy/sparse/construct.py 2007-12-13 16:45:21 UTC (rev 3640)
+++ trunk/scipy/sparse/construct.py 2007-12-13 20:25:41 UTC (rev 3641)
@@ -0,0 +1,186 @@
+""" Functions to construct sparse matrices
+"""
+
+
+__all__ = [ 'spdiags','speye','spidentity','spkron', 'lil_eye', 'lil_diags' ]
+
+import itertools
+
+import numpy
+from numpy import ones, clip, array, arange, intc
+
+from sparse import csr_matrix, csc_matrix, coo_matrix, \
+ dok_matrix, lil_matrix
+from sparse import isspmatrix, isspmatrix_csr, isspmatrix_csc
+import sparsetools
+
+def _spdiags_tosub(diag_num, a, b):
+ part1 = where(less(diag_num, a), abs(diag_num-a), 0)
+ part2 = where(greater(diag_num, b), abs(diag_num-b), 0)
+ return part1+part2
+
+# Note: sparsetools only offers diagonal -> CSC matrix conversion functions,
+# not to CSR
+def spdiags(diags, offsets, M, N):
+ """Return a sparse matrix in CSC format given its diagonals.
+
+ B = spdiags(diags, offsets, M, N)
+
+ Inputs:
+ diags -- rows contain diagonal values
+ offsets -- diagonals to set (0 is main)
+ M, N -- sparse matrix returned is M X N
+ """
+ # diags = array(transpose(diags), copy=True)
+ diags = array(diags, copy = True)
+ if diags.dtype.char not in 'fdFD':
+ diags = diags.astype('d')
+ if not hasattr(offsets, '__len__' ):
+ offsets = (offsets,)
+ offsets = array(offsets, copy=False, dtype=intc)
+ assert(len(offsets) == diags.shape[0])
+ indptr, rowind, data = sparsetools.spdiags(M, N, len(offsets), offsets, diags)
+ return csc_matrix((data, rowind, indptr), (M, N))
+
+
+
+def spidentity(n, dtype='d'):
+ """
+ spidentity( n ) returns the identity matrix of shape (n, n) stored
+ in CSC sparse matrix format.
+ """
+ return csc_matrix((ones(n,dtype=dtype),arange(n),arange(n+1)),(n,n))
+
+
+def speye(n, m, k = 0, dtype = 'd'):
+ """
+ speye(n, m) returns a (n x m) matrix stored
+ in CSC sparse matrix format, where the k-th diagonal is all ones,
+ and everything else is zeros.
+ """
+ diags = ones((1, n), dtype = dtype)
+ return spdiags(diags, k, n, m)
+
+def spkron(a,b):
+ """kronecker product of sparse matrices a and b
+
+ *Parameters*:
+ a,b : sparse matrices
+ E.g. csr_matrix, csc_matrix, coo_matrix, etc.
+
+ *Returns*:
+ coo_matrix
+ kronecker product in COOrdinate format
+
+ *Example*:
+ -------
+
+ >>> a = csr_matrix(array([[0,2],[5,0]]))
+ >>> b = csr_matrix(array([[1,2],[3,4]]))
+ >>> spkron(a,b).todense()
+ matrix([[ 0., 0., 2., 4.],
+ [ 0., 0., 6., 8.],
+ [ 5., 10., 0., 0.],
+ [ 15., 20., 0., 0.]])
+
+ """
+ if not isspmatrix(a) and isspmatrix(b):
+ raise ValueError,'expected sparse matrix'
+
+ a,b = a.tocoo(),b.tocoo()
+ output_shape = (a.shape[0]*b.shape[0],a.shape[1]*b.shape[1])
+
+ if a.nnz == 0 or b.nnz == 0:
+ # kronecker product is the zero matrix
+ return coo_matrix( output_shape )
+
+
+ # expand entries of a into blocks
+ row = a.row.repeat(b.nnz)
+ col = a.col.repeat(b.nnz)
+ data = a.data.repeat(b.nnz)
+
+ row *= b.shape[0]
+ col *= b.shape[1]
+
+ # increment block indices
+ row,col = row.reshape(-1,b.nnz),col.reshape(-1,b.nnz)
+ row += b.row
+ col += b.col
+ row,col = row.reshape(-1),col.reshape(-1)
+
+ # compute block entries
+ data = data.reshape(-1,b.nnz) * b.data
+ data = data.reshape(-1)
+
+ return coo_matrix((data,(row,col)), dims=output_shape)
+
+
+
+
+def lil_eye((r,c), k=0, dtype=float):
+ """Generate a lil_matrix of dimensions (r,c) with the k-th
+ diagonal set to 1.
+
+ :Parameters:
+ r,c : int
+ Row and column-dimensions of the output.
+ k : int
+ Diagonal offset. In the output matrix,
+ out[m,m+k] == 1 for all m.
+ dtype : dtype
+ Data-type of the output array.
+
+ """
+ out = lil_matrix((r,c),dtype=dtype)
+ for c in xrange(clip(k,0,c),clip(r+k,0,c)):
+ out.rows[c-k].append(c)
+ out.data[c-k].append(1)
+ return out
+
+def lil_diags(diags,offsets,(m,n),dtype=float):
+ """Generate a lil_matrix with the given diagonals.
+
+ :Parameters:
+ diags : list of list of values e.g. [[1,2,3],[4,5]]
+ Values to be placed on each indicated diagonal.
+ offsets : list of ints
+ Diagonal offsets. This indicates the diagonal on which
+ the given values should be placed.
+ (r,c) : tuple of ints
+ Row and column dimensions of the output.
+ dtype : dtype
+ Output data-type.
+
+ Example:
+ -------
+
+ >>> lil_diags([[1,2,3],[4,5],[6]],[0,1,2],(3,3)).todense()
+ matrix([[ 1., 4., 6.],
+ [ 0., 2., 5.],
+ [ 0., 0., 3.]])
+
+ """
+ offsets_unsorted = list(offsets)
+ diags_unsorted = list(diags)
+ if len(diags) != len(offsets):
+ raise ValueError("Number of diagonals provided should "
+ "agree with offsets.")
+
+ sort_indices = numpy.argsort(offsets_unsorted)
+ diags = [diags_unsorted[k] for k in sort_indices]
+ offsets = [offsets_unsorted[k] for k in sort_indices]
+
+ for i,k in enumerate(offsets):
+ if len(diags[i]) < m-abs(k):
+ raise ValueError("Not enough values specified to fill "
+ "diagonal %s." % k)
+
+ out = lil_matrix((m,n),dtype=dtype)
+ for k,diag in itertools.izip(offsets,diags):
+ for ix,c in enumerate(xrange(clip(k,0,n),clip(m+k,0,n))):
+ out.rows[c-k].append(c)
+ out.data[c-k].append(diag[ix])
+ return out
+
+
Modified: trunk/scipy/sparse/info.py
===================================================================
--- trunk/scipy/sparse/info.py 2007-12-13 16:45:21 UTC (rev 3640)
+++ trunk/scipy/sparse/info.py 2007-12-13 20:25:41 UTC (rev 3641)
@@ -12,7 +12,7 @@
(2) csr_matrix: Compressed Sparse Row format
(3) lil_matrix: List of Lists format
(4) dok_matrix: Dictionary of Keys format
- (5) coo_matrix: COOrdinate format (IJV triplets)
+ (5) coo_matrix: COOrdinate format (aka IJV, triplet format)
To construct a matrix efficiently, use either lil_matrix (recommended) or
dok_matrix. The lil_matrix class supports basic slicing and fancy
@@ -74,10 +74,11 @@
Further Details:
CSR column indices are not necessarily sorted. Likewise for CSC row
- indices. Use the .ensure_sorted_indices() method when sorted indices
- are necessary. Note that there is no expectation for sorted indices
- in the sparsetools module. Furthermore some sparsetools functions
- produce matrices with unsorted indices even when sorted input is given.
+ indices. Use the .sorted_indices() and .sort_indices() methods when
+ sorted indices are necessary. Note that there is no expectation for
+ sorted indices in the sparsetools module. Furthermore some sparsetools
+ functions produce matrices with unsorted indices even when sorted
+ input is given.
"""
postpone_import = 1
Modified: trunk/scipy/sparse/sparse.py
===================================================================
--- trunk/scipy/sparse/sparse.py 2007-12-13 16:45:21 UTC (rev 3640)
+++ trunk/scipy/sparse/sparse.py 2007-12-13 20:25:41 UTC (rev 3641)
@@ -6,11 +6,9 @@
__all__ = ['spmatrix','csc_matrix','csr_matrix','coo_matrix',
- 'lil_matrix','dok_matrix',
- 'spdiags','speye','spidentity','spkron','extract_diagonal',
- 'isspmatrix','issparse','isspmatrix_csc','isspmatrix_csr',
- 'isspmatrix_lil','isspmatrix_dok', 'isspmatrix_coo',
- 'lil_eye', 'lil_diags' ]
+ 'lil_matrix','dok_matrix' ]
+__all__ += [ 'isspmatrix','issparse','isspmatrix_csc','isspmatrix_csr',
+ 'isspmatrix_lil','isspmatrix_dok', 'isspmatrix_coo' ]
import warnings
@@ -24,6 +22,7 @@
from scipy.sparse.sparsetools import csrtodense, \
cootocsr, cootocsc, csctocsr, csrtocsc
+
import sparsetools
import itertools, operator, copy
from bisect import bisect_left
@@ -2593,6 +2592,9 @@
# unsymmetric sparse skyline
# variable block row
+
+
+
def _isinstance(x, _class):
##
# This makes scipy.sparse.sparse.csc_matrix == __main__.csc_matrix.
@@ -2653,6 +2655,10 @@
else:
return True
+def issequence(t):
+ return isinstance(t, (list, tuple))\
+ or (isinstance(t, ndarray) and (t.ndim == 1))
+
def getdtype(dtype, a=None, default=None):
"""Function used to simplify argument processing. If 'dtype' is not
specified (is None), returns a.dtype; otherwise returns a numpy.dtype
@@ -2676,191 +2682,3 @@
return newdtype
-def _spdiags_tosub(diag_num, a, b):
- part1 = where(less(diag_num, a), abs(diag_num-a), 0)
- part2 = where(greater(diag_num, b), abs(diag_num-b), 0)
- return part1+part2
-
-# Note: sparsetools only offers diagonal -> CSC matrix conversion functions,
-# not to CSR
-def spdiags(diags, offsets, M, N):
- """Return a sparse matrix in CSC format given its diagonals.
-
- B = spdiags(diags, offsets, M, N)
-
- Inputs:
- diags -- rows contain diagonal values
- offsets -- diagonals to set (0 is main)
- M, N -- sparse matrix returned is M X N
- """
- # diags = array(transpose(diags), copy=True)
- diags = numpy.array(diags, copy = True)
- if diags.dtype.char not in 'fdFD':
- diags = diags.astype('d')
- if not hasattr(offsets, '__len__' ):
- offsets = (offsets,)
- offsets = array(offsets, copy=False, dtype=numpy.intc)
- assert(len(offsets) == diags.shape[0])
- indptr, rowind, data = sparsetools.spdiags(M, N, len(offsets), offsets, diags)
- return csc_matrix((data, rowind, indptr), (M, N))
-
-def extract_diagonal(A):
- """
- extract_diagonal(A) returns the main diagonal of A.
- """
- #TODO extract k-th diagonal
- if isspmatrix_csr(A) or isspmatrix_csc(A):
- fn = getattr(sparsetools, "extract_" + A.format + "_diagonal")
- y = empty( min(A.shape), dtype=upcast(A.dtype) )
- fn(A.shape[0],A.shape[1],A.indptr,A.indices,A.data,y)
- return y
- elif isspmatrix(A):
- return extract_diagonal(A.tocsr())
- else:
- raise ValueError,'expected sparse matrix'
-
-
-
-
-def spidentity(n, dtype='d'):
- """
- spidentity( n ) returns the identity matrix of shape (n, n) stored
- in CSC sparse matrix format.
- """
- return csc_matrix((ones(n,dtype=dtype),arange(n),arange(n+1)),(n,n))
-
-
-def speye(n, m, k = 0, dtype = 'd'):
- """
- speye(n, m) returns a (n x m) matrix stored
- in CSC sparse matrix format, where the k-th diagonal is all ones,
- and everything else is zeros.
- """
- diags = ones((1, n), dtype = dtype)
- return spdiags(diags, k, n, m)
-
-def spkron(a,b):
- """kronecker product of sparse matrices a and b
-
- *Parameters*:
- a,b : sparse matrices
- E.g. csr_matrix, csc_matrix, coo_matrix, etc.
-
- *Returns*:
- coo_matrix
- kronecker product in COOrdinate format
-
- *Example*:
- -------
-
- >>> a = csr_matrix(array([[0,2],[5,0]]))
- >>> b = csr_matrix(array([[1,2],[3,4]]))
- >>> spkron(a,b).todense()
- matrix([[ 0., 0., 2., 4.],
- [ 0., 0., 6., 8.],
- [ 5., 10., 0., 0.],
- [ 15., 20., 0., 0.]])
-
- """
- if not isspmatrix(a) and isspmatrix(b):
- raise ValueError,'expected sparse matrix'
-
- a,b = a.tocoo(),b.tocoo()
- output_shape = (a.shape[0]*b.shape[0],a.shape[1]*b.shape[1])
-
- if a.nnz == 0 or b.nnz == 0:
- # kronecker product is the zero matrix
- return coo_matrix( output_shape )
-
-
- # expand entries of a into blocks
- row = a.row.repeat(b.nnz)
- col = a.col.repeat(b.nnz)
- data = a.data.repeat(b.nnz)
-
- row *= b.shape[0]
- col *= b.shape[1]
-
- # increment block indices
- row,col = row.reshape(-1,b.nnz),col.reshape(-1,b.nnz)
- row += b.row
- col += b.col
- row,col = row.reshape(-1),col.reshape(-1)
-
- # compute block entries
- data = data.reshape(-1,b.nnz) * b.data
- data = data.reshape(-1)
-
- return coo_matrix((data,(row,col)), dims=output_shape)
-
-
-
-
-def lil_eye((r,c), k=0, dtype=float):
- """Generate a lil_matrix of dimensions (r,c) with the k-th
- diagonal set to 1.
-
- :Parameters:
- r,c : int
- Row and column-dimensions of the output.
- k : int
- Diagonal offset. In the output matrix,
- out[m,m+k] == 1 for all m.
- dtype : dtype
- Data-type of the output array.
-
- """
- out = lil_matrix((r,c),dtype=dtype)
- for c in xrange(clip(k,0,c),clip(r+k,0,c)):
- out.rows[c-k].append(c)
- out.data[c-k].append(1)
- return out
-
-def lil_diags(diags,offsets,(m,n),dtype=float):
- """Generate a lil_matrix with the given diagonals.
-
- :Parameters:
- diags : list of list of values e.g. [[1,2,3],[4,5]]
- Values to be placed on each indicated diagonal.
- offsets : list of ints
- Diagonal offsets. This indicates the diagonal on which
- the given values should be placed.
- (r,c) : tuple of ints
- Row and column dimensions of the output.
- dtype : dtype
- Output data-type.
-
- Example:
- -------
-
- >>> lil_diags([[1,2,3],[4,5],[6]],[0,1,2],(3,3)).todense()
- matrix([[ 1., 4., 6.],
- [ 0., 2., 5.],
- [ 0., 0., 3.]])
-
- """
- offsets_unsorted = list(offsets)
- diags_unsorted = list(diags)
- if len(diags) != len(offsets):
- raise ValueError("Number of diagonals provided should "
- "agree with offsets.")
-
- sort_indices = numpy.argsort(offsets_unsorted)
- diags = [diags_unsorted[k] for k in sort_indices]
- offsets = [offsets_unsorted[k] for k in sort_indices]
-
- for i,k in enumerate(offsets):
- if len(diags[i]) < m-abs(k):
- raise ValueError("Not enough values specified to fill "
- "diagonal %s." % k)
-
- out = lil_matrix((m,n),dtype=dtype)
- for k,diag in itertools.izip(offsets,diags):
- for ix,c in enumerate(xrange(clip(k,0,n),clip(m+k,0,n))):
- out.rows[c-k].append(c)
- out.data[c-k].append(diag[ix])
- return out
-
-def issequence(t):
- return isinstance(t, (list, tuple))\
- or (isinstance(t, ndarray) and (t.ndim == 1))
Added: trunk/scipy/sparse/spfuncs.py
===================================================================
--- trunk/scipy/sparse/spfuncs.py 2007-12-13 16:45:21 UTC (rev 3640)
+++ trunk/scipy/sparse/spfuncs.py 2007-12-13 20:25:41 UTC (rev 3641)
@@ -0,0 +1,27 @@
+""" Functions that operate on sparse matrices
+"""
+
+from numpy import empty
+
+from sparse import isspmatrix_csr, isspmatrix_csc, isspmatrix
+from sparse import upcast
+
+import sparsetools
+
+def extract_diagonal(A):
+ """
+ extract_diagonal(A) returns the main diagonal of A.
+ """
+ #TODO extract k-th diagonal
+ if isspmatrix_csr(A) or isspmatrix_csc(A):
+ fn = getattr(sparsetools, "extract_" + A.format + "_diagonal")
+ y = empty( min(A.shape), dtype=upcast(A.dtype) )
+ fn(A.shape[0],A.shape[1],A.indptr,A.indices,A.data,y)
+ return y
+ elif isspmatrix(A):
+ return extract_diagonal(A.tocsr())
+ else:
+ raise ValueError,'expected sparse matrix'
+
+
+
More information about the Scipy-svn
mailing list