Merge branch 'develop' into documentation

This commit is contained in:
fjosw 2021-11-09 10:28:34 +00:00
commit 06a99c8be7
4 changed files with 11 additions and 72 deletions

View file

@ -1,6 +1,3 @@
#!/usr/bin/env python
# coding: utf-8
import numpy as np import numpy as np
from autograd import jacobian from autograd import jacobian
import autograd.numpy as anp # Thinly-wrapped numpy import autograd.numpy as anp # Thinly-wrapped numpy

View file

@ -1,6 +1,3 @@
#!/usr/bin/env python
# coding: utf-8
import numpy as np import numpy as np
from .obs import Obs from .obs import Obs

View file

@ -1,24 +1,24 @@
#!/usr/bin/env python
# coding: utf-8
import numpy as np import numpy as np
import scipy.linalg import scipy.linalg
from .obs import Obs from .obs import Obs
from .linalg import svd, eig, pinv from .linalg import svd, eig
def matrix_pencil_method(corrs, k=1, p=None, **kwargs): def matrix_pencil_method(corrs, k=1, p=None, **kwargs):
""" Matrix pencil method to extract k energy levels from data """Matrix pencil method to extract k energy levels from data
Implementation of the matrix pencil method based on Implementation of the matrix pencil method based on
eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990) eq. (2.17) of Y. Hua, T. K. Sarkar, IEEE Trans. Acoust. 38, 814-824 (1990)
Parameters Parameters
---------- ----------
data -- can be a list of Obs for the analysis of a single correlator, or a list of lists 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. of Obs if several correlators are to analyzed at once.
k -- Number of states to extract (default 1). k : int
p -- matrix pencil parameter which filters noise. The optimal value is expected between 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 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). to len(data)/2 but could possibly suppress more noise (default len(data)//2).
""" """
@ -56,55 +56,3 @@ def matrix_pencil_method(corrs, k=1, p=None, **kwargs):
# Return the sorted logarithms of the real eigenvalues as Obs # Return the sorted logarithms of the real eigenvalues as Obs
energy_levels = np.log(np.abs(eig(z, **kwargs))) energy_levels = np.log(np.abs(eig(z, **kwargs)))
return sorted(energy_levels, key=lambda x: abs(x.value)) return sorted(energy_levels, key=lambda x: abs(x.value))
def matrix_pencil_method_old(data, p, noise_level=None, verbose=1, **kwargs):
""" Older impleentation of the matrix pencil method with pencil p on given data to
extract energy levels.
Parameters
----------
data -- lists of Obs, where the nth entry is considered to be the correlation function
at x0=n+offset.
p -- matrix pencil parameter which corresponds to the number of energy levels to extract.
higher values for p can help decreasing noise.
noise_level -- If this argument is not None an additional prefiltering via singular
value decomposition is performed in which all singular values below 10^(-noise_level)
times the largest singular value are discarded. This increases the computation time.
verbose -- if larger than zero details about the noise filtering are printed to stdout
(default 1)
"""
n_data = len(data)
if n_data <= p:
raise Exception('The pencil p has to be smaller than the number of data samples.')
matrix = scipy.linalg.hankel(data[:n_data - p], data[n_data - p - 1:]) @ np.identity(p + 1)
if noise_level is not None:
u, s, vh = svd(matrix)
s_values = np.vectorize(lambda x: x.value)(s)
if verbose > 0:
print('Singular values: ', s_values)
digit = np.argwhere(s_values / s_values[0] < 10.0**(-noise_level))
if digit.size == 0:
digit = len(s_values)
else:
digit = int(digit[0])
if verbose > 0:
print('Consider only', digit, 'out of', len(s), 'singular values')
new_matrix = u[:, :digit] * s[:digit] @ vh[:digit, :]
y1 = new_matrix[:, :-1]
y2 = new_matrix[:, 1:]
else:
y1 = matrix[:, :-1]
y2 = matrix[:, 1:]
# MoorePenrose pseudoinverse
pinv_y1 = pinv(y1)
e = eig((pinv_y1 @ y2), **kwargs)
energy_levels = -np.log(np.abs(e))
return sorted(energy_levels, key=lambda x: abs(x.value))

View file

@ -1,6 +1,3 @@
#!/usr/bin/env python
# coding: utf-8
import warnings import warnings
import pickle import pickle
import numpy as np import numpy as np