pyerrors.mpm
1import numpy as np 2import scipy.linalg 3from .obs import Obs 4from .linalg import svd, eig 5 6 7def matrix_pencil_method(corrs, k=1, p=None, **kwargs): 8 """Matrix pencil method to extract k energy levels from data 9 10 Implementation of the matrix pencil method based on 11 eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990) 12 13 Parameters 14 ---------- 15 data : list 16 can be a list of Obs for the analysis of a single correlator, or a list of lists 17 of Obs if several correlators are to analyzed at once. 18 k : int 19 Number of states to extract (default 1). 20 p : int 21 matrix pencil parameter which filters noise. The optimal value is expected between 22 len(data)/3 and 2*len(data)/3. The computation is more expensive the closer p is 23 to len(data)/2 but could possibly suppress more noise (default len(data)//2). 24 25 Returns 26 ------- 27 energy_levels : list[Obs] 28 Extracted energy levels 29 """ 30 if isinstance(corrs[0], Obs): 31 data = [corrs] 32 else: 33 data = corrs 34 35 lengths = [len(d) for d in data] 36 if lengths.count(lengths[0]) != len(lengths): 37 raise Exception('All datasets have to have the same length.') 38 39 data_sets = len(data) 40 n_data = len(data[0]) 41 42 if p is None: 43 p = max(n_data // 2, k) 44 if n_data <= p: 45 raise Exception('The pencil p has to be smaller than the number of data samples.') 46 if p < k or n_data - p < k: 47 raise Exception('Cannot extract', k, 'energy levels with p=', p, 'and N-p=', n_data - p) 48 49 # Construct the hankel matrices 50 matrix = [] 51 for n in range(data_sets): 52 matrix.append(scipy.linalg.hankel(data[n][:n_data - p], data[n][n_data - p - 1:])) 53 matrix = np.array(matrix) 54 # Construct y1 and y2 55 y1 = np.concatenate(matrix[:, :, :p]) 56 y2 = np.concatenate(matrix[:, :, 1:]) 57 # Apply SVD to y2 58 u, s, vh = svd(y2, **kwargs) 59 # Construct z from y1 and SVD of y2, setting all singular values beyond the kth to zero 60 z = np.diag(1. / s[:k]) @ u[:, :k].T @ y1 @ vh.T[:, :k] 61 # Return the sorted logarithms of the real eigenvalues as Obs 62 energy_levels = np.log(np.abs(eig(z, **kwargs))) 63 return sorted(energy_levels, key=lambda x: abs(x.value))
def
matrix_pencil_method(corrs, k=1, p=None, **kwargs):
8def matrix_pencil_method(corrs, k=1, p=None, **kwargs): 9 """Matrix pencil method to extract k energy levels from data 10 11 Implementation of the matrix pencil method based on 12 eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990) 13 14 Parameters 15 ---------- 16 data : list 17 can be a list of Obs for the analysis of a single correlator, or a list of lists 18 of Obs if several correlators are to analyzed at once. 19 k : int 20 Number of states to extract (default 1). 21 p : int 22 matrix pencil parameter which filters noise. The optimal value is expected between 23 len(data)/3 and 2*len(data)/3. The computation is more expensive the closer p is 24 to len(data)/2 but could possibly suppress more noise (default len(data)//2). 25 26 Returns 27 ------- 28 energy_levels : list[Obs] 29 Extracted energy levels 30 """ 31 if isinstance(corrs[0], Obs): 32 data = [corrs] 33 else: 34 data = corrs 35 36 lengths = [len(d) for d in data] 37 if lengths.count(lengths[0]) != len(lengths): 38 raise Exception('All datasets have to have the same length.') 39 40 data_sets = len(data) 41 n_data = len(data[0]) 42 43 if p is None: 44 p = max(n_data // 2, k) 45 if n_data <= p: 46 raise Exception('The pencil p has to be smaller than the number of data samples.') 47 if p < k or n_data - p < k: 48 raise Exception('Cannot extract', k, 'energy levels with p=', p, 'and N-p=', n_data - p) 49 50 # Construct the hankel matrices 51 matrix = [] 52 for n in range(data_sets): 53 matrix.append(scipy.linalg.hankel(data[n][:n_data - p], data[n][n_data - p - 1:])) 54 matrix = np.array(matrix) 55 # Construct y1 and y2 56 y1 = np.concatenate(matrix[:, :, :p]) 57 y2 = np.concatenate(matrix[:, :, 1:]) 58 # Apply SVD to y2 59 u, s, vh = svd(y2, **kwargs) 60 # Construct z from y1 and SVD of y2, setting all singular values beyond the kth to zero 61 z = np.diag(1. / s[:k]) @ u[:, :k].T @ y1 @ vh.T[:, :k] 62 # Return the sorted logarithms of the real eigenvalues as Obs 63 energy_levels = np.log(np.abs(eig(z, **kwargs))) 64 return sorted(energy_levels, key=lambda x: abs(x.value))
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