diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py index 354a3fd8..18c6ec66 100644 --- a/pyerrors/input/hadrons.py +++ b/pyerrors/input/hadrons.py @@ -6,6 +6,7 @@ import numpy as np from ..obs import Obs, CObs from ..correlators import Corr from ..dirac import epsilon_tensor_rank4 +from .misc import fit_t0 def _get_files(path, filestem, idl): @@ -121,6 +122,72 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=Non return corr +def _extract_real_arrays(path, files, tree, keys): + corr_data = {} + for key in keys: + corr_data[key] = [] + for hd5_file in files: + h5file = h5py.File(path + '/' + hd5_file, "r") + for key in keys: + if not tree + '/' + key in h5file: + raise Exception("Entry '" + key + "' not contained in the files.") + raw_data = h5file[tree + '/' + key + '/data'] + real_data = raw_data[:].astype(np.double) + corr_data[key].append(real_data) + h5file.close() + for key in keys: + corr_data[key] = np.array(corr_data[key]) + return corr_data + + +def extract_t0_hd5(path, filestem, ens_id, obs='Clover energy density', fit_range=5, idl=None, **kwargs): + r'''Read hadrons FlowObservables hdf5 file and extract t0 + + Parameters + ----------------- + path : str + path to the files to read + filestem : str + namestem of the files to read + ens_id : str + name of the ensemble, required for internal bookkeeping + obs : str + label of the observable from which t0 should be extracted. + Options: 'Clover energy density' and 'Plaquette energy density' + fit_range : int + Number of data points left and right of the zero + crossing to be included in the linear fit. (Default: 5) + idl : range + If specified only configurations in the given range are read in. + plot_fit : bool + If true, the fit for the extraction of t0 is shown together with the data. + ''' + + files, idx = _get_files(path, filestem, idl) + tree = "FlowObservables" + + h5file = h5py.File(path + '/' + files[0], "r") + obs_key = None + for key in h5file[tree].keys(): + if obs == h5file[tree][key].attrs["description"][0].decode(): + obs_key = key + break + h5file.close() + if obs_key is None: + raise Exception(f"Observable {obs} not found.") + + corr_data = _extract_real_arrays(path, files, tree, ["FlowObservables_0", obs_key]) + + if not np.allclose(corr_data["FlowObservables_0"][0], corr_data["FlowObservables_0"][:]): + raise Exception("Not all flow times were equal.") + + t2E_dict = {} + for t2, dat in zip(corr_data["FlowObservables_0"][0], corr_data[obs_key].T): + t2E_dict[t2] = Obs([dat], [ens_id], idl=[idx]) - 0.3 + + return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit')) + + def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None): """Read hadrons DistillationContraction hdf5 files in given directory structure diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index f56c1ac5..b6cbccf7 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -3,7 +3,58 @@ import fnmatch import re import struct import numpy as np # Thinly-wrapped numpy +import matplotlib.pyplot as plt +from matplotlib import gridspec from ..obs import Obs +from ..fits import fit_lin + + +def fit_t0(t2E_dict, fit_range, plot_fit=False): + zero_crossing = np.argmax(np.array( + [o.value for o in t2E_dict.values()]) > 0.0) + + x = list(t2E_dict.keys())[zero_crossing - fit_range: + zero_crossing + fit_range] + y = list(t2E_dict.values())[zero_crossing - fit_range: + zero_crossing + fit_range] + [o.gamma_method() for o in y] + + fit_result = fit_lin(x, y) + + if plot_fit is True: + plt.figure() + gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) + ax0 = plt.subplot(gs[0]) + xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] + ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] + [o.gamma_method() for o in ymore] + ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') + xplot = np.linspace(np.min(x), np.max(x)) + yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] + [yi.gamma_method() for yi in yplot] + ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) + retval = (-fit_result[0] / fit_result[1]) + retval.gamma_method() + ylim = ax0.get_ylim() + ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) + ax0.set_ylim(ylim) + ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') + xlim = ax0.get_xlim() + + fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] + residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) + ax1 = plt.subplot(gs[1]) + ax1.plot(x, residuals, 'ko', ls='none', markersize=5) + ax1.tick_params(direction='out') + ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) + ax1.axhline(y=0.0, ls='--', color='k') + ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') + ax1.set_xlim(xlim) + ax1.set_ylabel('Residuals') + ax1.set_xlabel(r'$t/a^2$') + + plt.draw() + return -fit_result[0] / fit_result[1] def read_pbp(path, prefix, **kwargs): diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index d0e3e71d..3712f3c5 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -3,12 +3,10 @@ import fnmatch import struct import warnings import numpy as np # Thinly-wrapped numpy -import matplotlib.pyplot as plt -from matplotlib import gridspec from ..obs import Obs -from ..fits import fit_lin from ..obs import CObs from ..correlators import Corr +from .misc import fit_t0 from .utils import sort_names @@ -428,51 +426,7 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfi new_obs = Obs(samples, rep_names, idl=idl) t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3 - zero_crossing = np.argmax(np.array( - [o.value for o in t2E_dict.values()]) > 0.0) - - x = list(t2E_dict.keys())[zero_crossing - fit_range: - zero_crossing + fit_range] - y = list(t2E_dict.values())[zero_crossing - fit_range: - zero_crossing + fit_range] - [o.gamma_method() for o in y] - - fit_result = fit_lin(x, y) - - if kwargs.get('plot_fit'): - plt.figure() - gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) - ax0 = plt.subplot(gs[0]) - xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] - ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2] - [o.gamma_method() for o in ymore] - ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x') - xplot = np.linspace(np.min(x), np.max(x)) - yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot] - [yi.gamma_method() for yi in yplot] - ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot]) - retval = (-fit_result[0] / fit_result[1]) - retval.gamma_method() - ylim = ax0.get_ylim() - ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4) - ax0.set_ylim(ylim) - ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $') - xlim = ax0.get_xlim() - - fit_res = [fit_result[0] + fit_result[1] * xi for xi in x] - residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y]) - ax1 = plt.subplot(gs[1]) - ax1.plot(x, residuals, 'ko', ls='none', markersize=5) - ax1.tick_params(direction='out') - ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) - ax1.axhline(y=0.0, ls='--', color='k') - ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k') - ax1.set_xlim(xlim) - ax1.set_ylabel('Residuals') - ax1.set_xlabel(r'$t/a^2$') - - plt.draw() - return -fit_result[0] / fit_result[1] + return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit')) def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):