[Scipy-svn] r6267 - in trunk/scipy/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Mar 23 22:39:33 EDT 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