[Scipy-svn] r6284 - trunk/scipy/linalg

scipy-svn@scip... scipy-svn@scip...
Mon Mar 29 23:39:33 CDT 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