Rewritten large parts of Qtop routines, making it compatible with sfqcd

This commit is contained in:
Simon Kuberski 2022-02-08 15:44:56 +01:00
parent 65568c84a4
commit aadf9cf32f

View file

@ -217,8 +217,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
return result
def extract_t0(path, prefix, dtr_read, xmin,
spatial_extent, fit_range=5, **kwargs):
def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
"""Extract t0 from given .ms.dat files. Returns t0 as Obs.
It is assumed that all boundary effects have
@ -226,7 +225,6 @@ def extract_t0(path, prefix, dtr_read, xmin,
The data around the zero crossing of t^2<E> - 0.3
is fitted with a linear function
from which the exact root is extracted.
Only works with openQCD
It is assumed that one measurement is performed for each config.
If this is not the case, the resulting idl, as well as the handling
@ -510,8 +508,289 @@ def _read_array_openQCD2(fp):
return {'d': d, 'n': n, 'size': size, 'arr': arr}
def read_qtop(path, prefix, c, dtr_cnfg=1, version="1.2", **kwargs):
"""Read qtop format from given folder structure.
def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
"""Read the topologial charge based on openQCD gradient flow measurements.
Parameters
----------
path : str
path of the measurement files
prefix : str
prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
Ignored if file names are passed explicitly via keyword files.
c : double
Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
dtr_cnfg : int
(optional) parameter that specifies the number of measurements
between two configs.
If it is not set, the distance between two measurements
in the file is assumed to be the distance between two configurations.
steps : int
(optional) Distance between two configurations in units of trajectories /
cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
version : str
Either openQCD or sfqcd, depending on the data.
L : int
spatial length of the lattice in L/a.
HAS to be set if version != sfqcd, since openQCD does not provide
this in the header
r_start : list
list which contains the first config to be read for each replicum.
r_stop : list
list which contains the last config to be read for each replicum.
files : list
specify the exact files that need to be read
from path, practical if e.g. only one replicum is needed
names : list
Alternative labeling for replicas/ensembles.
Has to have the appropriate length.
Zeuthen_flow : bool
(optional) If True, the Zeuthen flow is used for Qtop. Only possible
for version=='sfqcd' If False, the Wilson flow is used.
integer_charge : bool
If True, the charge is rounded towards the nearest integer on each config.
"""
known_versions = ["openQCD", "sfqcd"]
if version not in known_versions:
raise Exception("Unknown openQCD version.")
if "steps" in kwargs:
steps = kwargs.get("steps")
if version == "sfqcd":
if "L" in kwargs:
supposed_L = kwargs.get("L")
else:
supposed_L = None
postfix = ".gfms.dat"
else:
if "L" not in kwargs:
raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
else:
L = kwargs.get("L")
postfix = ".ms.dat"
if "files" in kwargs:
files = kwargs.get("files")
postfix = ''
else:
found = []
files = []
for (dirpath, dirnames, filenames) in os.walk(path + "/"):
# print(filenames)
found.extend(filenames)
break
for f in found:
if fnmatch.fnmatch(f, prefix + "*" + postfix):
files.append(f)
if 'r_start' in kwargs:
r_start = kwargs.get('r_start')
if len(r_start) != len(files):
raise Exception('r_start does not match number of replicas')
r_start = [o if o else None for o in r_start]
else:
r_start = [None] * len(files)
if 'r_stop' in kwargs:
r_stop = kwargs.get('r_stop')
if len(r_stop) != len(files):
raise Exception('r_stop does not match number of replicas')
else:
r_stop = [None] * len(files)
rep_names = []
zeuthen = kwargs.get('Zeuthen_flow', False)
if zeuthen and version not in ['sfqcd']:
raise Exception('Zeuthen flow can only be used for version==sfqcd')
r_start_index = []
r_stop_index = []
deltas = []
configlist = []
for rep, file in enumerate(files):
with open(path + "/" + file, "rb") as fp:
Q = []
traj_list = []
if version in ['sfqcd']:
if zeuthen:
obspos = 0
else:
obspos = 8
t = fp.read(12)
header = struct.unpack('<iii', t)
zthfl = header[0] # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow)
ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
tmax = header[2] # lattice T/a
t = fp.read(12)
Ls = struct.unpack('<iii', t)
if(Ls[0] == Ls[1] and Ls[1] == Ls[2]):
L = Ls[0]
if not (supposed_L == L) and supposed_L:
raise Exception("It seems the length given in the header and by you contradict each other")
else:
raise Exception("Found more than one spatial length in header!")
t = fp.read(16)
header2 = struct.unpack('<dd', t)
tol = header2[0]
cmax = header2[1] # highest value of c used
if c > cmax:
raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
if(zthfl == 2):
nfl = 2 # number of flows
else:
nfl = 1
iobs = 8 * nfl # number of flow observables calculated
while 0 < 1:
t = fp.read(4)
if(len(t) < 4):
break
traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done
for j in range(ncs + 1):
for i in range(iobs):
t = fp.read(8 * tmax)
if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow
Q.append(struct.unpack('d' * tmax, t))
else:
t = fp.read(12)
header = struct.unpack('<iii', t)
# step size in integration steps "dnms"
dn = header[0]
# number of measurements, so "ntot"/dn
nn = header[1]
# lattice T/a
tmax = header[2]
t = fp.read(8)
eps = struct.unpack('d', t)[0]
while 0 < 1:
t = fp.read(4)
if(len(t) < 4):
break
traj_list.append(struct.unpack('i', t)[0])
# Wsl
t = fp.read(8 * tmax * (nn + 1))
# Ysl
t = fp.read(8 * tmax * (nn + 1))
# Qsl, which is asked for in this method
t = fp.read(8 * tmax * (nn + 1))
# unpack the array of Qtops,
# on each timeslice t=0,...,tmax-1 and the
# measurement number in = 0...nn (see README.qcd1)
tmpd = struct.unpack('d' * tmax * (nn + 1), t)
Q.append(tmpd)
if len(np.unique(np.diff(traj_list))) != 1:
raise Exception("Irregularities in stepsize found")
else:
if 'steps' in kwargs:
if steps != traj_list[1] - traj_list[0]:
raise Exception("steps and the found stepsize are not the same")
else:
steps = traj_list[1] - traj_list[0]
configlist.append([tr // steps / dtr_cnfg for tr in traj_list])
if configlist[-1][0] > 1:
offset = configlist[-1][0] - 1
warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
offset, offset * steps))
configlist[-1] = [item - offset for item in configlist[-1]]
if r_start[rep] is None:
r_start_index.append(0)
else:
try:
r_start_index.append(configlist[-1].index(r_start[rep]))
except ValueError:
raise Exception('Config %d not in file with range [%d, %d]' % (
r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
if r_stop[rep] is None:
r_stop_index.append(len(configlist[-1]) - 1)
else:
try:
r_stop_index.append(configlist[-1].index(r_stop[rep]))
except ValueError:
raise Exception('Config %d not in file with range [%d, %d]' % (
r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
if version in ['sfqcd']:
cstepsize = cmax / ncs
index_aim = round(c / cstepsize)
else:
t_aim = (c * L) ** 2 / 8
index_aim = round(t_aim / eps / dn)
Q_sum = []
for i, item in enumerate(Q):
Q_sum.append([sum(item[current:current + tmax])
for current in range(0, len(item), tmax)])
Q_top = []
if version in ['sfqcd']:
for i in range(len(Q_sum) // (ncs + 1)):
Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
else:
for i in range(len(Q) // dtr_cnfg):
Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
if len(Q_top) != len(traj_list) // dtr_cnfg:
raise Exception("qtops and traj_list dont have the same length")
if kwargs.get('integer_charge', False):
Q_top = [round(q) for q in Q_top]
truncated_file = file[:-len(postfix)]
if "names" not in kwargs:
try:
idx = truncated_file.index('r')
except Exception:
if "names" not in kwargs:
raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
ens_name = truncated_file[:idx]
rep_names.append(ens_name + '|' + truncated_file[idx:])
else:
names = kwargs.get("names")
rep_names = names
deltas.append(Q_top)
idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]), 1) for rep in range(len(deltas))]
deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep]] for nrep in range(len(deltas))]
result = Obs(deltas, rep_names, idl=idl)
return result
def qtop_projection(qtop, target=0):
"""Returns the projection to the topological charge sector defined by target.
Parameters
----------
path : Obs
Topological charge, rounded to nearest integer on each config.
target : int
Specifies the topological sector to be reweighted to (default 0)
"""
if qtop.reweighted:
raise Exception('You can not use a reweighted observable for reweighting!')
proj_qtop = []
for n in qtop.deltas:
proj_qtop.append(np.array([1 if int(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
reto.is_merged = qtop.is_merged
return reto
def read_qtop_sector(path, prefix, c, target=0, **kwargs):
"""Constructs reweighting factors to a specified topological sector.
Parameters
----------
@ -521,18 +800,19 @@ def read_qtop(path, prefix, c, dtr_cnfg=1, version="1.2", **kwargs):
prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
c : double
Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
target : int
Specifies the topological sector to be reweighted to (default 0)
dtr_cnfg : int
(optional) parameter that specifies the number of trajectories
between two configs.
if it is not set, the distance between two measurements
in the file is assumed to be
the distance between two configurations.
in the file is assumed to be the distance between two configurations.
steps : int
(optional) (maybe only necessary for openQCD2.0)
nt step size, guessed if not given
(optional) Distance between two configurations in units of trajectories /
cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
version : str
version string of the openQCD (sfqcd) version used to create
the ensemble
the ensemble. Default is 2.0. May also be set to sfqcd.
L : int
spatial length of the lattice in L/a.
HAS to be set if version != sfqcd, since openQCD does not provide
@ -548,200 +828,17 @@ def read_qtop(path, prefix, c, dtr_cnfg=1, version="1.2", **kwargs):
names : list
Alternative labeling for replicas/ensembles.
Has to have the appropriate length
"""
known_versions = ["1.0", "1.2", "1.4", "1.6", "2.0", "sfqcd"]
if version not in known_versions:
raise Exception("Unknown openQCD version.")
if "steps" in kwargs:
steps = kwargs.get("steps")
if version == "sfqcd":
if "L" in kwargs:
supposed_L = kwargs.get("L")
else:
if "L" not in kwargs:
raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
else:
L = kwargs.get("L")
r_start = 1
if "r_start" in kwargs:
r_start = kwargs.get("r_start")
if "r_stop" in kwargs:
r_stop = kwargs.get("r_stop")
if "files" in kwargs:
files = kwargs.get("files")
else:
found = []
files = []
for (dirpath, dirnames, filenames) in os.walk(path + "/"):
# print(filenames)
found.extend(filenames)
break
for f in found:
if fnmatch.fnmatch(f, prefix + "*" + ".ms.dat"):
files.append(f)
print(files)
rep_names = []
deltas = []
idl = []
for rep, file in enumerate(files):
with open(path + "/" + file, "rb") as fp:
t = fp.read(12)
header = struct.unpack('<iii', t)
# step size in integration steps "dnms"
dn = header[0]
# number of measurements, so "ntot"/dn
nn = header[1]
# lattice T/a
tmax = header[2]
if version == "sfqcd":
t = fp.read(12)
Ls = struct.unpack('<iii', t)
if(Ls[0] == Ls[1] and Ls[1] == Ls[2]):
L = Ls[0]
if not (supposed_L == L):
raise Exception("It seems the length given in the header and by you contradict each other")
else:
raise Exception("Found more than one spatial length in header!")
print('dnms:', dn)
print('nn:', nn)
print('tmax:', tmax)
t = fp.read(8)
eps = struct.unpack('d', t)[0]
print('eps:', eps)
Q = []
ncs = []
while 0 < 1:
t = fp.read(4)
if(len(t) < 4):
break
ncs.append(struct.unpack('i', t)[0])
# Wsl
t = fp.read(8 * tmax * (nn + 1))
# Ysl
t = fp.read(8 * tmax * (nn + 1))
# Qsl, which is asked for in this method
t = fp.read(8 * tmax * (nn + 1))
# unpack the array of Qtops,
# on each timeslice t=0,...,tmax-1 and the
# measurement number in = 0...nn (see README.qcd1)
tmpd = struct.unpack('d' * tmax * (nn + 1), t)
Q.append(tmpd)
if not len(set([ncs[i] - ncs[i - 1] for i in range(1, len(ncs))])):
raise Exception("Irregularities in stepsize found")
else:
if 'steps' in kwargs:
if steps != ncs[1] - ncs[0]:
raise Exception("steps and the found stepsize are not the same")
else:
steps = ncs[1] - ncs[0]
print(len(Q))
print('max_t:', dn * (nn) * eps)
t_aim = (c * L) ** 2 / 8
print('t_aim:', t_aim)
index_aim = round(t_aim / eps / dn)
print('index_aim:', index_aim)
Q_sum = []
for i, item in enumerate(Q):
Q_sum.append([sum(item[current:current + tmax])
for current in range(0, len(item), tmax)])
print(len(Q_sum))
print(len(Q_sum[0]))
Q_round = []
for i in range(len(Q) // dtr_cnfg):
Q_round.append(round(Q_sum[dtr_cnfg * i][index_aim]))
if len(Q_round) != len(ncs) // dtr_cnfg:
raise Exception("qtops and ncs dont have the same length")
truncated_file = file[:-7]
print(truncated_file)
idl_start = 1
if "r_start" in kwargs:
Q_round = Q_round[r_start[rep]:]
idl_start = r_start[rep]
if "r_stop" in kwargs:
Q_round = Q_round[:r_stop[rep]]
idl_stop = idl_start + len(Q_round)
# keyword "names" prevails over "ens_name"
if "names" not in kwargs:
try:
idx = truncated_file.index('r')
except Exception:
if "names" not in kwargs:
raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
if "ens_name" in kwargs:
ens_name = kwargs.get("ens_name")
else:
ens_name = truncated_file[:idx]
rep_names.append(ens_name + '|' + truncated_file[idx:])
else:
names = kwargs.get("names")
rep_names = names
deltas.append(np.array(Q_round))
idl.append(range(idl_start, idl_stop))
result = Obs(deltas, rep_names, idl=idl)
return result
def read_qtop_sector(target=0, **kwargs):
"""Constructs reweighting factors to a specified topological sector.
Parameters
----------
target : int
Specifies the topological sector to be reweighted to (default 0)
q_top : Obs
Alternatively takes args of read_qtop method as kwargs
Zeuthen_flow : bool
(optional) If True, the Zeuthen flow is used for Qtop. Only possible
for version=='sfqcd' If False, the Wilson flow is used.
integer_charge : bool
If True, the charge is rounded towards the nearest integer on each config.
"""
if not isinstance(target, int):
raise Exception("'target' has to be an integer.")
if "q_top" in kwargs:
qtop = kwargs.get("q_top")
else:
if "path" in kwargs:
path = kwargs.get("path")
del kwargs["path"]
else:
raise Exception("If you are not providing q_top, please provide path")
if "prefix" in kwargs:
prefix = kwargs.get("prefix")
del kwargs["prefix"]
else:
raise Exception("If you are not providing q_top, please provide prefix")
if "c" in kwargs:
c = kwargs.get("c")
del kwargs["c"]
else:
raise Exception("If you are not providing q_top, please provide c")
if "version" in kwargs:
version = kwargs.get("version")
del kwargs["version"]
else:
version = "1.2"
if "dtr_cnfg" in kwargs:
dtr_cnfg = kwargs.get("dtr_cnfg")
del kwargs["dtr_cnfg"]
else:
dtr_cnfg = 1
qtop = read_qtop(path, prefix, c, dtr_cnfg=dtr_cnfg,
version=version, **kwargs)
names = qtop.names
print(names)
print(qtop.deltas.keys())
proj_qtop = []
for n in qtop.deltas:
proj_qtop.append(np.array([1 if int(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
kwargs['integer_charge'] = True
qtop = read_qtop(path, prefix, c, **kwargs)
result = Obs(proj_qtop, qtop.names)
return result
return qtop_projection(qtop, target=target)