Merge pull request #157 from fjosw/fix/calc_rho_large_t

_calculate_drho for large time separations

.github/workflows/flake8.yml

@@ -21,6 +21,6 @@ jobs:
       - name: flake8 Lint
         uses: py-actions/flake8@v2
         with:
-          ignore: "E501"
+          ignore: "E501,W503"
           exclude: "__init__.py, input/__init__.py"
           path: "pyerrors"

CONTRIBUTING.md

@@ -33,6 +33,6 @@ pytest --cov=pyerrors --cov-report html
 ```
 The linter `flake8` is executed with the command
 ```
-flake8 --ignore=E501 --exclude=__init__.py pyerrors
+flake8 --ignore=E501,W503 --exclude=__init__.py pyerrors
 ```
 Please make sure that all tests are passed for a new pull requests.

pyerrors/obs.py

@@ -289,7 +289,10 @@ class Obs:
             self.e_n_dtauint[e_name][0] = 0.0
 
             def _compute_drho(i):
-                tmp = self.e_rho[e_name][i + 1:w_max] + np.concatenate([self.e_rho[e_name][i - 1::-1], self.e_rho[e_name][1:w_max - 2 * i]]) - 2 * self.e_rho[e_name][i] * self.e_rho[e_name][1:w_max - i]
+                tmp = (self.e_rho[e_name][i + 1:w_max]
+                       + np.concatenate([self.e_rho[e_name][i - 1:None if i - w_max // 2 < 0 else 2 * (i - w_max // 2):-1],
+                                         self.e_rho[e_name][1:max(1, w_max - 2 * i)]])
+                       - 2 * self.e_rho[e_name][i] * self.e_rho[e_name][1:w_max - i])
                 self.e_drho[e_name][i] = np.sqrt(np.sum(tmp ** 2) / e_N)
 
             if self.tau_exp[e_name] > 0:
@@ -320,8 +323,8 @@ class Obs:
                     # Standard automatic windowing procedure
                     tau = self.S[e_name] / np.log((2 * self.e_n_tauint[e_name][gapsize::gapsize] + 1) / (2 * self.e_n_tauint[e_name][gapsize::gapsize] - 1))
                     g_w = np.exp(- np.arange(1, len(tau) + 1) / tau) - tau / np.sqrt(np.arange(1, len(tau) + 1) * e_N)
-                    for n in range(1, w_max):
-                        if g_w[n - 1] < 0 or n >= w_max - 1:
+                    for n in range(1, w_max // gapsize):
+                        if g_w[n - 1] < 0 or n >= w_max // gapsize - 1:
                             _compute_drho(gapsize * n)
                             n *= gapsize
                             self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 1) / e_N) / (1 + 1 / e_N)  # Bias correction hep-lat/0306017 eq. (49)

tests/obs_test.py

@@ -1146,3 +1146,13 @@ def test_non_overlapping_operations_different_lengths():
 
         assert np.isclose(res1.value, res2.value)
         assert np.isclose(res1.dvalue, res2.dvalue, rtol=0.01)
+
+
+def test_nan_obs():
+    o = pe.pseudo_Obs(1, .1, 'test')
+    no = np.nan * o
+    no.gamma_method()
+
+    o.idl['test'] = [1, 5] + list(range(7, 2002, 2))
+    no = np.NaN * o
+    no.gamma_method()