[Scipy-svn] r5171 - in trunk/scipy/sparse: . linalg linalg/tests tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Nov 22 17:47:23 EST 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