lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-05-16 20:43:41 +02:00
Code Issues Projects Releases Packages Wiki Activity

correlator class cleaned up

Browse source
This commit is contained in:
Fabian Joswig 2021-10-12 10:13:58 +01:00
parent c5c673bdb6
commit 311b0bf91a
1 changed files with 189 additions and 264 deletions
Download patch file Download diff file Expand all files Collapse all files

437
pyerrors/correlators.py
View file

@ -1,15 +1,12 @@
import numpy as np import numpy as np
import autograd.numpy as anp import autograd.numpy as anp
#from scipy.special.orthogonal import _IntegerType import scipy.linalg
from .pyerrors import * from .pyerrors import Obs, dump_object
from .fits import standard_fit from .fits import standard_fit
from .linalg import * from .linalg import eigh, mat_mat_op
from .roots import find_root from .roots import find_root
from matplotlib import pyplot as plt import matplotlib.pyplot as plt
from matplotlib.ticker import NullFormatter # useful for `logit` scale
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
#import PySimpleGUI as sg
import matplotlib
class Corr: class Corr:
"""The class for a correlator (time dependent sequence of pe.Obs). """The class for a correlator (time dependent sequence of pe.Obs).
@ -25,22 +22,22 @@ class Corr:
""" """
def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): 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)):
raise TypeError('Corr__init__ expects a list of timeslices.') raise TypeError('Corr__init__ expects a list of timeslices.')
# data_input can have multiple shapes. The simplest one is a list of Obs. # data_input can have multiple shapes. The simplest one is a list of Obs.
#We check, if this is the case # We check, if this is the case
if all([isinstance(item, Obs) for item in data_input]): if all([isinstance(item, Obs) for item in data_input]):
self.content=[np.asarray([item]) for item in data_input] self.content = [np.asarray([item]) for item in data_input]
#Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices. # Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices.
self.N = 1 # number of smearings self.N = 1 # number of smearings
#data_input in the form [np.array(Obs,NxN)] # data_input in the form [np.array(Obs,NxN)]
elif all([isinstance(item,np.ndarray) or item==None for item in data_input]) and any([isinstance(item,np.ndarray)for item in data_input]): elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]):
self.content = data_input self.content = data_input
noNull=[a for a in self.content if not (a is None)] #To check if the matrices are correct for all undefined elements noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements
self.N = noNull[0].shape[0] self.N = noNull[0].shape[0]
# The checks are now identical to the case above # The checks are now identical to the case above
if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
@ -48,18 +45,18 @@ class Corr:
if (not all([item.shape == noNull[0].shape for item in noNull])): if (not all([item.shape == noNull[0].shape for item in noNull])):
raise Exception("Items in data_input are not of identical shape." + str(noNull)) raise Exception("Items in data_input are not of identical shape." + str(noNull))
else: # In case its a list of something else. else: # In case its a list of something else.
raise Exception ("data_input contains item of wrong type") raise Exception("data_input contains item of wrong type")
self.tag = None self.tag = None
#We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value. # We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value.
#An undefined timeslice is represented by the None object # An undefined timeslice is represented by the None object
self.content = [None] * padding_front + self.content + [None] * padding_back self.content = [None] * padding_front + self.content + [None] * padding_back
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. # The attribute "range" [start,end] marks a range of two timeslices.
#This is useful for keeping track of plateaus and fitranges. # 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. # The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant.
self.prange = prange self.prange = prange
self.gamma_method() self.gamma_method()
@ -72,64 +69,64 @@ class Corr:
else: else:
for i in range(self.N): for i in range(self.N):
for j in range(self.N): for j in range(self.N):
item[i,j].gamma_method() item[i, j].gamma_method()
#We need to project the Correlator with a Vector to get a single value at each timeslice. # We need to project the Correlator with a Vector to get a single value at each timeslice.
#The method can use one or two vectors. # The method can use one or two vectors.
#If two are specified it returns v1@G@v2 (the order might be very important.) # If two are specified it returns v1@G@v2 (the order might be very important.)
#By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to # By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
def projected(self, vector_l=None, vector_r=None): def projected(self, vector_l=None, vector_r=None):
if self.N == 1: if self.N == 1:
raise Exception("Trying to project a Corr, that already has N=1.") raise Exception("Trying to project a Corr, that already has N=1.")
#This Exception is in no way necessary. One could just return self # This Exception is in no way necessary. One could just return self
#But there is no scenario, where a user would want that to happen and the error message might be more informative. # But there is no scenario, where a user would want that to happen and the error message might be more informative.
self.gamma_method() self.gamma_method()
if vector_l is None: if vector_l is None:
vector_l,vector_r=np.asarray([1.] + (self.N - 1) * [0.]),np.asarray([1.] + (self.N - 1) * [0.]) vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.])
elif(vector_r is None): elif(vector_r is None):
vector_r=vector_l vector_r = vector_l
if not vector_l.shape == vector_r.shape == (self.N,): if not vector_l.shape == vector_r.shape == (self.N,):
raise Exception("Vectors are of wrong shape!") raise Exception("Vectors are of wrong shape!")
#We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized. # We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized.
if (not(0.95 < vector_r@vector_r < 1.05)) or (not(0.95 < vector_l@vector_l < 1.05)): if (not (0.95 < vector_r @ vector_r < 1.05)) or (not (0.95 < vector_l @ vector_l < 1.05)):
print("Vectors are normalized before projection!") print("Vectors are normalized before projection!")
vector_l,vector_r = vector_l / np.sqrt((vector_l@vector_l)), vector_r / np.sqrt(vector_r@vector_r) vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
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): def sum(self):
return np.sqrt(self.N)*self.projected(np.ones(self.N)) 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:
raise Exception("Trying to pick smearing from projected Corr") raise Exception("Trying to pick smearing from projected Corr")
newcontent=[None if(item is None) else item[i, j] for item in self.content] newcontent = [None if(item is None) else item[i, j] for item in self.content]
return Corr(newcontent) return Corr(newcontent)
# Obs and Matplotlib do not play nicely
#Obs and Matplotlib do not play nicely # We often want to retrieve x,y,y_err as lists to pass them to something like pyplot.errorbar
#We often want to retrieve x,y,y_err as lists to pass them to something like pyplot.errorbar
def plottable(self): def plottable(self):
if self.N != 1: if self.N != 1:
raise Exception("Can only make Corr[N=1] plottable") #We could also autoproject to the groundstate or expect vectors, but this is supposed to be a super simple function. raise Exception("Can only make Corr[N=1] plottable") # We could also autoproject to the groundstate or expect vectors, but this is supposed to be a super simple function.
x_list = [x for x in range(self.T) if (not self.content[x] is None)] x_list = [x for x in range(self.T) if (not self.content[x] is None)]
y_list = [y[0].value for y in self.content if (not y is None)] y_list = [y[0].value for y in self.content if (y is not None)]
y_err_list = [y[0].dvalue for y in self.content if (not y is None)] y_err_list = [y[0].dvalue for y in self.content if (y is not None)]
return x_list, y_list, y_err_list return x_list, y_list, y_err_list
#symmetric returns a Corr, that has been symmetrized. # symmetric returns a Corr, that has been symmetrized.
#A symmetry checker is still to be implemented # A symmetry checker is still to be implemented
#The method will not delete any redundant timeslices (Bad for memory, Great for convenience) # The method will not delete any redundant timeslices (Bad for memory, Great for convenience)
def symmetric(self): def symmetric(self):
if self.T%2 != 0: if self.T % 2 != 0:
raise Exception("Can not symmetrize odd T") raise Exception("Can not symmetrize odd T")
newcontent = [self.content[0]] newcontent = [self.content[0]]
@ -142,13 +139,12 @@ class Corr:
raise Exception("Corr could not be symmetrized: No redundant values") raise Exception("Corr could not be symmetrized: No redundant values")
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
def anti_symmetric(self): def anti_symmetric(self):
if self.T%2 != 0: if self.T % 2 != 0:
raise Exception("Can not symmetrize odd T") raise Exception("Can not symmetrize odd T")
newcontent=[self.content[0]] newcontent = [self.content[0]]
for t in range(1, self.T): for t in range(1, self.T):
if (self.content[t] is None) or (self.content[self.T - t] is None): if (self.content[t] is None) or (self.content[self.T - t] is None):
newcontent.append(None) newcontent.append(None)
@ -158,18 +154,16 @@ class Corr:
raise Exception("Corr could not be symmetrized: No redundant values") raise Exception("Corr could not be symmetrized: No redundant values")
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
# This method will symmetrice the matrices and therefore make them positive definit.
#This method will symmetrice the matrices and therefore make them positive definit.
def smearing_symmetric(self): def smearing_symmetric(self):
if self.N > 1: if self.N > 1:
transposed = [None if (G is None) else G.T for G in self.content] transposed = [None if (G is None) else G.T for G in self.content]
return 0.5 * (Corr(transposed)+self) return 0.5 * (Corr(transposed) + self)
if self.N == 1: if self.N == 1:
raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.") raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.")
# 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, state=1):
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")
@ -178,39 +172,35 @@ class Corr:
G0[i, j] = self.content[t0][i, j].value G0[i, j] = self.content[t0][i, j].value
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.argsort(sp_val)[-state]] #we only want the eigenvector belonging to the selected state 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):
def Eigenvalue(self,t0,state=1): G = self.smearing_symmetric()
G=self.smearing_symmetric() G0 = G.content[t0]
G0=G.content[t0]
L = mat_mat_op(anp.linalg.cholesky, G0) L = mat_mat_op(anp.linalg.cholesky, G0)
Li = mat_mat_op(anp.linalg.inv, L) Li = mat_mat_op(anp.linalg.inv, L)
LT=L.T LT = L.T
LTi=mat_mat_op(anp.linalg.inv, LT) LTi = mat_mat_op(anp.linalg.inv, LT)
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
Gt=G.content[t] Gt = G.content[t]
M=Li@Gt@LTi M = Li @ Gt @ LTi
eigenvalues = eigh(M)[0] eigenvalues = eigh(M)[0]
#print(eigenvalues) eigenvalue = eigenvalues[-state]
eigenvalue=eigenvalues[-state]
newcontent.append(eigenvalue) newcontent.append(eigenvalue)
return Corr(newcontent) 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)))
def deriv(self, symmetric=True): # Defaults to symmetric derivative
def deriv(self, symmetric=True): #Defaults to symmetric derivative
if not symmetric: if not symmetric:
newcontent = [] newcontent = []
for t in range(self.T - 1): for t in range(self.T - 1):
if (self.content[t] is None) or (self.content[t+1] is None): if (self.content[t] is None) or (self.content[t + 1] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t + 1] - self.content[t]) newcontent.append(self.content[t + 1] - self.content[t])
@ -219,8 +209,8 @@ class Corr:
return Corr(newcontent, padding_back=1) return Corr(newcontent, padding_back=1)
if symmetric: if symmetric:
newcontent = [] newcontent = []
for t in range(1, self.T-1): for t in range(1, self.T - 1):
if (self.content[t-1] is None) or (self.content[t+1] is None): if (self.content[t - 1] is None) or (self.content[t + 1] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
@ -228,11 +218,10 @@ class Corr:
raise Exception('Derivative is undefined at all timeslices') raise Exception('Derivative is undefined at all timeslices')
return Corr(newcontent, padding_back=1, padding_front=1) return Corr(newcontent, padding_back=1, padding_front=1)
def second_deriv(self): def second_deriv(self):
newcontent = [] newcontent = []
for t in range(1, self.T-1): for t in range(1, self.T - 1):
if (self.content[t-1] is None) or (self.content[t+1] is None): if (self.content[t - 1] is None) or (self.content[t + 1] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
@ -240,7 +229,6 @@ class Corr:
raise Exception("Derivative is undefined at all timeslices") raise Exception("Derivative is undefined at all timeslices")
return Corr(newcontent, padding_back=1, padding_front=1) return Corr(newcontent, padding_back=1, padding_front=1)
def m_eff(self, variant='log', guess=1.0): def m_eff(self, variant='log', guess=1.0):
"""Returns the effective mass of the correlator as correlator object """Returns the effective mass of the correlator as correlator object
@ -252,7 +240,7 @@ class Corr:
""" """
if self.N != 1: if self.N != 1:
raise Exception('Correlator must be projected before getting m_eff') raise Exception('Correlator must be projected before getting m_eff')
if variant is 'log': if variant == 'log':
newcontent = [] newcontent = []
for t in range(self.T - 1): for t in range(self.T - 1):
if (self.content[t] is None) or (self.content[t + 1] is None): if (self.content[t] is None) or (self.content[t + 1] is None):
@ -264,43 +252,43 @@ class Corr:
return np.log(Corr(newcontent, padding_back=1)) return np.log(Corr(newcontent, padding_back=1))
elif variant is 'periodic': elif variant == 'periodic':
newcontent = [] newcontent = []
for t in range(self.T - 1): for t in range(self.T - 1):
if (self.content[t] is None) or (self.content[t + 1] is None): if (self.content[t] is None) or (self.content[t + 1] is None):
newcontent.append(None) newcontent.append(None)
else: else:
func = lambda x, d : anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d func = lambda x, d: anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d
newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess))) newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess)))
if(all([x is None for x in newcontent])): if(all([x is None for x in newcontent])):
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=[] elif variant == 'arccosh':
for t in range(1,self.T-1): newcontent = []
if (self.content[t] is None) or (self.content[t+1] is None)or (self.content[t-1] is None): 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) newcontent.append(None)
else: else:
newcontent.append((self.content[t+1]+self.content[t-1])/(2*self.content[t])) newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
if(all([x is None for x in newcontent])): if(all([x is None for x in newcontent])):
raise Exception("m_eff is undefined at all timeslices") raise Exception("m_eff is undefined at all timeslices")
return np.arccosh(Corr(newcontent,padding_back=1,padding_front=1)) return np.arccosh(Corr(newcontent, padding_back=1, padding_front=1))
else: else:
raise Exception('Unkown variant.') raise Exception('Unkown variant.')
#We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range. # We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range.
def fit(self, function, fitrange=None, silent=False, **kwargs): def fit(self, function, fitrange=None, silent=False, **kwargs):
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: # The default behaviour is:
#1 use explicit fitrange # 1 use explicit fitrange
# if none is provided, use the range of the corr # 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 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:
if self.prange: if self.prange:
fitrange = self.prange fitrange = self.prange
@ -311,15 +299,13 @@ class Corr:
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]
result = standard_fit(xs, ys, function, silent=silent, **kwargs) result = standard_fit(xs, ys, function, silent=silent, **kwargs)
if isinstance(result, list): if isinstance(result, list):
[item.gamma_method() for item in result if isinstance(item,Obs)] [item.gamma_method() for item in result if isinstance(item, Obs)]
elif isinstance(result, dict): elif isinstance(result, dict):
[item.gamma_method() for item in result['fit_parameters'] if isinstance(item,Obs)] [item.gamma_method() for item in result['fit_parameters'] if isinstance(item, Obs)]
else: else:
raise Exception('Unexpected fit result.') raise Exception('Unexpected fit result.')
return result return result
#we want to quickly get a plateau
def plateau(self, plateau_range=None, method="fit"): def plateau(self, plateau_range=None, method="fit"):
if not plateau_range: if not plateau_range:
if self.prange: if self.prange:
@ -333,9 +319,9 @@ class Corr:
if method == "fit": if method == "fit":
def const_func(a, t): def const_func(a, t):
return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt
return self.fit(const_func,plateau_range)[0] return self.fit(const_func, plateau_range)[0]
elif method in ["avg","average","mean"]: elif method in ["avg", "average", "mean"]:
returnvalue= np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1]+1] if not item is None]) returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
returnvalue.gamma_method() returnvalue.gamma_method()
return returnvalue return returnvalue
@ -343,11 +329,11 @@ class Corr:
raise Exception("Unsupported plateau method: " + method) raise Exception("Unsupported plateau method: " + method)
def set_prange(self, prange): def set_prange(self, prange):
if not len(prange)==2: if not len(prange) == 2:
raise Exception("prange must be a list or array with two values") raise Exception("prange must be a list or array with two values")
if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))): if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
raise Exception("Start and end point must be integers") 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] ): 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") raise Exception("Start and end point must define a range in the interval 0,T")
self.prange = prange self.prange = prange
@ -364,7 +350,7 @@ class Corr:
logscale -- Sets y-axis to logscale logscale -- Sets y-axis to logscale
save -- path to file in which the figure should be saved save -- path to file in which the figure should be saved
""" """
if self.N!=1: if self.N != 1:
raise Exception("Correlator must be projected before plotting") raise Exception("Correlator must be projected before plotting")
if x_range is None: if x_range is None:
x_range = [0, self.T] x_range = [0, self.T]
@ -372,7 +358,7 @@ class Corr:
fig = plt.figure() fig = plt.figure()
ax1 = fig.add_subplot(111) ax1 = fig.add_subplot(111)
x,y,y_err=self.plottable() x, y, y_err = self.plottable()
ax1.errorbar(x, y, y_err, label=self.tag) ax1.errorbar(x, y, y_err, label=self.tag)
if logscale: if logscale:
ax1.set_yscale('log') ax1.set_yscale('log')
@ -380,8 +366,8 @@ class Corr:
# we generate ylim instead of using autoscaling. # we generate ylim instead of using autoscaling.
if y_range is None: if y_range is None:
try: try:
y_min=min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]:x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
y_max=max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]:x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1]] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
except: except:
pass pass
@ -390,7 +376,7 @@ class Corr:
if comp: if comp:
if isinstance(comp, Corr) or isinstance(comp, list): if isinstance(comp, Corr) or isinstance(comp, list):
for corr in comp if isinstance(comp, list) else [comp]: for corr in comp if isinstance(comp, list) else [comp]:
x,y,y_err=corr.plottable() x, y, y_err = corr.plottable()
plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor']) plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
else: else:
raise Exception('comp must be a correlator or a list of correlators.') raise Exception('comp must be a correlator or a list of correlators.')
@ -408,8 +394,8 @@ class Corr:
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)
ax1.plot(x_samples, ax1.plot(x_samples,
fit_res['fit_function']([o.value for o in fit_res['fit_parameters']], x_samples) fit_res['fit_function']([o.value for o in fit_res['fit_parameters']], x_samples),
, ls='-', marker=',', lw=2) ls='-', marker=',', lw=2)
ax1.set_xlabel(r'$x_0 / a$') ax1.set_xlabel(r'$x_0 / a$')
if ylabel: if ylabel:
@ -418,7 +404,7 @@ class Corr:
handles, labels = ax1.get_legend_handles_labels() handles, labels = ax1.get_legend_handles_labels()
if labels: if labels:
legend = ax1.legend() ax1.legend()
plt.draw() plt.draw()
if save: if save:
@ -429,8 +415,8 @@ class Corr:
return return
def dump(self,filename): def dump(self, filename):
dump_object(self,filename) dump_object(self, filename)
return return
def print(self, range=[0, None]): def print(self, range=[0, None]):
@ -449,251 +435,190 @@ class Corr:
content_string += '\t' + element.__repr__()[4:-1] content_string += '\t' + element.__repr__()[4:-1]
content_string += '\n' content_string += '\n'
return content_string return content_string
def __str__(self): def __str__(self):
return self.__repr__() return self.__repr__()
#return ("Corr[T="+str(self.T)+" , N="+str(self.N)+" , content="+str([o[0] for o in [o for o in self.content]])+"]")
#We define the basic operations, that can be performed with correlators. # We define the basic operations, that can be performed with correlators.
#While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
#This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
#One could try and tell Obs to check if the y in __mul__ is a Corr and # One could try and tell Obs to check if the y in __mul__ is a Corr and
def __add__(self, y): def __add__(self, y):
if isinstance(y, Corr): if isinstance(y, Corr):
if ((self.N!=y.N) or (self.T!=y.T) ): if ((self.N != y.N) or (self.T != y.T)):
raise Exception("Addition of Corrs with different shape") raise Exception("Addition of Corrs with different shape")
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None) or (y.content[t] is None): if (self.content[t] is None) or (y.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]+y.content[t]) newcontent.append(self.content[t] + y.content[t])
return Corr(newcontent) return Corr(newcontent)
elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y,float): elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y, float):
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None): if (self.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]+y) newcontent.append(self.content[t] + y)
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
else: else:
raise TypeError("Corr + wrong type") raise TypeError("Corr + wrong type")
def __mul__(self,y): def __mul__(self, y):
if isinstance(y,Corr): if isinstance(y, Corr):
if not((self.N==1 or y.N==1 or self.N==y.N) and self.T==y.T): if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None) or (y.content[t] is None): if (self.content[t] is None) or (y.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]*y.content[t]) newcontent.append(self.content[t] * y.content[t])
return Corr(newcontent) return Corr(newcontent)
elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y,float): elif isinstance(y, Obs) or isinstance(y, int) or isinstance(y, float):
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None): if (self.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]*y) newcontent.append(self.content[t] * y)
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
else: else:
raise TypeError("Corr * wrong type") raise TypeError("Corr * wrong type")
def __truediv__(self,y): def __truediv__(self, y):
if isinstance(y,Corr): if isinstance(y, Corr):
if not((self.N==1 or y.N==1 or self.N==y.N) and self.T==y.T): if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T")
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None) or (y.content[t] is None): if (self.content[t] is None) or (y.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]/y.content[t]) newcontent.append(self.content[t] / y.content[t])
#Here we set the entire timeslice to undefined, if one of the smearings has encountered an division by zero. # Here we set the entire timeslice to undefined, if one of the smearings has encountered an division by zero.
#While this might throw away perfectly good values in other smearings, we will never have to check, if all values in our matrix are defined # While this might throw away perfectly good values in other smearings, we will never have to check, if all values in our matrix are defined
for t in range(self.T): for t in range(self.T):
if newcontent[t] is None: if newcontent[t] is None:
continue continue
if np.isnan(np.sum(newcontent[t]).value): if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None newcontent[t] = None
if all([item is None for item in newcontent]): if all([item is None for item in newcontent]):
raise Exception("Division returns completely undefined correlator") raise Exception("Division returns completely undefined correlator")
return Corr(newcontent) return Corr(newcontent)
elif isinstance(y, Obs): elif isinstance(y, Obs):
if y.value==0: if y.value == 0:
raise Exception("Division by Zero will return undefined correlator") raise Exception('Division by zero will return undefined correlator')
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None): if (self.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]/y) newcontent.append(self.content[t] / y)
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
elif isinstance(y, int) or isinstance(y,float): elif isinstance(y, int) or isinstance(y, float):
if y==0: if y == 0:
raise Exception("Division by Zero will return undefined correlator") raise Exception('Division by zero will return undefined correlator')
newcontent=[] newcontent = []
for t in range(self.T): for t in range(self.T):
if (self.content[t] is None): if (self.content[t] is None):
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(self.content[t]/y) newcontent.append(self.content[t] / y)
return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None) return Corr(newcontent, prange=self.prange)
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,prange=self.prange) return Corr(newcontent, prange=self.prange)
def __sub__(self,y): def __sub__(self, y):
return self +(-y) return self + (-y)
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,prange=self.prange) return Corr(newcontent, prange=self.prange)
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,prange= self.prange if hasattr(self,"prange") else None) return Corr(newcontent, prange=self.prange)
#The numpy functions: # The numpy functions:
def sqrt(self): def sqrt(self):
return self**0.5 return self**0.5
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,prange= self.prange if hasattr(self,"prange") else None) return Corr(newcontent, prange=self.prange)
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,prange= self.prange if hasattr(self,"prange") else None) return Corr(newcontent, prange=self.prange)
def _apply_func_to_corr(self, func):
newcontent = [None if (item is None) else func(item) for item in self.content]
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t] = None
if all([item is None for item in newcontent]):
raise Exception('Operation returns undefined correlator')
return Corr(newcontent)
def sin(self): def sin(self):
newcontent=[None if (item is None) else np.sin(item) for item in self.content] return self._apply_func_to_corr(np.sin)
return Corr(newcontent)
def cos(self): def cos(self):
newcontent=[None if (item is None) else np.cos(item) for item in self.content] return self._apply_func_to_corr(np.cos)
return Corr(newcontent)
def tan(self): def tan(self):
newcontent=[None if (item is None) else np.tan(item) for item in self.content] return self._apply_func_to_corr(np.tan)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def sinh(self): def sinh(self):
newcontent=[None if (item is None) else np.sinh(item) for item in self.content] return self._apply_func_to_corr(np.sinh)
return Corr(newcontent)
def cosh(self): def cosh(self):
newcontent=[None if (item is None) else np.cosh(item) for item in self.content] return self._apply_func_to_corr(np.cosh)
return Corr(newcontent)
def tanh(self): def tanh(self):
newcontent=[None if (item is None) else np.tanh(item) for item in self.content] return self._apply_func_to_corr(np.tanh)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arcsin(self): def arcsin(self):
newcontent=[None if (item is None) else np.arcsin(item) for item in self.content] return self._apply_func_to_corr(np.arcsin)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arccos(self): def arccos(self):
newcontent=[None if (item is None) else np.arccos(item) for item in self.content] return self._apply_func_to_corr(np.arccos)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arctan(self): def arctan(self):
newcontent=[None if (item is None) else np.arctan(item) for item in self.content] return self._apply_func_to_corr(np.arctan)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arcsinh(self): def arcsinh(self):
newcontent=[None if (item is None) else np.arcsinh(item) for item in self.content] return self._apply_func_to_corr(np.arcsinh)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arccosh(self): def arccosh(self):
newcontent=[None if (item is None) else np.arccosh(item) for item in self.content] return self._apply_func_to_corr(np.arccosh)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
def arctanh(self): def arctanh(self):
newcontent=[None if (item is None) else np.arctanh(item) for item in self.content] return self._apply_func_to_corr(np.arctanh)
for t in range(self.T):
if newcontent[t] is None:
continue
if np.isnan(np.sum(newcontent[t]).value):
newcontent[t]=None
if all([item is None for item in newcontent]):
raise Exception("Operation returns completely undefined correlator")
return Corr(newcontent)
# Right hand side operations (require tweak in main module to work)
def __rsub__(self, y):
return -self + y
#right hand side operations (require tweak in main module to work)
def __rsub__(self,y):
return -self+y
def __rmul__(self, y): def __rmul__(self, y):
return self * y return self * y
def __radd__(self,y):
def __radd__(self, y):
return self + y return self + y