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

scipy-svn@scip... scipy-svn@scip...
Tue Mar 30 20:42:00 CDT 2010


Author: warren.weckesser
Date: 2010-03-30 20:42:00 -0500 (Tue, 30 Mar 2010)
New Revision: 6286

Added:
   trunk/scipy/linalg/tests/test_special_matrices.py
Modified:
   trunk/scipy/linalg/basic.py
   trunk/scipy/linalg/special_matrices.py
   trunk/scipy/linalg/tests/test_basic.py
Log:
ENH: Added the function hadamard() to linalg (ticket #675).
Cleaned up unused and duplicated imports in basic.py.
Moved tests of special matrices to their own file.
Added a test for kron.

Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py	2010-03-30 04:39:41 UTC (rev 6285)
+++ trunk/scipy/linalg/basic.py	2010-03-31 01:42:00 UTC (rev 6286)
@@ -7,29 +7,30 @@
 #
 # w/ additions by Travis Oliphant, March 2002
 
-__all__ = ['solve', 'inv', 'det', 'lstsq', 'norm', 'pinv', 'pinv2',
-           'tri','tril', 'triu', 'toeplitz', 'circulant', 'hankel', 'block_diag',
-           'lu_solve', 'cho_solve', 'solve_banded', 'LinAlgError', 'kron',
-           'all_mat', 'cholesky_banded', 'solveh_banded']
+__all__ = ['solve', 'inv', 'det', 'lstsq', 'pinv', 'pinv2',
+           'cholesky_banded', 'solveh_banded', 'lu_solve', 'cho_solve',
+           'solve_banded',
+           #  From special_matrices:
+           'tri','tril', 'triu', 'toeplitz', 'circulant', 'hankel', 'kron',
+           'hadamard', 'block_diag', 'all_mat',
+           #  From misc:
+           'LinAlgError', 'norm',
+            ]
 
+from numpy import asarray, zeros, sum, conjugate, dot, transpose, \
+        asarray_chkfinite,  single
+import numpy
+
 #from blas import get_blas_funcs
 from flinalg import get_flinalg_funcs
 from lapack import get_lapack_funcs
-from numpy import asarray, zeros, sum, greater_equal, subtract, arange,\
-     conjugate, dot, transpose
-import numpy
-from numpy import asarray_chkfinite, atleast_2d, outer, concatenate, reshape, single
-from numpy import matrix as Matrix
 from misc import LinAlgError, norm
 from special_matrices import tri, tril, triu, toeplitz, circulant, hankel, \
-        kron, block_diag, all_mat
+        hadamard, kron, block_diag, all_mat
 from scipy.linalg import calc_lwork
-
-from misc import LinAlgError, norm
-from special_matrices import tri, tril, triu, toeplitz, circulant, hankel, \
-    block_diag, kron
 import decomp
 
+
 def lu_solve((lu, piv), b, trans=0, overwrite_b=0):
     """Solve an equation system, a x = b, given the LU factorization of a
 

Modified: trunk/scipy/linalg/special_matrices.py
===================================================================
--- trunk/scipy/linalg/special_matrices.py	2010-03-30 04:39:41 UTC (rev 6285)
+++ trunk/scipy/linalg/special_matrices.py	2010-03-31 01:42:00 UTC (rev 6286)
@@ -1,3 +1,5 @@
+
+import math
 import numpy as np
 
 #-----------------------------------------------------------------------------
@@ -265,6 +267,59 @@
     # `vals[indx]` is the Hankel matrix.
     return vals[indx]
 
+def hadamard(n, dtype=int):
+    """Construct a Hadamard matrix.
+    
+    `hadamard(n)` constructs an n-by-n Hadamard matrix, using Sylvester's
+    construction.  `n` must be a power of 2.
+
+    Parameters
+    ----------
+    n : int
+        The order of the matrix.  `n` must be a power of 2.
+    dtype : numpy dtype
+        The data type of the array to be constructed.
+        
+    Returns
+    -------
+    H : ndarray with shape (n, n)
+        The Hadamard matrix.
+
+    Examples
+    --------
+    >>> hadamard(2, dtype=complex)
+    array([[ 1.+0.j,  1.+0.j],
+           [ 1.+0.j, -1.-0.j]])
+    >>> hadamard(4)
+    array([[ 1,  1,  1,  1],
+           [ 1, -1,  1, -1],
+           [ 1,  1, -1, -1],
+           [ 1, -1, -1,  1]])
+
+    Notes
+    -----
+    .. versionadded:: 0.8.0
+
+    """
+    
+    # This function is a slightly modified version of the
+    # function contributed by Ivo in ticket #675.
+
+    if n < 1:
+        lg2 = 0
+    else:
+        lg2 = int(math.log(n, 2))
+    if 2 ** lg2 != n:
+        raise ValueError("n must be an positive integer, and n must be power of 2")
+
+    H = np.array([[1]], dtype=dtype)
+
+    # Sylvester's construction
+    for i in range(0, lg2): 
+        H = np.vstack((np.hstack((H, H)), np.hstack((H, -H))))
+
+    return H
+
 def all_mat(*args):
     return map(np.matrix,args)
 
@@ -297,12 +352,12 @@
 
     """
     if not a.flags['CONTIGUOUS']:
-        a = reshape(a, a.shape)
+        a = np.reshape(a, a.shape)
     if not b.flags['CONTIGUOUS']:
-        b = reshape(b, b.shape)
-    o = outer(a,b)
-    o=o.reshape(a.shape + b.shape)
-    return concatenate(concatenate(o, axis=1), axis=1)
+        b = np.reshape(b, b.shape)
+    o = np.outer(a,b)
+    o = o.reshape(a.shape + b.shape)
+    return np.concatenate(np.concatenate(o, axis=1), axis=1)
 
 def block_diag(*arrs):
     """Create a block diagonal matrix from the provided arrays.

Modified: trunk/scipy/linalg/tests/test_basic.py
===================================================================
--- trunk/scipy/linalg/tests/test_basic.py	2010-03-30 04:39:41 UTC (rev 6285)
+++ trunk/scipy/linalg/tests/test_basic.py	2010-03-31 01:42:00 UTC (rev 6286)
@@ -19,23 +19,17 @@
   python tests/test_basic.py
 """
 
-from numpy import arange, add, array, dot, zeros, identity, conjugate, \
-     transpose, eye, all, copy
+from numpy import arange, array, dot, zeros, identity, conjugate, transpose
 import numpy.linalg as linalg
 
 from numpy.testing import *
 
-from scipy.linalg import solve,inv,det,lstsq, toeplitz, hankel, circulant, \
-     tri, triu, tril, pinv, pinv2, solve_banded, block_diag, norm
+from scipy.linalg import solve, inv, det, lstsq, pinv, pinv2, solve_banded, norm
 
 
 def random(size):
     return rand(*size)
 
-def get_mat(n):
-    data = arange(n)
-    data = add.outer(data,data)
-    return data
 
 class TestSolveBanded(TestCase):
 
@@ -305,182 +299,9 @@
             #XXX: check definition of res
             assert_array_almost_equal(x,direct_lstsq(a,b,1))
 
-class TestTri(TestCase):
-    def test_basic(self):
-        assert_equal(tri(4),array([[1,0,0,0],
-                                   [1,1,0,0],
-                                   [1,1,1,0],
-                                   [1,1,1,1]]))
-        assert_equal(tri(4,dtype='f'),array([[1,0,0,0],
-                                                [1,1,0,0],
-                                                [1,1,1,0],
-                                                [1,1,1,1]],'f'))
-    def test_diag(self):
-        assert_equal(tri(4,k=1),array([[1,1,0,0],
-                                       [1,1,1,0],
-                                       [1,1,1,1],
-                                       [1,1,1,1]]))
-        assert_equal(tri(4,k=-1),array([[0,0,0,0],
-                                        [1,0,0,0],
-                                        [1,1,0,0],
-                                        [1,1,1,0]]))
-    def test_2d(self):
-        assert_equal(tri(4,3),array([[1,0,0],
-                                     [1,1,0],
-                                     [1,1,1],
-                                     [1,1,1]]))
-        assert_equal(tri(3,4),array([[1,0,0,0],
-                                     [1,1,0,0],
-                                     [1,1,1,0]]))
-    def test_diag2d(self):
-        assert_equal(tri(3,4,k=2),array([[1,1,1,0],
-                                         [1,1,1,1],
-                                         [1,1,1,1]]))
-        assert_equal(tri(4,3,k=-2),array([[0,0,0],
-                                          [0,0,0],
-                                          [1,0,0],
-                                          [1,1,0]]))
 
-class TestTril(TestCase):
-    def test_basic(self):
-        a = (100*get_mat(5)).astype('l')
-        b = a.copy()
-        for k in range(5):
-            for l in range(k+1,5):
-                b[k,l] = 0
-        assert_equal(tril(a),b)
 
-    def test_diag(self):
-        a = (100*get_mat(5)).astype('f')
-        b = a.copy()
-        for k in range(5):
-            for l in range(k+3,5):
-                b[k,l] = 0
-        assert_equal(tril(a,k=2),b)
-        b = a.copy()
-        for k in range(5):
-            for l in range(max((k-1,0)),5):
-                b[k,l] = 0
-        assert_equal(tril(a,k=-2),b)
 
-class TestTriu(TestCase):
-    def test_basic(self):
-        a = (100*get_mat(5)).astype('l')
-        b = a.copy()
-        for k in range(5):
-            for l in range(k+1,5):
-                b[l,k] = 0
-        assert_equal(triu(a),b)
-
-    def test_diag(self):
-        a = (100*get_mat(5)).astype('f')
-        b = a.copy()
-        for k in range(5):
-            for l in range(max((k-1,0)),5):
-                b[l,k] = 0
-        assert_equal(triu(a,k=2),b)
-        b = a.copy()
-        for k in range(5):
-            for l in range(k+3,5):
-                b[l,k] = 0
-        assert_equal(triu(a,k=-2),b)
-
-
-class TestToeplitz(TestCase):
-    
-    def test_basic(self):
-        y = toeplitz([1,2,3])
-        assert_array_equal(y,[[1,2,3],[2,1,2],[3,2,1]])
-        y = toeplitz([1,2,3],[1,4,5])
-        assert_array_equal(y,[[1,4,5],[2,1,4],[3,2,1]])
-        
-    def test_complex_01(self):
-        data = (1.0 + arange(3.0)) * (1.0 + 1.0j)
-        x = copy(data)
-        t = toeplitz(x)
-        # Calling toeplitz should not change x.
-        assert_array_equal(x, data)
-        # According to the docstring, x should be the first column of t.
-        col0 = t[:,0]
-        assert_array_equal(col0, data)
-        assert_array_equal(t[0,1:], data[1:].conj())
-
-    def test_scalar_00(self):
-        """Scalar arguments still produce a 2D array."""
-        t = toeplitz(10)
-        assert_array_equal(t, [[10]])
-        t = toeplitz(10, 20)
-        assert_array_equal(t, [[10]])
-        
-    def test_scalar_01(self):
-        c = array([1,2,3])
-        t = toeplitz(c, 1)
-        assert_array_equal(t, [[1],[2],[3]])
-
-    def test_scalar_02(self):
-        c = array([1,2,3])
-        t = toeplitz(c, array(1))
-        assert_array_equal(t, [[1],[2],[3]])
-
-    def test_scalar_03(self):
-        c = array([1,2,3])
-        t = toeplitz(c, array([1]))
-        assert_array_equal(t, [[1],[2],[3]])
-
-    def test_scalar_04(self):
-        r = array([10,2,3])
-        t = toeplitz(1, r)
-        assert_array_equal(t, [[1,2,3]])
-
-
-class TestHankel(TestCase):
-    def test_basic(self):
-        y = hankel([1,2,3])
-        assert_array_equal(y,[[1,2,3],[2,3,0],[3,0,0]])
-        y = hankel([1,2,3],[3,4,5])
-        assert_array_equal(y,[[1,2,3],[2,3,4],[3,4,5]])
-
-
-class TestCirculant(TestCase):
-    def test_basic(self):
-        y = circulant([1,2,3])
-        assert_array_equal(y,[[1,3,2],[2,1,3],[3,2,1]])
-
-
-class TestBlockDiag:
-    def test_basic(self):
-        x = block_diag(eye(2), [[1,2], [3,4], [5,6]], [[1, 2, 3]])
-        assert all(x == [[1, 0, 0, 0, 0, 0, 0],
-                         [0, 1, 0, 0, 0, 0, 0],
-                         [0, 0, 1, 2, 0, 0, 0],
-                         [0, 0, 3, 4, 0, 0, 0],
-                         [0, 0, 5, 6, 0, 0, 0],
-                         [0, 0, 0, 0, 1, 2, 3]])
-
-    def test_dtype(self):
-        x = block_diag([[1.5]])
-        assert_equal(x.dtype, float)
-
-        x = block_diag([[True]])
-        assert_equal(x.dtype, bool)
-        
-    def test_scalar_and_1d_args(self):
-        a = block_diag(1)
-        assert_equal(a.shape, (1,1))
-        assert_array_equal(a, [[1]])
-        
-        a = block_diag([2,3], 4)
-        assert_array_equal(a, [[2, 3, 0], [0, 0, 4]])
-
-    def test_bad_arg(self):
-        assert_raises(ValueError, block_diag, [[[1]]])
-
-    def test_no_args(self):
-        a = block_diag()
-        assert_equal(a.ndim, 2)
-        assert_equal(a.nbytes, 0)
-
-
 class TestPinv(TestCase):
 
     def test_simple(self):

Added: trunk/scipy/linalg/tests/test_special_matrices.py
===================================================================
--- trunk/scipy/linalg/tests/test_special_matrices.py	                        (rev 0)
+++ trunk/scipy/linalg/tests/test_special_matrices.py	2010-03-31 01:42:00 UTC (rev 6286)
@@ -0,0 +1,229 @@
+"""Tests for functions in special_matrices.py."""
+
+from numpy import arange, add, array, eye, all, copy
+from numpy.testing import *
+
+from scipy.linalg import toeplitz, hankel, circulant, hadamard, tri, triu, tril, \
+                            kron, block_diag
+
+
+def get_mat(n):
+    data = arange(n)
+    data = add.outer(data,data)
+    return data
+
+
+class TestTri(TestCase):
+    def test_basic(self):
+        assert_equal(tri(4),array([[1,0,0,0],
+                                   [1,1,0,0],
+                                   [1,1,1,0],
+                                   [1,1,1,1]]))
+        assert_equal(tri(4,dtype='f'),array([[1,0,0,0],
+                                                [1,1,0,0],
+                                                [1,1,1,0],
+                                                [1,1,1,1]],'f'))
+    def test_diag(self):
+        assert_equal(tri(4,k=1),array([[1,1,0,0],
+                                       [1,1,1,0],
+                                       [1,1,1,1],
+                                       [1,1,1,1]]))
+        assert_equal(tri(4,k=-1),array([[0,0,0,0],
+                                        [1,0,0,0],
+                                        [1,1,0,0],
+                                        [1,1,1,0]]))
+    def test_2d(self):
+        assert_equal(tri(4,3),array([[1,0,0],
+                                     [1,1,0],
+                                     [1,1,1],
+                                     [1,1,1]]))
+        assert_equal(tri(3,4),array([[1,0,0,0],
+                                     [1,1,0,0],
+                                     [1,1,1,0]]))
+    def test_diag2d(self):
+        assert_equal(tri(3,4,k=2),array([[1,1,1,0],
+                                         [1,1,1,1],
+                                         [1,1,1,1]]))
+        assert_equal(tri(4,3,k=-2),array([[0,0,0],
+                                          [0,0,0],
+                                          [1,0,0],
+                                          [1,1,0]]))
+
+class TestTril(TestCase):
+    def test_basic(self):
+        a = (100*get_mat(5)).astype('l')
+        b = a.copy()
+        for k in range(5):
+            for l in range(k+1,5):
+                b[k,l] = 0
+        assert_equal(tril(a),b)
+
+    def test_diag(self):
+        a = (100*get_mat(5)).astype('f')
+        b = a.copy()
+        for k in range(5):
+            for l in range(k+3,5):
+                b[k,l] = 0
+        assert_equal(tril(a,k=2),b)
+        b = a.copy()
+        for k in range(5):
+            for l in range(max((k-1,0)),5):
+                b[k,l] = 0
+        assert_equal(tril(a,k=-2),b)
+
+
+class TestTriu(TestCase):
+    def test_basic(self):
+        a = (100*get_mat(5)).astype('l')
+        b = a.copy()
+        for k in range(5):
+            for l in range(k+1,5):
+                b[l,k] = 0
+        assert_equal(triu(a),b)
+
+    def test_diag(self):
+        a = (100*get_mat(5)).astype('f')
+        b = a.copy()
+        for k in range(5):
+            for l in range(max((k-1,0)),5):
+                b[l,k] = 0
+        assert_equal(triu(a,k=2),b)
+        b = a.copy()
+        for k in range(5):
+            for l in range(k+3,5):
+                b[l,k] = 0
+        assert_equal(triu(a,k=-2),b)
+
+
+class TestToeplitz(TestCase):
+    
+    def test_basic(self):
+        y = toeplitz([1,2,3])
+        assert_array_equal(y,[[1,2,3],[2,1,2],[3,2,1]])
+        y = toeplitz([1,2,3],[1,4,5])
+        assert_array_equal(y,[[1,4,5],[2,1,4],[3,2,1]])
+        
+    def test_complex_01(self):
+        data = (1.0 + arange(3.0)) * (1.0 + 1.0j)
+        x = copy(data)
+        t = toeplitz(x)
+        # Calling toeplitz should not change x.
+        assert_array_equal(x, data)
+        # According to the docstring, x should be the first column of t.
+        col0 = t[:,0]
+        assert_array_equal(col0, data)
+        assert_array_equal(t[0,1:], data[1:].conj())
+
+    def test_scalar_00(self):
+        """Scalar arguments still produce a 2D array."""
+        t = toeplitz(10)
+        assert_array_equal(t, [[10]])
+        t = toeplitz(10, 20)
+        assert_array_equal(t, [[10]])
+        
+    def test_scalar_01(self):
+        c = array([1,2,3])
+        t = toeplitz(c, 1)
+        assert_array_equal(t, [[1],[2],[3]])
+
+    def test_scalar_02(self):
+        c = array([1,2,3])
+        t = toeplitz(c, array(1))
+        assert_array_equal(t, [[1],[2],[3]])
+
+    def test_scalar_03(self):
+        c = array([1,2,3])
+        t = toeplitz(c, array([1]))
+        assert_array_equal(t, [[1],[2],[3]])
+
+    def test_scalar_04(self):
+        r = array([10,2,3])
+        t = toeplitz(1, r)
+        assert_array_equal(t, [[1,2,3]])
+
+
+class TestHankel(TestCase):
+    def test_basic(self):
+        y = hankel([1,2,3])
+        assert_array_equal(y, [[1,2,3], [2,3,0], [3,0,0]])
+        y = hankel([1,2,3], [3,4,5])
+        assert_array_equal(y, [[1,2,3], [2,3,4], [3,4,5]])
+
+
+class TestCirculant(TestCase):
+    def test_basic(self):
+        y = circulant([1,2,3])
+        assert_array_equal(y, [[1,3,2], [2,1,3], [3,2,1]])
+
+
+class TestHadamard(TestCase):
+
+    def test_basic(self):
+
+        y = hadamard(1)
+        assert_array_equal(y, [[1]])
+
+        y = hadamard(2, dtype=float)
+        assert_array_equal(y, [[1.0, 1.0], [1.0, -1.0]])
+
+        y = hadamard(4)        
+        assert_array_equal(y, [[1,1,1,1], [1,-1,1,-1], [1,1,-1,-1], [1,-1,-1,1]])
+
+        assert_raises(ValueError, hadamard, 0)
+        assert_raises(ValueError, hadamard, 5)
+
+
+class TestBlockDiag:
+    def test_basic(self):
+        x = block_diag(eye(2), [[1,2], [3,4], [5,6]], [[1, 2, 3]])
+        assert all(x == [[1, 0, 0, 0, 0, 0, 0],
+                         [0, 1, 0, 0, 0, 0, 0],
+                         [0, 0, 1, 2, 0, 0, 0],
+                         [0, 0, 3, 4, 0, 0, 0],
+                         [0, 0, 5, 6, 0, 0, 0],
+                         [0, 0, 0, 0, 1, 2, 3]])
+
+    def test_dtype(self):
+        x = block_diag([[1.5]])
+        assert_equal(x.dtype, float)
+
+        x = block_diag([[True]])
+        assert_equal(x.dtype, bool)
+        
+    def test_scalar_and_1d_args(self):
+        a = block_diag(1)
+        assert_equal(a.shape, (1,1))
+        assert_array_equal(a, [[1]])
+        
+        a = block_diag([2,3], 4)
+        assert_array_equal(a, [[2, 3, 0], [0, 0, 4]])
+
+    def test_bad_arg(self):
+        assert_raises(ValueError, block_diag, [[[1]]])
+
+    def test_no_args(self):
+        a = block_diag()
+        assert_equal(a.ndim, 2)
+        assert_equal(a.nbytes, 0)
+
+
+class TestKron:
+
+    def test_basic(self):
+
+        a = kron(array([[1, 2], [3, 4]]), array([[1, 1, 1]]))
+        assert_array_equal(a, array([[1, 1, 1, 2, 2, 2],
+                                     [3, 3, 3, 4, 4, 4]]))
+
+        m1 = array([[1, 2], [3, 4]])
+        m2 = array([[10], [11]])
+        a = kron(m1, m2)
+        expected = array([[ 10, 20 ],
+                          [ 11, 22 ],
+                          [ 30, 40 ],
+                          [ 33, 44 ]])
+        assert_array_equal(a, expected)
+
+
+if __name__ == "__main__":
+    run_module_suite()



More information about the Scipy-svn mailing list