[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