diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 4ad8c582..1f47de87 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -457,25 +457,16 @@ class Corr: return (self + T_partner) / 2 - def deriv(self, symmetric=True): + def deriv(self, variant="symmetric"): """Return the first derivative of the correlator with respect to x0. Parameters ---------- - symmetric : bool - decides whether symmetric of simple finite differences are used. Default: True + variant : str + decides which definition of the finite differences derivative is used. + Available choice: symmetric, forward, backward, improved, default: symmetric """ - if not symmetric: - newcontent = [] - for t in range(self.T - 1): - if (self.content[t] is None) or (self.content[t + 1] is None): - newcontent.append(None) - else: - newcontent.append(self.content[t + 1] - self.content[t]) - if(all([x is None for x in newcontent])): - raise Exception("Derivative is undefined at all timeslices") - return Corr(newcontent, padding=[0, 1]) - if symmetric: + if variant == "symmetric": newcontent = [] for t in range(1, self.T - 1): if (self.content[t - 1] is None) or (self.content[t + 1] is None): @@ -485,18 +476,70 @@ class Corr: if(all([x is None for x in newcontent])): raise Exception('Derivative is undefined at all timeslices') return Corr(newcontent, padding=[1, 1]) + elif variant == "forward": + newcontent = [] + for t in range(self.T - 1): + if (self.content[t] is None) or (self.content[t + 1] is None): + newcontent.append(None) + else: + newcontent.append(self.content[t + 1] - self.content[t]) + if(all([x is None for x in newcontent])): + raise Exception("Derivative is undefined at all timeslices") + return Corr(newcontent, padding=[0, 1]) + elif variant == "backward": + newcontent = [] + for t in range(1, self.T): + if (self.content[t - 1] is None) or (self.content[t] is None): + newcontent.append(None) + else: + newcontent.append(self.content[t] - self.content[t - 1]) + if(all([x is None for x in newcontent])): + raise Exception("Derivative is undefined at all timeslices") + return Corr(newcontent, padding=[1, 0]) + elif variant == "improved": + newcontent = [] + for t in range(2, self.T - 2): + if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): + newcontent.append(None) + else: + newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) + if(all([x is None for x in newcontent])): + raise Exception('Derivative is undefined at all timeslices') + return Corr(newcontent, padding=[2, 2]) + else: + raise Exception("Unknown variant.") - def second_deriv(self): - """Return the second derivative of the correlator with respect to x0.""" - newcontent = [] - for t in range(1, self.T - 1): - if (self.content[t - 1] is None) or (self.content[t + 1] is None): - newcontent.append(None) - else: - newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) - if(all([x is None for x in newcontent])): - raise Exception("Derivative is undefined at all timeslices") - return Corr(newcontent, padding=[1, 1]) + def second_deriv(self, variant="symmetric"): + """Return the second derivative of the correlator with respect to x0. + + Parameters + ---------- + variant : str + decides which definition of the finite differences derivative is used. + Available choice: symmetric, improved, default: symmetric + """ + if variant == "symmetric": + newcontent = [] + for t in range(1, self.T - 1): + if (self.content[t - 1] is None) or (self.content[t + 1] is None): + newcontent.append(None) + else: + newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) + if(all([x is None for x in newcontent])): + raise Exception("Derivative is undefined at all timeslices") + return Corr(newcontent, padding=[1, 1]) + elif variant == "improved": + newcontent = [] + for t in range(2, self.T - 2): + if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): + newcontent.append(None) + else: + newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2])) + if(all([x is None for x in newcontent])): + raise Exception("Derivative is undefined at all timeslices") + return Corr(newcontent, padding=[2, 2]) + else: + raise Exception("Unknown variant.") def m_eff(self, variant='log', guess=1.0): """Returns the effective mass of the correlator as correlator object diff --git a/tests/correlators_test.py b/tests/correlators_test.py index fb91b1b0..8cec767c 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -45,7 +45,7 @@ def test_modify_correlator(): exponent = np.random.normal(3, 5) corr_content.append(pe.pseudo_Obs(2 + 10 ** exponent, 10 ** (exponent - 1), 't')) - corr = pe.correlators.Corr(corr_content) + corr = pe.Corr(corr_content) with pytest.warns(RuntimeWarning): corr.symmetric() @@ -53,14 +53,40 @@ def test_modify_correlator(): corr.anti_symmetric() for pad in [0, 2]: - corr = pe.correlators.Corr(corr_content, padding=[pad, pad]) + corr = pe.Corr(corr_content, padding=[pad, pad]) corr.roll(np.random.randint(100)) - corr.deriv(symmetric=True) - corr.deriv(symmetric=False) + corr.deriv(variant="forward") + corr.deriv(variant="symmetric") + corr.deriv(variant="improved") corr.deriv().deriv() - corr.second_deriv() + corr.second_deriv(variant="symmetric") + corr.second_deriv(variant="improved") corr.second_deriv().second_deriv() + for i, e in enumerate(corr.content): + corr.content[i] = None + + for func in [pe.Corr.deriv, pe.Corr.second_deriv]: + for variant in ["symmetric", "improved", "forward", "gibberish", None]: + with pytest.raises(Exception): + func(corr, variant=variant) + + +def test_deriv(): + corr_content = [] + for t in range(24): + exponent = 1.2 + corr_content.append(pe.pseudo_Obs(2 + t ** exponent, 0.2, 't')) + + corr = pe.Corr(corr_content) + + forward = corr.deriv(variant="forward") + backward = corr.deriv(variant="backward") + sym = corr.deriv(variant="symmetric") + assert np.all([o == 0 for o in (0.5 * (forward + backward) - sym)[1:-1]]) + 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]]) + def test_m_eff(): for padding in [0, 4]: