From c9f65298b467cfe389bac7cccb21a4030ee51112 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 7 Apr 2022 17:14:26 +0100 Subject: [PATCH 1/4] tests: test added that checks that correlated and uncorrelated fits give the same result if the input data is uncorrelated. --- tests/fits_test.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/fits_test.py b/tests/fits_test.py index b6236fb9..374550e4 100644 --- a/tests/fits_test.py +++ b/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) From a5394521e7caec93dc60944bb91a19848e6b3714 Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Wed, 13 Apr 2022 18:20:13 +0200 Subject: [PATCH 2/4] Fix: read_mesons and read_dSdm now return Obs with the correct idl based on the information in the bdio file. --- pyerrors/input/bdio.py | 88 ++++++++++++++++++++++++++++++++---------- 1 file changed, 68 insertions(+), 20 deletions(-) diff --git a/pyerrors/input/bdio.py b/pyerrors/input/bdio.py index ef0e7f68..75ee676b 100644 --- a/pyerrors/input/bdio.py +++ b/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 - bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so) - stop -- stops reading at given configuration number (default None) - alternative_ensemble_name -- Manually overwrite ensemble name + file_path : str + path to the bdio file + 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 @@ -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: - if cnfg_no >= kwargs.get('stop') - 1: + cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID=')) + 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') - 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,11 +467,20 @@ 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) + 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): tmp_corr = [] for t in range(d0 - 2): - tmp_corr.append(Obs([np.asarray(data[c])[:, t]], [ensemble_name])) + deltas = [np.asarray(data[c])[:, t][index] for index in indices] + tmp_corr.append(Obs([deltas], [ensemble_name], idl=[idl_target])) result[(corr_name[c], corr_source[c]) + tuple(sorted(corr_kappa[c]))] = tmp_corr # Check that all data entries have the same number of configurations @@ -480,10 +504,24 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): Parameters ---------- - file_path -- path to the bdio file - bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so) - stop -- stops reading at given configuration number (default None) + file_path : str + path to the bdio file + 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 +567,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 +624,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: - if cnfg_no >= kwargs.get('stop') - 1: + cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID=')) + 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') - if cnfg_no == 0: + if len(idl) == 1: no_corrs = len(corr_name) data = [] for c in range(no_corrs): @@ -612,9 +651,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: From ee631c3ee028ae7071c6f080866c67837f62328c Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Tue, 19 Apr 2022 12:27:49 +0200 Subject: [PATCH 3/4] Removed sorting of kappas and speed improvement in read_mesons --- pyerrors/input/bdio.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/pyerrors/input/bdio.py b/pyerrors/input/bdio.py index 75ee676b..f54dfbcd 100644 --- a/pyerrors/input/bdio.py +++ b/pyerrors/input/bdio.py @@ -436,7 +436,7 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs): if cnfg_no > kwargs.get('stop'): break 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 len(idl) == 1: no_corrs = len(corr_name) data = [] @@ -470,18 +470,26 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs): 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) + + 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): - deltas = [np.asarray(data[c])[:, t][index] for index in indices] + if indices: + 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(sorted(corr_kappa[c]))] = tmp_corr + 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: @@ -629,7 +637,7 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): if cnfg_no > kwargs.get('stop'): break 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 len(idl) == 1: no_corrs = len(corr_name) data = [] From 9b356c9e97197006d9a0d74a2ad9eb099385de72 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 19 Apr 2022 13:13:45 +0100 Subject: [PATCH 4/4] docs: clarification on Corr.fit fitrange added. --- pyerrors/correlators.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index ac3f8b07..6ddf3204 100644 --- a/pyerrors/correlators.py +++ b/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")