mirror of
https://github.com/fjosw/pyerrors.git
synced 2026-08-04 11:31:21 +02:00
64 lines
2.2 KiB
Python
64 lines
2.2 KiB
Python
import numpy as np
|
|
import scipy.linalg
|
|
|
|
from .linalg import eig, svd
|
|
from .obs import Obs
|
|
|
|
|
|
def matrix_pencil_method(corrs, k=1, p=None, **kwargs):
|
|
"""Matrix pencil method to extract k energy levels from data
|
|
|
|
Implementation of the matrix pencil method based on
|
|
eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990)
|
|
|
|
Parameters
|
|
----------
|
|
data : list
|
|
can be a list of Obs for the analysis of a single correlator, or a list of lists
|
|
of Obs if several correlators are to analyzed at once.
|
|
k : int
|
|
Number of states to extract (default 1).
|
|
p : int
|
|
matrix pencil parameter which filters noise. The optimal value is expected between
|
|
len(data)/3 and 2*len(data)/3. The computation is more expensive the closer p is
|
|
to len(data)/2 but could possibly suppress more noise (default len(data)//2).
|
|
|
|
Returns
|
|
-------
|
|
energy_levels : list[Obs]
|
|
Extracted energy levels
|
|
"""
|
|
if isinstance(corrs[0], Obs):
|
|
data = [corrs]
|
|
else:
|
|
data = corrs
|
|
|
|
lengths = [len(d) for d in data]
|
|
if lengths.count(lengths[0]) != len(lengths):
|
|
raise ValueError('All datasets have to have the same length.')
|
|
|
|
data_sets = len(data)
|
|
n_data = len(data[0])
|
|
|
|
if p is None:
|
|
p = max(n_data // 2, k)
|
|
if n_data <= p:
|
|
raise ValueError('The pencil p has to be smaller than the number of data samples.')
|
|
if p < k or n_data - p < k:
|
|
raise ValueError(f'Cannot extract {k} energy levels with p={p} and N-p={n_data - p}')
|
|
|
|
# Construct the hankel matrices
|
|
matrix = []
|
|
for n in range(data_sets):
|
|
matrix.append(scipy.linalg.hankel(data[n][:n_data - p], data[n][n_data - p - 1:]))
|
|
matrix = np.array(matrix)
|
|
# Construct y1 and y2
|
|
y1 = np.concatenate(matrix[:, :, :p])
|
|
y2 = np.concatenate(matrix[:, :, 1:])
|
|
# Apply SVD to y2
|
|
u, s, vh = svd(y2, **kwargs)
|
|
# Construct z from y1 and SVD of y2, setting all singular values beyond the kth to zero
|
|
z = np.diag(1. / s[:k]) @ u[:, :k].T @ y1 @ vh.T[:, :k]
|
|
# Return the sorted logarithms of the real eigenvalues as Obs
|
|
energy_levels = np.log(np.abs(eig(z, **kwargs)))
|
|
return sorted(energy_levels, key=lambda x: abs(x.value))
|