mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-05-14 19:43:41 +02:00
commit
090a26919b
2 changed files with 126 additions and 26 deletions
|
@ -1,7 +1,9 @@
|
||||||
import numpy as np
|
import numpy as np
|
||||||
import autograd.numpy as anp
|
import autograd.numpy as anp
|
||||||
|
#from scipy.special.orthogonal import _IntegerType
|
||||||
from .pyerrors import *
|
from .pyerrors import *
|
||||||
from .fits import standard_fit
|
from .fits import standard_fit
|
||||||
|
from .linalg import *
|
||||||
from .roots import find_root
|
from .roots import find_root
|
||||||
from matplotlib import pyplot as plt
|
from matplotlib import pyplot as plt
|
||||||
from matplotlib.ticker import NullFormatter # useful for `logit` scale
|
from matplotlib.ticker import NullFormatter # useful for `logit` scale
|
||||||
|
@ -22,7 +24,7 @@ class Corr:
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
|
||||||
def __init__(self, data_input, padding_front=0, padding_back=0):
|
def __init__(self, data_input, padding_front=0, padding_back=0,prange=None):
|
||||||
#All data_input should be a list of things at different timeslices. This needs to be verified
|
#All data_input should be a list of things at different timeslices. This needs to be verified
|
||||||
|
|
||||||
if not (isinstance(data_input, list)):
|
if not (isinstance(data_input, list)):
|
||||||
|
@ -56,6 +58,12 @@ class Corr:
|
||||||
self.T = len(self.content) #for convenience: will be used a lot
|
self.T = len(self.content) #for convenience: will be used a lot
|
||||||
|
|
||||||
|
|
||||||
|
#The attribute "range" [start,end] marks a range of two timeslices.
|
||||||
|
#This is useful for keeping track of plateaus and fitranges.
|
||||||
|
#The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant.
|
||||||
|
if not range is None:
|
||||||
|
self.prange=prange
|
||||||
|
|
||||||
self.gamma_method()
|
self.gamma_method()
|
||||||
|
|
||||||
|
|
||||||
|
@ -98,6 +106,8 @@ class Corr:
|
||||||
newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content]
|
newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content]
|
||||||
return Corr(newcontent)
|
return Corr(newcontent)
|
||||||
|
|
||||||
|
def sum(self):
|
||||||
|
return np.sqrt(self.N)*self.projected(np.ones(self.N))
|
||||||
#For purposes of debugging and verification, one might want to see a single smearing level. smearing will return a Corr at the specified i,j. where both are integers 0<=i,j<N.
|
#For purposes of debugging and verification, one might want to see a single smearing level. smearing will return a Corr at the specified i,j. where both are integers 0<=i,j<N.
|
||||||
def smearing(self, i, j):
|
def smearing(self, i, j):
|
||||||
if self.N == 1:
|
if self.N == 1:
|
||||||
|
@ -133,7 +143,8 @@ class Corr:
|
||||||
newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
|
newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
|
||||||
if(all([x is None for x in newcontent])):
|
if(all([x is None for x in newcontent])):
|
||||||
raise Exception("Corr could not be symmetrized: No redundant values")
|
raise Exception("Corr could not be symmetrized: No redundant values")
|
||||||
return Corr(newcontent)
|
|
||||||
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
|
|
||||||
def anti_symmetric(self):
|
def anti_symmetric(self):
|
||||||
|
|
||||||
|
@ -148,7 +159,7 @@ class Corr:
|
||||||
newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
|
newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
|
||||||
if(all([x is None for x in newcontent])):
|
if(all([x is None for x in newcontent])):
|
||||||
raise Exception("Corr could not be symmetrized: No redundant values")
|
raise Exception("Corr could not be symmetrized: No redundant values")
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
|
|
||||||
|
|
||||||
#This method will symmetrice the matrices and therefore make them positive definit.
|
#This method will symmetrice the matrices and therefore make them positive definit.
|
||||||
|
@ -161,7 +172,7 @@ class Corr:
|
||||||
|
|
||||||
|
|
||||||
#We also include a simple GEVP method based on Scipy.linalg
|
#We also include a simple GEVP method based on Scipy.linalg
|
||||||
def GEVP(self, t0, ts):
|
def GEVP(self, t0, ts,state=1):
|
||||||
if (self.content[t0] is None) or (self.content[ts] is None):
|
if (self.content[t0] is None) or (self.content[ts] is None):
|
||||||
raise Exception("Corr not defined at t0/ts")
|
raise Exception("Corr not defined at t0/ts")
|
||||||
G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
|
G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
|
||||||
|
@ -171,11 +182,37 @@ class Corr:
|
||||||
Gt[i, j] = self.content[ts][i, j].value
|
Gt[i, j] = self.content[ts][i, j].value
|
||||||
|
|
||||||
sp_val,sp_vec = scipy.linalg.eig(Gt, G0)
|
sp_val,sp_vec = scipy.linalg.eig(Gt, G0)
|
||||||
sp_vec = sp_vec[:,np.argmax(sp_val)] #we only want the eigenvector belonging to the biggest eigenvalue.
|
sp_vec=sp_vec[:,np.argsort(sp_val)[-state]] #we only want the eigenvector belonging to the selected state
|
||||||
sp_vec = sp_vec/np.sqrt(sp_vec@sp_vec)
|
sp_vec = sp_vec/np.sqrt(sp_vec@sp_vec)
|
||||||
return sp_vec
|
return sp_vec
|
||||||
|
|
||||||
|
|
||||||
|
def Eigenvalue(self,t0,state=1):
|
||||||
|
G=self.smearing_symmetric()
|
||||||
|
G0=G.content[t0]
|
||||||
|
L = mat_mat_op(anp.linalg.cholesky, G0)
|
||||||
|
Li = mat_mat_op(anp.linalg.inv, L)
|
||||||
|
LT=L.T
|
||||||
|
LTi=mat_mat_op(anp.linalg.inv, LT)
|
||||||
|
newcontent=[]
|
||||||
|
for t in range(self.T):
|
||||||
|
Gt=G.content[t]
|
||||||
|
M=Li@Gt@LTi
|
||||||
|
eigenvalues = eigh(M)[0]
|
||||||
|
#print(eigenvalues)
|
||||||
|
eigenvalue=eigenvalues[-state]
|
||||||
|
newcontent.append(eigenvalue)
|
||||||
|
return Corr(newcontent)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def roll(self, dt):
|
def roll(self, dt):
|
||||||
return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
|
return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
|
||||||
|
|
||||||
|
@ -250,6 +287,17 @@ class Corr:
|
||||||
raise Exception('m_eff is undefined at all timeslices')
|
raise Exception('m_eff is undefined at all timeslices')
|
||||||
|
|
||||||
return Corr(newcontent, padding_back=1)
|
return Corr(newcontent, padding_back=1)
|
||||||
|
elif variant is 'arccosh':
|
||||||
|
newcontent=[]
|
||||||
|
for t in range(1,self.T-1):
|
||||||
|
if (self.content[t] is None) or (self.content[t+1] is None)or (self.content[t-1] is None):
|
||||||
|
newcontent.append(None)
|
||||||
|
else:
|
||||||
|
newcontent.append((self.content[t+1]+self.content[t-1])/(2*self.content[t]))
|
||||||
|
if(all([x is None for x in newcontent])):
|
||||||
|
raise Exception("m_eff is undefined at all timeslices")
|
||||||
|
return np.arccosh(Corr(newcontent,padding_back=1,padding_front=1))
|
||||||
|
|
||||||
else:
|
else:
|
||||||
raise Exception('Unkown variant.')
|
raise Exception('Unkown variant.')
|
||||||
|
|
||||||
|
@ -258,8 +306,17 @@ class Corr:
|
||||||
if self.N != 1:
|
if self.N != 1:
|
||||||
raise Exception("Correlator must be projected before fitting")
|
raise Exception("Correlator must be projected before fitting")
|
||||||
|
|
||||||
|
#The default behaviour is:
|
||||||
|
#1 use explicit fitrange
|
||||||
|
# if none is provided, use the range of the corr
|
||||||
|
# if this is also not set, use the whole length of the corr (This could come with a warning!)
|
||||||
|
|
||||||
|
|
||||||
if fitrange is None:
|
if fitrange is None:
|
||||||
fitrange=[0, self.T]
|
if hasattr(self,"prange"):
|
||||||
|
fitrange=self.prange
|
||||||
|
else:
|
||||||
|
fitrange=[0, self.T]
|
||||||
|
|
||||||
xs = [x for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
|
xs = [x for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
|
||||||
ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
|
ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None]
|
||||||
|
@ -273,7 +330,12 @@ class Corr:
|
||||||
return result
|
return result
|
||||||
|
|
||||||
#we want to quickly get a plateau
|
#we want to quickly get a plateau
|
||||||
def plateau(self, plateau_range, method="fit"):
|
def plateau(self, plateau_range=None, method="fit"):
|
||||||
|
if not plateau_range:
|
||||||
|
if hasattr(self,"prange"):
|
||||||
|
plateau_range=self.prange
|
||||||
|
else:
|
||||||
|
raise Exception("no plateau range provided")
|
||||||
if self.N != 1:
|
if self.N != 1:
|
||||||
raise Exception("Correlator must be projected before getting a plateau.")
|
raise Exception("Correlator must be projected before getting a plateau.")
|
||||||
if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1])])):
|
if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1])])):
|
||||||
|
@ -288,7 +350,26 @@ class Corr:
|
||||||
return returnvalue
|
return returnvalue
|
||||||
|
|
||||||
else:
|
else:
|
||||||
raise Exception("Unsupported plateau method: " + method)
|
raise Exception("Unsupported plateau method: " + method)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def set_prange(self,prange):
|
||||||
|
if not len(prange)==2:
|
||||||
|
raise Exception("range must be a list or array with two values")
|
||||||
|
if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))):
|
||||||
|
raise Exception("start and end point must be integers")
|
||||||
|
if not (0<=prange[0]<=self.T and 0<=prange[1]<=self.T and prange[0]<prange[1] ):
|
||||||
|
raise Exception("start and end point must define a range in the interval 0,T")
|
||||||
|
|
||||||
|
self.prange=prange
|
||||||
|
return
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#quick and dirty plotting function to view Correlator inside Jupyter
|
#quick and dirty plotting function to view Correlator inside Jupyter
|
||||||
#If one would not want to import pyplot, this could easily be replaced by a call to pe.plot_corrs
|
#If one would not want to import pyplot, this could easily be replaced by a call to pe.plot_corrs
|
||||||
|
@ -337,6 +418,9 @@ class Corr:
|
||||||
ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
|
ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
|
||||||
else:
|
else:
|
||||||
raise Exception('plateau must be an Obs')
|
raise Exception('plateau must be an Obs')
|
||||||
|
if hasattr(self,"prange"):
|
||||||
|
ax1.axvline(self.prange[0],0,1)
|
||||||
|
ax1.axvline(self.prange[1],0,1)
|
||||||
|
|
||||||
if fit_res:
|
if fit_res:
|
||||||
x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
|
x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
|
||||||
|
@ -369,6 +453,10 @@ class Corr:
|
||||||
def print(self, range=[0, None]):
|
def print(self, range=[0, None]):
|
||||||
print(self.__repr__(range))
|
print(self.__repr__(range))
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def __repr__(self, range=[0, None]):
|
def __repr__(self, range=[0, None]):
|
||||||
if range[1]:
|
if range[1]:
|
||||||
range[1] += 1
|
range[1] += 1
|
||||||
|
@ -410,7 +498,7 @@ class Corr:
|
||||||
newcontent.append(None)
|
newcontent.append(None)
|
||||||
else:
|
else:
|
||||||
newcontent.append(self.content[t]+y)
|
newcontent.append(self.content[t]+y)
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
else:
|
else:
|
||||||
raise TypeError("Corr + wrong type")
|
raise TypeError("Corr + wrong type")
|
||||||
|
|
||||||
|
@ -433,7 +521,7 @@ class Corr:
|
||||||
newcontent.append(None)
|
newcontent.append(None)
|
||||||
else:
|
else:
|
||||||
newcontent.append(self.content[t]*y)
|
newcontent.append(self.content[t]*y)
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
else:
|
else:
|
||||||
raise TypeError("Corr * wrong type")
|
raise TypeError("Corr * wrong type")
|
||||||
|
|
||||||
|
@ -471,7 +559,11 @@ class Corr:
|
||||||
newcontent.append(None)
|
newcontent.append(None)
|
||||||
else:
|
else:
|
||||||
newcontent.append(self.content[t]/y)
|
newcontent.append(self.content[t]/y)
|
||||||
return Corr(newcontent)
|
if hasattr(self,"prange"):
|
||||||
|
newrange=self.prange
|
||||||
|
else:
|
||||||
|
newrange=None
|
||||||
|
return Corr(newcontent,prange=newrange)
|
||||||
|
|
||||||
elif isinstance(y, int) or isinstance(y,float):
|
elif isinstance(y, int) or isinstance(y,float):
|
||||||
if y==0:
|
if y==0:
|
||||||
|
@ -482,13 +574,17 @@ class Corr:
|
||||||
newcontent.append(None)
|
newcontent.append(None)
|
||||||
else:
|
else:
|
||||||
newcontent.append(self.content[t]/y)
|
newcontent.append(self.content[t]/y)
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
else:
|
else:
|
||||||
raise TypeError("Corr / wrong type")
|
raise TypeError("Corr / wrong type")
|
||||||
|
|
||||||
def __neg__(self):
|
def __neg__(self):
|
||||||
newcontent=[None if (item is None) else -1.*item for item in self.content]
|
newcontent=[None if (item is None) else -1.*item for item in self.content]
|
||||||
return Corr(newcontent)
|
if hasattr(self,"prange"):
|
||||||
|
newrange=self.prange
|
||||||
|
else:
|
||||||
|
newrange=None
|
||||||
|
return Corr(newcontent,prange=newrange)
|
||||||
|
|
||||||
def __sub__(self,y):
|
def __sub__(self,y):
|
||||||
return self +(-y)
|
return self +(-y)
|
||||||
|
@ -496,13 +592,13 @@ class Corr:
|
||||||
def __pow__(self, y):
|
def __pow__(self, y):
|
||||||
if isinstance(y, Obs) or isinstance(y,int) or isinstance(y,float):
|
if isinstance(y, Obs) or isinstance(y,int) or isinstance(y,float):
|
||||||
newcontent=[None if (item is None) else item**y for item in self.content]
|
newcontent=[None if (item is None) else item**y for item in self.content]
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
else:
|
else:
|
||||||
raise TypeError("type of exponent not supported")
|
raise TypeError("type of exponent not supported")
|
||||||
|
|
||||||
def __abs__(self):
|
def __abs__(self):
|
||||||
newcontent=[None if (item is None) else np.abs(item) for item in self.content]
|
newcontent=[None if (item is None) else np.abs(item) for item in self.content]
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
|
|
||||||
#The numpy functions:
|
#The numpy functions:
|
||||||
def sqrt(self):
|
def sqrt(self):
|
||||||
|
@ -510,11 +606,11 @@ class Corr:
|
||||||
|
|
||||||
def log(self):
|
def log(self):
|
||||||
newcontent=[None if (item is None) else np.log(item) for item in self.content]
|
newcontent=[None if (item is None) else np.log(item) for item in self.content]
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
|
|
||||||
def exp(self):
|
def exp(self):
|
||||||
newcontent=[None if (item is None) else np.exp(item) for item in self.content]
|
newcontent=[None if (item is None) else np.exp(item) for item in self.content]
|
||||||
return Corr(newcontent)
|
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
|
||||||
|
|
||||||
def sin(self):
|
def sin(self):
|
||||||
newcontent=[None if (item is None) else np.sin(item) for item in self.content]
|
newcontent=[None if (item is None) else np.sin(item) for item in self.content]
|
||||||
|
|
|
@ -14,7 +14,7 @@ from autograd import elementwise_grad as egrad
|
||||||
from .pyerrors import Obs, derived_observable, covariance, pseudo_Obs
|
from .pyerrors import Obs, derived_observable, covariance, pseudo_Obs
|
||||||
|
|
||||||
|
|
||||||
def standard_fit(x, y, func, silent=False, **kwargs):
|
def standard_fit(x, y, func,n_parms="auto", silent=False, **kwargs):
|
||||||
"""Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
|
"""Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
|
||||||
|
|
||||||
x has to be a list of floats.
|
x has to be a list of floats.
|
||||||
|
@ -68,15 +68,19 @@ def standard_fit(x, y, func, silent=False, **kwargs):
|
||||||
if not callable(func):
|
if not callable(func):
|
||||||
raise TypeError('func has to be a function.')
|
raise TypeError('func has to be a function.')
|
||||||
|
|
||||||
for i in range(25):
|
if n_parms=="auto":
|
||||||
try:
|
for i in range(25):
|
||||||
func(np.arange(i), x.T[0])
|
try:
|
||||||
except:
|
func(np.arange(i), x.T[0])
|
||||||
pass
|
except:
|
||||||
else:
|
pass
|
||||||
break
|
else:
|
||||||
|
break
|
||||||
|
|
||||||
|
|
||||||
n_parms = i
|
n_parms = i
|
||||||
|
|
||||||
|
|
||||||
if not silent:
|
if not silent:
|
||||||
print('Fit with', n_parms, 'parameters')
|
print('Fit with', n_parms, 'parameters')
|
||||||
|
|
||||||
|
|
Loading…
Add table
Add a link
Reference in a new issue