mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-03-16 23:30:24 +01:00
t0 extractions for new Hadrons module (#171)
* feat: first version of read flow observables for new hadrons module. * feat: refactored t0 fit in seperate function and added extract_t0_hd5 to hadrons submodule.
This commit is contained in:
Fabian Joswig
GitHub
parent
be061fdc23
commit
1184a0fe76
3 changed files with 120 additions and 48 deletions
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue