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

scipy-svn@scip... scipy-svn@scip...
Tue Mar 23 21:39:33 CDT 2010


Author: warren.weckesser
Date: 2010-03-23 21:39:31 -0500 (Tue, 23 Mar 2010)
New Revision: 6267

Modified:
   trunk/scipy/linalg/basic.py
   trunk/scipy/linalg/tests/test_basic.py
Log:
BUG+ENH: Updated linalg.toeplitz to fix problems reported in ticket #667.  Updated linalg.hankel to have behavior consistent with toeplitz.  Added a new function, linalg.circulant.

Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py	2010-03-21 21:16:09 UTC (rev 6266)
+++ trunk/scipy/linalg/basic.py	2010-03-24 02:39:31 UTC (rev 6267)
@@ -7,16 +7,16 @@
 #
 # w/ additions by Travis Oliphant, March 2002
 
-__all__ = ['solve','inv','det','lstsq','norm','pinv','pinv2',
-           'tri','tril','triu','toeplitz','hankel','block_diag','lu_solve',
-           'cho_solve','solve_banded','LinAlgError','kron',
+__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']
 
 #from blas import get_blas_funcs
 from flinalg import get_flinalg_funcs
 from lapack import get_lapack_funcs
-from numpy import asarray,zeros,sum,newaxis,greater_equal,subtract,arange,\
-     conjugate,ravel,r_,mgrid,take,ones,dot,transpose
+from numpy import asarray, zeros, sum, greater_equal, subtract, arange,\
+     conjugate, dot, transpose
 import numpy
 from numpy import asarray_chkfinite, outer, concatenate, reshape, single
 from numpy import matrix as Matrix
@@ -537,6 +537,7 @@
 feps = numpy.finfo(single).eps
 
 _array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
+
 def pinv2(a, cond=None, rcond=None):
     """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
 
@@ -695,24 +696,31 @@
     out = (1-tri(m.shape[0], m.shape[1], k-1, m.dtype.char))*m
     return out
 
-def toeplitz(c,r=None):
+
+def toeplitz(c, r=None):
     """Construct a Toeplitz matrix.
 
-    The Toepliz matrix has constant diagonals, c as its first column,
-    and r as its first row (if not given, r == c is assumed).
+    The Toepliz matrix has constant diagonals, with c as its first column
+    and r as its first row.  If r is not given, r == conjugate(c) is
+    assumed.
 
     Parameters
     ----------
-    c : array
-        First column of the matrix
-    r : array
-        First row of the matrix. If None, r == c is assumed.
+    c : array-like, 1D
+        First column of the matrix.  Whatever the actual shape of `c`, it
+        will be converted to a 1D array.
+    r : array-like, 1D
+        First row of the matrix. If None, `r = conjugate(c)` is assumed; in
+        this case, if `c[0]` is real, the result is a Hermitian matrix.
+        `r[0]` is ignored; the first row of the returned matrix is
+        `[c[0], r[1:]]`.  Whatever the actual shape of `r`, it will be
+        converted to a 1D array.
 
     Returns
     -------
     A : array, shape (len(c), len(r))
-        Constructed Toeplitz matrix.
-        dtype is the same as (c[0] + r[0]).dtype
+        The Toeplitz matrix.
+        dtype is the same as `(c[0] + r[0]).dtype`.
 
     Examples
     --------
@@ -721,54 +729,105 @@
     array([[1, 4, 5, 6],
            [2, 1, 4, 5],
            [3, 2, 1, 4]])
+    >>> toeplitz([1.0, 2+3j, 4-1j])
+    array([[ 1.+0.j,  2.-3.j,  4.+1.j],
+           [ 2.+3.j,  1.+0.j,  2.-3.j],
+           [ 4.-1.j,  2.+3.j,  1.+0.j]])
 
     See also
     --------
+    circulant : circulant matrix
     hankel : Hankel matrix
 
+    Notes
+    -----
+    The behavior when `c` or `r` is a scalar, or when `c` is complex and
+    `r` is None, was changed in version 0.8.0.  The behavior in previous
+    versions was undocumented and is no longer supported. 
     """
-    isscalar = numpy.isscalar
-    if isscalar(c) or isscalar(r):
-        return c
+    c = numpy.asarray(c).ravel()
     if r is None:
-        r = c
-        r[0] = conjugate(r[0])
-        c = conjugate(c)
-    r,c = map(asarray_chkfinite,(r,c))
-    r,c = map(ravel,(r,c))
-    rN,cN = map(len,(r,c))
-    if r[0] != c[0]:
-        print "Warning: column and row values don't agree; column value used."
-    vals = r_[r[rN-1:0:-1], c]
-    cols = mgrid[0:cN]
-    rows = mgrid[rN:0:-1]
-    indx = cols[:,newaxis]*ones((1,rN),dtype=int) + \
-           rows[newaxis,:]*ones((cN,1),dtype=int) - 1
-    return take(vals, indx, 0)
+        r = c.conjugate()
+    else:
+        r = numpy.asarray(r).ravel()
+    # Form a 1D array of values to be used in the matrix, containing a reversed
+    # copy of r[1:], followed by c.
+    vals = numpy.concatenate((r[-1:0:-1], c))
+    a, b = numpy.ogrid[0:len(c), len(r)-1:-1:-1]
+    indx = a + b
+    # `indx` is a 2D array of indices into the 1D array `vals`, arranged so that
+    # `vals[indx]` is the Toeplitz matrix.
+    return vals[indx]
 
+def circulant(c):
+    """Construct a circulant matrix.
 
-def hankel(c,r=None):
+    Parameters
+    ----------
+    c : array-like, 1D
+        First column of the matrix.
+
+    Returns
+    -------
+    A : array, shape (len(c), len(c))
+        A circulant matrix whose first column is `c`.
+
+    Examples
+    --------
+    >>> from scipy.linalg import circulant
+    >>> circulant([1, 2, 3])
+    array([[1, 3, 2],
+           [2, 1, 3],
+           [3, 2, 1]])
+
+    See also
+    --------
+    toeplitz : Toeplitz matrix
+    hankel : Hankel matrix
+    
+    Notes
+    -----
+    .. versionadded:: 0.8.0
+
+    """
+    c = numpy.asarray(c).ravel()
+    a, b = numpy.ogrid[0:len(c), 0:-len(c):-1]
+    indx = a + b
+    # `indx` is a 2D array of indices into `c`, arranged so that `c[indx]` is
+    # the circulant matrix.
+    return c[indx]
+
+def hankel(c, r=None):
     """Construct a Hankel matrix.
 
-    The Hankel matrix has constant anti-diagonals, c as its first column,
-    and r as its last row (if not given, r == 0 os assumed).
+    The Hankel matrix has constant anti-diagonals, with `c` as its
+    first column and `r` as its last row.  If `r` is not given, then
+    `r = zeros_like(c)` is assumed.
 
     Parameters
     ----------
-    c : array
-        First column of the matrix
-    r : array
-        Last row of the matrix. If None, r == 0 is assumed.
+    c : array-like, 1D
+        First column of the matrix.  Whatever the actual shape of `c`, it
+        will be converted to a 1D array.
+    r : array-like, 1D
+        Last row of the matrix. If None, `r` == 0 is assumed.
+        `r[0]` is ignored; the last row of the returned matrix is
+        `[c[0], r[1:]]`.  Whatever the actual shape of `r`, it will be
+        converted to a 1D array.
 
     Returns
     -------
     A : array, shape (len(c), len(r))
-        Constructed Hankel matrix.
-        dtype is the same as (c[0] + r[0]).dtype
+        The Hankel matrix.
+        dtype is the same as `(c[0] + r[0]).dtype`.
 
     Examples
     --------
     >>> from scipy.linalg import hankel
+    >>> hankel([1, 17, 99])
+    array([[ 1, 17, 99],
+           [17, 99,  0],
+           [99,  0,  0]])
     >>> hankel([1,2,3,4], [4,7,7,8,9])
     array([[1, 2, 3, 4, 7],
            [2, 3, 4, 7, 7],
@@ -778,24 +837,22 @@
     See also
     --------
     toeplitz : Toeplitz matrix
+    circulant : circulant matrix
 
     """
-    isscalar = numpy.isscalar
-    if isscalar(c) or isscalar(r):
-        return c
+    c = numpy.asarray(c).ravel()
     if r is None:
-        r = zeros(len(c))
-    elif r[0] != c[-1]:
-        print "Warning: column and row values don't agree; column value used."
-    r,c = map(asarray_chkfinite,(r,c))
-    r,c = map(ravel,(r,c))
-    rN,cN = map(len,(r,c))
-    vals = r_[c, r[1:rN]]
-    cols = mgrid[1:cN+1]
-    rows = mgrid[0:rN]
-    indx = cols[:,newaxis]*ones((1,rN),dtype=int) + \
-           rows[newaxis,:]*ones((cN,1),dtype=int) - 1
-    return take(vals, indx, 0)
+        r = numpy.zeros_like(c)
+    else:
+        r = numpy.asarray(r).ravel()
+    # Form a 1D array of values to be used in the matrix, containing `c`
+    # followed by r[1:].
+    vals = numpy.concatenate((c, r[1:]))
+    a, b = numpy.ogrid[0:len(c), 0:len(r)]
+    indx = a + b
+    # `indx` is a 2D array of indices into the 1D array `vals`, arranged so that
+    # `vals[indx]` is the Hankel matrix.
+    return vals[indx]
 
 def all_mat(*args):
     return map(Matrix,args)

Modified: trunk/scipy/linalg/tests/test_basic.py
===================================================================
--- trunk/scipy/linalg/tests/test_basic.py	2010-03-21 21:16:09 UTC (rev 6266)
+++ trunk/scipy/linalg/tests/test_basic.py	2010-03-24 02:39:31 UTC (rev 6267)
@@ -20,13 +20,13 @@
 """
 
 from numpy import arange, add, array, dot, zeros, identity, conjugate, \
-     transpose, eye, all
+     transpose, eye, all, copy
 import numpy.linalg as linalg
 
 from numpy.testing import *
 
-from scipy.linalg import solve,inv,det,lstsq, toeplitz, hankel, tri, triu, \
-     tril, pinv, pinv2, solve_banded, block_diag, norm
+from scipy.linalg import solve,inv,det,lstsq, toeplitz, hankel, circulant, \
+     tri, triu, tril, pinv, pinv2, solve_banded, block_diag, norm
 
 
 def random(size):
@@ -385,13 +385,54 @@
                 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])
@@ -399,6 +440,13 @@
         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]])



More information about the Scipy-svn mailing list