[Scipy-svn] r3910 - in trunk/scipy/sparse: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Feb 9 15:02:39 CST 2008


Author: wnbell
Date: 2008-02-09 15:02:27 -0600 (Sat, 09 Feb 2008)
New Revision: 3910

Modified:
   trunk/scipy/sparse/construct.py
   trunk/scipy/sparse/tests/test_base.py
   trunk/scipy/sparse/tests/test_construct.py
Log:
added sparse.kronsum and sparse.kron
deprecated sparse.spkron


Modified: trunk/scipy/sparse/construct.py
===================================================================
--- trunk/scipy/sparse/construct.py	2008-02-09 07:26:50 UTC (rev 3909)
+++ trunk/scipy/sparse/construct.py	2008-02-09 21:02:27 UTC (rev 3910)
@@ -2,7 +2,7 @@
 """
 
 
-__all__ = [ 'spdiags','speye','spidentity','spkron', 'bmat', 'lil_eye', 'lil_diags' ]
+__all__ = [ 'spdiags','speye','spidentity', 'spkron', 'kron', 'kronsum', 'bmat', 'lil_eye', 'lil_diags' ]
 
 from itertools import izip
 from warnings import warn
@@ -81,7 +81,7 @@
     diags = ones((1, m), dtype = dtype)
     return spdiags(diags, k, m, n).asformat(format)
 
-def spkron(A, B, format=None):
+def kron(A, B, format=None):
     """kronecker product of sparse matrices A and B
 
     Parameters
@@ -100,13 +100,13 @@
 
     >>> A = csr_matrix(array([[0,2],[5,0]]))
     >>> B = csr_matrix(array([[1,2],[3,4]]))
-    >>> spkron(A,B).todense()
+    >>> kron(A,B).todense()
     matrix([[ 0,  0,  2,  4],
             [ 0,  0,  6,  8],
             [ 5, 10,  0,  0],
             [15, 20,  0,  0]])
 
-    >>> spkron(A,[[1,2],[3,4]]).todense()
+    >>> kron(A,[[1,2],[3,4]]).todense()
     matrix([[ 0,  0,  2,  4],
             [ 0,  0,  6,  8],
             [ 5, 10,  0,  0],
@@ -159,78 +159,49 @@
 
         return coo_matrix((data,(row,col)), shape=output_shape).asformat(format)
 
+def kronsum(A, B, format=None):
+    """kronecker sum of sparse matrices A and B
 
+    Kronecker sum of two sparse matrices is a sum of two Kronecker 
+    products kron(I_n,A) + kron(B,I_m) where A has shape (m,m)
+    and B has shape (n,n) and I_m and I_n are identity matrices 
+    of shape (m,m) and (n,n) respectively.
 
-
-def lil_eye((r,c), k=0, dtype='d'):
-    """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.
+        A,B    : squared dense or sparse matrices
+        format : format of the result (e.g. "csr")
+            -  By default (format=None) an appropriate sparse matrix 
+               format is returned.  This choice is subject to change.
 
-    """
-    warn("lil_eye is deprecated. use speye(... , format='lil') instead", \
-            DeprecationWarning)
-    return speye(r,c,k,dtype=dtype,format='lil')
+    Returns
+    =======
+        kronecker sum in a sparse matrix format
 
+    Examples
+    ========
 
+    
+    """
+    A = coo_matrix(A)
+    B = coo_matrix(B)
 
-#TODO remove this function
-def lil_diags(diags,offsets,(m,n),dtype='d'):
-    """Generate a lil_matrix with the given diagonals.
+    if A.shape[0] != A.shape[1]:
+        raise ValueError('A is not square')
+    
+    if B.shape[0] != B.shape[1]:
+        raise ValueError('B is not square')
 
-    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.
+    dtype = upcast(A.dtype,B.dtype)
 
-    Example
-    =======
+    L = kron(spidentity(B.shape[0],dtype=dtype), A, format=format) 
+    R = kron(B, spidentity(A.shape[0],dtype=dtype), format=format)
 
-    >>> lil_diags([[1,2,3],[4,5],[6]],[0,1,2],(3,3)).todense()
-    matrix([[ 1.,  4.,  6.],
-            [ 0.,  2.,  5.],
-            [ 0.,  0.,  3.]])
+    return (L+R).asformat(format) #since L + R is not always same format
 
-    """
-    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 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 bmat( blocks, format=None, dtype=None ):
     """
     Build a sparse matrix from sparse sub-blocks
@@ -332,3 +303,79 @@
     return coo_matrix( (data, (row, col)), shape=shape ).asformat(format)
 
 
+
+#################################
+# Deprecated functions
+################################
+from numpy import deprecate
+
+spkron = deprecate(kron, oldname='spkron', newname='kron')
+
+def lil_eye((r,c), k=0, dtype='d'):
+    """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.
+
+    """
+    warn("lil_eye is deprecated. use speye(... , format='lil') instead", \
+            DeprecationWarning)
+    return speye(r,c,k,dtype=dtype,format='lil')
+
+
+
+#TODO remove this function
+def lil_diags(diags,offsets,(m,n),dtype='d'):
+    """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 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/tests/test_base.py
===================================================================
--- trunk/scipy/sparse/tests/test_base.py	2008-02-09 07:26:50 UTC (rev 3909)
+++ trunk/scipy/sparse/tests/test_base.py	2008-02-09 21:02:27 UTC (rev 3910)
@@ -17,14 +17,15 @@
 
 import numpy
 from numpy import arange, zeros, array, dot, ones, matrix, asmatrix, \
-        asarray, vstack, ndarray, kron, transpose, diag
+        asarray, vstack, ndarray, transpose, diag
 
 import random
 from scipy.testing import *
 
+import scipy.sparse as sparse
 from scipy.sparse import csc_matrix, csr_matrix, dok_matrix, \
         coo_matrix, lil_matrix, dia_matrix, bsr_matrix, \
-        speye, spkron, SparseEfficiencyWarning
+        speye, SparseEfficiencyWarning
 from scipy.sparse.sputils import supported_dtypes
 from scipy.splinalg import splu
 
@@ -84,12 +85,12 @@
         mats.append( [[0,1],[0,2],[0,3]] )
         mats.append( [[0,0,1],[0,0,2],[0,3,0]] )
 
-        mats.append( kron(mats[0],[[1,2]]) )
-        mats.append( kron(mats[0],[[1],[2]]) )
-        mats.append( kron(mats[1],[[1,2],[3,4]]) )
-        mats.append( kron(mats[2],[[1,2],[3,4]]) )
-        mats.append( kron(mats[3],[[1,2],[3,4]]) )
-        mats.append( kron(mats[3],[[1,2,3,4]]) )
+        mats.append( numpy.kron(mats[0],[[1,2]]) )
+        mats.append( numpy.kron(mats[0],[[1],[2]]) )
+        mats.append( numpy.kron(mats[1],[[1,2],[3,4]]) )
+        mats.append( numpy.kron(mats[2],[[1,2],[3,4]]) )
+        mats.append( numpy.kron(mats[3],[[1,2],[3,4]]) )
+        mats.append( numpy.kron(mats[3],[[1,2,3,4]]) )
 
         for m in mats:
             assert_equal(self.spmatrix(m).diagonal(),diag(m))
@@ -345,7 +346,7 @@
             assert_equal( result, dot(a,b) )
 
     def test_formatconversions(self):
-        A = spkron([[1,0,1],[0,1,1],[1,0,0]], [[1,1],[0,1]] )
+        A = sparse.kron([[1,0,1],[0,1,1],[1,0,0]], [[1,1],[0,1]] )
         D = A.todense()
         A = self.spmatrix(A)
 
@@ -371,7 +372,7 @@
     def test_tocompressedblock(self):
         x = array([[1,0,2,0],[0,0,0,0],[0,0,4,5]])
         y = array([[0,1,2],[3,0,5]])
-        A = kron(x,y)
+        A = numpy.kron(x,y)
         Asp = self.spmatrix(A)
         for format in ['bsr']:
             fn = getattr(Asp, 'to' + format )
@@ -1277,7 +1278,7 @@
         data[3] = array([[ 0,  5, 10],
                          [15,  0, 25]])
 
-        A = kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
+        A = numpy.kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
         Asp = bsr_matrix((data,indices,indptr),shape=(6,12))
         assert_equal(Asp.todense(),A)
         
@@ -1296,7 +1297,7 @@
         assert_equal(bsr_matrix(A,blocksize=(2,2)).todense(),A)
         assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
 
-        A = kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
+        A = numpy.kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
         assert_equal(bsr_matrix(A).todense(),A)
         assert_equal(bsr_matrix(A,shape=(6,12)).todense(),A)
         assert_equal(bsr_matrix(A,blocksize=(1,1)).todense(),A)
@@ -1306,11 +1307,11 @@
         assert_equal(bsr_matrix(A,blocksize=(3,12)).todense(),A)
         assert_equal(bsr_matrix(A,blocksize=(6,12)).todense(),A)
         
-        A = kron( [[1,0,2,0],[0,1,0,0],[0,0,0,0]], [[0,1,2],[3,0,5]] )
+        A = numpy.kron( [[1,0,2,0],[0,1,0,0],[0,0,0,0]], [[0,1,2],[3,0,5]] )
         assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
         
     def test_eliminate_zeros(self):
-        data = kron([1, 0, 0, 0, 2, 0, 3, 0], [[1,1],[1,1]]).T
+        data = numpy.kron([1, 0, 0, 0, 2, 0, 3, 0], [[1,1],[1,1]]).T
         data = data.reshape(-1,2,2)
         indices = array( [1, 2, 3, 4, 5, 6, 7, 8] )
         indptr  = array( [0, 3, 8] )

Modified: trunk/scipy/sparse/tests/test_construct.py
===================================================================
--- trunk/scipy/sparse/tests/test_construct.py	2008-02-09 07:26:50 UTC (rev 3909)
+++ trunk/scipy/sparse/tests/test_construct.py	2008-02-09 21:02:27 UTC (rev 3910)
@@ -1,6 +1,7 @@
 """test sparse matrix construction functions"""
 
-from numpy import array, matrix, kron
+import numpy
+from numpy import array, matrix
 from scipy.testing import *
 
 
@@ -77,7 +78,7 @@
         b = array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype='d')
         assert_array_equal(a.toarray(), b)
 
-    def test_spkron(self):
+    def test_kron(self):
         cases = []
 
         cases.append(array([[ 0]]))
@@ -96,11 +97,30 @@
         
         for a in cases:
             for b in cases:
-                result = spkron(csr_matrix(a),csr_matrix(b)).todense()
-                expected = kron(a,b)
+                result   = kron(csr_matrix(a),csr_matrix(b)).todense()
+                expected = numpy.kron(a,b)
+                assert_array_equal(result,expected)
 
+    def test_kronsum(self):
+        cases = []
+
+        cases.append(array([[ 0]]))
+        cases.append(array([[-1]]))
+        cases.append(array([[ 4]]))
+        cases.append(array([[10]]))
+        cases.append(array([[1,2],[3,4]]))
+        cases.append(array([[0,2],[5,0]]))
+        cases.append(array([[0,2,-6],[8,0,14],[0,3,0]]))
+        cases.append(array([[1,0,0],[0,5,-1],[4,-2,8]]))
+        
+        for a in cases:
+            for b in cases:
+                result   = kronsum(csr_matrix(a),csr_matrix(b)).todense()
+                expected = numpy.kron(numpy.eye(len(b)), a) + \
+                        numpy.kron(b, numpy.eye(len(a)))
                 assert_array_equal(result,expected)
 
+
     def test_bmat(self):
 
         A = coo_matrix([[1,2],[3,4]])



More information about the Scipy-svn mailing list