[Scipy-svn] r6286 - in trunk/scipy/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Mar 30 21:42:00 EDT 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