lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-05-14 11:33:42 +02:00
Code Issues Projects Releases Packages Wiki Activity

Merge 4089238ddd into dcb95265ac

Browse source
This commit is contained in:
Fabian Joswig 2025-04-28 20:18:24 +00:00 • committed by GitHub
commit 67b2bf237f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
12 changed files with 94 additions and 68 deletions
Download patch file Download diff file Expand all files Collapse all files

4
pyerrors/__init__.py
View file

@ -1,4 +1,4 @@
r''' r"""
# What is pyerrors? # What is pyerrors?
`pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data.
It is based on the gamma method [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017). Some of its features are: It is based on the gamma method [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017). Some of its features are:
@ -476,7 +476,7 @@ The array `cdata` contains information about the contribution of auxiliary obser
A JSON schema that may be used to verify the correctness of a file with respect to the format definition is stored in ./examples/json_schema.json. The schema is a self-descriptive format definition and contains an exemplary file. A JSON schema that may be used to verify the correctness of a file with respect to the format definition is stored in ./examples/json_schema.json. The schema is a self-descriptive format definition and contains an exemplary file.
Julia I/O routines for the json.gz format, compatible with [ADerrors.jl](https://gitlab.ift.uam-csic.es/alberto/aderrors.jl), can be found [here](https://github.com/fjosw/ADjson.jl). Julia I/O routines for the json.gz format, compatible with [ADerrors.jl](https://gitlab.ift.uam-csic.es/alberto/aderrors.jl), can be found [here](https://github.com/fjosw/ADjson.jl).
''' """
from .obs import * from .obs import *
from .correlators import * from .correlators import *
from .fits import * from .fits import *

42
pyerrors/correlators.py
View file

@ -42,7 +42,7 @@ class Corr:
__slots__ = ["content", "N", "T", "tag", "prange"] __slots__ = ["content", "N", "T", "tag", "prange"]
def __init__(self, data_input, padding=[0, 0], prange=None): def __init__(self, data_input, padding=None, prange=None):
""" Initialize a Corr object. """ Initialize a Corr object.
Parameters Parameters
@ -58,6 +58,8 @@ class Corr:
region identified for this correlator. region identified for this correlator.
""" """
if padding is None:
padding = [0, 0]
if isinstance(data_input, np.ndarray): if isinstance(data_input, np.ndarray):
if data_input.ndim == 1: if data_input.ndim == 1:
data_input = list(data_input) data_input = list(data_input)
@ -105,7 +107,7 @@ class Corr:
self.N = noNull[0].shape[0] self.N = noNull[0].shape[0]
if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]:
raise ValueError("Smearing matrices are not NxN.") raise ValueError("Smearing matrices are not NxN.")
if (not all([item.shape == noNull[0].shape for item in noNull])): if not all([item.shape == noNull[0].shape for item in noNull]):
raise ValueError("Items in data_input are not of identical shape." + str(noNull)) raise ValueError("Items in data_input are not of identical shape." + str(noNull))
else: else:
raise TypeError("'data_input' contains item of wrong type.") raise TypeError("'data_input' contains item of wrong type.")
@ -236,7 +238,7 @@ class Corr:
newcontent.append(None) newcontent.append(None)
else: else:
newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
if (all([x is None for x in newcontent])): if all([x is None for x in newcontent]):
raise ValueError("Corr could not be symmetrized: No redundant values") raise ValueError("Corr could not be symmetrized: No redundant values")
return Corr(newcontent, prange=self.prange) return Corr(newcontent, prange=self.prange)
@ -300,7 +302,7 @@ class Corr:
return 0.5 * (Corr(transposed) + self) return 0.5 * (Corr(transposed) + self)
def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **kwargs): def GEVP(self, t0, ts=None, sort="Eigenvalue", vector_obs=False, **kwargs):
r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors. r"""Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
@ -333,12 +335,12 @@ class Corr:
Method used to solve the GEVP. Method used to solve the GEVP.
- "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False) - "eigh": Use scipy.linalg.eigh to solve the GEVP. (default for vector_obs=False)
- "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True. - "cholesky": Use manually implemented solution via the Cholesky decomposition. Automatically chosen if vector_obs==True.
''' """
if self.N == 1: if self.N == 1:
raise ValueError("GEVP methods only works on correlator matrices and not single correlators.") raise ValueError("GEVP methods only works on correlator matrices and not single correlators.")
if ts is not None: if ts is not None:
if (ts <= t0): if ts <= t0:
raise ValueError("ts has to be larger than t0.") raise ValueError("ts has to be larger than t0.")
if "sorted_list" in kwargs: if "sorted_list" in kwargs:
@ -786,7 +788,7 @@ class Corr:
raise ValueError('Unknown variant.') raise ValueError('Unknown variant.')
def fit(self, function, fitrange=None, silent=False, **kwargs): def fit(self, function, fitrange=None, silent=False, **kwargs):
r'''Fits function to the data r"""Fits function to the data
Parameters Parameters
---------- ----------
@ -799,7 +801,7 @@ class Corr:
If not specified, self.prange or all timeslices are used. If not specified, self.prange or all timeslices are used.
silent : bool silent : bool
Decides whether output is printed to the standard output. Decides whether output is printed to the standard output.
''' """
if self.N != 1: if self.N != 1:
raise ValueError("Correlator must be projected before fitting") raise ValueError("Correlator must be projected before fitting")
@ -878,6 +880,8 @@ class Corr:
comp : Corr or list of Corr comp : Corr or list of Corr
Correlator or list of correlators which are plotted for comparison. Correlator or list of correlators which are plotted for comparison.
The tags of these correlators are used as labels if available. The tags of these correlators are used as labels if available.
y_range : list
list of two values, determining the range of the y-axis e.g. [0, 12].
logscale : bool logscale : bool
Sets y-axis to logscale. Sets y-axis to logscale.
plateau : Obs plateau : Obs
@ -1093,7 +1097,7 @@ class Corr:
def __add__(self, y): def __add__(self, y):
if isinstance(y, Corr): if isinstance(y, Corr):
if ((self.N != y.N) or (self.T != y.T)): if (self.N != y.N) or (self.T != y.T):
raise ValueError("Addition of Corrs with different shape") raise ValueError("Addition of Corrs with different shape")
newcontent = [] newcontent = []
for t in range(self.T): for t in range(self.T):
@ -1338,21 +1342,21 @@ class Corr:
@property @property
def real(self): def real(self):
def return_real(obs_OR_cobs): def return_real(obs_or_cobs):
if isinstance(obs_OR_cobs.flatten()[0], CObs): if isinstance(obs_or_cobs.flatten()[0], CObs):
return np.vectorize(lambda x: x.real)(obs_OR_cobs) return np.vectorize(lambda x: x.real)(obs_or_cobs)
else: else:
return obs_OR_cobs return obs_or_cobs
return self._apply_func_to_corr(return_real) return self._apply_func_to_corr(return_real)
@property @property
def imag(self): def imag(self):
def return_imag(obs_OR_cobs): def return_imag(obs_or_cobs):
if isinstance(obs_OR_cobs.flatten()[0], CObs): if isinstance(obs_or_cobs.flatten()[0], CObs):
return np.vectorize(lambda x: x.imag)(obs_OR_cobs) return np.vectorize(lambda x: x.imag)(obs_or_cobs)
else: else:
return obs_OR_cobs * 0 # So it stays the right type return obs_or_cobs * 0 # So it stays the right type
return self._apply_func_to_corr(return_imag) return self._apply_func_to_corr(return_imag)
@ -1396,7 +1400,7 @@ class Corr:
if basematrix is None: if basematrix is None:
basematrix = self basematrix = self
if Ntrunc >= basematrix.N: if Ntrunc >= basematrix.N:
raise ValueError('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) raise ValueError('Cannot truncate using Ntrunc <= %d' % basematrix.N)
if basematrix.N != self.N: if basematrix.N != self.N:
raise ValueError('basematrix and targetmatrix have to be of the same size.') raise ValueError('basematrix and targetmatrix have to be of the same size.')
@ -1495,7 +1499,7 @@ def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
def matmul(*operands): def matmul(*operands):
return np.linalg.multi_dot(operands) return np.linalg.multi_dot(operands)
N = Gt.shape[0] N = Gt.shape[0]
output = [[] for j in range(N)] output = [[] for _ in range(N)]
if chol_inv is None: if chol_inv is None:
chol = cholesky(G0) # This will automatically report if the matrix is not pos-def chol = cholesky(G0) # This will automatically report if the matrix is not pos-def
chol_inv = inv(chol) chol_inv = inv(chol)

8
pyerrors/dirac.py
View file

@ -32,8 +32,8 @@ def epsilon_tensor(i, j, k):
elem : int elem : int
Element (i,j,k) of the epsilon tensor of rank 3 Element (i,j,k) of the epsilon tensor of rank 3
""" """
test_set = set((i, j, k)) test_set = {i, j, k}
if not (test_set <= set((1, 2, 3)) or test_set <= set((0, 1, 2))): if not (test_set <= {1, 2, 3} or test_set <= {0, 1, 2}):
raise ValueError("Unexpected input", i, j, k) raise ValueError("Unexpected input", i, j, k)
return (i - j) * (j - k) * (k - i) / 2 return (i - j) * (j - k) * (k - i) / 2
@ -50,8 +50,8 @@ def epsilon_tensor_rank4(i, j, k, o):
elem : int elem : int
Element (i,j,k,o) of the epsilon tensor of rank 4 Element (i,j,k,o) of the epsilon tensor of rank 4
""" """
test_set = set((i, j, k, o)) test_set = {i, j, k, o}
if not (test_set <= set((1, 2, 3, 4)) or test_set <= set((0, 1, 2, 3))): if not (test_set <= {1, 2, 3, 4} or test_set <= {0, 1, 2, 3}):
raise ValueError("Unexpected input", i, j, k, o) raise ValueError("Unexpected input", i, j, k, o)
return (i - j) * (j - k) * (k - i) * (i - o) * (j - o) * (o - k) / 12 return (i - j) * (j - k) * (k - i) * (i - o) * (j - o) * (o - k) / 12

8
pyerrors/fits.py
View file

@ -71,7 +71,7 @@ class Fit_result(Sequence):
def least_squares(x, y, func, priors=None, silent=False, **kwargs): def least_squares(x, y, func, priors=None, silent=False, **kwargs):
r'''Performs a non-linear fit to y = func(x). r"""Performs a non-linear fit to y = func(x).
``` ```
Parameters Parameters
@ -224,7 +224,7 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
chisquare/d.o.f.: 0.5388013574561786 # random chisquare/d.o.f.: 0.5388013574561786 # random
fit parameters [1.11897846 0.96361162 0.92325319] # random fit parameters [1.11897846 0.96361162 0.92325319] # random
''' """
output = Fit_result() output = Fit_result()
if (isinstance(x, dict) and isinstance(y, dict) and isinstance(func, dict)): if (isinstance(x, dict) and isinstance(y, dict) and isinstance(func, dict)):
@ -506,7 +506,7 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
def total_least_squares(x, y, func, silent=False, **kwargs): def total_least_squares(x, y, func, silent=False, **kwargs):
r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. r"""Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
Parameters Parameters
---------- ----------
@ -555,7 +555,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs):
------- -------
output : Fit_result output : Fit_result
Parameters and information on the fitted result. Parameters and information on the fitted result.
''' """
output = Fit_result() output = Fit_result()

4
pyerrors/input/__init__.py
View file

@ -1,10 +1,10 @@
r''' r"""
`pyerrors` includes an `input` submodule in which input routines and parsers for the output of various numerical programs are contained. `pyerrors` includes an `input` submodule in which input routines and parsers for the output of various numerical programs are contained.
# Jackknife samples # Jackknife samples
For comparison with other analysis workflows `pyerrors` can also generate jackknife samples from an `Obs` object or import jackknife samples into an `Obs` object. For comparison with other analysis workflows `pyerrors` can also generate jackknife samples from an `Obs` object or import jackknife samples into an `Obs` object.
See `pyerrors.obs.Obs.export_jackknife` and `pyerrors.obs.import_jackknife` for details. See `pyerrors.obs.Obs.export_jackknife` and `pyerrors.obs.import_jackknife` for details.
''' """
from . import bdio as bdio from . import bdio as bdio
from . import dobs as dobs from . import dobs as dobs
from . import hadrons as hadrons from . import hadrons as hadrons

20
pyerrors/input/dobs.py
View file

@ -85,7 +85,7 @@ def _dict_to_xmlstring_spaces(d, space=' '):
return o return o
def create_pobs_string(obsl, name, spec='', origin='', symbol=[], enstag=None): def create_pobs_string(obsl, name, spec='', origin='', symbol=None, enstag=None):
"""Export a list of Obs or structures containing Obs to an xml string """Export a list of Obs or structures containing Obs to an xml string
according to the Zeuthen pobs format. according to the Zeuthen pobs format.
@ -113,6 +113,8 @@ def create_pobs_string(obsl, name, spec='', origin='', symbol=[], enstag=None):
XML formatted string of the input data XML formatted string of the input data
""" """
if symbol is None:
symbol = []
od = {} od = {}
ename = obsl[0].e_names[0] ename = obsl[0].e_names[0]
names = list(obsl[0].deltas.keys()) names = list(obsl[0].deltas.keys())
@ -176,7 +178,7 @@ def create_pobs_string(obsl, name, spec='', origin='', symbol=[], enstag=None):
return rs return rs
def write_pobs(obsl, fname, name, spec='', origin='', symbol=[], enstag=None, gz=True): def write_pobs(obsl, fname, name, spec='', origin='', symbol=None, enstag=None, gz=True):
"""Export a list of Obs or structures containing Obs to a .xml.gz file """Export a list of Obs or structures containing Obs to a .xml.gz file
according to the Zeuthen pobs format. according to the Zeuthen pobs format.
@ -206,6 +208,8 @@ def write_pobs(obsl, fname, name, spec='', origin='', symbol=[], enstag=None, gz
------- -------
None None
""" """
if symbol is None:
symbol = []
pobsstring = create_pobs_string(obsl, name, spec, origin, symbol, enstag) pobsstring = create_pobs_string(obsl, name, spec, origin, symbol, enstag)
if not fname.endswith('.xml') and not fname.endswith('.gz'): if not fname.endswith('.xml') and not fname.endswith('.gz'):
@ -309,7 +313,7 @@ def read_pobs(fname, full_output=False, gz=True, separator_insertion=None):
full_output : bool full_output : bool
If True, a dict containing auxiliary information and the data is returned. If True, a dict containing auxiliary information and the data is returned.
If False, only the data is returned as list. If False, only the data is returned as list.
separatior_insertion: str or int separator_insertion: str or int
str: replace all occurences of "separator_insertion" within the replica names str: replace all occurences of "separator_insertion" within the replica names
by "|%s" % (separator_insertion) when constructing the names of the replica. by "|%s" % (separator_insertion) when constructing the names of the replica.
int: Insert the separator "|" at the position given by separator_insertion. int: Insert the separator "|" at the position given by separator_insertion.
@ -409,7 +413,7 @@ def import_dobs_string(content, full_output=False, separator_insertion=True):
full_output : bool full_output : bool
If True, a dict containing auxiliary information and the data is returned. If True, a dict containing auxiliary information and the data is returned.
If False, only the data is returned as list. If False, only the data is returned as list.
separatior_insertion: str, int or bool separator_insertion: str, int or bool
str: replace all occurences of "separator_insertion" within the replica names str: replace all occurences of "separator_insertion" within the replica names
by "|%s" % (separator_insertion) when constructing the names of the replica. by "|%s" % (separator_insertion) when constructing the names of the replica.
int: Insert the separator "|" at the position given by separator_insertion. int: Insert the separator "|" at the position given by separator_insertion.
@ -678,7 +682,7 @@ def _dobsdict_to_xmlstring_spaces(d, space=' '):
return o return o
def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags=None): def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=None, who=None, enstags=None):
"""Generate the string for the export of a list of Obs or structures containing Obs """Generate the string for the export of a list of Obs or structures containing Obs
to a .xml.gz file according to the Zeuthen dobs format. to a .xml.gz file according to the Zeuthen dobs format.
@ -709,6 +713,8 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
xml_str : str xml_str : str
XML string generated from the data XML string generated from the data
""" """
if symbol is None:
symbol = []
if enstags is None: if enstags is None:
enstags = {} enstags = {}
od = {} od = {}
@ -867,7 +873,7 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
return rs return rs
def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags=None, gz=True): def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=None, who=None, enstags=None, gz=True):
"""Export a list of Obs or structures containing Obs to a .xml.gz file """Export a list of Obs or structures containing Obs to a .xml.gz file
according to the Zeuthen dobs format. according to the Zeuthen dobs format.
@ -901,6 +907,8 @@ def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=[], who=No
------- -------
None None
""" """
if symbol is None:
symbol = []
if enstags is None: if enstags is None:
enstags = {} enstags = {}

12
pyerrors/input/hadrons.py
View file

@ -245,7 +245,7 @@ def extract_t0_hd5(path, filestem, ens_id, obs='Clover energy density', fit_rang
return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit')) return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))
def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None): def read_DistillationContraction_hd5(path, ens_id, diagrams=None, idl=None):
"""Read hadrons DistillationContraction hdf5 files in given directory structure """Read hadrons DistillationContraction hdf5 files in given directory structure
Parameters Parameters
@ -265,6 +265,8 @@ def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None
extracted DistillationContration data extracted DistillationContration data
""" """
if diagrams is None:
diagrams = ["direct"]
res_dict = {} res_dict = {}
directories, idx = _get_files(path, "data", idl) directories, idx = _get_files(path, "data", idl)
@ -486,7 +488,7 @@ def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
return result_dict return result_dict
def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]): def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=None):
"""Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
Parameters Parameters
@ -508,6 +510,8 @@ def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
extracted fourquark matrizes extracted fourquark matrizes
""" """
if vertices is None:
vertices = ["VA", "AV"]
files, idx = _get_files(path, filestem, idl) files, idx = _get_files(path, filestem, idl)
mom_in = None mom_in = None
@ -596,7 +600,7 @@ def _get_lorentz_names(name):
assert len(name) == 2 assert len(name) == 2
if 'S' in name or 'P' in name: if 'S' in name or 'P' in name:
if not set(name) <= set(['S', 'P']): if not set(name) <= {'S', 'P'}:
raise Exception("'" + name + "' is not a Lorentz scalar") raise Exception("'" + name + "' is not a Lorentz scalar")
g_names = {'S': 'Identity', g_names = {'S': 'Identity',
@ -605,7 +609,7 @@ def _get_lorentz_names(name):
res.append((g_names[name[0]], g_names[name[1]])) res.append((g_names[name[0]], g_names[name[1]]))
else: else:
if not set(name) <= set(['V', 'A']): if not set(name) <= {'V', 'A'}:
raise Exception("'" + name + "' is not a Lorentz scalar") raise Exception("'" + name + "' is not a Lorentz scalar")
for ind in lorentz_index: for ind in lorentz_index:

6
pyerrors/input/openQCD.py
View file

@ -596,7 +596,9 @@ def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):
return arr return arr
def _find_files(path, prefix, postfix, ext, known_files=[]): def _find_files(path, prefix, postfix, ext, known_files=None):
if known_files is None:
known_files = []
found = [] found = []
files = [] files = []
@ -1268,7 +1270,7 @@ def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
idl_wanted = True idl_wanted = True
if 'idl' in kwargs: if 'idl' in kwargs:
idl_wanted = (cnfg in expected_idl[repnum]) idl_wanted = (cnfg in expected_idl[repnum])
left_idl = left_idl - set([cnfg]) left_idl = left_idl - {cnfg}
if idl_wanted: if idl_wanted:
cnfgs[repnum].append(cnfg) cnfgs[repnum].append(cnfg)

14
pyerrors/input/sfcf.py
View file

@ -177,7 +177,7 @@ def read_sfcf_multi(path, prefix, name_list, quarks_list=['.*'], corr_type_list=
# Exclude folders with different names # Exclude folders with different names
for exc in ls: for exc in ls:
if not fnmatch.fnmatch(exc, prefix + '*'): if not fnmatch.fnmatch(exc, prefix + '*'):
ls = list(set(ls) - set([exc])) ls = list(set(ls) - {exc})
if not appended: if not appended:
ls = sort_names(ls) ls = sort_names(ls)
@ -344,7 +344,7 @@ def read_sfcf_multi(path, prefix, name_list, quarks_list=['.*'], corr_type_list=
name_ls = ls name_ls = ls
for exc in name_ls: for exc in name_ls:
if not fnmatch.fnmatch(exc, prefix + '*.' + name): if not fnmatch.fnmatch(exc, prefix + '*.' + name):
name_ls = list(set(name_ls) - set([exc])) name_ls = list(set(name_ls) - {exc})
name_ls = sort_names(name_ls) name_ls = sort_names(name_ls)
pattern = intern[name]['spec'][quarks][off][w][w2]['pattern'] pattern = intern[name]['spec'][quarks][off][w][w2]['pattern']
deltas = [] deltas = []
@ -461,7 +461,9 @@ def _extract_corr_type(corr_type):
return b2b, single return b2b, single
def _find_files(rep_path, prefix, compact, files=[]): def _find_files(rep_path, prefix, compact, files=None):
if files is None:
files = []
sub_ls = [] sub_ls = []
if not files == []: if not files == []:
files.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) files.sort(key=lambda x: int(re.findall(r'\d+', x)[-1]))
@ -475,12 +477,12 @@ def _find_files(rep_path, prefix, compact, files=[]):
if compact: if compact:
for exc in sub_ls: for exc in sub_ls:
if not fnmatch.fnmatch(exc, prefix + '*'): if not fnmatch.fnmatch(exc, prefix + '*'):
sub_ls = list(set(sub_ls) - set([exc])) sub_ls = list(set(sub_ls) - {exc})
sub_ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) sub_ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1]))
else: else:
for exc in sub_ls: for exc in sub_ls:
if not fnmatch.fnmatch(exc, 'cfg*'): if not fnmatch.fnmatch(exc, 'cfg*'):
sub_ls = list(set(sub_ls) - set([exc])) sub_ls = list(set(sub_ls) - {exc})
sub_ls.sort(key=lambda x: int(x[3:])) sub_ls.sort(key=lambda x: int(x[3:]))
files = sub_ls files = sub_ls
if len(files) == 0: if len(files) == 0:
@ -666,7 +668,7 @@ def _get_appended_rep_names(ls, prefix, name, ens_name=None, rep_sep='r'):
new_names = [] new_names = []
for exc in ls: for exc in ls:
if not fnmatch.fnmatch(exc, prefix + '*.' + name): if not fnmatch.fnmatch(exc, prefix + '*.' + name):
ls = list(set(ls) - set([exc])) ls = list(set(ls) - {exc})
ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1]))
for entry in ls: for entry in ls:
myentry = entry[:-len(name) - 1] myentry = entry[:-len(name) - 1]

2
pyerrors/input/utils.py
View file

@ -112,7 +112,7 @@ def check_params(path, param_hash, prefix, param_prefix="parameters_"):
# Exclude folders with different names # Exclude folders with different names
for exc in ls: for exc in ls:
if not fnmatch.fnmatch(exc, prefix + '*'): if not fnmatch.fnmatch(exc, prefix + '*'):
ls = list(set(ls) - set([exc])) ls = list(set(ls) - {exc})
ls = sort_names(ls) ls = sort_names(ls)
nums = {} nums = {}

2
pyerrors/misc.py
View file

@ -174,7 +174,7 @@ def gen_correlated_data(means, cov, name, tau=0.5, samples=1000):
return [Obs([dat], [name]) for dat in corr_data.T] return [Obs([dat], [name]) for dat in corr_data.T]
def _assert_equal_properties(ol, otype=Obs): def _assert_equal_properties(ol):
otype = type(ol[0]) otype = type(ol[0])
for o in ol[1:]: for o in ol[1:]:
if not isinstance(o, otype): if not isinstance(o, otype):

40
pyerrors/obs.py
View file

@ -114,7 +114,7 @@ class Obs:
else: else:
self.idl[name] = list(idx) self.idl[name] = list(idx)
else: else:
raise TypeError('incompatible type for idl[%s].' % (name)) raise TypeError('incompatible type for idl[%s].' % name)
else: else:
for name, sample in sorted(zip(names, samples)): for name, sample in sorted(zip(names, samples)):
self.idl[name] = range(1, len(sample) + 1) self.idl[name] = range(1, len(sample) + 1)
@ -390,7 +390,7 @@ class Obs:
if self.tag is not None: if self.tag is not None:
print("Description:", self.tag) print("Description:", self.tag)
if not hasattr(self, 'e_dvalue'): if not hasattr(self, 'e_dvalue'):
print('Result\t %3.8e' % (self.value)) print('Result\t %3.8e' % self.value)
else: else:
if self.value == 0.0: if self.value == 0.0:
percentage = np.nan percentage = np.nan
@ -448,7 +448,7 @@ class Obs:
my_string_list.append(my_string) my_string_list.append(my_string)
print('\n'.join(my_string_list)) print('\n'.join(my_string_list))
def reweight(self, weight): def reweight(self, weight, all_configs=False):
"""Reweight the obs with given rewighting factors. """Reweight the obs with given rewighting factors.
Parameters Parameters
@ -461,7 +461,7 @@ class Obs:
the reweighting factor on all configurations in weight.idl and not the reweighting factor on all configurations in weight.idl and not
on the configurations in obs[i].idl. Default False. on the configurations in obs[i].idl. Default False.
""" """
return reweight(weight, [self])[0] return reweight(weight, [self], all_configs=all_configs)[0]
def is_zero_within_error(self, sigma=1): def is_zero_within_error(self, sigma=1):
"""Checks whether the observable is zero within 'sigma' standard errors. """Checks whether the observable is zero within 'sigma' standard errors.
@ -1080,7 +1080,7 @@ def _expand_deltas(deltas, idx, shape, gapsize):
are found in idx, the data is expanded. are found in idx, the data is expanded.
""" """
if isinstance(idx, range): if isinstance(idx, range):
if (idx.step == gapsize): if idx.step == gapsize:
return deltas return deltas
ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize) ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize)
for i in range(shape): for i in range(shape):
@ -1190,6 +1190,10 @@ def derived_observable(func, data, array_mode=False, **kwargs):
of func. Use cautiously, supplying the wrong derivative will of func. Use cautiously, supplying the wrong derivative will
not be intercepted. not be intercepted.
array_mode: bool
If True, the function is applied to the full array of data.
Default: False
Notes Notes
----- -----
For simple mathematical operations it can be practical to use anonymous For simple mathematical operations it can be practical to use anonymous
@ -1212,7 +1216,7 @@ def derived_observable(func, data, array_mode=False, **kwargs):
for name in o.cov_names: for name in o.cov_names:
if name in allcov: if name in allcov:
if not np.allclose(allcov[name], o.covobs[name].cov): if not np.allclose(allcov[name], o.covobs[name].cov):
raise Exception('Inconsistent covariance matrices for %s!' % (name)) raise Exception('Inconsistent covariance matrices for %s!' % name)
else: else:
allcov[name] = o.covobs[name].cov allcov[name] = o.covobs[name].cov
@ -1262,7 +1266,7 @@ def derived_observable(func, data, array_mode=False, **kwargs):
for mc_name in obs.mc_names: for mc_name in obs.mc_names:
mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')] mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')] new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d): if 0 < len(mc_idl_d) < len(new_mc_idl_d):
scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d]) scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
return scalef_d return scalef_d
@ -1388,7 +1392,7 @@ def _reduce_deltas(deltas, idx_old, idx_new):
return np.array(deltas)[indices] return np.array(deltas)[indices]
def reweight(weight, obs, **kwargs): def reweight(weight, obs, all_configs=False):
"""Reweight a list of observables. """Reweight a list of observables.
Parameters Parameters
@ -1421,7 +1425,7 @@ def reweight(weight, obs, **kwargs):
new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name])) new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
if kwargs.get('all_configs'): if all_configs:
new_weight = weight new_weight = weight
else: else:
new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
@ -1478,8 +1482,8 @@ def correlate(obs_a, obs_b):
return o return o
def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): def covariance(obs, visualize=False, correlation=False, smooth=None):
r'''Calculates the error covariance matrix of a set of observables. r"""Calculates the error covariance matrix of a set of observables.
WARNING: This function should be used with care, especially for observables with support on multiple WARNING: This function should be used with care, especially for observables with support on multiple
ensembles with differing autocorrelations. See the notes below for details. ensembles with differing autocorrelations. See the notes below for details.
@ -1510,7 +1514,7 @@ def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements. For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
$$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$ $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors). This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
''' """
length = len(obs) length = len(obs)
@ -1564,7 +1568,7 @@ def invert_corr_cov_cholesky(corr, inverrdiag):
if condn > 0.1 / np.finfo(float).eps: if condn > 0.1 / np.finfo(float).eps:
raise ValueError(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})") raise ValueError(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
if condn > 1e13: if condn > 1e13:
warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning) warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % condn, RuntimeWarning)
chol = np.linalg.cholesky(corr) chol = np.linalg.cholesky(corr)
chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True) chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True)
@ -1721,6 +1725,8 @@ def import_jackknife(jacks, name, idl=None):
the N jackknife samples as first to Nth entry. the N jackknife samples as first to Nth entry.
name : str name : str
name of the ensemble the samples are defined on. name of the ensemble the samples are defined on.
idl : list, optional
list of ranges or lists on which the samples are defined
""" """
length = len(jacks) - 1 length = len(jacks) - 1
prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
@ -1800,7 +1806,7 @@ def cov_Obs(means, cov, name, grad=None):
Parameters Parameters
---------- ----------
mean : list of floats or float means : list of floats or float
N mean value(s) of the new Obs N mean value(s) of the new Obs
cov : list or array cov : list or array
2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
@ -1832,7 +1838,7 @@ def cov_Obs(means, cov, name, grad=None):
for i in range(len(means)): for i in range(len(means)):
ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad))) ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
if ol[0].covobs[name].N != len(means): if ol[0].covobs[name].N != len(means):
raise ValueError('You have to provide %d mean values!' % (ol[0].N)) raise ValueError('You have to provide %d mean values!' % ol[0].N)
if len(ol) == 1: if len(ol) == 1:
return ol[0] return ol[0]
return ol return ol
@ -1854,14 +1860,14 @@ def _determine_gap(o, e_content, e_name):
def _check_lists_equal(idl): def _check_lists_equal(idl):
''' """
Use groupby to efficiently check whether all elements of idl are identical. Use groupby to efficiently check whether all elements of idl are identical.
Returns True if all elements are equal, otherwise False. Returns True if all elements are equal, otherwise False.
Parameters Parameters
---------- ----------
idl : list of lists, ranges or np.ndarrays idl : list of lists, ranges or np.ndarrays
''' """
g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl]) g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl])
if next(g, True) and not next(g, False): if next(g, True) and not next(g, False):
return True return True