# Source code for spectrum.linalg

"""
.. topic:: Linear algebra tools

.. autosummary::

csvd
corrmtx
pascal

.. codeauthor:: Thomas Cokelaer 2011

"""
import numpy

import spectrum.errors as errors

__all__ = ["pascal", "csvd", "corrmtx"]

[docs]def pascal(n):
"""Return Pascal matrix

:param int n: size of the matrix

.. doctest::

>>> from spectrum import pascal
>>> pascal(6)
array([[   1.,    1.,    1.,    1.,    1.,    1.],
[   1.,    2.,    3.,    4.,    5.,    6.],
[   1.,    3.,    6.,   10.,   15.,   21.],
[   1.,    4.,   10.,   20.,   35.,   56.],
[   1.,    5.,   15.,   35.,   70.,  126.],
[   1.,    6.,   21.,   56.,  126.,  252.]])

.. todo:: use the symmetric property to improve computational time if needed
"""
errors.is_positive_integer(n)
result = numpy.zeros((n, n))

#fill the first row and column
for i in range(0, n):
result[i, 0] = 1
result[0, i] = 1
if n > 1:
for i in range(1, n):
for j in range(1, n):
result[i, j] = result[i-1, j] + result[i, j-1]
return result

[docs]def corrmtx(x_input, m, method='autocorrelation'):
r"""Correlation matrix

This function is used by PSD estimator functions. It generates
the correlation matrix from a correlation data set and a maximum lag.

:param array x: autocorrelation samples (1D)
:param int m: the maximum lag

Depending on the choice of the method, the correlation matrix has different
sizes, but the number of rows is always m+1.

Method can be :

* 'autocorrelation': (default) X is the (n+m)-by-(m+1) rectangular Toeplitz
matrix derived using prewindowed and postwindowed data.
* 'prewindowed': X is the n-by-(m+1) rectangular Toeplitz matrix derived
using prewindowed data only.
* 'postwindowed': X is the n-by-(m+1) rectangular Toeplitz matrix that
derived using postwindowed data only.
* 'covariance': X is the (n-m)-by-(m+1) rectangular Toeplitz matrix
derived using nonwindowed data.
* 'modified': X is the 2(n-m)-by-(m+1) modified rectangular Toeplitz
matrix that generates an autocorrelation estimate for the length n data
vector x, derived using forward and backward prediction error estimates.

:return:
* the autocorrelation matrix
* R, the (m+1)-by-(m+1) autocorrelation matrix estimate R= X'*X.

.. rubric:: Algorithm details:

The **autocorrelation** matrix is a :math:(N+p) \times (p+1) rectangular Toeplitz
data matrix:

.. math:: X_p = \begin{pmatrix}L_p\\T_p\\Up\end{pmatrix}

where the lower triangular :math:p \times (p+1) matrix :math:L_p is

.. math:: L_p =
\begin{pmatrix}
x[1]     &  \cdots     & 0     & 0        \\
\vdots    &  \ddots     & \vdots & \vdots    \\
x[p]     &  \cdots     & x[1]  & 0
\end{pmatrix}

where the rectangular :math:(N-p) \times (p+1) matrix :math:T_p is

.. math:: T_p =
\begin{pmatrix}
x[p+1]     &  \cdots    & x[1]        \\
x[p+2]     &  \cdots    & x[2]        \\
\vdots    &  \ddots     & \vdots      \\
x[N-1]     &  \cdots    & x[N-p-1]    \\
x[N]     &  \cdots      & x[N-p]
\end{pmatrix}

and where the upper triangular :math:p \times (p+1) matrix :math:U_p is

.. math:: U_p =
\begin{pmatrix}
0         &  x[N]      & \cdots     & x[N-p+1]        \\
\vdots    &  \vdots    & \ddots & \vdots    \\
0         &  0         & \cdots  & x[N]
\end{pmatrix}

From this definition, the prewindowed matrix is

.. math:: X_p = \begin{pmatrix}L_p\\T_p\end{pmatrix}

the postwindowed matrix is

.. math:: X_p = \begin{pmatrix}T_p\\U_p\end{pmatrix}

the covariance matrix is:

.. math:: X_p = \begin{pmatrix}T_p\end{pmatrix}

and the modified covariance matrix is:

.. math:: X_p = \begin{pmatrix}T_p\\T_p^*\end{pmatrix}

"""

valid_methods = ['autocorrelation', 'prewindowed', 'postwindowed',
'covariance', 'modified']
if method not in valid_methods:
raise ValueError("Invalid method. Try one of %s" % valid_methods)

from scipy.linalg import toeplitz

# create the relevant matrices that will be useful to create
# the correlation matrices
N = len(x_input)

# FIXME:do we need a copy ?
if isinstance(x_input, list):
x = numpy.array(x_input)
else:
x = x_input.copy()

if x.dtype == complex:
complex_type = True
else:
complex_type = False

# Compute the Lp, Up and Tp matrices according to the requested method
if method in ['autocorrelation', 'prewindowed']:
Lp = toeplitz(x[0:m], [0]*(m+1))
Tp = toeplitz(x[m:N], x[m::-1])
if method in ['autocorrelation', 'postwindowed']:
Up = toeplitz([0]*(m+1), numpy.insert(x[N:N-m-1:-1],0,0))

# Create the output matrix

if method == 'autocorrelation':
if complex_type == True:
C = numpy.zeros((N+m, m+1), dtype=complex)
else:
C = numpy.zeros((N+m, m+1))
for i in range(0, m):
C[i] = Lp[i]
for i in range(m, N):
C[i] = Tp[i-m]
for i in range(N, N+m):
C[i] = Up[i-N]
elif method == 'prewindowed':
if complex_type == True:
C = numpy.zeros((N, m+1), dtype=complex)
else:
C = numpy.zeros((N, m+1))

for i in range(0, m):
C[i] = Lp[i]
for i in range(m, N):
C[i] = Tp[i-m]
elif method == 'postwindowed':
if complex_type == True:
C = numpy.zeros((N, m+1), dtype=complex)
else:
C = numpy.zeros((N, m+1))
for i in range(0, N-m):
C[i] = Tp[i]
for i in range(N-m, N):
C[i] = Up[i-N+m]
elif method == 'covariance':
return Tp
elif method == 'modified':
if complex_type == True:
C = numpy.zeros((2*(N-m), m+1), dtype=complex)
else:
C = numpy.zeros((2*(N-m), m+1))
for i in range(0, N-m):
C[i] = Tp[i]
Tp = numpy.fliplr(Tp.conj())
for i in range(N-m, 2*(N-m)):
C[i] = Tp[i-N+m]

return C
#autocorrelation n+m by m+1 matrix (Lp, Tp, Up)

#Rx = toeplitz(r)

[docs]def csvd(A):
"""SVD decomposition using numpy.linalg.svd

:param A: a M by N matrix

:return:
* U, a M by M matrix
* S the N eigen values
* V a N by N matrix

See :func:numpy.linalg.svd for a detailed documentation.

Should return the same as in [Marple]_ , CSVD routine.

::

U, S, V = numpy.linalg.svd(A)
U, S, V = cvsd(A)

"""
U, S, V = numpy.linalg.svd(A)
return U, S, V