Merge branch 'develop' into fix/covariance_idl

pyerrors/correlators.py

@@ -621,7 +621,7 @@ class Corr:
             raise Exception('Unknown variant.')
 
     def fit(self, function, fitrange=None, silent=False, **kwargs):
-        """Fits function to the data
+        r'''Fits function to the data
 
         Parameters
         ----------
@@ -629,10 +629,12 @@ class Corr:
             function to fit to the data. See fits.least_squares for details.
         fitrange : list
             Two element list containing the timeslices on which the fit is supposed to start and stop.
+            Caution: This range is inclusive as opposed to standard python indexing.
+            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
             If not specified, self.prange or all timeslices are used.
         silent : bool
             Decides whether output is printed to the standard output.
-        """
+        '''
         if self.N != 1:
             raise Exception("Correlator must be projected before fitting")
 

pyerrors/input/bdio.py

@@ -302,11 +302,24 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
 
     Parameters
     ----------
-    file_path -- path to the bdio file
+    file_path : str
-    bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
+        path to the bdio file
-    stop -- stops reading at given configuration number (default None)
+    bdio_path : str
-    alternative_ensemble_name -- Manually overwrite ensemble name
+        path to the shared bdio library libbdio.so (default ./libbdio.so)
+    start : int
+        The first configuration to be read (default 1)
+    stop : int
+        The last configuration to be read (default None)
+    step : int
+        Fixed step size between two measurements (default 1)
+    alternative_ensemble_name : str
+        Manually overwrite ensemble name
     """
+
+    start = kwargs.get('start', 1)
+    stop = kwargs.get('stop', None)
+    step = kwargs.get('step', 1)
+
     bdio = ctypes.cdll.LoadLibrary(bdio_path)
 
     bdio_open = bdio.bdio_open
@@ -353,9 +366,9 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
     prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
     prop_source = []  # Contains propagator source positions
     # Check noise type for multiple replica?
-    cnfg_no = -1
     corr_no = -1
     data = []
+    idl = []
 
     fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
 
@@ -418,18 +431,20 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
                 prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
                 prop_source.append(_get_kwd(tmp_string, 'x0='))
             if ruinfo == 4:
-                if 'stop' in kwargs:
+                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
-                    if cnfg_no >= kwargs.get('stop') - 1:
+                if stop:
+                    if cnfg_no > kwargs.get('stop'):
                         break
-                cnfg_no += 1
+                idl.append(cnfg_no)
-                print('\r%s %i' % ('Reading configuration', cnfg_no + 1), end='\r')
+                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
-                if cnfg_no == 0:
+                if len(idl) == 1:
                     no_corrs = len(corr_name)
                     data = []
                     for c in range(no_corrs):
                         data.append([])
 
                 corr_no = 0
+
     bdio_close(fbdio)
 
     print('\nEnsemble: ', ensemble_name)
@@ -441,7 +456,7 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
     print('Number of time values: ', d0)
     print('Number of random sources: ', d1)
     print('Number of corrs: ', len(corr_name))
-    print('Number of configurations: ', cnfg_no + 1)
+    print('Number of configurations: ', len(idl))
 
     corr_kappa = []  # Contains kappa values for both propagators of given correlation function
     corr_source = []
@@ -452,12 +467,29 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
         else:
             corr_source.append(int(prop_source[int(item[0])]))
 
+    if stop is None:
+        stop = idl[-1]
+    idl_target = range(start, stop + 1, step)
+
+    if set(idl) != set(idl_target):
+        try:
+            indices = [idl.index(i) for i in idl_target]
+        except ValueError as err:
+            raise Exception('Configurations in file do no match target list!', err)
+    else:
+        indices = None
+
     result = {}
     for c in range(no_corrs):
         tmp_corr = []
+        tmp_data = np.asarray(data[c])
         for t in range(d0 - 2):
-            tmp_corr.append(Obs([np.asarray(data[c])[:, t]], [ensemble_name]))
+            if indices:
-        result[(corr_name[c], corr_source[c]) + tuple(sorted(corr_kappa[c]))] = tmp_corr
+                deltas = [tmp_data[:, t][index] for index in indices]
+            else:
+                deltas = tmp_data[:, t]
+            tmp_corr.append(Obs([deltas], [ensemble_name], idl=[idl_target]))
+        result[(corr_name[c], corr_source[c]) + tuple(corr_kappa[c])] = tmp_corr
 
     # Check that all data entries have the same number of configurations
     if len(set([o[0].N for o in list(result.values())])) != 1:
@@ -480,10 +512,24 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
 
     Parameters
     ----------
-    file_path -- path to the bdio file
+    file_path : str
-    bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
+        path to the bdio file
-    stop -- stops reading at given configuration number (default None)
+    bdio_path : str
+        path to the shared bdio library libbdio.so (default ./libbdio.so)
+    start : int
+        The first configuration to be read (default 1)
+    stop : int
+        The last configuration to be read (default None)
+    step : int
+        Fixed step size between two measurements (default 1)
+    alternative_ensemble_name : str
+        Manually overwrite ensemble name
     """
+
+    start = kwargs.get('start', 1)
+    stop = kwargs.get('stop', None)
+    step = kwargs.get('step', 1)
+
     bdio = ctypes.cdll.LoadLibrary(bdio_path)
 
     bdio_open = bdio.bdio_open
@@ -529,9 +575,9 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
     # d1 = 0  # nnoise
     prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
     # Check noise type for multiple replica?
-    cnfg_no = -1
     corr_no = -1
     data = []
+    idl = []
 
     fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
 
@@ -586,12 +632,13 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
             if ruinfo == 2:
                 prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
             if ruinfo == 4:
-                if 'stop' in kwargs:
+                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
-                    if cnfg_no >= kwargs.get('stop') - 1:
+                if stop:
+                    if cnfg_no > kwargs.get('stop'):
                         break
-                cnfg_no += 1
+                idl.append(cnfg_no)
-                print('\r%s %i' % ('Reading configuration', cnfg_no + 1), end='\r')
+                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
-                if cnfg_no == 0:
+                if len(idl) == 1:
                     no_corrs = len(corr_name)
                     data = []
                     for c in range(no_corrs):
@@ -612,9 +659,18 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
     for item in corr_props:
         corr_kappa.append(float(prop_kappa[int(item)]))
 
+    if stop is None:
+        stop = idl[-1]
+    idl_target = range(start, stop + 1, step)
+    try:
+        indices = [idl.index(i) for i in idl_target]
+    except ValueError as err:
+        raise Exception('Configurations in file do no match target list!', err)
+
     result = {}
     for c in range(no_corrs):
-        result[(corr_name[c], str(corr_kappa[c]))] = Obs([np.asarray(data[c])], [ensemble_name])
+        deltas = [np.asarray(data[c])[index] for index in indices]
+        result[(corr_name[c], str(corr_kappa[c]))] = Obs([deltas], [ensemble_name], idl=[idl_target])
 
     # Check that all data entries have the same number of configurations
     if len(set([o.N for o in list(result.values())])) != 1:

tests/fits_test.py

@@ -157,6 +157,28 @@ def test_correlated_fit():
             assert(diff.is_zero_within_error(sigma=5))
 
 
+def test_fit_corr_independent():
+    dim = 50
+    x = np.arange(dim)
+    y = 0.84 * np.exp(-0.12 * x) + np.random.normal(0.0, 0.1, dim)
+    yerr = [0.1] * dim
+
+    oy = []
+    for i, item in enumerate(x):
+        oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i)))
+
+    def func(a, x):
+        y = a[0] * anp.exp(-a[1] * x)
+        return y
+
+    out = pe.least_squares(x, oy, func)
+    out_corr = pe.least_squares(x, oy, func, correlated_fit=True)
+
+    assert np.isclose(out.chisquare, out_corr.chisquare)
+    assert (out[0] - out_corr[0]).is_zero(atol=1e-5)
+    assert (out[1] - out_corr[1]).is_zero(atol=1e-5)
+
+
 def test_total_least_squares():
     dim = 10 + int(30 * np.random.rand())
     x = np.arange(dim) + np.random.normal(0.0, 0.15, dim)