[Scipy-svn] r6284 - trunk/scipy/linalg
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Mar 30 00:39:33 EDT 2010
Author: cdavid
Date: 2010-03-29 23:39:32 -0500 (Mon, 29 Mar 2010)
New Revision: 6284
Added:
trunk/scipy/linalg/misc.py
trunk/scipy/linalg/special_matrices.py
Modified:
trunk/scipy/linalg/basic.py
trunk/scipy/linalg/decomp.py
Log:
PY3K: split linalg decomp/basic to avoid circular imports.
2to3 and python 3 are stricted in terms of circular imports, and that's
bad practice in general.
Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py 2010-03-30 04:39:23 UTC (rev 6283)
+++ trunk/scipy/linalg/basic.py 2010-03-30 04:39:32 UTC (rev 6284)
@@ -20,9 +20,14 @@
import numpy
from numpy import asarray_chkfinite, atleast_2d, outer, concatenate, reshape, single
from numpy import matrix as Matrix
-from numpy.linalg import LinAlgError
+from misc import LinAlgError, norm
+from special_matrices import tri, tril, triu, toeplitz, circulant, hankel, \
+ 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):
@@ -386,13 +391,6 @@
return inv_a
-### Norm
-
-def norm(a, ord=None):
- # Differs from numpy only in non-finite handling
- return numpy.linalg.norm(asarray_chkfinite(a), ord=ord)
-norm.__doc__ = numpy.linalg.norm.__doc__
-
### Determinant
def det(a, overwrite_a=0):
@@ -588,376 +586,3 @@
psigma[i,i] = 1.0/conjugate(s[i])
#XXX: use lapack/blas routines for dot
return transpose(conjugate(dot(dot(u,psigma),vh)))
-
-#-----------------------------------------------------------------------------
-# matrix construction functions
-#-----------------------------------------------------------------------------
-
-def tri(N, M=None, k=0, dtype=None):
- """Construct (N, M) matrix filled with ones at and below the k-th diagonal.
-
- The matrix has A[i,j] == 1 for i <= j + k
-
- Parameters
- ----------
- N : integer
- M : integer
- Size of the matrix. If M is None, M == N is assumed.
- k : integer
- Number of subdiagonal below which matrix is filled with ones.
- k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
- dtype : dtype
- Data type of the matrix.
-
- Returns
- -------
- A : array, shape (N, M)
-
- Examples
- --------
- >>> from scipy.linalg import tri
- >>> tri(3, 5, 2, dtype=int)
- array([[1, 1, 1, 0, 0],
- [1, 1, 1, 1, 0],
- [1, 1, 1, 1, 1]])
- >>> tri(3, 5, -1, dtype=int)
- array([[0, 0, 0, 0, 0],
- [1, 0, 0, 0, 0],
- [1, 1, 0, 0, 0]])
-
- """
- if M is None: M = N
- if type(M) == type('d'):
- #pearu: any objections to remove this feature?
- # As tri(N,'d') is equivalent to tri(N,dtype='d')
- dtype = M
- M = N
- m = greater_equal(subtract.outer(arange(N), arange(M)),-k)
- if dtype is None:
- return m
- else:
- return m.astype(dtype)
-
-def tril(m, k=0):
- """Construct a copy of a matrix with elements above the k-th diagonal zeroed.
-
- Parameters
- ----------
- m : array
- Matrix whose elements to return
- k : integer
- Diagonal above which to zero elements.
- k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
-
- Returns
- -------
- A : array, shape m.shape, dtype m.dtype
-
- Examples
- --------
- >>> from scipy.linalg import tril
- >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
- array([[ 0, 0, 0],
- [ 4, 0, 0],
- [ 7, 8, 0],
- [10, 11, 12]])
-
- """
- m = asarray(m)
- out = tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype.char)*m
- return out
-
-def triu(m, k=0):
- """Construct a copy of a matrix with elements below the k-th diagonal zeroed.
-
- Parameters
- ----------
- m : array
- Matrix whose elements to return
- k : integer
- Diagonal below which to zero elements.
- k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
-
- Returns
- -------
- A : array, shape m.shape, dtype m.dtype
-
- Examples
- --------
- >>> from scipy.linalg import tril
- >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
- array([[ 1, 2, 3],
- [ 4, 5, 6],
- [ 0, 8, 9],
- [ 0, 0, 12]])
-
- """
- m = asarray(m)
- out = (1-tri(m.shape[0], m.shape[1], k-1, m.dtype.char))*m
- return out
-
-
-def toeplitz(c, r=None):
- """Construct a Toeplitz matrix.
-
- 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-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))
- The Toeplitz matrix.
- dtype is the same as `(c[0] + r[0]).dtype`.
-
- Examples
- --------
- >>> from scipy.linalg import toeplitz
- >>> toeplitz([1,2,3], [1,4,5,6])
- 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.
- """
- c = numpy.asarray(c).ravel()
- if r is None:
- 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.
-
- 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, 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-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))
- 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],
- [3, 4, 7, 7, 8],
- [4, 7, 7, 8, 9]])
-
- See also
- --------
- toeplitz : Toeplitz matrix
- circulant : circulant matrix
-
- """
- c = numpy.asarray(c).ravel()
- if r is None:
- 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)
-
-def kron(a,b):
- """Kronecker product of a and b.
-
- The result is the block matrix::
-
- a[0,0]*b a[0,1]*b ... a[0,-1]*b
- a[1,0]*b a[1,1]*b ... a[1,-1]*b
- ...
- a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b
-
- Parameters
- ----------
- a : array, shape (M, N)
- b : array, shape (P, Q)
-
- Returns
- -------
- A : array, shape (M*P, N*Q)
- Kronecker product of a and b
-
- Examples
- --------
- >>> from scipy import kron, array
- >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
- array([[1, 1, 1, 2, 2, 2],
- [3, 3, 3, 4, 4, 4]])
-
- """
- if not a.flags['CONTIGUOUS']:
- a = 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)
-
-def block_diag(*arrs):
- """Create a block diagonal matrix from the provided arrays.
-
- Given the inputs `A`, `B` and `C`, the output will have these
- arrays arranged on the diagonal::
-
- [[A, 0, 0],
- [0, B, 0],
- [0, 0, C]]
-
- If all the input arrays are square, the output is known as a
- block diagonal matrix.
-
- Parameters
- ----------
- A, B, C, ... : array-like, up to 2D
- Input arrays. A 1D array or array-like sequence with length n is
- treated as a 2D array with shape (1,n).
-
- Returns
- -------
- D : ndarray
- Array with `A`, `B`, `C`, ... on the diagonal. `D` has the
- same dtype as `A`.
-
- References
- ----------
- .. [1] Wikipedia, "Block matrix",
- http://en.wikipedia.org/wiki/Block_diagonal_matrix
-
- Examples
- --------
- >>> A = [[1, 0],
- ... [0, 1]]
- >>> B = [[3, 4, 5],
- ... [6, 7, 8]]
- >>> C = [[7]]
- >>> print(block_diag(A, B, C))
- [[1 0 0 0 0 0]
- [0 1 0 0 0 0]
- [0 0 3 4 5 0]
- [0 0 6 7 8 0]
- [0 0 0 0 0 7]]
- >>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
- array([[ 1., 0., 0., 0., 0.],
- [ 0., 2., 3., 0., 0.],
- [ 0., 0., 0., 4., 5.],
- [ 0., 0., 0., 6., 7.]])
-
- """
- if arrs == ():
- arrs = ([],)
- arrs = [atleast_2d(a) for a in arrs]
-
- bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2]
- if bad_args:
- raise ValueError("arguments in the following positions have dimension "
- "greater than 2: %s" % bad_args)
-
- shapes = numpy.array([a.shape for a in arrs])
- out = zeros(sum(shapes, axis=0), dtype=arrs[0].dtype)
-
- r, c = 0, 0
- for i, (rr, cc) in enumerate(shapes):
- out[r:r + rr, c:c + cc] = arrs[i]
- r += rr
- c += cc
- return out
Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py 2010-03-30 04:39:23 UTC (rev 6283)
+++ trunk/scipy/linalg/decomp.py 2010-03-30 04:39:32 UTC (rev 6284)
@@ -15,8 +15,9 @@
'schur','rsf2csf','lu_factor','cho_factor','cho_solve','orth',
'hessenberg']
-from basic import LinAlgError
-import basic
+from misc import LinAlgError
+import misc
+import special_matrices
from warnings import warn
from lapack import get_lapack_funcs, find_best_lapack_type
@@ -1196,9 +1197,9 @@
% -info)
if not econ or M<N:
- R = basic.triu(qr)
+ R = special_matrices.triu(qr)
else:
- R = basic.triu(qr[0:N,0:N])
+ R = special_matrices.triu(qr[0:N,0:N])
if mode=='r':
return R
@@ -1276,7 +1277,7 @@
'illegal value in %-th argument of internal geqrf'%(-info)
gemm, = get_blas_funcs(('gemm',),(qr,))
t = qr.dtype.char
- R = basic.triu(qr)
+ R = special_matrices.triu(qr)
Q = numpy.identity(M,dtype=t)
ident = numpy.identity(M,dtype=t)
zeros = numpy.zeros
@@ -1336,7 +1337,7 @@
'illegal value in %-th argument of internal geqrf'%(-info)
gemm, = get_blas_funcs(('gemm',),(rq,))
t = rq.dtype.char
- R = basic.triu(rq)
+ R = special_matrices.triu(rq)
Q = numpy.identity(M,dtype=t)
ident = numpy.identity(M,dtype=t)
zeros = numpy.zeros
@@ -1490,7 +1491,7 @@
if abs(T[m,m-1]) > eps*(abs(T[m-1,m-1]) + abs(T[m,m])):
k = slice(m-1,m+1)
mu = eigvals(T[k,k]) - T[m,m]
- r = basic.norm([mu[0], T[m,m-1]])
+ r = misc.norm([mu[0], T[m,m-1]])
c = mu[0] / r
s = T[m,m-1] / r
G = r_[arr([[conj(c),s]],dtype=t),arr([[-s,c]],dtype=t)]
Added: trunk/scipy/linalg/misc.py
===================================================================
--- trunk/scipy/linalg/misc.py (rev 0)
+++ trunk/scipy/linalg/misc.py 2010-03-30 04:39:32 UTC (rev 6284)
@@ -0,0 +1,10 @@
+import numpy as np
+from numpy.linalg import LinAlgError
+
+### Norm
+
+def norm(a, ord=None):
+ # Differs from numpy only in non-finite handling
+ return np.linalg.norm(np.asarray_chkfinite(a), ord=ord)
+norm.__doc__ = np.linalg.norm.__doc__
+
Added: trunk/scipy/linalg/special_matrices.py
===================================================================
--- trunk/scipy/linalg/special_matrices.py (rev 0)
+++ trunk/scipy/linalg/special_matrices.py 2010-03-30 04:39:32 UTC (rev 6284)
@@ -0,0 +1,374 @@
+import numpy as np
+
+#-----------------------------------------------------------------------------
+# matrix construction functions
+#-----------------------------------------------------------------------------
+
+def tri(N, M=None, k=0, dtype=None):
+ """Construct (N, M) matrix filled with ones at and below the k-th diagonal.
+
+ The matrix has A[i,j] == 1 for i <= j + k
+
+ Parameters
+ ----------
+ N : integer
+ M : integer
+ Size of the matrix. If M is None, M == N is assumed.
+ k : integer
+ Number of subdiagonal below which matrix is filled with ones.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+ dtype : dtype
+ Data type of the matrix.
+
+ Returns
+ -------
+ A : array, shape (N, M)
+
+ Examples
+ --------
+ >>> from scipy.linalg import tri
+ >>> tri(3, 5, 2, dtype=int)
+ array([[1, 1, 1, 0, 0],
+ [1, 1, 1, 1, 0],
+ [1, 1, 1, 1, 1]])
+ >>> tri(3, 5, -1, dtype=int)
+ array([[0, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0],
+ [1, 1, 0, 0, 0]])
+
+ """
+ if M is None: M = N
+ if type(M) == type('d'):
+ #pearu: any objections to remove this feature?
+ # As tri(N,'d') is equivalent to tri(N,dtype='d')
+ dtype = M
+ M = N
+ m = np.greater_equal(np.subtract.outer(np.arange(N), np.arange(M)),-k)
+ if dtype is None:
+ return m
+ else:
+ return m.astype(dtype)
+
+def tril(m, k=0):
+ """Construct a copy of a matrix with elements above the k-th diagonal zeroed.
+
+ Parameters
+ ----------
+ m : array
+ Matrix whose elements to return
+ k : integer
+ Diagonal above which to zero elements.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+ Returns
+ -------
+ A : array, shape m.shape, dtype m.dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import tril
+ >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+ array([[ 0, 0, 0],
+ [ 4, 0, 0],
+ [ 7, 8, 0],
+ [10, 11, 12]])
+
+ """
+ m = np.asarray(m)
+ out = tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype.char)*m
+ return out
+
+def triu(m, k=0):
+ """Construct a copy of a matrix with elements below the k-th diagonal zeroed.
+
+ Parameters
+ ----------
+ m : array
+ Matrix whose elements to return
+ k : integer
+ Diagonal below which to zero elements.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+ Returns
+ -------
+ A : array, shape m.shape, dtype m.dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import tril
+ >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+ array([[ 1, 2, 3],
+ [ 4, 5, 6],
+ [ 0, 8, 9],
+ [ 0, 0, 12]])
+
+ """
+ m = np.asarray(m)
+ out = (1-tri(m.shape[0], m.shape[1], k-1, m.dtype.char))*m
+ return out
+
+
+def toeplitz(c, r=None):
+ """Construct a Toeplitz matrix.
+
+ 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-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))
+ The Toeplitz matrix.
+ dtype is the same as `(c[0] + r[0]).dtype`.
+
+ Examples
+ --------
+ >>> from scipy.linalg import toeplitz
+ >>> toeplitz([1,2,3], [1,4,5,6])
+ 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.
+ """
+ c = np.asarray(c).ravel()
+ if r is None:
+ r = c.conjugate()
+ else:
+ r = np.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 = np.concatenate((r[-1:0:-1], c))
+ a, b = np.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.
+
+ 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 = np.asarray(c).ravel()
+ a, b = np.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, 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-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))
+ 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],
+ [3, 4, 7, 7, 8],
+ [4, 7, 7, 8, 9]])
+
+ See also
+ --------
+ toeplitz : Toeplitz matrix
+ circulant : circulant matrix
+
+ """
+ c = np.asarray(c).ravel()
+ if r is None:
+ r = np.zeros_like(c)
+ else:
+ r = np.asarray(r).ravel()
+ # Form a 1D array of values to be used in the matrix, containing `c`
+ # followed by r[1:].
+ vals = np.concatenate((c, r[1:]))
+ a, b = np.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(np.matrix,args)
+
+def kron(a,b):
+ """Kronecker product of a and b.
+
+ The result is the block matrix::
+
+ a[0,0]*b a[0,1]*b ... a[0,-1]*b
+ a[1,0]*b a[1,1]*b ... a[1,-1]*b
+ ...
+ a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ b : array, shape (P, Q)
+
+ Returns
+ -------
+ A : array, shape (M*P, N*Q)
+ Kronecker product of a and b
+
+ Examples
+ --------
+ >>> from scipy import kron, array
+ >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
+ array([[1, 1, 1, 2, 2, 2],
+ [3, 3, 3, 4, 4, 4]])
+
+ """
+ if not a.flags['CONTIGUOUS']:
+ a = 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)
+
+def block_diag(*arrs):
+ """Create a block diagonal matrix from the provided arrays.
+
+ Given the inputs `A`, `B` and `C`, the output will have these
+ arrays arranged on the diagonal::
+
+ [[A, 0, 0],
+ [0, B, 0],
+ [0, 0, C]]
+
+ If all the input arrays are square, the output is known as a
+ block diagonal matrix.
+
+ Parameters
+ ----------
+ A, B, C, ... : array-like, up to 2D
+ Input arrays. A 1D array or array-like sequence with length n is
+ treated as a 2D array with shape (1,n).
+
+ Returns
+ -------
+ D : ndarray
+ Array with `A`, `B`, `C`, ... on the diagonal. `D` has the
+ same dtype as `A`.
+
+ References
+ ----------
+ .. [1] Wikipedia, "Block matrix",
+ http://en.wikipedia.org/wiki/Block_diagonal_matrix
+
+ Examples
+ --------
+ >>> A = [[1, 0],
+ ... [0, 1]]
+ >>> B = [[3, 4, 5],
+ ... [6, 7, 8]]
+ >>> C = [[7]]
+ >>> print(block_diag(A, B, C))
+ [[1 0 0 0 0 0]
+ [0 1 0 0 0 0]
+ [0 0 3 4 5 0]
+ [0 0 6 7 8 0]
+ [0 0 0 0 0 7]]
+ >>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
+ array([[ 1., 0., 0., 0., 0.],
+ [ 0., 2., 3., 0., 0.],
+ [ 0., 0., 0., 4., 5.],
+ [ 0., 0., 0., 6., 7.]])
+
+ """
+ if arrs == ():
+ arrs = ([],)
+ arrs = [np.atleast_2d(a) for a in arrs]
+
+ bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2]
+ if bad_args:
+ raise ValueError("arguments in the following positions have dimension "
+ "greater than 2: %s" % bad_args)
+
+ shapes = np.array([a.shape for a in arrs])
+ out = np.zeros(np.sum(shapes, axis=0), dtype=arrs[0].dtype)
+
+ r, c = 0, 0
+ for i, (rr, cc) in enumerate(shapes):
+ out[r:r + rr, c:c + cc] = arrs[i]
+ r += rr
+ c += cc
+ return out
More information about the Scipy-svn
mailing list