Merge pull request #127 from s-kuberski/feat/deriv

feat: added log-derivatives and symmetric effective mass

pyerrors/correlators.py

@@ -523,7 +523,7 @@ class Corr:
         ----------
         variant : str
             decides which definition of the finite differences derivative is used.
-            Available choice: symmetric, forward, backward, improved, default: symmetric
+            Available choice: symmetric, forward, backward, improved, log, default: symmetric
         """
         if self.N != 1:
             raise Exception("deriv only implemented for one-dimensional correlators.")
@@ -567,6 +567,17 @@ class Corr:
             if (all([x is None for x in newcontent])):
                 raise Exception('Derivative is undefined at all timeslices')
             return Corr(newcontent, padding=[2, 2])
+        elif variant == 'log':
+            newcontent = []
+            for t in range(self.T):
+                if (self.content[t] is None) or (self.content[t] <= 0):
+                    newcontent.append(None)
+                else:
+                    newcontent.append(np.log(self.content[t]))
+            if (all([x is None for x in newcontent])):
+                raise Exception("Log is undefined at all timeslices")
+            logcorr = Corr(newcontent)
+            return self * logcorr.deriv('symmetric')
         else:
             raise Exception("Unknown variant.")
 
@@ -577,7 +588,7 @@ class Corr:
         ----------
         variant : str
             decides which definition of the finite differences derivative is used.
-            Available choice: symmetric, improved, default: symmetric
+            Available choice: symmetric, improved, log, default: symmetric
         """
         if self.N != 1:
             raise Exception("second_deriv only implemented for one-dimensional correlators.")
@@ -601,6 +612,17 @@ class Corr:
             if (all([x is None for x in newcontent])):
                 raise Exception("Derivative is undefined at all timeslices")
             return Corr(newcontent, padding=[2, 2])
+        elif variant == 'log':
+            newcontent = []
+            for t in range(self.T):
+                if (self.content[t] is None) or (self.content[t] <= 0):
+                    newcontent.append(None)
+                else:
+                    newcontent.append(np.log(self.content[t]))
+            if (all([x is None for x in newcontent])):
+                raise Exception("Log is undefined at all timeslices")
+            logcorr = Corr(newcontent)
+            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
         else:
             raise Exception("Unknown variant.")
 
@@ -615,6 +637,7 @@ class Corr:
             sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
             See, e.g., arXiv:1205.5380
             arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
         guess : float
             guess for the root finder, only relevant for the root variant
         """
@@ -634,6 +657,20 @@ class Corr:
 
             return np.log(Corr(newcontent, padding=[0, 1]))
 
+        elif variant == 'logsym':
+            newcontent = []
+            for t in range(1, self.T - 1):
+                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+                    newcontent.append(None)
+                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
+                    newcontent.append(None)
+                else:
+                    newcontent.append(self.content[t - 1] / self.content[t + 1])
+            if (all([x is None for x in newcontent])):
+                raise Exception('m_eff is undefined at all timeslices')
+
+            return np.log(Corr(newcontent, padding=[1, 1])) / 2
+
         elif variant in ['periodic', 'cosh', 'sinh']:
             if variant in ['periodic', 'cosh']:
                 func = anp.cosh

tests/correlators_test.py

@@ -90,14 +90,38 @@ def test_deriv():
     assert np.all([o == 0 for o in (corr.deriv('forward').deriv('backward') - corr.second_deriv())[1:-1]])
     assert np.all([o == 0 for o in (corr.deriv('backward').deriv('forward') - corr.second_deriv())[1:-1]])
 
+    corr_content = []
+    exponent = -0.05
+    for t in range(24):
+        corr_content.append(pe.pseudo_Obs(np.exp(t * exponent), np.exp(t * exponent) * 0.02, 't'))
+
+    corr = pe.Corr(corr_content)
+
+    for o in [(corr.deriv('log') / corr / exponent - 1)[10], (corr.second_deriv('log') / corr / exponent**2 - 1)[12]]:
+        o.gamma_method()
+        assert (o.is_zero_within_error() and np.isclose(0.0, o.value, 1e-12, 1e-12))
+
 
 def test_m_eff():
     for padding in [0, 4]:
         my_corr = pe.correlators.Corr([pe.pseudo_Obs(10, 0.1, 't'), pe.pseudo_Obs(9, 0.05, 't'), pe.pseudo_Obs(9, 0.1, 't'), pe.pseudo_Obs(10, 0.05, 't')], padding=[padding, padding])
         my_corr.m_eff('log')
+        my_corr.m_eff('logsym')
         my_corr.m_eff('cosh')
         my_corr.m_eff('arccosh')
 
+    corr_content = []
+    exponent = -2.2
+    for t in range(24):
+        corr_content.append(pe.pseudo_Obs(np.exp(t * exponent), np.exp(t * exponent) * 0.02, 't'))
+
+    corr = pe.Corr(corr_content)
+
+    for variant in ['log', 'logsym']:
+        o = (corr.m_eff(variant) / exponent + 1)[7]
+        o.gamma_method()
+        assert (o.is_zero_within_error() and np.isclose(0.0, o.value, 1e-12, 1e-12))
+
     with pytest.warns(RuntimeWarning):
         my_corr.m_eff('sinh')
 
@@ -112,6 +136,8 @@ def test_m_eff_negative_values():
         assert m_eff_log[padding + 1] is None
         m_eff_cosh = my_corr.m_eff('cosh')
         assert m_eff_cosh[padding + 1] is None
+        with pytest.raises(Exception):
+            my_corr.m_eff('logsym')
 
 
 def test_reweighting():

tests/fits_test.py

@@ -235,8 +235,8 @@ def test_fit_corr_independent():
         assert np.isclose(out.chisquare, out_corr.chisquare)
         assert np.isclose(out.dof, out_corr.dof)
         assert np.isclose(out.chisquare_by_dof, out_corr.chisquare_by_dof)
-        assert (out[0] - out_corr[0]).is_zero(atol=1e-5)
+        assert (out[0] - out_corr[0]).is_zero(atol=1e-4)
-        assert (out[1] - out_corr[1]).is_zero(atol=1e-5)
+        assert (out[1] - out_corr[1]).is_zero(atol=1e-4)
 
 
 def test_linear_fit_guesses():