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

scipy-svn@scip... scipy-svn@scip...
Sat Nov 22 16:47:23 CST 2008


Author: wnbell
Date: 2008-11-22 16:47:13 -0600 (Sat, 22 Nov 2008)
New Revision: 5171

Modified:
   trunk/scipy/sparse/base.py
   trunk/scipy/sparse/linalg/interface.py
   trunk/scipy/sparse/linalg/tests/test_interface.py
   trunk/scipy/sparse/tests/test_base.py
Log:
LinearOperator now wraps user-defined matvec and matmat routines
corrected sparse __mul__ handling of arguments with small shape


Modified: trunk/scipy/sparse/base.py
===================================================================
--- trunk/scipy/sparse/base.py	2008-11-22 20:18:12 UTC (rev 5170)
+++ trunk/scipy/sparse/base.py	2008-11-22 22:47:13 UTC (rev 5171)
@@ -291,8 +291,9 @@
             # If it's a list or whatever, treat it like a matrix
             other = np.asanyarray(other)
 
-        if isdense(other) and np.asarray(other).squeeze().ndim <= 1:
-            ##
+        other = np.asanyarray(other)
+    
+        if other.ndim == 1 or other.ndim == 2 and other.shape[1] == 1:
             # dense row or column vector
             if other.shape != (N,) and other.shape != (N,1):
                 raise ValueError('dimension mismatch')
@@ -308,7 +309,7 @@
 
             return result
 
-        elif len(other.shape) == 2:
+        elif other.ndim == 2:
             ##
             # dense 2D array or matrix ("multivector")
 

Modified: trunk/scipy/sparse/linalg/interface.py
===================================================================
--- trunk/scipy/sparse/linalg/interface.py	2008-11-22 20:18:12 UTC (rev 5170)
+++ trunk/scipy/sparse/linalg/interface.py	2008-11-22 22:47:13 UTC (rev 5171)
@@ -1,5 +1,4 @@
-import numpy
-from numpy import matrix, ndarray, asarray, dot, atleast_2d, hstack
+import numpy as np
 from scipy.sparse.sputils import isshape
 from scipy.sparse import isspmatrix
 
@@ -20,13 +19,12 @@
     shape : tuple
         Matrix dimensions (M,N)
     matvec : callable f(v)
-        Returns returns A * v
+        Returns returns A * v.
 
     Optional Parameters
     -------------------
     rmatvec : callable f(v)
-        Returns A^H * v, where A^H represents the Hermitian
-        (conjugate) transpose of A.
+        Returns A^H * v, where A^H is the conjugate transpose of A.
     matmat : callable f(V)
         Returns A * V, where V is a dense matrix with dimensions (N,K).
     dtype : dtype
@@ -36,6 +34,12 @@
     --------
     aslinearoperator : Construct LinearOperators
 
+    Notes
+    -----
+    The user-defined matvec() function must properly handle the case
+    where v has shape (N,) as well as the (N,1) case.  The shape of
+    the return type is handled internally by LinearOperator.
+
     Examples
     --------
     >>> from scipy.sparse.linalg import LinearOperator
@@ -52,7 +56,7 @@
     array([ 2.,  3.])
 
     """
-    def __init__( self, shape, matvec, rmatvec=None, matmat=None, dtype=None ):
+    def __init__(self, shape, matvec, rmatvec=None, matmat=None, dtype=None):
 
         shape = tuple(shape)
 
@@ -60,7 +64,7 @@
             raise ValueError('invalid shape')
 
         self.shape  = shape
-        self.matvec = matvec
+        self._matvec = matvec
 
         if rmatvec is None:
             def rmatvec(v):
@@ -69,26 +73,118 @@
         else:
             self.rmatvec = rmatvec
 
-        if matmat is None:
+        if matmat is not None:
             # matvec each column of V
-            def matmat(V):
-                V = asarray(V)
-                return hstack( [ matvec(col.reshape(-1,1)) for col in V.T ] )
-            self.matmat = matmat
-        else:
-            self.matmat = matmat
+            self._matmat = matmat
 
         if dtype is not None:
-            self.dtype = numpy.dtype(dtype)
+            self.dtype = np.dtype(dtype)
 
+
+    def _matmat(self, X):
+        """Default matrix-matrix multiplication handler.  Falls back on
+        the user-defined matvec() routine, which is always provided.
+        """
+        
+        return np.hstack( [ self.matvec(col.reshape(-1,1)) for col in X.T ] )
+
+
+    def matvec(self, x):
+        """Matrix-vector multiplication
+
+        Performs the operation y=A*x where A is an MxN linear 
+        operator and x is a column vector or rank-1 array.
+
+        Parameters
+        ----------
+        x : {matrix, ndarray}
+            An array with shape (N,) or (N,1).
+
+        Returns
+        -------
+        y : {matrix, ndarray}
+            A matrix or ndarray with shape (M,) or (M,1) depending 
+            on the type and shape of the x argument.
+
+        Notes
+        -----
+        This matvec wraps the user-specified matvec routine to ensure that
+        y has the correct shape and type.
+
+        """
+
+        x = np.asanyarray(x)
+        
+        M,N = self.shape
+            
+        if x.shape != (N,) and x.shape != (N,1):
+            raise ValueError('dimension mismatch')
+
+        y = self._matvec(x)
+
+        if x.ndim == 2:
+            # If 'x' is a column vector, reshape the result
+            y = y.reshape(-1,1)
+
+        if isinstance(x, np.matrix):
+            y = np.asmatrix(y)
+
+        return y
+
+
+    def matmat(self, X):
+        """Matrix-matrix multiplication
+
+        Performs the operation y=A*X where A is an MxN linear 
+        operator and X dense N*K matrix or ndarray.
+
+        Parameters
+        ----------
+        X : {matrix, ndarray}
+            An array with shape (N,K).
+
+        Returns
+        -------
+        Y : {matrix, ndarray}
+            A matrix or ndarray with shape (M,K) depending on
+            the type of the X argument.
+
+        Notes
+        -----
+        This matmat wraps any user-specified matmat routine to ensure that
+        y has the correct type.
+
+        """
+
+        X = np.asanyarray(X)
+        
+        if X.ndim != 2:
+            raise ValueError('expected rank-2 ndarray or matrix')
+        
+        M,N = self.shape
+
+        if X.shape[0] != N:
+            raise ValueError('dimension mismatch')
+
+        Y = self._matmat(X)
+
+        if isinstance(Y, np.matrix):
+            Y = np.asmatrix(Y)
+
+        return Y
+        
+
     def __mul__(self,x):
-        x = numpy.asarray(x)
+        x = np.asarray(x)
 
-        if numpy.rank(x.squeeze()) == 1:
+        if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1:
             return self.matvec(x)
+        elif x.ndim == 2:
+            return self.matmat(x)
         else:
-            return self.matmat(x)
+            raise ValueError('expected rank-1 or rank-2 array or matrix')
 
+
     def __repr__(self):
         M,N = self.shape
         if hasattr(self,'dtype'):
@@ -121,18 +217,18 @@
     if isinstance(A, LinearOperator):
         return A
 
-    elif isinstance(A, ndarray) or isinstance(A, matrix):
-        if len(A.shape) > 2:
+    elif isinstance(A, np.ndarray) or isinstance(A, np.matrix):
+        if A.ndim > 2:
             raise ValueError('array must have rank <= 2')
 
-        A = atleast_2d(asarray(A))
+        A = np.atleast_2d(np.asarray(A))
 
         def matvec(v):
-            return dot(A, v)
+            return np.dot(A, v)
         def rmatvec(v):
-            return dot(A.conj().transpose(), v)
+            return np.dot(A.conj().transpose(), v)
         def matmat(V):
-            return dot(A, V)
+            return np.dot(A, V)
         return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
                               matmat=matmat, dtype=A.dtype)
 
@@ -147,13 +243,13 @@
                               matmat=matmat, dtype=A.dtype)
 
     else:
-        if hasattr(A,'shape') and hasattr(A,'matvec'):
+        if hasattr(A, 'shape') and hasattr(A, 'matvec'):
             rmatvec = None
             dtype = None
 
-            if hasattr(A,'rmatvec'):
+            if hasattr(A, 'rmatvec'):
                 rmatvec = A.rmatvec
-            if hasattr(A,'dtype'):
+            if hasattr(A, 'dtype'):
                 dtype = A.dtype
             return LinearOperator(A.shape, A.matvec,
                                   rmatvec=rmatvec, dtype=dtype)

Modified: trunk/scipy/sparse/linalg/tests/test_interface.py
===================================================================
--- trunk/scipy/sparse/linalg/tests/test_interface.py	2008-11-22 20:18:12 UTC (rev 5170)
+++ trunk/scipy/sparse/linalg/tests/test_interface.py	2008-11-22 22:47:13 UTC (rev 5171)
@@ -1,62 +1,95 @@
-#!/usr/bin/env python
-""" Test functions for the sparse.linalg.interface module
+"""Test functions for the sparse.linalg.interface module
 """
 
 from numpy.testing import *
 
-import numpy
-from numpy import array, matrix, dtype
-from scipy.sparse import csr_matrix
+import numpy as np
+import scipy.sparse as sparse 
 
 from scipy.sparse.linalg.interface import *
 
 
-class TestInterface(TestCase):
-    def test_aslinearoperator(self):
-        cases = []
+class TestLinearOperator(TestCase):
+    def test_matvec(self):
+        def matvec(x):
+            # note, this matvec does not preserve type or shape
+            y = np.array([ 1*x[0] + 2*x[1] + 3*x[2],
+                           4*x[0] + 5*x[1] + 6*x[2]])
+            return y
 
-        cases.append( matrix([[1,2,3],[4,5,6]]) )
-        cases.append( array([[1,2,3],[4,5,6]]) )
-        cases.append( csr_matrix([[1,2,3],[4,5,6]]) )
+        A = LinearOperator((2,3), matvec)
 
-        class matlike:
-            def __init__(self):
-                self.dtype = dtype('int')
-                self.shape = (2,3)
-            def matvec(self,x):
-                y = array([ 1*x[0] + 2*x[1] + 3*x[2],
-                            4*x[0] + 5*x[1] + 6*x[2]])
-                if len(x.shape) == 2:
-                    y = y.reshape(-1,1)
-                return y
+        assert_equal(A.matvec(np.array([1,2,3])),       [14,32])
+        assert_equal(A.matvec(np.array([[1],[2],[3]])), [[14],[32]])
+        assert_equal(A * np.array([1,2,3]),             [14,32])
+        assert_equal(A * np.array([[1],[2],[3]]),       [[14],[32]])
+        
+        assert_equal(A.matvec(np.matrix([[1],[2],[3]])), [[14],[32]])
+        assert_equal(A * np.matrix([[1],[2],[3]]),       [[14],[32]])
 
-            def rmatvec(self,x):
-                return array([ 1*x[0] + 4*x[1],
-                               2*x[0] + 5*x[1],
-                               3*x[0] + 6*x[1]])
+        assert( isinstance(A.matvec(np.array([1,2,3])),       np.ndarray) )
+        assert( isinstance(A.matvec(np.array([[1],[2],[3]])), np.ndarray) )
+        assert( isinstance(A * np.array([1,2,3]),             np.ndarray) )
+        assert( isinstance(A * np.array([[1],[2],[3]]),       np.ndarray) )
 
-        cases.append( matlike() )
+        assert( isinstance(A.matvec(np.matrix([[1],[2],[3]])), np.ndarray) )
+        assert( isinstance(A * np.matrix([[1],[2],[3]]),       np.ndarray) )
 
+        assert_raises(ValueError, A.matvec, np.array([1,2]))
+        assert_raises(ValueError, A.matvec, np.array([1,2,3,4]))
+        assert_raises(ValueError, A.matvec, np.array([[1],[2]]))
+        assert_raises(ValueError, A.matvec, np.array([[1],[2],[3],[4]]))
 
-        for M in cases:
+
+class TestAsLinearOperator(TestCase):
+    def setUp(self):
+        self.cases = []
+
+        def make_cases(dtype):
+            self.cases.append( np.matrix([[1,2,3],[4,5,6]], dtype=dtype) )
+            self.cases.append( np.array([[1,2,3],[4,5,6]], dtype=dtype) )
+            self.cases.append( sparse.csr_matrix([[1,2,3],[4,5,6]], dtype=dtype) )
+
+            class matlike:
+                def __init__(self, dtype):
+                    self.dtype = np.dtype(dtype)
+                    self.shape = (2,3)
+                def matvec(self,x):
+                    y = np.array([ 1*x[0] + 2*x[1] + 3*x[2],
+                                   4*x[0] + 5*x[1] + 6*x[2]], dtype=self.dtype)
+                    if len(x.shape) == 2:
+                        y = y.reshape(-1,1)
+                    return y
+
+                def rmatvec(self,x):
+                    return np.array([ 1*x[0] + 4*x[1],
+                                      2*x[0] + 5*x[1],
+                                      3*x[0] + 6*x[1]], dtype=self.dtype)
+            self.cases.append( matlike('int') )
+
+        make_cases('int32')
+        make_cases('float32')
+        make_cases('float64')
+
+    def test_basic(self):
+
+        for M in self.cases:
             A = aslinearoperator(M)
             M,N = A.shape
 
-            assert_equal(A.matvec(array([1,2,3])),       [14,32])
-            assert_equal(A.matvec(array([[1],[2],[3]])), [[14],[32]])
+            assert_equal(A.matvec(np.array([1,2,3])),       [14,32])
+            assert_equal(A.matvec(np.array([[1],[2],[3]])), [[14],[32]])
 
-            assert_equal(A * array([1,2,3]),       [14,32])
-            assert_equal(A * array([[1],[2],[3]]), [[14],[32]])
+            assert_equal(A * np.array([1,2,3]),       [14,32])
+            assert_equal(A * np.array([[1],[2],[3]]), [[14],[32]])
 
-            assert_equal(A.rmatvec(array([1,2])),     [9,12,15])
-            assert_equal(A.rmatvec(array([[1],[2]])), [[9],[12],[15]])
+            assert_equal(A.rmatvec(np.array([1,2])),     [9,12,15])
+            assert_equal(A.rmatvec(np.array([[1],[2]])), [[9],[12],[15]])
 
-            assert_equal(A.matmat(array([[1,4],[2,5],[3,6]])), [[14,32],[32,77]] )
+            assert_equal(A.matmat(np.array([[1,4],[2,5],[3,6]])), [[14,32],[32,77]] )
 
-            assert_equal(A * array([[1,4],[2,5],[3,6]]), [[14,32],[32,77]] )
+            assert_equal(A * np.array([[1,4],[2,5],[3,6]]), [[14,32],[32,77]] )
 
             if hasattr(M,'dtype'):
                 assert_equal(A.dtype, M.dtype)
 
-if __name__ == "__main__":
-    nose.run(argv=['', __file__])

Modified: trunk/scipy/sparse/tests/test_base.py
===================================================================
--- trunk/scipy/sparse/tests/test_base.py	2008-11-22 20:18:12 UTC (rev 5170)
+++ trunk/scipy/sparse/tests/test_base.py	2008-11-22 22:47:13 UTC (rev 5171)
@@ -26,12 +26,12 @@
 import scipy.sparse as sparse
 from scipy.sparse import csc_matrix, csr_matrix, dok_matrix, \
         coo_matrix, lil_matrix, dia_matrix, bsr_matrix, \
-        eye, SparseEfficiencyWarning
+        eye, isspmatrix, SparseEfficiencyWarning
 from scipy.sparse.sputils import supported_dtypes
 from scipy.sparse.linalg import splu
 
 
-warnings.simplefilter('ignore',SparseEfficiencyWarning)
+warnings.simplefilter('ignore', SparseEfficiencyWarning)
 
 
 #TODO check that spmatrix( ... , copy=X ) is respected
@@ -347,7 +347,17 @@
         assert_array_almost_equal([1,2,3,4]*M, dot([1,2,3,4], M.toarray()))
         row = matrix([[1,2,3,4]])
         assert_array_almost_equal(row*M, row*M.todense())
+    
+    def test_small_multiplication(self):
+        """test that A*x works for x with shape () (1,) and (1,1)
+        """
+        A = self.spmatrix([[1],[2],[3]])
 
+        assert(isspmatrix(A * array(1)))
+        assert_equal((A * array(1)).todense(), [[1],[2],[3]])
+        assert_equal(A * array([1]), array([1,2,3]))
+        assert_equal(A * array([[1]]), array([[1],[2],[3]]))
+
     def test_matvec(self):
         M = self.spmatrix(matrix([[3,0,0],[0,1,0],[2,0,3.0],[2,3,0]]))
         col = matrix([1,2,3]).T



More information about the Scipy-svn mailing list