1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
"""Test functions for the sparse.linalg.interface module
"""
from numpy.testing import *
import numpy as np
import scipy.sparse as sparse
from scipy.sparse.linalg.interface import *
class TestLinearOperator(TestCase):
def setUp(self):
self.matvecs = []
# these matvecs do not preserve type or shape
def matvec1(x):
return np.array([ 1*x[0] + 2*x[1] + 3*x[2],
4*x[0] + 5*x[1] + 6*x[2]])
def matvec2(x):
return np.matrix(matvec1(x).reshape(2,1))
self.matvecs.append(matvec1)
self.matvecs.append(matvec2)
def test_matvec(self):
for matvec in self.matvecs:
A = LinearOperator((2,3), matvec)
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]])
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) )
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]]))
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(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.rmatvec(np.array([1,2])), [9,12,15])
assert_equal(A.rmatvec(np.array([[1],[2]])), [[9],[12],[15]])
assert_equal(A.matmat(np.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)
|
