From 014c0d12ce850f429ca893bb6fbe1b0b17063b39 Mon Sep 17 00:00:00 2001
From: Jan Neuendorf <jp333n@gmx.de>
Date: Tue, 5 Oct 2021 12:19:54 +0200
Subject: [PATCH 1/2] The GEVP method now has the parameter "state", deciding
 which eigenvector to return. Correlators get theoptional attribute "range",
 which is preserved under certain operations. Fits and plateaus will default
 to this range.

---
 pyerrors/correlators.py | 92 +++++++++++++++++++++++++++++++++--------
 1 file changed, 75 insertions(+), 17 deletions(-)

diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py
index dc5f181e..a7f16cb8 100644
--- a/pyerrors/correlators.py
+++ b/pyerrors/correlators.py
@@ -1,5 +1,6 @@
 import numpy as np
 import autograd.numpy as anp
+#from scipy.special.orthogonal import _IntegerType
 from .pyerrors import *
 from .fits import standard_fit
 from .roots import find_root
@@ -22,7 +23,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,range=None):
         #All data_input should be a list of things at different timeslices. This needs to be verified
 
         if not (isinstance(data_input, list)):
@@ -56,6 +57,12 @@ class Corr:
         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.range=range
+
         self.gamma_method()
 
 
@@ -98,6 +105,8 @@ class Corr:
         newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content]
         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.
     def smearing(self, i, j):
         if self.N == 1:
@@ -133,7 +142,8 @@ class Corr:
                 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
         if(all([x is None for x in newcontent])):
             raise Exception("Corr could not be symmetrized: No redundant values")
-        return Corr(newcontent)
+        
+        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
 
     def anti_symmetric(self):
 
@@ -148,7 +158,7 @@ class Corr:
                 newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
         if(all([x is None for x in newcontent])):
             raise Exception("Corr could not be symmetrized: No redundant values")
-        return Corr(newcontent)
+        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
 
 
     #This method will symmetrice the matrices and therefore make them positive definit.
@@ -161,7 +171,7 @@ class Corr:
 
 
     #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):
             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")
@@ -171,7 +181,7 @@ class Corr:
                 Gt[i, j] = self.content[ts][i, j].value
 
         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)
         return sp_vec
 
@@ -258,8 +268,17 @@ class Corr:
         if self.N != 1:
             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:
-            fitrange=[0, self.T]
+            if hasattr(self,"range"):
+                fitrange=self.range
+            else:
+                fitrange=[0, self.T] 
 
         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]
@@ -273,7 +292,12 @@ class Corr:
         return result
 
     #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,"range"):
+                plateau_range=self.range 
+            else: 
+                raise Exception("no plateau range provided")
         if self.N != 1:
             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])])):
@@ -288,7 +312,26 @@ class Corr:
             return returnvalue
 
         else:
-            raise Exception("Unsupported plateau method: " + method)
+            raise Exception("Unsupported plateau method: " + method) 
+
+
+
+    def set_range(self,range): 
+        if not len(range)==2:
+            raise Exception("range must be a list or array with two values")
+        if not ((isinstance(range[0],int)) and (isinstance(range[1],int))):
+            raise Exception("start and end point must be integers")
+        if not (0<=range[0]<=self.T and 0<=range[1]<=self.T and range[0]<range[1]  ):
+            raise Exception("start and end point must define a range in the interval 0,T")
+
+        self.range=range 
+        return 
+
+
+
+
+
+
 
     #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
@@ -337,6 +380,9 @@ class Corr:
                 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
             else:
                 raise Exception('plateau must be an Obs')
+        if hasattr(self,"range"):
+            ax1.axvline(self.range[0],0,1)
+            ax1.axvline(self.range[1],0,1)
 
         if fit_res:
             x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
@@ -367,6 +413,10 @@ class Corr:
     def print(self, range=[0, None]):
         print(self.__repr__(range))
 
+
+
+    
+
     def __repr__(self, range=[0, None]):
         if range[1]:
             range[1] += 1
@@ -408,7 +458,7 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]+y)
-            return Corr(newcontent)
+            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
         else:
             raise TypeError("Corr + wrong type")
 
@@ -431,7 +481,7 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]*y)
-            return Corr(newcontent)
+            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
         else:
             raise TypeError("Corr * wrong type")
 
@@ -469,7 +519,11 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]/y)
-            return Corr(newcontent)
+            if hasattr(self,"range"):
+                newrange=self.range
+            else:
+                newrange=None
+            return Corr(newcontent,range=newrange)
 
         elif isinstance(y, int) or isinstance(y,float):
             if y==0:
@@ -480,13 +534,17 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]/y)
-            return Corr(newcontent)
+            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
         else:
             raise TypeError("Corr / wrong type")
 
     def __neg__(self):
         newcontent=[None if (item is None) else -1.*item for item in self.content]
-        return Corr(newcontent)
+        if hasattr(self,"range"):
+            newrange=self.range
+        else:
+            newrange=None
+        return Corr(newcontent,range=newrange)
 
     def __sub__(self,y):
         return self +(-y)
@@ -494,13 +552,13 @@ class Corr:
     def __pow__(self, y):
         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]
-            return Corr(newcontent)
+            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
         else:
             raise TypeError("type of exponent not supported")
 
     def __abs__(self):
         newcontent=[None if (item is None) else np.abs(item) for item in self.content]
-        return Corr(newcontent)
+        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
 
 #The numpy functions:
     def sqrt(self):
@@ -508,11 +566,11 @@ class Corr:
 
     def log(self):
         newcontent=[None if (item is None) else np.log(item) for item in self.content]
-        return Corr(newcontent)
+        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
 
     def exp(self):
         newcontent=[None if (item is None) else np.exp(item) for item in self.content]
-        return Corr(newcontent)
+        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
 
     def sin(self):
         newcontent=[None if (item is None) else np.sin(item) for item in self.content]

From 7c03cff42fc802e5a8cb0e9b4c613ad0c20274ca Mon Sep 17 00:00:00 2001
From: JanNeuendorf <jp333n@gmx.de>
Date: Wed, 6 Oct 2021 13:16:04 +0200
Subject: [PATCH 2/2] I changed the name from range to prange (hope i did not
 miss something) There is now an Eigenvalue method using the cholesky method
 you put in the example. I do not really use that, but it seems like a logical
 inclusion. I gave the standard fit an option to seht the number of fit
 parameters by hand, because it does not work automatically  if the function
 is not consistent between calls.

---
 pyerrors/correlators.py | 96 ++++++++++++++++++++++++++++-------------
 pyerrors/fits.py        | 22 ++++++----
 2 files changed, 80 insertions(+), 38 deletions(-)

diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py
index a7f16cb8..aedf0e54 100644
--- a/pyerrors/correlators.py
+++ b/pyerrors/correlators.py
@@ -3,6 +3,7 @@ import autograd.numpy as anp
 #from scipy.special.orthogonal import _IntegerType
 from .pyerrors import *
 from .fits import standard_fit
+from .linalg import *
 from .roots import find_root
 from matplotlib import pyplot as plt
 from matplotlib.ticker import NullFormatter  # useful for `logit` scale
@@ -23,7 +24,7 @@ class Corr:
 
     """
 
-    def __init__(self, data_input, padding_front=0, padding_back=0,range=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
 
         if not (isinstance(data_input, list)):
@@ -61,7 +62,7 @@ class Corr:
         #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.range=range
+            self.prange=prange
 
         self.gamma_method()
 
@@ -143,7 +144,7 @@ class Corr:
         if(all([x is None for x in newcontent])):
             raise Exception("Corr could not be symmetrized: No redundant values")
         
-        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
 
     def anti_symmetric(self):
 
@@ -158,7 +159,7 @@ class Corr:
                 newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
         if(all([x is None for x in newcontent])):
             raise Exception("Corr could not be symmetrized: No redundant values")
-        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
 
 
     #This method will symmetrice the matrices and therefore make them positive definit.
@@ -186,6 +187,32 @@ class Corr:
         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):
         return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 
@@ -260,6 +287,17 @@ class Corr:
                 raise Exception('m_eff is undefined at all timeslices')
 
             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:
             raise Exception('Unkown variant.')
 
@@ -275,8 +313,8 @@ class Corr:
 
 
         if fitrange is None:
-            if hasattr(self,"range"):
-                fitrange=self.range
+            if hasattr(self,"prange"):
+                fitrange=self.prange
             else:
                 fitrange=[0, self.T] 
 
@@ -294,8 +332,8 @@ class Corr:
     #we want to quickly get a plateau
     def plateau(self, plateau_range=None, method="fit"):
         if not plateau_range:
-            if hasattr(self,"range"):
-                plateau_range=self.range 
+            if hasattr(self,"prange"):
+                plateau_range=self.prange 
             else: 
                 raise Exception("no plateau range provided")
         if self.N != 1:
@@ -316,15 +354,15 @@ class Corr:
 
 
 
-    def set_range(self,range): 
-        if not len(range)==2:
+    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(range[0],int)) and (isinstance(range[1],int))):
+        if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))):
             raise Exception("start and end point must be integers")
-        if not (0<=range[0]<=self.T and 0<=range[1]<=self.T and range[0]<range[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")
 
-        self.range=range 
+        self.prange=prange 
         return 
 
 
@@ -380,9 +418,9 @@ class Corr:
                 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
             else:
                 raise Exception('plateau must be an Obs')
-        if hasattr(self,"range"):
-            ax1.axvline(self.range[0],0,1)
-            ax1.axvline(self.range[1],0,1)
+        if hasattr(self,"prange"):
+            ax1.axvline(self.prange[0],0,1)
+            ax1.axvline(self.prange[1],0,1)
 
         if fit_res:
             x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
@@ -458,7 +496,7 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]+y)
-            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
         else:
             raise TypeError("Corr + wrong type")
 
@@ -481,7 +519,7 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]*y)
-            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
         else:
             raise TypeError("Corr * wrong type")
 
@@ -519,11 +557,11 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]/y)
-            if hasattr(self,"range"):
-                newrange=self.range
+            if hasattr(self,"prange"):
+                newrange=self.prange
             else:
                 newrange=None
-            return Corr(newcontent,range=newrange)
+            return Corr(newcontent,prange=newrange)
 
         elif isinstance(y, int) or isinstance(y,float):
             if y==0:
@@ -534,17 +572,17 @@ class Corr:
                     newcontent.append(None)
                 else:
                     newcontent.append(self.content[t]/y)
-            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
         else:
             raise TypeError("Corr / wrong type")
 
     def __neg__(self):
         newcontent=[None if (item is None) else -1.*item for item in self.content]
-        if hasattr(self,"range"):
-            newrange=self.range
+        if hasattr(self,"prange"):
+            newrange=self.prange
         else:
             newrange=None
-        return Corr(newcontent,range=newrange)
+        return Corr(newcontent,prange=newrange)
 
     def __sub__(self,y):
         return self +(-y)
@@ -552,13 +590,13 @@ class Corr:
     def __pow__(self, y):
         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]
-            return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+            return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
         else:
             raise TypeError("type of exponent not supported")
 
     def __abs__(self):
         newcontent=[None if (item is None) else np.abs(item) for item in self.content]
-        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
 
 #The numpy functions:
     def sqrt(self):
@@ -566,11 +604,11 @@ class Corr:
 
     def log(self):
         newcontent=[None if (item is None) else np.log(item) for item in self.content]
-        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
 
     def exp(self):
         newcontent=[None if (item is None) else np.exp(item) for item in self.content]
-        return Corr(newcontent,range= self.range if hasattr(self,"range") else None)
+        return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None)
 
     def sin(self):
         newcontent=[None if (item is None) else np.sin(item) for item in self.content]
diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index c0f394ca..1d5ba31d 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -14,7 +14,7 @@ from autograd import elementwise_grad as egrad
 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.
 
     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):
         raise TypeError('func has to be a function.')
 
-    for i in range(25):
-        try:
-            func(np.arange(i), x.T[0])
-        except:
-            pass
-        else:
-            break
+    if n_parms=="auto":
+        for i in range(25):
+            try:
+                func(np.arange(i), x.T[0])
+            except:
+                pass
+            else:
+                break
+        
 
-    n_parms = i
+        n_parms = i
+
+        
     if not silent:
         print('Fit with', n_parms, 'parameters')