From b14f430f186e38daca8817f598e1f939907e6282 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 09:37:44 +0100 Subject: [PATCH 01/76] prange set to default value None --- pyerrors/correlators.py | 86 ++++++++++++++--------------------------- 1 file changed, 29 insertions(+), 57 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 9a332c6c..5f238328 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -24,7 +24,7 @@ class Corr: """ - def __init__(self, data_input, padding_front=0, padding_back=0,prange=None): + def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): #All data_input should be a list of things at different timeslices. This needs to be verified if not (isinstance(data_input, list)): @@ -57,16 +57,13 @@ class Corr: self.content = [None] * padding_front + self.content + [None] * padding_back self.T = len(self.content) #for convenience: will be used a lot - - #The attribute "range" [start,end] marks a range of two timeslices. + #The attribute "range" [start,end] marks a range of two timeslices. #This is useful for keeping track of plateaus and fitranges. - #The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. - if not prange is None: - self.prange=prange + #The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. + self.prange = prange self.gamma_method() - def gamma_method(self): for item in self.content: if not(item is None): @@ -143,8 +140,7 @@ class Corr: newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) if(all([x is None for x in newcontent])): raise Exception("Corr could not be symmetrized: No redundant values") - - return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None) + return Corr(newcontent, prange=self.prange) def anti_symmetric(self): @@ -159,7 +155,7 @@ class Corr: newcontent.append(0.5 * (self.content[t] - self.content[self.T - t])) if(all([x is None for x in newcontent])): raise Exception("Corr could not be symmetrized: No redundant values") - return Corr(newcontent,prange= self.prange if hasattr(self,"prange") else None) + return Corr(newcontent, prange=self.prange) #This method will symmetrice the matrices and therefore make them positive definit. @@ -192,12 +188,12 @@ class Corr: G0=G.content[t0] L = mat_mat_op(anp.linalg.cholesky, G0) Li = mat_mat_op(anp.linalg.inv, L) - LT=L.T + LT=L.T LTi=mat_mat_op(anp.linalg.inv, LT) newcontent=[] for t in range(self.T): Gt=G.content[t] - M=Li@Gt@LTi + M=Li@Gt@LTi eigenvalues = eigh(M)[0] #print(eigenvalues) eigenvalue=eigenvalues[-state] @@ -205,14 +201,6 @@ class Corr: return Corr(newcontent) - - - - - - - - def roll(self, dt): return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) @@ -307,16 +295,16 @@ class Corr: raise Exception("Correlator must be projected before fitting") #The default behaviour is: - #1 use explicit fitrange + #1 use explicit fitrange # if none is provided, use the range of the corr # if this is also not set, use the whole length of the corr (This could come with a warning!) if fitrange is None: - if hasattr(self,"prange"): - fitrange=self.prange + if self.prange: + fitrange = self.prange else: - fitrange=[0, self.T] + fitrange = [0, self.T] xs = [x for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None] ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None] @@ -329,12 +317,13 @@ class Corr: raise Exception('Unexpected fit result.') return result + #we want to quickly get a plateau def plateau(self, plateau_range=None, method="fit"): if not plateau_range: - if hasattr(self,"prange"): - plateau_range=self.prange - else: + if self.prange: + plateau_range = self.prange + else: raise Exception("no plateau range provided") if self.N != 1: raise Exception("Correlator must be projected before getting a plateau.") @@ -350,11 +339,10 @@ class Corr: return returnvalue else: - raise Exception("Unsupported plateau method: " + method) + raise Exception("Unsupported plateau method: " + method) - - def set_prange(self,prange): + def set_prange(self, prange): if not len(prange)==2: raise Exception("range must be a list or array with two values") if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))): @@ -362,18 +350,10 @@ class Corr: if not (0<=prange[0]<=self.T and 0<=prange[1]<=self.T and prange[0] Date: Mon, 11 Oct 2021 11:44:48 +0100 Subject: [PATCH 02/76] input submodule refactored --- pyerrors/input/__init__.py | 4 +- pyerrors/input/input.py | 733 ------------------------------------- pyerrors/input/misc.py | 126 +++++++ pyerrors/input/openQCD.py | 337 +++++++++++++++++ pyerrors/input/sfcf.py | 278 ++++++++++++++ 5 files changed, 744 insertions(+), 734 deletions(-) delete mode 100644 pyerrors/input/input.py create mode 100644 pyerrors/input/misc.py create mode 100644 pyerrors/input/openQCD.py create mode 100644 pyerrors/input/sfcf.py diff --git a/pyerrors/input/__init__.py b/pyerrors/input/__init__.py index 5c7496e9..24766e77 100644 --- a/pyerrors/input/__init__.py +++ b/pyerrors/input/__init__.py @@ -1,3 +1,5 @@ -from .input import * from . import bdio from . import hadrons +from . import sfcf +from . import openQCD +from . import misc diff --git a/pyerrors/input/input.py b/pyerrors/input/input.py deleted file mode 100644 index 98d57175..00000000 --- a/pyerrors/input/input.py +++ /dev/null @@ -1,733 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 - -import sys -import os -import fnmatch -import re -import struct -import autograd.numpy as np # Thinly-wrapped numpy -from ..pyerrors import Obs -from ..fits import fit_lin - - -def read_sfcf(path, prefix, name, **kwargs): - """Read sfcf C format from given folder structure. - - Keyword arguments - ----------------- - im -- if True, read imaginary instead of real part of the correlation function. - single -- if True, read a boundary-to-boundary correlation function with a single value - b2b -- if True, read a time-dependent boundary-to-boundary correlation function - names -- Alternative labeling for replicas/ensembles. Has to have the appropriate length - """ - if kwargs.get('im'): - im = 1 - part = 'imaginary' - else: - im = 0 - part = 'real' - - if kwargs.get('single'): - b2b = 1 - single = 1 - else: - b2b = 0 - single = 0 - - if kwargs.get('b2b'): - b2b = 1 - - read = 0 - T = 0 - start = 0 - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(dirnames) - break - if not ls: - print('Error, directory not found') - sys.exit() - for exc in ls: - if fnmatch.fnmatch(exc, prefix + '*'): - ls = list(set(ls) - set(exc)) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - replica = len(ls) - print('Read', part, 'part of', name, 'from', prefix, ',', replica, 'replica') - if 'names' in kwargs: - new_names = kwargs.get('names') - if len(new_names) != replica: - raise Exception('Names does not have the required length', replica) - else: - new_names = ls - print(replica, 'replica') - for i, item in enumerate(ls): - print(item) - sub_ls = [] - for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): - sub_ls.extend(dirnames) - break - for exc in sub_ls: - if fnmatch.fnmatch(exc, 'cfg*'): - sub_ls = list(set(sub_ls) - set(exc)) - sub_ls.sort(key=lambda x: int(x[3:])) - no_cfg = len(sub_ls) - print(no_cfg, 'configurations') - - if i == 0: - with open(path + '/' + item + '/' + sub_ls[0] + '/' + name) as fp: - for k, line in enumerate(fp): - if read == 1 and not line.strip() and k > start + 1: - break - if read == 1 and k >= start: - T += 1 - if '[correlator]' in line: - read = 1 - start = k + 7 + b2b - T -= b2b - - deltas = [] - for j in range(T): - deltas.append([]) - - sublength = len(sub_ls) - for j in range(T): - deltas[j].append(np.zeros(sublength)) - - for cnfg, subitem in enumerate(sub_ls): - with open(path + '/' + item + '/' + subitem + '/'+name) as fp: - for k, line in enumerate(fp): - if(k >= start and k < start + T): - floats = list(map(float, line.split())) - deltas[k-start][i][cnfg] = floats[1 + im - single] - - result = [] - for t in range(T): - result.append(Obs(deltas[t], new_names)) - - return result - - -def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwargs): - """Read sfcf c format from given folder structure. - - Arguments - ----------------- - quarks -- Label of the quarks used in the sfcf input file - noffset -- Offset of the source (only relevant when wavefunctions are used) - wf -- ID of wave function - wf2 -- ID of the second wavefunction (only relevant for boundary-to-boundary correlation functions) - - Keyword arguments - ----------------- - im -- if True, read imaginary instead of real part of the correlation function. - b2b -- if True, read a time-dependent boundary-to-boundary correlation function - names -- Alternative labeling for replicas/ensembles. Has to have the appropriate length - """ - - if kwargs.get('im'): - im = 1 - part = 'imaginary' - else: - im = 0 - part = 'real' - - if kwargs.get('b2b'): - b2b = 1 - else: - b2b = 0 - - read = 0 - T = 0 - start = 0 - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(dirnames) - break - if not ls: - print('Error, directory not found') - sys.exit() - # Exclude folders with different names - for exc in ls: - if not fnmatch.fnmatch(exc, prefix+'*'): - ls = list(set(ls) - set([exc])) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. - replica = len(ls) - if 'names' in kwargs: - new_names = kwargs.get('names') - if len(new_names) != replica: - raise Exception('Names does not have the required length', replica) - else: - new_names = ls - print('Read', part, 'part of', name, 'from', prefix[:-1], ',', replica, 'replica') - for i, item in enumerate(ls): - sub_ls = [] - for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): - sub_ls.extend(filenames) - break - for exc in sub_ls: - if not fnmatch.fnmatch(exc, prefix+'*'): - sub_ls = list(set(sub_ls) - set([exc])) - sub_ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) - - first_cfg = int(re.findall(r'\d+', sub_ls[0])[-1]) - - last_cfg = len(sub_ls) + first_cfg - 1 - - for cfg in range(1, len(sub_ls)): - if int(re.findall(r'\d+', sub_ls[cfg])[-1]) != first_cfg + cfg: - last_cfg = cfg + first_cfg - 1 - break - - no_cfg = last_cfg - first_cfg + 1 - print(item, ':', no_cfg, 'evenly spaced configurations (', first_cfg, '-', last_cfg, ') ,', len(sub_ls) - no_cfg, 'configs omitted\n') - - if i == 0: - pattern = 'name ' + name + '\nquarks ' + quarks + '\noffset ' + str(noffset) + '\nwf ' + str(wf) - if b2b: - pattern += '\nwf_2 ' + str(wf2) - - with open(path+'/'+item+'/'+sub_ls[0], 'r') as file: - content = file.read() - match = re.search(pattern, content) - if match: - start_read = content.count('\n', 0, match.start()) + 5 + b2b - end_match = re.search('\n\s*\n', content[match.start():]) - T = content[match.start():].count('\n', 0, end_match.start()) - 4 - b2b - assert T > 0 - print(T, 'entries, starting to read in line', start_read) - else: - raise Exception('Correlator with pattern\n' + pattern + '\nnot found.') - - deltas = [] - for j in range(T): - deltas.append([]) - - sublength = no_cfg - for j in range(T): - deltas[j].append(np.zeros(sublength)) - - for cfg in range(no_cfg): - with open(path+'/'+item+'/'+sub_ls[cfg]) as fp: - for k, line in enumerate(fp): - if k == start_read - 5 - b2b: - if line.strip() != 'name ' + name: - raise Exception('Wrong format', sub_ls[cfg]) - if(k >= start_read and k < start_read + T): - floats = list(map(float, line.split())) - deltas[k-start_read][i][cfg] = floats[-2:][im] - - result = [] - for t in range(T): - result.append(Obs(deltas[t], new_names)) - return result - - -def read_qtop(path, prefix, **kwargs): - """Read qtop format from given folder structure. - - Keyword arguments - ----------------- - target -- specifies the topological sector to be reweighted to (default 0) - full -- if true read the charge instead of the reweighting factor. - """ - - if 'target' in kwargs: - target = kwargs.get('target') - else: - target = 0 - - if kwargs.get('full'): - full = 1 - else: - full = 0 - - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(filenames) - break - - if not ls: - print('Error, directory not found') - sys.exit() - - # Exclude files with different names - for exc in ls: - if not fnmatch.fnmatch(exc, prefix+'*'): - ls = list(set(ls) - set([exc])) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. - replica = len(ls) - print('Read Q_top from', prefix[:-1], ',', replica, 'replica') - - deltas = [] - - for rep in range(replica): - tmp = [] - with open(path+'/'+ls[rep]) as fp: - for k, line in enumerate(fp): - floats = list(map(float, line.split())) - if full == 1: - tmp.append(floats[1]) - else: - if int(floats[1]) == target: - tmp.append(1.0) - else: - tmp.append(0.0) - - deltas.append(np.array(tmp)) - - result = Obs(deltas, [(w.split('.'))[0] for w in ls]) - - return result - -def parse_array_openQCD2(d, n, size, wa, quadrupel=False): - arr = [] - if d == 2: - tot = 0 - for i in range(n[d-1]-1): - if quadrupel: - tmp = wa[tot:n[d-1]] - tmp2 = [] - for i in range(len(tmp)): - if i % 2 == 0: - tmp2.append(tmp[i]) - arr.append(tmp2) - else: - arr.append(np.asarray(wa[tot:n[d-1]])) - return arr - - -# mimic the read_array routine of openQCD-2.0. -# fp is the opened file handle -# returns the dict array -# at this point we only parse a 2d array -# d = 2 -# n = [nfct[irw], 2*nsrc[irw]] -def read_array_openQCD2(fp): - t = fp.read(4) - d = struct.unpack('i', t)[0] - t = fp.read(4*d) - n = struct.unpack('%di' % (d), t) - t = fp.read(4) - size = struct.unpack('i', t)[0] - if size == 4: - types = 'i' - elif size == 8: - types = 'd' - elif size == 16: - types = 'dd' - else: - print('Type not known!') - m = n[0] - for i in range(1,d): - m *= n[i] - - t = fp.read(m*size) - tmp = struct.unpack('%d%s' % (m, types), t) - - arr = parse_array_openQCD2(d, n, size, tmp, quadrupel=True) - return {'d': d, 'n': n, 'size': size, 'arr': arr} - -def read_rwms(path, prefix, names=None, **kwargs): - """Read rwms format from given folder structure. Returns a list of length nrw - - Keyword arguments - ----------------- - new_format -- if True, the array of the associated numbers of Hasenbusch factors is extracted (v>=openQCD1.6) - r_start -- list which contains the first config to be read for each replicum - r_stop -- list which contains the last config to be read for each replicum - oqcd_ver_string -- openQCD version - postfix -- postfix of the file to read, e.g. '.ms1' for openQCD-files - """ - #oqcd versions implemented in this method - known_oqcd_versions = ['1.4','1.6','2.0'] - if 'oqcd_ver_string' in kwargs: - ver_str = kwargs.get('oqcd_ver_string') - if not (ver_str in known_oqcd_versions): - raise Exception('Unknown openQCD version defined!') - else: #Set defaults for openQCD Version to be version 1.4, emulate the old behaviour of this method - # Deprecate this kwarg in version 2.0. - if 'new_format': - ver_str = '1.6' - else: - ver_str = '1.4' - print("Working with openQCD version " + ver_str) - if 'postfix' in kwargs: - postfix = kwargs.get('postfix') - else: - postfix = '' - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(filenames) - break - - if not ls: - print('Error, directory not found') - sys.exit() - - # Exclude files with different names - for exc in ls: - if not fnmatch.fnmatch(exc, prefix + '*'+postfix+'.dat'): - ls = list(set(ls) - set([exc])) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - #ls = fnames - #print(ls) - replica = len(ls) - - if 'r_start' in kwargs: - r_start = kwargs.get('r_start') - if len(r_start) != replica: - raise Exception('r_start does not match number of replicas') - # Adjust Configuration numbering to python index - r_start = [o - 1 if o else None for o in r_start] - else: - r_start = [None] * replica - - if 'r_stop' in kwargs: - r_stop = kwargs.get('r_stop') - if len(r_stop) != replica: - raise Exception('r_stop does not match number of replicas') - else: - r_stop = [None] * replica - - print('Read reweighting factors from', prefix[:-1], ',', replica, 'replica', end='') - - print_err = 0 - if 'print_err' in kwargs: - print_err = 1 - print() - - deltas = [] - - for rep in range(replica): - tmp_array = [] - with open(path+ '/' + ls[rep], 'rb') as fp: - - #header - t = fp.read(4) # number of reweighting factors - if rep == 0: - nrw = struct.unpack('i', t)[0] - if ver_str == '2.0': - nrw = int(nrw/2) - for k in range(nrw): - deltas.append([]) - else: - if ((nrw != struct.unpack('i', t)[0] and (not ver_str == '2.0')) or (nrw != struct.unpack('i', t)[0]/2 and ver_str == '2.0')):# little weird if-clause due to the /2 operation needed. - print('Error: different number of reweighting factors for replicum', rep) - sys.exit() - - for k in range(nrw): - tmp_array.append([]) - - # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files - nfct = [] - if ver_str in ['1.6','2.0']: - for i in range(nrw): - t = fp.read(4) - nfct.append(struct.unpack('i', t)[0]) - #print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting - else: - for i in range(nrw): - nfct.append(1) - - nsrc = [] - for i in range(nrw): - t = fp.read(4) - nsrc.append(struct.unpack('i', t)[0]) - if (ver_str == '2.0'): - if not struct.unpack('i', fp.read(4))[0] == 0: - print('something is wrong!') - - #body - while 0 < 1: - t = fp.read(4) - if len(t) < 4: - break - if print_err: - config_no = struct.unpack('i', t) - for i in range(nrw): - if(ver_str == '2.0'): - tmpd = read_array_openQCD2(fp) - tmpd = read_array_openQCD2(fp) - tmp_rw = tmpd['arr'] - tmp_nfct = 1.0 - for j in range(tmpd['n'][0]): - tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) - if print_err: - print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw[j]))), np.std(np.exp(-np.asarray(tmp_rw[j])))) - print('Sources:', np.exp(-np.asarray(tmp_rw[j]))) - print('Partial factor:', tmp_nfct) - elif(ver_str=='1.6' or ver_str=='1.4'): - tmp_nfct = 1.0 - for j in range(nfct[i]): - t = fp.read(8 * nsrc[i]) - t = fp.read(8 * nsrc[i]) - tmp_rw = struct.unpack('d' * nsrc[i], t) - tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) - if print_err: - print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw))), np.std(np.exp(-np.asarray(tmp_rw)))) - print('Sources:', np.exp(-np.asarray(tmp_rw))) - print('Partial factor:', tmp_nfct) - tmp_array[i].append(tmp_nfct) - - for k in range(nrw): - deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]]) - - print(',', nrw, 'reweighting factors with', nsrc, 'sources') - result = [] - for t in range(nrw): - if names == None: - result.append(Obs(deltas[t], [w.split(".")[0] for w in ls])) - else: - print(names) - result.append(Obs(deltas[t], names)) - return result - - -def read_pbp(path, prefix, **kwargs): - """Read pbp format from given folder structure. Returns a list of length nrw - - Keyword arguments - ----------------- - r_start -- list which contains the first config to be read for each replicum - r_stop -- list which contains the last config to be read for each replicum - - """ - - extract_nfct = 1 - - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(filenames) - break - - if not ls: - print('Error, directory not found') - sys.exit() - - # Exclude files with different names - for exc in ls: - if not fnmatch.fnmatch(exc, prefix + '*.dat'): - ls = list(set(ls) - set([exc])) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - replica = len(ls) - - if 'r_start' in kwargs: - r_start = kwargs.get('r_start') - if len(r_start) != replica: - raise Exception('r_start does not match number of replicas') - # Adjust Configuration numbering to python index - r_start = [o - 1 if o else None for o in r_start] - else: - r_start = [None] * replica - - if 'r_stop' in kwargs: - r_stop = kwargs.get('r_stop') - if len(r_stop) != replica: - raise Exception('r_stop does not match number of replicas') - else: - r_stop = [None] * replica - - print('Read from', prefix[:-1], ',', replica, 'replica', end='') - - print_err = 0 - if 'print_err' in kwargs: - print_err = 1 - print() - - deltas = [] - - for rep in range(replica): - tmp_array = [] - with open(path+ '/' + ls[rep], 'rb') as fp: - - #header - t = fp.read(4) # number of reweighting factors - if rep == 0: - nrw = struct.unpack('i', t)[0] - for k in range(nrw): - deltas.append([]) - else: - if nrw != struct.unpack('i', t)[0]: - print('Error: different number of reweighting factors for replicum', rep) - sys.exit() - - for k in range(nrw): - tmp_array.append([]) - - # This block is necessary for openQCD1.6 ms1 files - nfct = [] - if extract_nfct == 1: - for i in range(nrw): - t = fp.read(4) - nfct.append(struct.unpack('i', t)[0]) - print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting - else: - for i in range(nrw): - nfct.append(1) - - nsrc = [] - for i in range(nrw): - t = fp.read(4) - nsrc.append(struct.unpack('i', t)[0]) - - #body - while 0 < 1: - t = fp.read(4) - if len(t) < 4: - break - if print_err: - config_no = struct.unpack('i', t) - for i in range(nrw): - tmp_nfct = 1.0 - for j in range(nfct[i]): - t = fp.read(8 * nsrc[i]) - t = fp.read(8 * nsrc[i]) - tmp_rw = struct.unpack('d' * nsrc[i], t) - tmp_nfct *= np.mean(np.asarray(tmp_rw)) - if print_err: - print(config_no, i, j, np.mean(np.asarray(tmp_rw)), np.std(np.asarray(tmp_rw))) - print('Sources:', np.asarray(tmp_rw)) - print('Partial factor:', tmp_nfct) - tmp_array[i].append(tmp_nfct) - - for k in range(nrw): - deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]]) - - print(',', nrw, ' with', nsrc, 'sources') - result = [] - for t in range(nrw): - result.append(Obs(deltas[t], [(w.split('.'))[0] for w in ls])) - - return result - - -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 sufficiently decayed at x0=xmin. - The data around the zero crossing of t^2 - 0.3 is fitted with a linear function - from which the exact root is extracted. - Only works with openQCD v 1.2. - - Parameters - ---------- - path -- Path to .ms.dat files - prefix -- Ensemble prefix - dtr_read -- Determines how many trajectories should be skipped when reading the ms.dat files. - Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. - xmin -- First timeslice where the boundary effects have sufficiently decayed. - spatial_extent -- spatial extent of the lattice, required for normalization. - fit_range -- Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5) - - Keyword arguments - ----------------- - r_start -- list which contains the first config to be read for each replicum. - r_stop -- list which contains the last config to be read for each replicum. - plaquette -- If true extract the plaquette estimate of t0 instead. - """ - - ls = [] - for (dirpath, dirnames, filenames) in os.walk(path): - ls.extend(filenames) - break - - if not ls: - print('Error, directory not found') - sys.exit() - - # Exclude files with different names - for exc in ls: - if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'): - ls = list(set(ls) - set([exc])) - if len(ls) > 1: - ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - replica = len(ls) - - if 'r_start' in kwargs: - r_start = kwargs.get('r_start') - if len(r_start) != replica: - raise Exception('r_start does not match number of replicas') - # Adjust Configuration numbering to python index - r_start = [o - 1 if o else None for o in r_start] - else: - r_start = [None] * replica - - if 'r_stop' in kwargs: - r_stop = kwargs.get('r_stop') - if len(r_stop) != replica: - raise Exception('r_stop does not match number of replicas') - else: - r_stop = [None] * replica - - print('Extract t0 from', prefix, ',', replica, 'replica') - - Ysum = [] - - for rep in range(replica): - - with open(path + '/' + ls[rep], 'rb') as fp: - # Read header - t = fp.read(12) - header = struct.unpack('iii', t) - if rep == 0: - dn = header[0] - nn = header[1] - tmax = header[2] - elif dn != header[0] or nn != header[1] or tmax != header[2]: - raise Exception('Replica parameters do not match.') - - t = fp.read(8) - if rep == 0: - eps = struct.unpack('d', t)[0] - print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) - elif eps != struct.unpack('d', t)[0]: - raise Exception('Values for eps do not match among replica.') - - Ysl = [] - - # Read body - while 0 < 1: - t = fp.read(4) - if(len(t) < 4): - break - nc = struct.unpack('i', t)[0] - - t = fp.read(8 * tmax * (nn + 1)) - if kwargs.get('plaquette'): - if nc % dtr_read == 0: - Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) - t = fp.read(8 * tmax * (nn + 1)) - if not kwargs.get('plaquette'): - if nc % dtr_read == 0: - Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) - t = fp.read(8 * tmax * (nn + 1)) - - Ysum.append([]) - for i, item in enumerate(Ysl): - Ysum[-1].append([np.mean(item[current + xmin:current + tmax - xmin]) for current in range(0, len(item), tmax)]) - - t2E_dict = {} - for n in range(nn + 1): - samples = [] - for nrep, rep in enumerate(Ysum): - samples.append([]) - for cnfg in rep: - samples[-1].append(cnfg[n]) - samples[-1] = samples[-1][r_start[nrep]:r_stop[nrep]] - new_obs = Obs(samples, [(w.split('.'))[0] for w in ls]) - 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) - return -fit_result[0] / fit_result[1] diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py new file mode 100644 index 00000000..28f8b1b2 --- /dev/null +++ b/pyerrors/input/misc.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python +# coding: utf-8 + +import os +import fnmatch +import re +import struct +import numpy as np # Thinly-wrapped numpy +from ..pyerrors import Obs + + +def read_pbp(path, prefix, **kwargs): + """Read pbp format from given folder structure. Returns a list of length nrw + + Keyword arguments + ----------------- + r_start -- list which contains the first config to be read for each replicum + r_stop -- list which contains the last config to be read for each replicum + + """ + + extract_nfct = 1 + + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(filenames) + break + + if not ls: + raise Exception('Error, directory not found') + + # Exclude files with different names + for exc in ls: + if not fnmatch.fnmatch(exc, prefix + '*.dat'): + ls = list(set(ls) - set([exc])) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) + replica = len(ls) + + if 'r_start' in kwargs: + r_start = kwargs.get('r_start') + if len(r_start) != replica: + raise Exception('r_start does not match number of replicas') + # Adjust Configuration numbering to python index + r_start = [o - 1 if o else None for o in r_start] + else: + r_start = [None] * replica + + if 'r_stop' in kwargs: + r_stop = kwargs.get('r_stop') + if len(r_stop) != replica: + raise Exception('r_stop does not match number of replicas') + else: + r_stop = [None] * replica + + print('Read from', prefix[:-1], ',', replica, 'replica', end='') + + print_err = 0 + if 'print_err' in kwargs: + print_err = 1 + print() + + deltas = [] + + for rep in range(replica): + tmp_array = [] + with open(path+ '/' + ls[rep], 'rb') as fp: + + #header + t = fp.read(4) # number of reweighting factors + if rep == 0: + nrw = struct.unpack('i', t)[0] + for k in range(nrw): + deltas.append([]) + else: + if nrw != struct.unpack('i', t)[0]: + raise Exception('Error: different number of reweighting factors for replicum', rep) + + for k in range(nrw): + tmp_array.append([]) + + # This block is necessary for openQCD1.6 ms1 files + nfct = [] + if extract_nfct == 1: + for i in range(nrw): + t = fp.read(4) + nfct.append(struct.unpack('i', t)[0]) + print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting + else: + for i in range(nrw): + nfct.append(1) + + nsrc = [] + for i in range(nrw): + t = fp.read(4) + nsrc.append(struct.unpack('i', t)[0]) + + #body + while 0 < 1: + t = fp.read(4) + if len(t) < 4: + break + if print_err: + config_no = struct.unpack('i', t) + for i in range(nrw): + tmp_nfct = 1.0 + for j in range(nfct[i]): + t = fp.read(8 * nsrc[i]) + t = fp.read(8 * nsrc[i]) + tmp_rw = struct.unpack('d' * nsrc[i], t) + tmp_nfct *= np.mean(np.asarray(tmp_rw)) + if print_err: + print(config_no, i, j, np.mean(np.asarray(tmp_rw)), np.std(np.asarray(tmp_rw))) + print('Sources:', np.asarray(tmp_rw)) + print('Partial factor:', tmp_nfct) + tmp_array[i].append(tmp_nfct) + + for k in range(nrw): + deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]]) + + print(',', nrw, ' with', nsrc, 'sources') + result = [] + for t in range(nrw): + result.append(Obs(deltas[t], [(w.split('.'))[0] for w in ls])) + + return result diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py new file mode 100644 index 00000000..046496ee --- /dev/null +++ b/pyerrors/input/openQCD.py @@ -0,0 +1,337 @@ +#!/usr/bin/env python +# coding: utf-8 + +import os +import fnmatch +import re +import struct +import numpy as np # Thinly-wrapped numpy +from ..pyerrors import Obs +from ..fits import fit_lin + + +def read_rwms(path, prefix, version='2.0', names=None, **kwargs): + """Read rwms format from given folder structure. Returns a list of length nrw + + Attributes + ----------------- + version -- version of openQCD, default 2.0 + + Keyword arguments + ----------------- + r_start -- list which contains the first config to be read for each replicum + r_stop -- list which contains the last config to be read for each replicum + postfix -- postfix of the file to read, e.g. '.ms1' for openQCD-files + """ + #oqcd versions implemented in this method + known_oqcd_versions = ['1.4','1.6','2.0'] + if not (version in known_oqcd_versions): + raise Exception('Unknown openQCD version defined!') + else: #Set defaults for openQCD Version to be version 1.4, emulate the old behaviour of this method + # Deprecate this kwarg in version 2.0. + print("Working with openQCD version " + version) + if 'postfix' in kwargs: + postfix = kwargs.get('postfix') + else: + postfix = '' + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(filenames) + break + + if not ls: + raise Exception('Error, directory not found') + + # Exclude files with different names + for exc in ls: + if not fnmatch.fnmatch(exc, prefix + '*'+postfix+'.dat'): + ls = list(set(ls) - set([exc])) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) + #ls = fnames + #print(ls) + replica = len(ls) + + if 'r_start' in kwargs: + r_start = kwargs.get('r_start') + if len(r_start) != replica: + raise Exception('r_start does not match number of replicas') + # Adjust Configuration numbering to python index + r_start = [o - 1 if o else None for o in r_start] + else: + r_start = [None] * replica + + if 'r_stop' in kwargs: + r_stop = kwargs.get('r_stop') + if len(r_stop) != replica: + raise Exception('r_stop does not match number of replicas') + else: + r_stop = [None] * replica + + print('Read reweighting factors from', prefix[:-1], ',', replica, 'replica', end='') + + print_err = 0 + if 'print_err' in kwargs: + print_err = 1 + print() + + deltas = [] + + for rep in range(replica): + tmp_array = [] + with open(path+ '/' + ls[rep], 'rb') as fp: + + #header + t = fp.read(4) # number of reweighting factors + if rep == 0: + nrw = struct.unpack('i', t)[0] + if version == '2.0': + nrw = int(nrw/2) + for k in range(nrw): + deltas.append([]) + else: + if ((nrw != struct.unpack('i', t)[0] and (not verion == '2.0')) or (nrw != struct.unpack('i', t)[0]/2 and version == '2.0')):# little weird if-clause due to the /2 operation needed. + raise Exception('Error: different number of reweighting factors for replicum', rep) + + for k in range(nrw): + tmp_array.append([]) + + # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files + nfct = [] + if version in ['1.6','2.0']: + for i in range(nrw): + t = fp.read(4) + nfct.append(struct.unpack('i', t)[0]) + #print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting + else: + for i in range(nrw): + nfct.append(1) + + nsrc = [] + for i in range(nrw): + t = fp.read(4) + nsrc.append(struct.unpack('i', t)[0]) + if version is '2.0': + if not struct.unpack('i', fp.read(4))[0] == 0: + print('something is wrong!') + + #body + while 0 < 1: + t = fp.read(4) + if len(t) < 4: + break + if print_err: + config_no = struct.unpack('i', t) + for i in range(nrw): + if(version == '2.0'): + tmpd = _read_array_openQCD2(fp) + tmpd = _read_array_openQCD2(fp) + tmp_rw = tmpd['arr'] + tmp_nfct = 1.0 + for j in range(tmpd['n'][0]): + tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j]))) + if print_err: + print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw[j]))), np.std(np.exp(-np.asarray(tmp_rw[j])))) + print('Sources:', np.exp(-np.asarray(tmp_rw[j]))) + print('Partial factor:', tmp_nfct) + elif version is '1.6' or version is '1.4': + tmp_nfct = 1.0 + for j in range(nfct[i]): + t = fp.read(8 * nsrc[i]) + t = fp.read(8 * nsrc[i]) + tmp_rw = struct.unpack('d' * nsrc[i], t) + tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw))) + if print_err: + print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw))), np.std(np.exp(-np.asarray(tmp_rw)))) + print('Sources:', np.exp(-np.asarray(tmp_rw))) + print('Partial factor:', tmp_nfct) + tmp_array[i].append(tmp_nfct) + + for k in range(nrw): + deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]]) + + print(',', nrw, 'reweighting factors with', nsrc, 'sources') + result = [] + for t in range(nrw): + if names == None: + result.append(Obs(deltas[t], [w.split(".")[0] for w in ls])) + else: + print(names) + result.append(Obs(deltas[t], names)) + return result + + + + +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 sufficiently decayed at x0=xmin. + The data around the zero crossing of t^2 - 0.3 is fitted with a linear function + from which the exact root is extracted. + Only works with openQCD v 1.2. + + Parameters + ---------- + path -- Path to .ms.dat files + prefix -- Ensemble prefix + dtr_read -- Determines how many trajectories should be skipped when reading the ms.dat files. + Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. + xmin -- First timeslice where the boundary effects have sufficiently decayed. + spatial_extent -- spatial extent of the lattice, required for normalization. + fit_range -- Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5) + + Keyword arguments + ----------------- + r_start -- list which contains the first config to be read for each replicum. + r_stop -- list which contains the last config to be read for each replicum. + plaquette -- If true extract the plaquette estimate of t0 instead. + """ + + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(filenames) + break + + if not ls: + raise Exception('Error, directory not found') + + # Exclude files with different names + for exc in ls: + if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'): + ls = list(set(ls) - set([exc])) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) + replica = len(ls) + + if 'r_start' in kwargs: + r_start = kwargs.get('r_start') + if len(r_start) != replica: + raise Exception('r_start does not match number of replicas') + # Adjust Configuration numbering to python index + r_start = [o - 1 if o else None for o in r_start] + else: + r_start = [None] * replica + + if 'r_stop' in kwargs: + r_stop = kwargs.get('r_stop') + if len(r_stop) != replica: + raise Exception('r_stop does not match number of replicas') + else: + r_stop = [None] * replica + + print('Extract t0 from', prefix, ',', replica, 'replica') + + Ysum = [] + + for rep in range(replica): + + with open(path + '/' + ls[rep], 'rb') as fp: + # Read header + t = fp.read(12) + header = struct.unpack('iii', t) + if rep == 0: + dn = header[0] + nn = header[1] + tmax = header[2] + elif dn != header[0] or nn != header[1] or tmax != header[2]: + raise Exception('Replica parameters do not match.') + + t = fp.read(8) + if rep == 0: + eps = struct.unpack('d', t)[0] + print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps) + elif eps != struct.unpack('d', t)[0]: + raise Exception('Values for eps do not match among replica.') + + Ysl = [] + + # Read body + while 0 < 1: + t = fp.read(4) + if(len(t) < 4): + break + nc = struct.unpack('i', t)[0] + + t = fp.read(8 * tmax * (nn + 1)) + if kwargs.get('plaquette'): + if nc % dtr_read == 0: + Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + t = fp.read(8 * tmax * (nn + 1)) + if not kwargs.get('plaquette'): + if nc % dtr_read == 0: + Ysl.append(struct.unpack('d' * tmax * (nn + 1), t)) + t = fp.read(8 * tmax * (nn + 1)) + + Ysum.append([]) + for i, item in enumerate(Ysl): + Ysum[-1].append([np.mean(item[current + xmin:current + tmax - xmin]) for current in range(0, len(item), tmax)]) + + t2E_dict = {} + for n in range(nn + 1): + samples = [] + for nrep, rep in enumerate(Ysum): + samples.append([]) + for cnfg in rep: + samples[-1].append(cnfg[n]) + samples[-1] = samples[-1][r_start[nrep]:r_stop[nrep]] + new_obs = Obs(samples, [(w.split('.'))[0] for w in ls]) + 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) + return -fit_result[0] / fit_result[1] + + +def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): + arr = [] + if d == 2: + tot = 0 + for i in range(n[d-1]-1): + if quadrupel: + tmp = wa[tot:n[d-1]] + tmp2 = [] + for i in range(len(tmp)): + if i % 2 == 0: + tmp2.append(tmp[i]) + arr.append(tmp2) + else: + arr.append(np.asarray(wa[tot:n[d-1]])) + return arr + + +# mimic the read_array routine of openQCD-2.0. +# fp is the opened file handle +# returns the dict array +# at this point we only parse a 2d array +# d = 2 +# n = [nfct[irw], 2*nsrc[irw]] +def _read_array_openQCD2(fp): + t = fp.read(4) + d = struct.unpack('i', t)[0] + t = fp.read(4*d) + n = struct.unpack('%di' % (d), t) + t = fp.read(4) + size = struct.unpack('i', t)[0] + if size == 4: + types = 'i' + elif size == 8: + types = 'd' + elif size == 16: + types = 'dd' + else: + print('Type not known!') + m = n[0] + for i in range(1,d): + m *= n[i] + + t = fp.read(m*size) + tmp = struct.unpack('%d%s' % (m, types), t) + + arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True) + return {'d': d, 'n': n, 'size': size, 'arr': arr} diff --git a/pyerrors/input/sfcf.py b/pyerrors/input/sfcf.py new file mode 100644 index 00000000..255f0784 --- /dev/null +++ b/pyerrors/input/sfcf.py @@ -0,0 +1,278 @@ +#!/usr/bin/env python +# coding: utf-8 + +import os +import fnmatch +import re +import numpy as np # Thinly-wrapped numpy +from ..pyerrors import Obs + + +def read_sfcf(path, prefix, name, **kwargs): + """Read sfcf C format from given folder structure. + + Keyword arguments + ----------------- + im -- if True, read imaginary instead of real part of the correlation function. + single -- if True, read a boundary-to-boundary correlation function with a single value + b2b -- if True, read a time-dependent boundary-to-boundary correlation function + names -- Alternative labeling for replicas/ensembles. Has to have the appropriate length + """ + if kwargs.get('im'): + im = 1 + part = 'imaginary' + else: + im = 0 + part = 'real' + + if kwargs.get('single'): + b2b = 1 + single = 1 + else: + b2b = 0 + single = 0 + + if kwargs.get('b2b'): + b2b = 1 + + read = 0 + T = 0 + start = 0 + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(dirnames) + break + if not ls: + raise Exception('Error, directory not found') + for exc in ls: + if fnmatch.fnmatch(exc, prefix + '*'): + ls = list(set(ls) - set(exc)) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) + replica = len(ls) + print('Read', part, 'part of', name, 'from', prefix, ',', replica, 'replica') + if 'names' in kwargs: + new_names = kwargs.get('names') + if len(new_names) != replica: + raise Exception('Names does not have the required length', replica) + else: + new_names = ls + print(replica, 'replica') + for i, item in enumerate(ls): + print(item) + sub_ls = [] + for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): + sub_ls.extend(dirnames) + break + for exc in sub_ls: + if fnmatch.fnmatch(exc, 'cfg*'): + sub_ls = list(set(sub_ls) - set(exc)) + sub_ls.sort(key=lambda x: int(x[3:])) + no_cfg = len(sub_ls) + print(no_cfg, 'configurations') + + if i == 0: + with open(path + '/' + item + '/' + sub_ls[0] + '/' + name) as fp: + for k, line in enumerate(fp): + if read == 1 and not line.strip() and k > start + 1: + break + if read == 1 and k >= start: + T += 1 + if '[correlator]' in line: + read = 1 + start = k + 7 + b2b + T -= b2b + + deltas = [] + for j in range(T): + deltas.append([]) + + sublength = len(sub_ls) + for j in range(T): + deltas[j].append(np.zeros(sublength)) + + for cnfg, subitem in enumerate(sub_ls): + with open(path + '/' + item + '/' + subitem + '/'+name) as fp: + for k, line in enumerate(fp): + if(k >= start and k < start + T): + floats = list(map(float, line.split())) + deltas[k-start][i][cnfg] = floats[1 + im - single] + + result = [] + for t in range(T): + result.append(Obs(deltas[t], new_names)) + + return result + + +def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwargs): + """Read sfcf c format from given folder structure. + + Arguments + ----------------- + quarks -- Label of the quarks used in the sfcf input file + noffset -- Offset of the source (only relevant when wavefunctions are used) + wf -- ID of wave function + wf2 -- ID of the second wavefunction (only relevant for boundary-to-boundary correlation functions) + + Keyword arguments + ----------------- + im -- if True, read imaginary instead of real part of the correlation function. + b2b -- if True, read a time-dependent boundary-to-boundary correlation function + names -- Alternative labeling for replicas/ensembles. Has to have the appropriate length + """ + + if kwargs.get('im'): + im = 1 + part = 'imaginary' + else: + im = 0 + part = 'real' + + if kwargs.get('b2b'): + b2b = 1 + else: + b2b = 0 + + read = 0 + T = 0 + start = 0 + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(dirnames) + break + if not ls: + raise Exception('Error, directory not found') + # Exclude folders with different names + for exc in ls: + if not fnmatch.fnmatch(exc, prefix+'*'): + ls = list(set(ls) - set([exc])) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. + replica = len(ls) + if 'names' in kwargs: + new_names = kwargs.get('names') + if len(new_names) != replica: + raise Exception('Names does not have the required length', replica) + else: + new_names = ls + print('Read', part, 'part of', name, 'from', prefix[:-1], ',', replica, 'replica') + for i, item in enumerate(ls): + sub_ls = [] + for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): + sub_ls.extend(filenames) + break + for exc in sub_ls: + if not fnmatch.fnmatch(exc, prefix+'*'): + sub_ls = list(set(sub_ls) - set([exc])) + sub_ls.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) + + first_cfg = int(re.findall(r'\d+', sub_ls[0])[-1]) + + last_cfg = len(sub_ls) + first_cfg - 1 + + for cfg in range(1, len(sub_ls)): + if int(re.findall(r'\d+', sub_ls[cfg])[-1]) != first_cfg + cfg: + last_cfg = cfg + first_cfg - 1 + break + + no_cfg = last_cfg - first_cfg + 1 + print(item, ':', no_cfg, 'evenly spaced configurations (', first_cfg, '-', last_cfg, ') ,', len(sub_ls) - no_cfg, 'configs omitted\n') + + if i == 0: + pattern = 'name ' + name + '\nquarks ' + quarks + '\noffset ' + str(noffset) + '\nwf ' + str(wf) + if b2b: + pattern += '\nwf_2 ' + str(wf2) + + with open(path+'/'+item+'/'+sub_ls[0], 'r') as file: + content = file.read() + match = re.search(pattern, content) + if match: + start_read = content.count('\n', 0, match.start()) + 5 + b2b + end_match = re.search('\n\s*\n', content[match.start():]) + T = content[match.start():].count('\n', 0, end_match.start()) - 4 - b2b + assert T > 0 + print(T, 'entries, starting to read in line', start_read) + else: + raise Exception('Correlator with pattern\n' + pattern + '\nnot found.') + + deltas = [] + for j in range(T): + deltas.append([]) + + sublength = no_cfg + for j in range(T): + deltas[j].append(np.zeros(sublength)) + + for cfg in range(no_cfg): + with open(path+'/'+item+'/'+sub_ls[cfg]) as fp: + for k, line in enumerate(fp): + if k == start_read - 5 - b2b: + if line.strip() != 'name ' + name: + raise Exception('Wrong format', sub_ls[cfg]) + if(k >= start_read and k < start_read + T): + floats = list(map(float, line.split())) + deltas[k-start_read][i][cfg] = floats[-2:][im] + + result = [] + for t in range(T): + result.append(Obs(deltas[t], new_names)) + return result + + +def read_qtop(path, prefix, **kwargs): + """Read qtop format from given folder structure. + + Keyword arguments + ----------------- + target -- specifies the topological sector to be reweighted to (default 0) + full -- if true read the charge instead of the reweighting factor. + """ + + if 'target' in kwargs: + target = kwargs.get('target') + else: + target = 0 + + if kwargs.get('full'): + full = 1 + else: + full = 0 + + ls = [] + for (dirpath, dirnames, filenames) in os.walk(path): + ls.extend(filenames) + break + + if not ls: + raise Exception('Error, directory not found') + + # Exclude files with different names + for exc in ls: + if not fnmatch.fnmatch(exc, prefix+'*'): + ls = list(set(ls) - set([exc])) + if len(ls) > 1: + ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. + replica = len(ls) + print('Read Q_top from', prefix[:-1], ',', replica, 'replica') + + deltas = [] + + for rep in range(replica): + tmp = [] + with open(path+'/'+ls[rep]) as fp: + for k, line in enumerate(fp): + floats = list(map(float, line.split())) + if full == 1: + tmp.append(floats[1]) + else: + if int(floats[1]) == target: + tmp.append(1.0) + else: + tmp.append(0.0) + + deltas.append(np.array(tmp)) + + result = Obs(deltas, [(w.split('.'))[0] for w in ls]) + + return result From a93944cf27c192a782b9d3eec179d21d3e20d62f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 11:46:00 +0100 Subject: [PATCH 03/76] bug fixed in input/openQCD --- pyerrors/input/openQCD.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 046496ee..f0756360 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -27,8 +27,6 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): known_oqcd_versions = ['1.4','1.6','2.0'] if not (version in known_oqcd_versions): raise Exception('Unknown openQCD version defined!') - else: #Set defaults for openQCD Version to be version 1.4, emulate the old behaviour of this method - # Deprecate this kwarg in version 2.0. print("Working with openQCD version " + version) if 'postfix' in kwargs: postfix = kwargs.get('postfix') From 950cd13c07ae2c694a8d0bcb018c951475d22b2e Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 11:48:17 +0100 Subject: [PATCH 04/76] removed plot_corrs --- pyerrors/pyerrors.py | 118 ------------------------------------------- 1 file changed, 118 deletions(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 45e8b82c..62238d66 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -1106,124 +1106,6 @@ def load_object(path): return pickle.load(file) -def plot_corrs(observables, **kwargs): - """Plot lists of Obs. - - Parameters - ---------- - observables -- list of lists of Obs, where the nth entry is considered to be the - correlation function - at x0=n e.g. [[f_A_0,f_A_1],[f_P_0,f_P_1]] or [f_A,f_P], where f_A and f_P are lists of Obs. - - Keyword arguments - ----------------- - xrange -- list of two values, determining the range of the x-axis e.g. [4, 8] - yrange -- list of two values, determining the range of the y-axis e.g. [0.2, 1.1] - prange -- list of two values, visualizing the width of the plateau e.g. [10, 15] - reference -- float valued variable which is shown as horizontal line for reference - plateau -- Obs which is shown as horizontal line with errorbar for reference - shift -- shift x by given value - label -- list of labels, has to have the same length as observables - exp -- plot exponential from fitting routine - """ - - if 'shift' in kwargs: - shift = kwargs.get('shift') - else: - shift = 0 - - if 'label' in kwargs: - label = kwargs.get('label') - if len(label) != len(observables): - raise Exception('label has to be a list with exactly one entry per entry of observables.') - else: - label = [] - for j in range(len(observables)): - label.append(str(j + 1)) - - - f = plt.figure() - for j in range(len(observables)): - T = len(observables[j]) - - x = np.arange(T) + shift - y = np.zeros(T) - y_err = np.zeros(T) - - for i in range(T): - y[i] = observables[j][i].value - y_err[i] = observables[j][i].dvalue - - plt.errorbar(x, y, yerr=y_err, ls='none', fmt='o', capsize=3, - markersize=5, lw=1, label=label[j]) - - if kwargs.get('logscale'): - plt.yscale('log') - - if 'xrange' in kwargs: - xrange = kwargs.get('xrange') - plt.xlim(xrange[0], xrange[1]) - visible_y = y[int(xrange[0] + 0.5):int(xrange[1] + 0.5)] - visible_y_err = y_err[int(xrange[0] + 0.5):int(xrange[1] + 0.5)] - y_start = np.min(visible_y - visible_y_err) - y_stop = np.max(visible_y + visible_y_err) - span = y_stop - y_start - if np.isfinite(y_start) and np.isfinite(y_stop): - plt.ylim(y_start - 0.1 * span, y_stop + 0.1 * span) - - if 'yrange' in kwargs: - yrange = kwargs.get('yrange') - plt.ylim(yrange[0], yrange[1]) - - if 'reference' in kwargs: - y_value = kwargs.get('reference') - plt.axhline(y=y_value, linewidth=2, color='k', alpha=0.25) - - if 'prange' in kwargs: - prange = kwargs.get('prange') - plt.axvline(x=prange[0] - 0.5, ls='--', c='k', lw=1, alpha=0.5, marker=',') - plt.axvline(x=prange[1] + 0.5, ls='--', c='k', lw=1, alpha=0.5, marker=',') - - if 'plateau' in kwargs: - plateau = kwargs.get('plateau') - if isinstance(plateau, Obs): - plt.axhline(y=plateau.value, linewidth=2, color='k', alpha=0.6, label='Plateau', marker=',', ls='--') - plt.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color='k') - elif isinstance(plateau, list): - for i in range(len(plateau)): - plt.axhline(y=plateau[i].value, linewidth=2, color='C' + str(i), alpha=0.6, label='Plateau' + str(i + 1), marker=',', ls='--') - plt.axhspan(plateau[i].value - plateau[i].dvalue, plateau[i].value + plateau[i].dvalue, - color='C' + str(i), alpha=0.25) - else: - raise Exception('Improper input for plateau.') - - if kwargs.get('exp'): - fit_result = kwargs.get('exp') - y_fit = fit_result[1].value * np.exp(-fit_result[0].value * x) - plt.plot(x, y_fit, color='k') - if not (fit_result[0].e_names == {} and fit_result[1].e_names == {}): - y_fit_err = np.sqrt((y_fit * fit_result[0].dvalue) ** 2 + 2 * covariance(fit_result[0], fit_result[1])* y_fit * - np.exp(-fit_result[0].value * x) + (np.exp(-fit_result[0].value * x) * fit_result[1].dvalue) ** 2) - plt.fill_between(x, y_fit + y_fit_err, y_fit - y_fit_err, color='k', alpha=0.1) - - plt.xlabel('$x_0/a$') - - if 'ylabel' in kwargs: - plt.ylabel(kwargs.get('ylabel')) - - if 'save' in kwargs: - lgd = plt.legend(loc=0) - else: - lgd = plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left') - plt.show() - - if 'save' in kwargs: - save = kwargs.get('save') - if not isinstance(save, str): - raise Exception('save has to be a string.') - f.savefig(save + '.pdf', bbox_extra_artists=(lgd,), bbox_inches='tight') - - def merge_obs(list_of_obs): """Combine all observables in list_of_obs into one new observable From d4eaf3f48aae957358520f68eba571658785beea Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:30:48 +0100 Subject: [PATCH 05/76] flake8 formatting improvements --- pyerrors/input/misc.py | 14 ++++++------- pyerrors/input/openQCD.py | 43 +++++++++++++++++---------------------- pyerrors/input/sfcf.py | 4 +--- 3 files changed, 27 insertions(+), 34 deletions(-) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 28f8b1b2..6dac12d9 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -53,7 +53,7 @@ def read_pbp(path, prefix, **kwargs): else: r_stop = [None] * replica - print('Read from', prefix[:-1], ',', replica, 'replica', end='') + print(r'Read from', prefix[:-1], ',', replica, 'replica', end='') print_err = 0 if 'print_err' in kwargs: @@ -64,10 +64,10 @@ def read_pbp(path, prefix, **kwargs): for rep in range(replica): tmp_array = [] - with open(path+ '/' + ls[rep], 'rb') as fp: + with open(path + '/' + ls[rep], 'rb') as fp: - #header - t = fp.read(4) # number of reweighting factors + # header + t = fp.read(4) # number of reweighting factors if rep == 0: nrw = struct.unpack('i', t)[0] for k in range(nrw): @@ -85,7 +85,7 @@ def read_pbp(path, prefix, **kwargs): for i in range(nrw): t = fp.read(4) nfct.append(struct.unpack('i', t)[0]) - print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting + print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting else: for i in range(nrw): nfct.append(1) @@ -95,7 +95,7 @@ def read_pbp(path, prefix, **kwargs): t = fp.read(4) nsrc.append(struct.unpack('i', t)[0]) - #body + # body while 0 < 1: t = fp.read(4) if len(t) < 4: @@ -118,7 +118,7 @@ def read_pbp(path, prefix, **kwargs): for k in range(nrw): deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]]) - print(',', nrw, ' with', nsrc, 'sources') + print(',', nrw, r' with', nsrc, 'sources') result = [] for t in range(nrw): result.append(Obs(deltas[t], [(w.split('.'))[0] for w in ls])) diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index f0756360..39135510 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -23,8 +23,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): r_stop -- list which contains the last config to be read for each replicum postfix -- postfix of the file to read, e.g. '.ms1' for openQCD-files """ - #oqcd versions implemented in this method - known_oqcd_versions = ['1.4','1.6','2.0'] + known_oqcd_versions = ['1.4', '1.6', '2.0'] if not (version in known_oqcd_versions): raise Exception('Unknown openQCD version defined!') print("Working with openQCD version " + version) @@ -42,12 +41,10 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): # Exclude files with different names for exc in ls: - if not fnmatch.fnmatch(exc, prefix + '*'+postfix+'.dat'): + if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'): ls = list(set(ls) - set([exc])) if len(ls) > 1: ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) - #ls = fnames - #print(ls) replica = len(ls) if 'r_start' in kwargs: @@ -77,18 +74,18 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for rep in range(replica): tmp_array = [] - with open(path+ '/' + ls[rep], 'rb') as fp: + with open(path + '/' + ls[rep], 'rb') as fp: - #header - t = fp.read(4) # number of reweighting factors + # header + t = fp.read(4) # number of reweighting factors if rep == 0: nrw = struct.unpack('i', t)[0] if version == '2.0': - nrw = int(nrw/2) + nrw = int(nrw / 2) for k in range(nrw): deltas.append([]) else: - if ((nrw != struct.unpack('i', t)[0] and (not verion == '2.0')) or (nrw != struct.unpack('i', t)[0]/2 and version == '2.0')):# little weird if-clause due to the /2 operation needed. + if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')): # little weird if-clause due to the /2 operation needed. raise Exception('Error: different number of reweighting factors for replicum', rep) for k in range(nrw): @@ -96,11 +93,11 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files nfct = [] - if version in ['1.6','2.0']: + if version in ['1.6', '2.0']: for i in range(nrw): t = fp.read(4) nfct.append(struct.unpack('i', t)[0]) - #print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting + # print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting else: for i in range(nrw): nfct.append(1) @@ -109,11 +106,11 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): for i in range(nrw): t = fp.read(4) nsrc.append(struct.unpack('i', t)[0]) - if version is '2.0': + if version == '2.0': if not struct.unpack('i', fp.read(4))[0] == 0: print('something is wrong!') - #body + # body while 0 < 1: t = fp.read(4) if len(t) < 4: @@ -132,7 +129,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): print(config_no, i, j, np.mean(np.exp(-np.asarray(tmp_rw[j]))), np.std(np.exp(-np.asarray(tmp_rw[j])))) print('Sources:', np.exp(-np.asarray(tmp_rw[j]))) print('Partial factor:', tmp_nfct) - elif version is '1.6' or version is '1.4': + elif version == '1.6' or version == '1.4': tmp_nfct = 1.0 for j in range(nfct[i]): t = fp.read(8 * nsrc[i]) @@ -151,7 +148,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): print(',', nrw, 'reweighting factors with', nsrc, 'sources') result = [] for t in range(nrw): - if names == None: + if names is None: result.append(Obs(deltas[t], [w.split(".")[0] for w in ls])) else: print(names) @@ -159,8 +156,6 @@ 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): """Extract t0 from given .ms.dat files. Returns t0 as Obs. @@ -290,16 +285,16 @@ def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): arr = [] if d == 2: tot = 0 - for i in range(n[d-1]-1): + for i in range(n[d - 1] - 1): if quadrupel: - tmp = wa[tot:n[d-1]] + tmp = wa[tot:n[d - 1]] tmp2 = [] for i in range(len(tmp)): if i % 2 == 0: tmp2.append(tmp[i]) arr.append(tmp2) else: - arr.append(np.asarray(wa[tot:n[d-1]])) + arr.append(np.asarray(wa[tot:n[d - 1]])) return arr @@ -312,7 +307,7 @@ def _parse_array_openQCD2(d, n, size, wa, quadrupel=False): def _read_array_openQCD2(fp): t = fp.read(4) d = struct.unpack('i', t)[0] - t = fp.read(4*d) + t = fp.read(4 * d) n = struct.unpack('%di' % (d), t) t = fp.read(4) size = struct.unpack('i', t)[0] @@ -325,10 +320,10 @@ def _read_array_openQCD2(fp): else: print('Type not known!') m = n[0] - for i in range(1,d): + for i in range(1, d): m *= n[i] - t = fp.read(m*size) + t = fp.read(m * size) tmp = struct.unpack('%d%s' % (m, types), t) arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True) diff --git a/pyerrors/input/sfcf.py b/pyerrors/input/sfcf.py index 255f0784..b732ecda 100644 --- a/pyerrors/input/sfcf.py +++ b/pyerrors/input/sfcf.py @@ -134,9 +134,7 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg else: b2b = 0 - read = 0 T = 0 - start = 0 ls = [] for (dirpath, dirnames, filenames) in os.walk(path): ls.extend(dirnames) @@ -189,7 +187,7 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg match = re.search(pattern, content) if match: start_read = content.count('\n', 0, match.start()) + 5 + b2b - end_match = re.search('\n\s*\n', content[match.start():]) + end_match = re.search(r'\n\s*\n', content[match.start():]) T = content[match.start():].count('\n', 0, end_match.start()) - 4 - b2b assert T > 0 print(T, 'entries, starting to read in line', start_read) From 3cbc455d106fec42e57ed2aa541fe548501383d8 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:41:56 +0100 Subject: [PATCH 06/76] flake8 ci added --- .github/workflows/flake8.yml | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 .github/workflows/flake8.yml diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml new file mode 100644 index 00000000..91d2770e --- /dev/null +++ b/.github/workflows/flake8.yml @@ -0,0 +1,21 @@ +name: flake8 Lint + +on: [push, pull_request] + +jobs: + flake8-lint: + runs-on: ubuntu-latest + name: Lint + steps: + - name: Check out source repository + uses: actions/checkout@v2 + - name: Set up Python environment + uses: actions/setup-python@v1 + with: + python-version: "3.8" + - name: flake8 Lint + uses: py-actions/flake8@v1 + with: + ignore: "E501" + exclude: "pyerrors/__init__.py" + path: "pyerrors" From ba9ec249443816db4b3283510b8fbd70401c31cd Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:42:51 +0100 Subject: [PATCH 07/76] flake8 badge added to README --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index eb288c4a..38dd0359 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,4 @@ +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) # pyerrors 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: From ff1efd30961bfcb9b89a68b61771cd06e930a5c3 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:51:11 +0100 Subject: [PATCH 08/76] formatting improvements --- .github/workflows/flake8.yml | 1 + pyerrors/input/hadrons.py | 4 +++- pyerrors/input/sfcf.py | 22 +++++++++++----------- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index 91d2770e..d6efe561 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -18,4 +18,5 @@ jobs: with: ignore: "E501" exclude: "pyerrors/__init__.py" + exclude: "pyerrors/input/__init__.py" path: "pyerrors" diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py index c70cc1e2..6254e912 100644 --- a/pyerrors/input/hadrons.py +++ b/pyerrors/input/hadrons.py @@ -32,8 +32,10 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', tree='meson'): if not files: raise Exception('No files starting with', filestem, 'in folder', path) + def get_cnfg_number(n): + return int(n[len(filestem) + 1:-3]) + # Sort according to configuration number - get_cnfg_number = lambda x: int(x[len(filestem) + 1:-3]) files.sort(key=get_cnfg_number) # Check that configurations are evenly spaced diff --git a/pyerrors/input/sfcf.py b/pyerrors/input/sfcf.py index b732ecda..7dd3544e 100644 --- a/pyerrors/input/sfcf.py +++ b/pyerrors/input/sfcf.py @@ -61,7 +61,7 @@ def read_sfcf(path, prefix, name, **kwargs): for i, item in enumerate(ls): print(item) sub_ls = [] - for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): + for (dirpath, dirnames, filenames) in os.walk(path + '/' + item): sub_ls.extend(dirnames) break for exc in sub_ls: @@ -92,11 +92,11 @@ def read_sfcf(path, prefix, name, **kwargs): deltas[j].append(np.zeros(sublength)) for cnfg, subitem in enumerate(sub_ls): - with open(path + '/' + item + '/' + subitem + '/'+name) as fp: + with open(path + '/' + item + '/' + subitem + '/' + name) as fp: for k, line in enumerate(fp): if(k >= start and k < start + T): floats = list(map(float, line.split())) - deltas[k-start][i][cnfg] = floats[1 + im - single] + deltas[k - start][i][cnfg] = floats[1 + im - single] result = [] for t in range(T): @@ -143,7 +143,7 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg raise Exception('Error, directory not found') # Exclude folders with different names for exc in ls: - if not fnmatch.fnmatch(exc, prefix+'*'): + if not fnmatch.fnmatch(exc, prefix + '*'): ls = list(set(ls) - set([exc])) if len(ls) > 1: ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. @@ -157,11 +157,11 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg print('Read', part, 'part of', name, 'from', prefix[:-1], ',', replica, 'replica') for i, item in enumerate(ls): sub_ls = [] - for (dirpath, dirnames, filenames) in os.walk(path+'/'+item): + for (dirpath, dirnames, filenames) in os.walk(path + '/' + item): sub_ls.extend(filenames) break 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.sort(key=lambda x: int(re.findall(r'\d+', x)[-1])) @@ -182,7 +182,7 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg if b2b: pattern += '\nwf_2 ' + str(wf2) - with open(path+'/'+item+'/'+sub_ls[0], 'r') as file: + with open(path + '/' + item + '/' + sub_ls[0], 'r') as file: content = file.read() match = re.search(pattern, content) if match: @@ -203,14 +203,14 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg deltas[j].append(np.zeros(sublength)) for cfg in range(no_cfg): - with open(path+'/'+item+'/'+sub_ls[cfg]) as fp: + with open(path + '/' + item + '/' + sub_ls[cfg]) as fp: for k, line in enumerate(fp): if k == start_read - 5 - b2b: if line.strip() != 'name ' + name: raise Exception('Wrong format', sub_ls[cfg]) if(k >= start_read and k < start_read + T): floats = list(map(float, line.split())) - deltas[k-start_read][i][cfg] = floats[-2:][im] + deltas[k - start_read][i][cfg] = floats[-2:][im] result = [] for t in range(T): @@ -247,7 +247,7 @@ def read_qtop(path, prefix, **kwargs): # Exclude files with different names for exc in ls: - if not fnmatch.fnmatch(exc, prefix+'*'): + if not fnmatch.fnmatch(exc, prefix + '*'): ls = list(set(ls) - set([exc])) if len(ls) > 1: ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0])) # New version, to cope with ids, etc. @@ -258,7 +258,7 @@ def read_qtop(path, prefix, **kwargs): for rep in range(replica): tmp = [] - with open(path+'/'+ls[rep]) as fp: + with open(path + '/' + ls[rep]) as fp: for k, line in enumerate(fp): floats = list(map(float, line.split())) if full == 1: From 50bdc70d5541518e212b4cd2ef9afc87e5687b1d Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:53:23 +0100 Subject: [PATCH 09/76] flake8.yml fixed --- .github/workflows/flake8.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index d6efe561..53710fbf 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,6 +17,5 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: "pyerrors/__init__.py" - exclude: "pyerrors/input/__init__.py" + exclude: ["pyerrors/__init__.py", "pyerrors/input/__init__.py"] path: "pyerrors" From d16eccb6990c34f29f073b7167e134366221bb31 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:55:11 +0100 Subject: [PATCH 10/76] fl fix --- .github/workflows/flake8.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index 53710fbf..59a7e756 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,5 +17,7 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: ["pyerrors/__init__.py", "pyerrors/input/__init__.py"] + exclude: + - "pyerrors/__init__.py" + - "pyerrors/input/__init__.py" path: "pyerrors" From da3bb97c9babdbabd8fee7c5df682b5ec119473b Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:57:16 +0100 Subject: [PATCH 11/76] array yaml --- .github/workflows/flake8.yml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index 59a7e756..364a2b39 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,7 +17,5 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: - - "pyerrors/__init__.py" - - "pyerrors/input/__init__.py" + exclude: '["pyerrors/__init__.py", "pyerrors/input/__init__.py"]' path: "pyerrors" From 1a295befe2a06e7654bfb05ea41bfe990c3cb1ab Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 13:58:40 +0100 Subject: [PATCH 12/76] comma separated list in yaml file --- .github/workflows/flake8.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index 364a2b39..e0936062 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,5 +17,5 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: '["pyerrors/__init__.py", "pyerrors/input/__init__.py"]' + exclude: "pyerrors/__init__.py", "pyerrors/input/__init__.py" path: "pyerrors" From 31882064e38161e7751719dd1c57e4bc8f5e31d3 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 14:00:48 +0100 Subject: [PATCH 13/76] comma separated list --- .github/workflows/flake8.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index e0936062..41f6d9e4 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,5 +17,5 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: "pyerrors/__init__.py", "pyerrors/input/__init__.py" + exclude: "pyerrors/__init__.py, pyerrors/input/__init__.py" path: "pyerrors" From 7cf6f9a85b22a8273ca4d805f4f9f1ceb1ce0fa5 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 11 Oct 2021 14:01:57 +0100 Subject: [PATCH 14/76] clean up flake8 file --- .github/workflows/flake8.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index 41f6d9e4..cf7cd178 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -17,5 +17,5 @@ jobs: uses: py-actions/flake8@v1 with: ignore: "E501" - exclude: "pyerrors/__init__.py, pyerrors/input/__init__.py" + exclude: "__init__.py, input/__init__.py" path: "pyerrors" From 311b0bf91aa7e8abac70314677d5f52a8fa25a9c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 10:13:58 +0100 Subject: [PATCH 15/76] correlator class cleaned up --- pyerrors/correlators.py | 453 +++++++++++++++++----------------------- 1 file changed, 189 insertions(+), 264 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index b49256b0..feb91147 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1,15 +1,12 @@ import numpy as np import autograd.numpy as anp -#from scipy.special.orthogonal import _IntegerType -from .pyerrors import * +import scipy.linalg +from .pyerrors import Obs, dump_object from .fits import standard_fit -from .linalg import * +from .linalg import eigh, mat_mat_op from .roots import find_root -from matplotlib import pyplot as plt -from matplotlib.ticker import NullFormatter # useful for `logit` scale -from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg -#import PySimpleGUI as sg -import matplotlib +import matplotlib.pyplot as plt + class Corr: """The class for a correlator (time dependent sequence of pe.Obs). @@ -25,41 +22,41 @@ class Corr: """ def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): - #All data_input should be a list of things at different timeslices. This needs to be verified + # All data_input should be a list of things at different timeslices. This needs to be verified if not (isinstance(data_input, list)): raise TypeError('Corr__init__ expects a list of timeslices.') # data_input can have multiple shapes. The simplest one is a list of Obs. - #We check, if this is the case + # We check, if this is the case if all([isinstance(item, Obs) for item in data_input]): - self.content=[np.asarray([item]) for item in data_input] - #Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices. - self.N = 1 # number of smearings + self.content = [np.asarray([item]) for item in data_input] + # Wrapping the Obs in an array ensures that the data structure is consistent with smearing matrices. + self.N = 1 # number of smearings - #data_input in the form [np.array(Obs,NxN)] - elif all([isinstance(item,np.ndarray) or item==None for item in data_input]) and any([isinstance(item,np.ndarray)for item in data_input]): + # data_input in the form [np.array(Obs,NxN)] + elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]): self.content = data_input - noNull=[a for a in self.content if not (a is None)] #To check if the matrices are correct for all undefined elements + noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements self.N = noNull[0].shape[0] # The checks are now identical to the case above if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: raise Exception("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 Exception("Items in data_input are not of identical shape." + str(noNull)) - else: # In case its a list of something else. - raise Exception ("data_input contains item of wrong type") + else: # In case its a list of something else. + raise Exception("data_input contains item of wrong type") self.tag = None - #We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value. - #An undefined timeslice is represented by the None object + # We now apply some padding to our list. In case that our list represents a correlator of length T but is not defined at every value. + # An undefined timeslice is represented by the None object self.content = [None] * padding_front + self.content + [None] * padding_back - self.T = len(self.content) #for convenience: will be used a lot + self.T = len(self.content) # for convenience: will be used a lot - #The attribute "range" [start,end] marks a range of two timeslices. - #This is useful for keeping track of plateaus and fitranges. - #The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. + # The attribute "range" [start,end] marks a range of two timeslices. + # This is useful for keeping track of plateaus and fitranges. + # The range can be inherited from other Corrs, if the operation should not alter a chosen range eg. multiplication with a constant. self.prange = prange self.gamma_method() @@ -72,64 +69,64 @@ class Corr: else: for i in range(self.N): for j in range(self.N): - item[i,j].gamma_method() + item[i, j].gamma_method() - #We need to project the Correlator with a Vector to get a single value at each timeslice. - #The method can use one or two vectors. - #If two are specified it returns v1@G@v2 (the order might be very important.) - #By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to + # We need to project the Correlator with a Vector to get a single value at each timeslice. + # The method can use one or two vectors. + # If two are specified it returns v1@G@v2 (the order might be very important.) + # By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to def projected(self, vector_l=None, vector_r=None): if self.N == 1: raise Exception("Trying to project a Corr, that already has N=1.") - #This Exception is in no way necessary. One could just return self - #But there is no scenario, where a user would want that to happen and the error message might be more informative. + # This Exception is in no way necessary. One could just return self + # But there is no scenario, where a user would want that to happen and the error message might be more informative. self.gamma_method() if vector_l is None: - vector_l,vector_r=np.asarray([1.] + (self.N - 1) * [0.]),np.asarray([1.] + (self.N - 1) * [0.]) + vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.]) elif(vector_r is None): - vector_r=vector_l + vector_r = vector_l if not vector_l.shape == vector_r.shape == (self.N,): raise Exception("Vectors are of wrong shape!") - #We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized. - if (not(0.95 < vector_r@vector_r < 1.05)) or (not(0.95 < vector_l@vector_l < 1.05)): + # We always normalize before projecting! But we only raise a warning, when it is clear, they where not meant to be normalized. + if (not (0.95 < vector_r @ vector_r < 1.05)) or (not (0.95 < vector_l @ vector_l < 1.05)): print("Vectors are normalized before projection!") - vector_l,vector_r = vector_l / np.sqrt((vector_l@vector_l)), vector_r / np.sqrt(vector_r@vector_r) + vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) - newcontent = [None if (item is None) else np.asarray([vector_l.T@item@vector_r]) for item in self.content] + newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] return Corr(newcontent) def sum(self): - return np.sqrt(self.N)*self.projected(np.ones(self.N)) - #For purposes of debugging and verification, one might want to see a single smearing level. smearing will return a Corr at the specified i,j. where both are integers 0<=i,j 1: transposed = [None if (G is None) else G.T for G in self.content] - return 0.5 * (Corr(transposed)+self) + return 0.5 * (Corr(transposed) + self) if self.N == 1: raise Exception("Trying to symmetrize a smearing matrix, that already has N=1.") - - #We also include a simple GEVP method based on Scipy.linalg - def GEVP(self, t0, ts,state=1): + # We also include a simple GEVP method based on Scipy.linalg + def GEVP(self, t0, ts, state=1): if (self.content[t0] is None) or (self.content[ts] is None): raise Exception("Corr not defined at t0/ts") G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") @@ -178,39 +172,35 @@ class Corr: G0[i, j] = self.content[t0][i, j].value Gt[i, j] = self.content[ts][i, j].value - sp_val,sp_vec = scipy.linalg.eig(Gt, G0) - sp_vec=sp_vec[:,np.argsort(sp_val)[-state]] #we only want the eigenvector belonging to the selected state - sp_vec = sp_vec/np.sqrt(sp_vec@sp_vec) + sp_val, sp_vec = scipy.linalg.eig(Gt, G0) + sp_vec = sp_vec[:, np.argsort(sp_val)[-state]] # We only want the eigenvector belonging to the selected state + sp_vec = sp_vec / np.sqrt(sp_vec @ sp_vec) return sp_vec - - def Eigenvalue(self,t0,state=1): - G=self.smearing_symmetric() - G0=G.content[t0] + def Eigenvalue(self, t0, state=1): + G = self.smearing_symmetric() + G0 = G.content[t0] L = mat_mat_op(anp.linalg.cholesky, G0) Li = mat_mat_op(anp.linalg.inv, L) - LT=L.T - LTi=mat_mat_op(anp.linalg.inv, LT) - newcontent=[] + LT = L.T + LTi = mat_mat_op(anp.linalg.inv, LT) + newcontent = [] for t in range(self.T): - Gt=G.content[t] - M=Li@Gt@LTi + Gt = G.content[t] + M = Li @ Gt @ LTi eigenvalues = eigh(M)[0] - #print(eigenvalues) - eigenvalue=eigenvalues[-state] + eigenvalue = eigenvalues[-state] newcontent.append(eigenvalue) return Corr(newcontent) - def roll(self, dt): return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) - - def deriv(self, symmetric=True): #Defaults to symmetric derivative + def deriv(self, symmetric=True): # Defaults to symmetric derivative if not symmetric: newcontent = [] for t in range(self.T - 1): - if (self.content[t] is None) or (self.content[t+1] is None): + if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append(self.content[t + 1] - self.content[t]) @@ -219,8 +209,8 @@ class Corr: return Corr(newcontent, padding_back=1) if symmetric: newcontent = [] - for t in range(1, self.T-1): - if (self.content[t-1] is None) or (self.content[t+1] is None): + for t in range(1, self.T - 1): + if (self.content[t - 1] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) @@ -228,11 +218,10 @@ class Corr: raise Exception('Derivative is undefined at all timeslices') return Corr(newcontent, padding_back=1, padding_front=1) - def second_deriv(self): newcontent = [] - for t in range(1, self.T-1): - if (self.content[t-1] is None) or (self.content[t+1] is None): + for t in range(1, self.T - 1): + if (self.content[t - 1] is None) or (self.content[t + 1] is None): newcontent.append(None) else: newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) @@ -240,7 +229,6 @@ class Corr: raise Exception("Derivative is undefined at all timeslices") return Corr(newcontent, padding_back=1, padding_front=1) - def m_eff(self, variant='log', guess=1.0): """Returns the effective mass of the correlator as correlator object @@ -252,7 +240,7 @@ class Corr: """ if self.N != 1: raise Exception('Correlator must be projected before getting m_eff') - if variant is 'log': + if variant == 'log': newcontent = [] for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): @@ -264,42 +252,42 @@ class Corr: return np.log(Corr(newcontent, padding_back=1)) - elif variant is 'periodic': + elif variant == 'periodic': newcontent = [] for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) else: - func = lambda x, d : anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d + func = lambda x, d: anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess))) if(all([x is None for x in newcontent])): raise Exception('m_eff is undefined at all timeslices') return Corr(newcontent, padding_back=1) - elif variant is 'arccosh': - newcontent=[] - for t in range(1,self.T-1): - if (self.content[t] is None) or (self.content[t+1] is None)or (self.content[t-1] is None): + + elif variant == 'arccosh': + newcontent = [] + for t in range(1, self.T - 1): + if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): newcontent.append(None) else: - newcontent.append((self.content[t+1]+self.content[t-1])/(2*self.content[t])) + newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) if(all([x is None for x in newcontent])): raise Exception("m_eff is undefined at all timeslices") - return np.arccosh(Corr(newcontent,padding_back=1,padding_front=1)) + return np.arccosh(Corr(newcontent, padding_back=1, padding_front=1)) else: raise Exception('Unkown variant.') - #We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range. + # We want to apply a pe.standard_fit directly to the Corr using an arbitrary function and range. def fit(self, function, fitrange=None, silent=False, **kwargs): if self.N != 1: raise Exception("Correlator must be projected before fitting") - #The default behaviour is: - #1 use explicit fitrange - # if none is provided, use the range of the corr - # if this is also not set, use the whole length of the corr (This could come with a warning!) - + # The default behaviour is: + # 1 use explicit fitrange + # if none is provided, use the range of the corr + # if this is also not set, use the whole length of the corr (This could come with a warning!) if fitrange is None: if self.prange: @@ -311,15 +299,13 @@ class Corr: ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1]) if not self.content[x] is None] result = standard_fit(xs, ys, function, silent=silent, **kwargs) if isinstance(result, list): - [item.gamma_method() for item in result if isinstance(item,Obs)] + [item.gamma_method() for item in result if isinstance(item, Obs)] elif isinstance(result, dict): - [item.gamma_method() for item in result['fit_parameters'] if isinstance(item,Obs)] + [item.gamma_method() for item in result['fit_parameters'] if isinstance(item, Obs)] else: raise Exception('Unexpected fit result.') return result - - #we want to quickly get a plateau def plateau(self, plateau_range=None, method="fit"): if not plateau_range: if self.prange: @@ -329,13 +315,13 @@ class Corr: if self.N != 1: raise Exception("Correlator must be projected before getting a plateau.") if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1])])): - raise Exception("plateau is undefined at all timeslices in plateaurange.") + raise Exception("plateau is undefined at all timeslices in plateaurange.") if method == "fit": def const_func(a, t): - return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt - return self.fit(const_func,plateau_range)[0] - elif method in ["avg","average","mean"]: - returnvalue= np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1]+1] if not item is None]) + return a[0] # At some point pe.standard fit had an issue with single parameter fits. Being careful does not hurt + return self.fit(const_func, plateau_range)[0] + elif method in ["avg", "average", "mean"]: + returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) returnvalue.gamma_method() return returnvalue @@ -343,11 +329,11 @@ class Corr: raise Exception("Unsupported plateau method: " + method) def set_prange(self, prange): - if not len(prange)==2: + if not len(prange) == 2: raise Exception("prange must be a list or array with two values") - if not ((isinstance(prange[0],int)) and (isinstance(prange[1],int))): + if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): raise Exception("Start and end point must be integers") - if not (0<=prange[0]<=self.T and 0<=prange[1]<=self.T and prange[0] Date: Tue, 12 Oct 2021 10:17:58 +0100 Subject: [PATCH 16/76] indents fixed --- pyerrors/fits.py | 4 ++-- pyerrors/pyerrors.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 47687932..1faf1d3b 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -146,7 +146,7 @@ def standard_fit(x, y, func, silent=False, **kwargs): result_dict['chisquare/expected_chisquare'] = chisquare / expected_chisquare if not silent: print('chisquare/expected_chisquare:', - result_dict['chisquare/expected_chisquare']) + result_dict['chisquare/expected_chisquare']) hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(fit_result.x)) @@ -305,7 +305,7 @@ def odr_fit(x, y, func, silent=False, **kwargs): result_dict['chisquare/expected_chisquare'] = odr_chisquare(np.concatenate((output.beta, output.xplus.ravel()))) / expected_chisquare if not silent: print('chisquare/expected_chisquare:', - result_dict['chisquare/expected_chisquare']) + result_dict['chisquare/expected_chisquare']) hess_inv = np.linalg.pinv(jacobian(jacobian(odr_chisquare))(np.concatenate((output.beta, output.xplus.ravel())))) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index d36a20fa..b7b52f82 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -1064,6 +1064,7 @@ def load_object(path): with open(path, 'rb') as file: return pickle.load(file) + def merge_obs(list_of_obs): """Combine all observables in list_of_obs into one new observable From 117a9257759b2c30d915e6f6017388ded6dc61e0 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 11:24:43 +0100 Subject: [PATCH 17/76] sinh effective mass implemented --- .github/workflows/flake8.yml | 2 +- pyerrors/correlators.py | 13 ++++++++++--- pyerrors/input/bdio.py | 24 +++++++++++------------- pyerrors/linalg.py | 9 ++++++--- 4 files changed, 28 insertions(+), 20 deletions(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index cf7cd178..a96c5633 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -16,6 +16,6 @@ jobs: - name: flake8 Lint uses: py-actions/flake8@v1 with: - ignore: "E501" + ignore: "E501,E722" exclude: "__init__.py, input/__init__.py" path: "pyerrors" diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index feb91147..98e20e6b 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -252,14 +252,21 @@ class Corr: return np.log(Corr(newcontent, padding_back=1)) - elif variant == 'periodic': + elif variant in ['periodic', 'cosh', 'sinh']: + if variant in ['periodic', 'cosh']: + func = anp.cosh + else: + func = anp.sinh + + def root_function(x, d): + return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d + newcontent = [] for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) else: - func = lambda x, d: anp.cosh(x * (t - self.T / 2)) / anp.cosh(x * (t + 1 - self.T / 2)) - d - newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], func, guess=guess))) + newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) if(all([x is None for x in newcontent])): raise Exception('m_eff is undefined at all timeslices') diff --git a/pyerrors/input/bdio.py b/pyerrors/input/bdio.py index 7b522fd9..8d14b440 100644 --- a/pyerrors/input/bdio.py +++ b/pyerrors/input/bdio.py @@ -65,7 +65,7 @@ def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): print('Reading of bdio file started') while 1 > 0: - record = bdio_seek_record(fbdio) + bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) if ruinfo == 7: @@ -75,13 +75,13 @@ def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): if ruinfo < 0: # EOF reached break - rlen = bdio_get_rlen(fbdio) + bdio_get_rlen(fbdio) def read_c_double(): d_buf = ctypes.c_double pd_buf = d_buf() ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf)) - iread = bdio_read_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio)) + bdio_read_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio)) return pd_buf.value mean = read_c_double() @@ -91,7 +91,7 @@ def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): d_buf = ctypes.c_size_t pd_buf = d_buf() ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf)) - iread = bdio_read_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio)) + bdio_read_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio)) return pd_buf.value neid = read_c_size_t() @@ -137,7 +137,7 @@ def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): d_buf = ctypes.c_double * np.sum(ndata) pd_buf = d_buf() ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf)) - iread = bdio_read_f64(ppd_buf, ctypes.c_size_t(8 * np.sum(ndata)), ctypes.c_void_p(fbdio)) + bdio_read_f64(ppd_buf, ctypes.c_size_t(8 * np.sum(ndata)), ctypes.c_void_p(fbdio)) delta = pd_buf[:] samples = np.split(np.asarray(delta) + mean, np.cumsum([a for su in vrep for a in su])[:-1]) @@ -212,8 +212,7 @@ def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs): fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_write), b_form) for obs in obs_list: - - mean = obs.value + # mean = obs.value neid = len(obs.e_names) vrep = [[obs.shape[o] for o in sl] for sl in list(obs.e_content.values())] vrep_write = [item for sublist in vrep for item in sublist] @@ -251,12 +250,12 @@ def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs): def write_c_double(double): pd_buf = ctypes.c_double(double) ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf)) - iwrite = bdio_write_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio)) + bdio_write_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio)) def write_c_size_t(int32): pd_buf = ctypes.c_size_t(int32) ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf)) - iwrite = bdio_write_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio)) + bdio_write_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio)) write_c_double(obs.value) write_c_size_t(neid) @@ -365,7 +364,7 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs): print('Reading of bdio file started') while 1 > 0: - record = bdio_seek_record(fbdio) + bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) if ruinfo < 0: # EOF reached @@ -530,7 +529,7 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): corr_type = [] # Contains correlator data type (important for reading out numerical data) corr_props = [] # Contains propagator types (Component of corr_kappa) d0 = 0 # tvals - d1 = 0 # nnoise + # d1 = 0 # nnoise prop_kappa = [] # Contains propagator kappas (Component of corr_kappa) # Check noise type for multiple replica? cnfg_no = -1 @@ -541,7 +540,7 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): print('Reading of bdio file started') while 1 > 0: - record = bdio_seek_record(fbdio) + bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) if ruinfo < 0: # EOF reached @@ -613,7 +612,6 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): print('Number of configurations: ', cnfg_no + 1) corr_kappa = [] # Contains kappa values for both propagators of given correlation function - corr_source = [] for item in corr_props: corr_kappa.append(float(prop_kappa[int(item)])) diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index f12fdeb8..f8dd6f8b 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -12,11 +12,14 @@ from functools import partial from autograd.extend import defvjp _dot = partial(anp.einsum, '...ij,...jk->...ik') + + # batched diag -_diag = lambda a: anp.eye(a.shape[-1]) * a +def _diag(a): + return anp.eye(a.shape[-1]) * a + + # batched diagonal, similar to matrix_diag in tensorflow - - def _matrix_diag(a): reps = anp.array(a.shape) reps[:-1] = 1 From 7af1768a853c46478ac3ec90940f02c256739991 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 13:57:41 +0100 Subject: [PATCH 18/76] CI set up --- .github/workflows/CI.yml | 29 ++++++ pyerrors/correlators.py | 10 +- tests/test_linalg.py | 56 +++++++++++ tests/test_pyerrors.py | 202 +++------------------------------------ 4 files changed, 102 insertions(+), 195 deletions(-) create mode 100644 .github/workflows/CI.yml create mode 100644 tests/test_linalg.py diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 00000000..83d85491 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,29 @@ +name: CI + +on: [push, pull_request] + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: true + matrix: + os: ["ubuntu-latest", "macos-latest"] + python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] + + steps: + - name: Checkout source + uses: actions/checkout@v2 + + - name: Setup python + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + architecture: x64 + + - name: Install + run: | + pip install -e . + pip install pytest + - name: Run tests + run: pytest diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 98e20e6b..33809267 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1,11 +1,11 @@ import numpy as np import autograd.numpy as anp +import matplotlib.pyplot as plt import scipy.linalg from .pyerrors import Obs, dump_object from .fits import standard_fit from .linalg import eigh, mat_mat_op from .roots import find_root -import matplotlib.pyplot as plt class Corr: @@ -24,7 +24,7 @@ class Corr: def __init__(self, data_input, padding_front=0, padding_back=0, prange=None): # All data_input should be a list of things at different timeslices. This needs to be verified - if not (isinstance(data_input, list)): + if not isinstance(data_input, list): raise TypeError('Corr__init__ expects a list of timeslices.') # data_input can have multiple shapes. The simplest one is a list of Obs. # We check, if this is the case @@ -115,9 +115,9 @@ class Corr: def plottable(self): if self.N != 1: raise Exception("Can only make Corr[N=1] plottable") # We could also autoproject to the groundstate or expect vectors, but this is supposed to be a super simple function. - x_list = [x for x in range(self.T) if (not self.content[x] is None)] - y_list = [y[0].value for y in self.content if (y is not None)] - y_err_list = [y[0].dvalue for y in self.content if (y is not None)] + x_list = [x for x in range(self.T) if not self.content[x] is None] + y_list = [y[0].value for y in self.content if y is not None] + y_err_list = [y[0].dvalue for y in self.content if y is not None] return x_list, y_list, y_err_list diff --git a/tests/test_linalg.py b/tests/test_linalg.py new file mode 100644 index 00000000..9d7f6377 --- /dev/null +++ b/tests/test_linalg.py @@ -0,0 +1,56 @@ +import sys +sys.path.append('..') +import autograd.numpy as np +import os +import random +import math +import string +import copy +import scipy.optimize +from scipy.odr import ODR, Model, Data, RealData +import pyerrors as pe +import pytest + +def test_matrix_functions(): + dim = 3 + int(4 * np.random.rand()) + print(dim) + matrix = [] + for i in range(dim): + row = [] + for j in range(dim): + row.append(pe.pseudo_Obs(np.random.rand(), 0.2 + 0.1 * np.random.rand(), 'e1')) + matrix.append(row) + matrix = np.array(matrix) @ np.identity(dim) + + # Check inverse of matrix + inv = pe.linalg.mat_mat_op(np.linalg.inv, matrix) + check_inv = matrix @ inv + + for (i, j), entry in np.ndenumerate(check_inv): + entry.gamma_method() + if(i == j): + assert math.isclose(entry.value, 1.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) + else: + assert math.isclose(entry.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) + assert math.isclose(entry.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + ' ' + str(entry.dvalue) + + # Check Cholesky decomposition + sym = np.dot(matrix, matrix.T) + cholesky = pe.linalg.mat_mat_op(np.linalg.cholesky, sym) + check = cholesky @ cholesky.T + + for (i, j), entry in np.ndenumerate(check): + diff = entry - sym[i, j] + diff.gamma_method() + assert math.isclose(diff.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + assert math.isclose(diff.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + + # Check eigh + e, v = pe.linalg.eigh(sym) + for i in range(dim): + tmp = sym @ v[:, i] - v[:, i] * e[i] + for j in range(dim): + tmp[j].gamma_method() + assert math.isclose(tmp[j].value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + assert math.isclose(tmp[j].dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index caf36760..0f61e8d2 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -1,18 +1,11 @@ -import sys -sys.path.append('..') import autograd.numpy as np import os import random -import math import string import copy -import scipy.optimize -from scipy.odr import ODR, Model, Data, RealData import pyerrors as pe import pytest -test_iterations = 100 - def test_dump(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) @@ -32,9 +25,9 @@ def test_comparison(): assert (value1 < value2) == (test_obs1 < test_obs2) -def test_man_grad(): - a = pe.pseudo_Obs(17,2.9,'e1') - b = pe.pseudo_Obs(4,0.8,'e1') +def test_function_overloading(): + a = pe.pseudo_Obs(17, 2.9, 'e1') + b = pe.pseudo_Obs(4, 0.8, 'e1') fs = [lambda x: x[0] + x[1], lambda x: x[1] + x[0], lambda x: x[0] - x[1], lambda x: x[1] - x[0], lambda x: x[0] * x[1], lambda x: x[1] * x[0], lambda x: x[0] / x[1], lambda x: x[1] / x[0], @@ -51,8 +44,8 @@ def test_man_grad(): def test_overloading_vectorization(): - a = np.array([5, 4, 8]) - b = pe.pseudo_Obs(4,0.8,'e1') + a = np.random.randint(0, 100, 10) + b = pe.pseudo_Obs(4, 0.8, 'e1') assert [o.value for o in a * b] == [o.value for o in b * a] assert [o.value for o in a + b] == [o.value for o in b + a] @@ -61,8 +54,7 @@ def test_overloading_vectorization(): assert [o.value for o in b / a] == [o.value for o in [b / p for p in a]] -@pytest.mark.parametrize("n", np.arange(test_iterations // 10)) -def test_covariance_is_variance(n): +def test_covariance_is_variance(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) test_obs = pe.pseudo_Obs(value, dvalue, 't') @@ -73,8 +65,7 @@ def test_covariance_is_variance(n): assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float).eps -@pytest.mark.parametrize("n", np.arange(test_iterations // 10)) -def test_fft(n): +def test_fft(): value = np.random.normal(5, 100) dvalue = np.abs(np.random.normal(0, 5)) test_obs1 = pe.pseudo_Obs(value, dvalue, 't', int(500 + 1000 * np.random.rand())) @@ -85,106 +76,7 @@ def test_fft(n): assert np.abs(test_obs1.dvalue - test_obs2.dvalue) <= 10 * max(test_obs1.dvalue, test_obs2.dvalue) * np.finfo(np.float).eps -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_standard_fit(n): - dim = 10 + int(30 * np.random.rand()) - x = np.arange(dim) - y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) - yerr = 0.1 + 0.1 * np.random.rand(dim) - - oy = [] - for i, item in enumerate(x): - oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) - - def f(x, a, b): - return a * np.exp(-b * x) - - popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=[o.dvalue for o in oy], absolute_sigma=True) - - def func(a, x): - y = a[0] * np.exp(-a[1] * x) - return y - - beta = pe.fits.standard_fit(x, oy, func) - - pe.Obs.e_tag_global = 5 - for i in range(2): - beta[i].gamma_method(e_tag=5, S=1.0) - assert math.isclose(beta[i].value, popt[i], abs_tol=1e-5) - assert math.isclose(pcov[i, i], beta[i].dvalue ** 2, abs_tol=1e-3) - assert math.isclose(pe.covariance(beta[0], beta[1]), pcov[0, 1], abs_tol=1e-3) - pe.Obs.e_tag_global = 0 - - chi2_pyerrors = np.sum(((f(x, *[o.value for o in beta]) - y) / yerr) ** 2) / (len(x) - 2) - chi2_scipy = np.sum(((f(x, *popt) - y) / yerr) ** 2) / (len(x) - 2) - assert math.isclose(chi2_pyerrors, chi2_scipy, abs_tol=1e-10) - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_odr_fit(n): - dim = 10 + int(30 * np.random.rand()) - x = np.arange(dim) + np.random.normal(0.0, 0.15, dim) - xerr = 0.1 + 0.1 * np.random.rand(dim) - y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) - yerr = 0.1 + 0.1 * np.random.rand(dim) - - ox = [] - for i, item in enumerate(x): - ox.append(pe.pseudo_Obs(x[i], xerr[i], str(i))) - - oy = [] - for i, item in enumerate(x): - oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) - - def f(x, a, b): - return a * np.exp(-b * x) - - def func(a, x): - y = a[0] * np.exp(-a[1] * x) - return y - - data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy]) - model = Model(func) - odr = ODR(data, model, [0,0], partol=np.finfo(np.float).eps) - odr.set_job(fit_type=0, deriv=1) - output = odr.run() - - beta = pe.fits.odr_fit(ox, oy, func) - - pe.Obs.e_tag_global = 5 - for i in range(2): - beta[i].gamma_method(e_tag=5, S=1.0) - assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5) - assert math.isclose(output.cov_beta[i,i], beta[i].dvalue**2, rel_tol=2.5e-1), str(output.cov_beta[i,i]) + ' ' + str(beta[i].dvalue**2) - assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0,1], rel_tol=2.5e-1) - pe.Obs.e_tag_global = 0 - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_odr_derivatives(n): - x = [] - y = [] - x_err = 0.01 - y_err = 0.01 - - for n in np.arange(1, 9, 2): - loc_xvalue = n + np.random.normal(0.0, x_err) - x.append(pe.pseudo_Obs(loc_xvalue, x_err, str(n))) - y.append(pe.pseudo_Obs((lambda x: x ** 2 - 1)(loc_xvalue) + - np.random.normal(0.0, y_err), y_err, str(n))) - - def func(a, x): - return a[0] + a[1] * x ** 2 - - fit1 = pe.fits.odr_fit(x, y, func) - - tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20) - assert np.abs(np.max(np.array(list(fit1[1].deltas.values())) - - np.array(list(tfit[1].deltas.values())))) < 10e-8 - - -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_covariance_symmetry(n): +def test_covariance_symmetry(): value1 = np.random.normal(5, 10) dvalue1 = np.abs(np.random.normal(0, 1)) test_obs1 = pe.pseudo_Obs(value1, dvalue1, 't') @@ -199,8 +91,7 @@ def test_covariance_symmetry(n): assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float).eps) -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_gamma_method(n): +def test_gamma_method(): # Construct pseudo Obs with random shape value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) @@ -213,28 +104,7 @@ def test_gamma_method(n): assert abs(test_obs.dvalue - dvalue) < 1e-10 * dvalue -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_overloading(n): - # Construct pseudo Obs with random shape - obs_list = [] - for i in range(5): - value = np.abs(np.random.normal(5, 2)) + 2.0 - dvalue = np.abs(np.random.normal(0, 0.1)) + 1e-5 - obs_list.append(pe.pseudo_Obs(value, dvalue, 't', 2000)) - - # Test if the error is processed correctly - def f(x): - return x[0] * x[1] + np.sin(x[2]) * np.exp(x[3] / x[1] / x[0]) - np.sqrt(2) / np.cosh(x[4] / x[0]) - - o_obs = f(obs_list) - d_obs = pe.derived_observable(f, obs_list) - - assert np.max(np.abs((o_obs.deltas['t'] - d_obs.deltas['t']) / o_obs.deltas['t'])) < 1e-7, str(obs_list) - assert np.abs((o_obs.value - d_obs.value) / o_obs.value) < 1e-10 - - -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_derived_observables(n): +def test_derived_observables(): # Construct pseudo Obs with random shape test_obs = pe.pseudo_Obs(2, 0.1 * (1 + np.random.rand()), 't', int(1000 * (1 + np.random.rand()))) @@ -257,8 +127,7 @@ def test_derived_observables(n): assert i_am_one.e_ddvalue['t'] <= 2 * np.finfo(np.float).eps -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_multi_ens_system(n): +def test_multi_ens_system(): names = [] for i in range(100 + int(np.random.rand() * 50)): tmp_string = '' @@ -276,8 +145,7 @@ def test_multi_ens_system(n): assert sorted(x for y in sorted(new_obs.e_content.values()) for x in y) == sorted(new_obs.names) -@pytest.mark.parametrize('n', np.arange(test_iterations)) -def test_overloaded_functions(n): +def test_overloaded_functions(): funcs = [np.exp, np.log, np.sin, np.cos, np.tan, np.sinh, np.cosh, np.arcsinh, np.arccosh] deriv = [np.exp, lambda x: 1 / x, np.cos, lambda x: -np.sin(x), lambda x: 1 / np.cos(x) ** 2, np.cosh, np.sinh, lambda x: 1 / np.sqrt(x ** 2 + 1), lambda x: 1 / np.sqrt(x ** 2 - 1)] val = 3 + 0.5 * np.random.rand() @@ -291,49 +159,3 @@ def test_overloaded_functions(n): assert np.max((ad_obs.deltas['t'] - fd_obs.deltas['t']) / ad_obs.deltas['t']) < 1e-8, item.__name__ assert np.abs((ad_obs.value - item(val)) / ad_obs.value) < 1e-10, item.__name__ assert np.abs(ad_obs.dvalue - dval * np.abs(deriv[i](val))) < 1e-6, item.__name__ - - -@pytest.mark.parametrize('n', np.arange(test_iterations // 10)) -def test_matrix_functions(n): - dim = 3 + int(4 * np.random.rand()) - print(dim) - matrix = [] - for i in range(dim): - row = [] - for j in range(dim): - row.append(pe.pseudo_Obs(np.random.rand(), 0.2 + 0.1 * np.random.rand(), 'e1')) - matrix.append(row) - matrix = np.array(matrix) @ np.identity(dim) - - # Check inverse of matrix - inv = pe.linalg.mat_mat_op(np.linalg.inv, matrix) - check_inv = matrix @ inv - - for (i, j), entry in np.ndenumerate(check_inv): - entry.gamma_method() - if(i == j): - assert math.isclose(entry.value, 1.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) - else: - assert math.isclose(entry.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) + ' ' + str(entry.value) - assert math.isclose(entry.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) + ' ' + str(entry.dvalue) - - # Check Cholesky decomposition - sym = np.dot(matrix, matrix.T) - cholesky = pe.linalg.mat_mat_op(np.linalg.cholesky, sym) - check = cholesky @ cholesky.T - - for (i, j), entry in np.ndenumerate(check): - diff = entry - sym[i, j] - diff.gamma_method() - assert math.isclose(diff.value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) - assert math.isclose(diff.dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) - - # Check eigh - e, v = pe.linalg.eigh(sym) - for i in range(dim): - tmp = sym @ v[:, i] - v[:, i] * e[i] - for j in range(dim): - tmp[j].gamma_method() - assert math.isclose(tmp[j].value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) - assert math.isclose(tmp[j].dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) - From 68580325d06a2649a67a20e78dfc5d59f0cbd8bb Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:06:36 +0100 Subject: [PATCH 19/76] CI extended --- .github/workflows/CI.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 83d85491..a0b0c751 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -9,7 +9,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] + python-version: ["3.5", "3.6", "3.8"] steps: - name: Checkout source @@ -21,9 +21,14 @@ jobs: python-version: ${{ matrix.python-version }} architecture: x64 + - name: Display Python version + run: python -c "import sys; print(sys.version)" + - name: Install run: | + python -m pip install --upgrade pip pip install -e . pip install pytest + - name: Run tests run: pytest From ef35dd868418d0e94ee14adfe73609b047db63fa Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:11:01 +0100 Subject: [PATCH 20/76] np.float changed to np.float64 --- tests/test_pyerrors.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 0f61e8d2..393cde8d 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -59,10 +59,10 @@ def test_covariance_is_variance(): dvalue = np.abs(np.random.normal(0, 1)) test_obs = pe.pseudo_Obs(value, dvalue, 't') test_obs.gamma_method() - assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float).eps + assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float64).eps test_obs = test_obs + pe.pseudo_Obs(value, dvalue, 'q', 200) test_obs.gamma_method(e_tag=0) - assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float).eps + assert np.abs(test_obs.dvalue ** 2 - pe.covariance(test_obs, test_obs)) <= 10 * np.finfo(np.float64).eps def test_fft(): @@ -72,8 +72,8 @@ def test_fft(): test_obs2 = copy.deepcopy(test_obs1) test_obs1.gamma_method() test_obs2.gamma_method(fft=False) - assert max(np.abs(test_obs1.e_rho[''] - test_obs2.e_rho[''])) <= 10 * np.finfo(np.float).eps - assert np.abs(test_obs1.dvalue - test_obs2.dvalue) <= 10 * max(test_obs1.dvalue, test_obs2.dvalue) * np.finfo(np.float).eps + assert max(np.abs(test_obs1.e_rho[''] - test_obs2.e_rho[''])) <= 10 * np.finfo(np.float64).eps + assert np.abs(test_obs1.dvalue - test_obs2.dvalue) <= 10 * max(test_obs1.dvalue, test_obs2.dvalue) * np.finfo(np.float64).eps def test_covariance_symmetry(): @@ -87,8 +87,8 @@ def test_covariance_symmetry(): test_obs2.gamma_method() cov_ab = pe.covariance(test_obs1, test_obs2) cov_ba = pe.covariance(test_obs2, test_obs1) - assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float).eps - assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float).eps) + assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps + assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps) def test_gamma_method(): @@ -115,16 +115,16 @@ def test_derived_observables(): d_Obs_fd.gamma_method() assert d_Obs_ad.value == d_Obs_fd.value - assert np.abs(4.0 * np.sin(4.0) - d_Obs_ad.value) < 1000 * np.finfo(np.float).eps * np.abs(d_Obs_ad.value) - assert np.abs(d_Obs_ad.dvalue-d_Obs_fd.dvalue) < 1000 * np.finfo(np.float).eps * d_Obs_ad.dvalue + assert np.abs(4.0 * np.sin(4.0) - d_Obs_ad.value) < 1000 * np.finfo(np.float64).eps * np.abs(d_Obs_ad.value) + assert np.abs(d_Obs_ad.dvalue-d_Obs_fd.dvalue) < 1000 * np.finfo(np.float64).eps * d_Obs_ad.dvalue i_am_one = pe.derived_observable(lambda x, **kwargs: x[0] / x[1], [d_Obs_ad, d_Obs_ad]) i_am_one.gamma_method(e_tag=1) assert i_am_one.value == 1.0 - assert i_am_one.dvalue < 2 * np.finfo(np.float).eps - assert i_am_one.e_dvalue['t'] <= 2 * np.finfo(np.float).eps - assert i_am_one.e_ddvalue['t'] <= 2 * np.finfo(np.float).eps + assert i_am_one.dvalue < 2 * np.finfo(np.float64).eps + assert i_am_one.e_dvalue['t'] <= 2 * np.finfo(np.float64).eps + assert i_am_one.e_ddvalue['t'] <= 2 * np.finfo(np.float64).eps def test_multi_ens_system(): From cddf6ddf6b60a69db6b663ecbd6128683cdc8a6f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:12:21 +0100 Subject: [PATCH 21/76] further instances of np.float removed --- pyerrors/fits.py | 6 +++--- pyerrors/misc.py | 2 +- pyerrors/pyerrors.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 1faf1d3b..137ede2f 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -256,7 +256,7 @@ def odr_fit(x, y, func, silent=False, **kwargs): data = RealData(x_f, y_f, sx=dx_f, sy=dy_f) model = Model(func) - odr = ODR(data, model, x0, partol=np.finfo(np.float).eps) + odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps) odr.set_job(fit_type=0, deriv=1) output = odr.run() @@ -610,7 +610,7 @@ def covariance_matrix(y): def error_band(x, func, beta): """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.""" cov = covariance_matrix(beta) - if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float).eps): + if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): print('Warning, Covariance matrix is not symmetric within floating point precision') print('cov - cov.T:') print(cov - cov.T) @@ -716,7 +716,7 @@ def fit_general(x, y, func, silent=False, **kwargs): model = Model(func) - odr = ODR(data, model, beta0, partol=np.finfo(np.float).eps) + odr = ODR(data, model, beta0, partol=np.finfo(np.float64).eps) odr.set_job(fit_type=fit_type, deriv=1) output = odr.run() if print_output and not silent: diff --git a/pyerrors/misc.py b/pyerrors/misc.py index b393e49f..a4d71a5f 100644 --- a/pyerrors/misc.py +++ b/pyerrors/misc.py @@ -71,7 +71,7 @@ def ks_test(obs=None): plt.ylabel('Cumulative probability') plt.title(str(bins) + ' Q values') - n = np.arange(1, bins + 1) / np.float(bins) + n = np.arange(1, bins + 1) / np.float64(bins) Xs = np.sort(Qs) plt.step(Xs, n) diffs = n - Xs diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index b7b52f82..947be666 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -297,7 +297,7 @@ class Obs: self.e_windowsize[e_name] = n break - if len(self.e_content[e_name]) > 1 and self.e_dvalue[e_name] > np.finfo(np.float).eps: + if len(self.e_content[e_name]) > 1 and self.e_dvalue[e_name] > np.finfo(np.float64).eps: e_mean = 0 for r_name in self.e_content[e_name]: e_mean += self.shape[r_name] * self.r_values[r_name] From 7e0534eab55d9ec411597abbf85c7e284018226f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:20:47 +0100 Subject: [PATCH 22/76] numpy and matpltolib version numbers added to CI --- .github/workflows/CI.yml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index a0b0c751..991c0201 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,8 +21,11 @@ jobs: python-version: ${{ matrix.python-version }} architecture: x64 - - name: Display Python version - run: python -c "import sys; print(sys.version)" + - name: Display Python version numbers + run: | + python -c "import sys; print(sys.version)" + python -c "import numpy; print(numpy.__version__)" + python -c "import matplotlib; print(matplotlib.__version__)" - name: Install run: | From 83fc7c46293f8f2696bfef93a1d5bb7231d9ad89 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:22:01 +0100 Subject: [PATCH 23/76] numpy and mpl version numbers moved --- .github/workflows/CI.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 991c0201..6e628c86 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,17 +21,17 @@ jobs: python-version: ${{ matrix.python-version }} architecture: x64 - - name: Display Python version numbers + - name: Display Python version number run: | python -c "import sys; print(sys.version)" - python -c "import numpy; print(numpy.__version__)" - python -c "import matplotlib; print(matplotlib.__version__)" - name: Install run: | python -m pip install --upgrade pip pip install -e . pip install pytest + python -c "import numpy; print(numpy.__version__)" + python -c "import matplotlib; print(matplotlib.__version__)" - name: Run tests run: pytest From 0efd8a3c6d1bd10b24e9802b4b12d7ab1cd88b29 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:28:32 +0100 Subject: [PATCH 24/76] CI python versions changed --- .github/workflows/CI.yml | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 6e628c86..4e0e9b3b 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -9,7 +9,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.5", "3.6", "3.8"] + python-version: ["3.6", "3.8", "3.10"] steps: - name: Checkout source diff --git a/README.md b/README.md index 325c8358..5bed276d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) # pyerrors 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: From 62c51998650748fb4136b1faa13ec1a19404dfa0 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:45:20 +0100 Subject: [PATCH 25/76] CI targets changed --- .github/workflows/CI.yml | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4e0e9b3b..fd65c11c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,6 +1,11 @@ name: CI -on: [push, pull_request] +on: + push: + branches: + - master + - develop + pull_request: jobs: test: @@ -9,7 +14,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.6", "3.8", "3.10"] + python-version: ["3.6", "3.8", "3.9"] steps: - name: Checkout source @@ -21,17 +26,11 @@ jobs: python-version: ${{ matrix.python-version }} architecture: x64 - - name: Display Python version number - run: | - python -c "import sys; print(sys.version)" - - name: Install run: | python -m pip install --upgrade pip pip install -e . pip install pytest - python -c "import numpy; print(numpy.__version__)" - python -c "import matplotlib; print(matplotlib.__version__)" - name: Run tests run: pytest From 0098e74aa96bceff8a9c8e42aa01969c6ef60849 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:48:37 +0100 Subject: [PATCH 26/76] target of flake8 Lint changed --- .github/workflows/flake8.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml index a96c5633..0766352d 100644 --- a/.github/workflows/flake8.yml +++ b/.github/workflows/flake8.yml @@ -1,6 +1,11 @@ name: flake8 Lint -on: [push, pull_request] +on: + push: + branches: + - master + - develop + pull_request: jobs: flake8-lint: From 7b607ede2aa8c7050709b2e49df37068ff308862 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 12 Oct 2021 14:57:00 +0100 Subject: [PATCH 27/76] macos removed from CI --- .github/workflows/CI.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index fd65c11c..05445774 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -9,12 +9,11 @@ on: jobs: test: - runs-on: ${{ matrix.os }} + runs-on: ubuntu-latest strategy: fail-fast: true matrix: - os: ["ubuntu-latest", "macos-latest"] - python-version: ["3.6", "3.8", "3.9"] + python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] steps: - name: Checkout source From cac85f53f8a85d1ebe438b2985744fc466f1211f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 13 Oct 2021 09:43:40 +0100 Subject: [PATCH 28/76] CI output more verbose --- .github/workflows/CI.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 05445774..7c4bdf41 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -8,7 +8,7 @@ on: pull_request: jobs: - test: + pytest: runs-on: ubuntu-latest strategy: fail-fast: true @@ -23,7 +23,6 @@ jobs: uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} - architecture: x64 - name: Install run: | @@ -32,4 +31,4 @@ jobs: pip install pytest - name: Run tests - run: pytest + run: pytest -v From 68c5e68149518475230339076b0e5e6604d94f84 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 13 Oct 2021 09:46:05 +0100 Subject: [PATCH 29/76] coverage report added to CI --- .github/workflows/CI.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7c4bdf41..d0608ae5 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,6 +29,7 @@ jobs: python -m pip install --upgrade pip pip install -e . pip install pytest + pip install pytest-cov - name: Run tests - run: pytest -v + run: pytest --cov=pyerrors -v From a5f9583229546ecf55524c4e20682be2907479ce Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 13 Oct 2021 09:47:47 +0100 Subject: [PATCH 30/76] gitignore extended --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index b35fbaaa..18c31bfa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,8 @@ __pycache__ *.pyc .ipynb_* +.coverage +.cache examples/B1k2_pcac_plateau.p examples/Untitled.* core.* From 1ad4f798e4e946ba6ba43c2788182ca04897309c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 14 Oct 2021 09:12:45 +0100 Subject: [PATCH 31/76] README shortened, python badge added --- README.md | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/README.md b/README.md index 5bed276d..721ee96c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml)[![](https://img.shields.io/badge/python-3.4+-blue.svg)] # pyerrors 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: @@ -16,18 +16,11 @@ There exist similar implementations of gamma method error analysis suites in - [Python 3](https://github.com/mbruno46/pyobs) ## Installation -pyerrors requires python versions >= 3.5.0 - Install the package for the local user: ```bash pip install . --user ``` -Run tests to verify the installation: -```bash -pytest . -``` - ## Usage The basic objects of a pyerrors analysis are instances of the class `Obs`. They can be initialized with an array of Monte Carlo data (e.g. `samples1`) and a name for the given ensemble (e.g. `'ensemble1'`). The `gamma_method` can then be used to compute the statistical error, taking into account autocorrelations. The `print` method outputs a human readable result. ```python From c6c02cdc57236345bdf584a2b08dc23b0e9b4ce0 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 14 Oct 2021 09:14:22 +0100 Subject: [PATCH 32/76] python badge fixed --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 721ee96c..31adec9d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml)[![](https://img.shields.io/badge/python-3.4+-blue.svg)] +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml)[![](https://img.shields.io/badge/python-3.5+-blue.svg)](https://www.python.org/downloads/) # pyerrors 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: From 8488df00fd5293723cb40ffbfcbc34d9c9678b5e Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 11:47:41 +0100 Subject: [PATCH 33/76] zero_within_error method added, warnings added to Corr methods symmetric and anti-symmetric --- README.md | 2 +- pyerrors/correlators.py | 6 ++++++ pyerrors/pyerrors.py | 3 +++ 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 31adec9d..7401d2b1 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml)[![](https://img.shields.io/badge/python-3.5+-blue.svg)](https://www.python.org/downloads/) +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.5+-blue.svg)](https://www.python.org/downloads/) # pyerrors 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: diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 33809267..10bc3b14 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -129,6 +129,9 @@ class Corr: if self.T % 2 != 0: raise Exception("Can not symmetrize odd T") + if np.argmax(np.abs(self.content)) != 0: + print('Warning: correlator does not seem to be symmetric around x0=0.') + newcontent = [self.content[0]] for t in range(1, self.T): if (self.content[t] is None) or (self.content[self.T - t] is None): @@ -144,6 +147,9 @@ class Corr: if self.T % 2 != 0: raise Exception("Can not symmetrize odd T") + if not all([o.zero_within_error() for o in self.content[0]]): + print('Warning: correlator does not seem to be anti-symmetric around x0=0.') + newcontent = [self.content[0]] for t in range(1, self.T): if (self.content[t] is None) or (self.content[self.T - t] is None): diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 947be666..9a0d02c0 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -340,6 +340,9 @@ class Obs: for e_name in self.e_names: print(e_name, ':', self.e_content[e_name]) + def zero_within_error(self): + return np.abs(self.value) <= self.dvalue + def plot_tauint(self, save=None): """Plot integrated autocorrelation time for each ensemble.""" if not self.e_names: From 4ca6f8ea496405686ad8b32c342d6e4f6722f559 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:11:06 +0100 Subject: [PATCH 34/76] Warning messages promoted to RuntimeWarnings --- pyerrors/correlators.py | 5 +++-- pyerrors/fits.py | 7 +++---- pyerrors/pyerrors.py | 5 +++-- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 10bc3b14..01a7d61e 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -1,3 +1,4 @@ +import warnings import numpy as np import autograd.numpy as anp import matplotlib.pyplot as plt @@ -130,7 +131,7 @@ class Corr: raise Exception("Can not symmetrize odd T") if np.argmax(np.abs(self.content)) != 0: - print('Warning: correlator does not seem to be symmetric around x0=0.') + warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) newcontent = [self.content[0]] for t in range(1, self.T): @@ -148,7 +149,7 @@ class Corr: raise Exception("Can not symmetrize odd T") if not all([o.zero_within_error() for o in self.content[0]]): - print('Warning: correlator does not seem to be anti-symmetric around x0=0.') + warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning) newcontent = [self.content[0]] for t in range(1, self.T): diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 137ede2f..da8aed6b 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # coding: utf-8 +import warnings import numpy as np import autograd.numpy as anp import scipy.optimize @@ -300,7 +301,7 @@ def odr_fit(x, y, func, silent=False, **kwargs): P_phi = A @ np.linalg.inv(A.T @ A) @ A.T expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) if expected_chisquare <= 0.0: - print('Warning, negative expected_chisquare.') + warnings.warn("Negative expected_chisquare.", RuntimeWarning) expected_chisquare = np.abs(expected_chisquare) result_dict['chisquare/expected_chisquare'] = odr_chisquare(np.concatenate((output.beta, output.xplus.ravel()))) / expected_chisquare if not silent: @@ -611,9 +612,7 @@ def error_band(x, func, beta): """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.""" cov = covariance_matrix(beta) if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): - print('Warning, Covariance matrix is not symmetric within floating point precision') - print('cov - cov.T:') - print(cov - cov.T) + warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning) deriv = [] for i, item in enumerate(x): diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 9a0d02c0..15d01f7f 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # coding: utf-8 +import warnings import pickle import numpy as np import autograd.numpy as anp # Thinly-wrapped numpy @@ -784,9 +785,9 @@ def correlate(obs_a, obs_b): raise Exception('Shapes of ensemble', name, 'do not fit') if obs_a.reweighted == 1: - print('Warning: The first observable is already reweighted.') + warnings.warn("The first observable is already reweighted.", RuntimeWarning) if obs_b.reweighted == 1: - print('Warning: The second observable is already reweighted.') + warnings.warn("The second observable is already reweighted.", RuntimeWarning) new_samples = [] for name in sorted(obs_a.names): From c169da550195817d14223dae3dc17fc18ee2cadc Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:15:04 +0100 Subject: [PATCH 35/76] Deprecation warning added to fit general --- pyerrors/fits.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index da8aed6b..1781df37 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -379,7 +379,7 @@ def prior_fit(x, y, func, priors, silent=False, **kwargs): result_dict['fit_function'] = func if Obs.e_tag_global < 4: - print('WARNING: e_tag_global is smaller than 4, this can cause problems when calculating errors from fits with priors') + warnings.warn("e_tag_global is smaller than 4, this can cause problems when calculating errors from fits with priors", RuntimeWarning) x = np.asarray(x) @@ -629,7 +629,6 @@ def error_band(x, func, beta): def fit_general(x, y, func, silent=False, **kwargs): """Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. - WARNING: In the current version the fits are performed with numerical derivatives. Plausibility of the results should be checked. To control the numerical differentiation the kwargs of numdifftools.step_generators.MaxStepGenerator can be used. @@ -650,9 +649,7 @@ def fit_general(x, y, func, silent=False, **kwargs): with many parameters. """ - if not silent: - print('WARNING: This function is deprecated and will be removed in future versions.') - print('New fit functions with exact error propagation are now available as alternative.') + warnings.warn("New fit functions with exact error propagation are now available as alternative.", DeprecationWarning) if not callable(func): raise TypeError('func has to be a function.') From 6ce591dda5d0ee2371cf0dda6dd4756cd8fd9218 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:23:42 +0100 Subject: [PATCH 36/76] Installation instructions updated --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 7401d2b1..06ef30a4 100644 --- a/README.md +++ b/README.md @@ -16,9 +16,9 @@ There exist similar implementations of gamma method error analysis suites in - [Python 3](https://github.com/mbruno46/pyobs) ## Installation -Install the package for the local user: +To install the current `develop` version of `pyerrors` run ```bash -pip install . --user +pip install git+https://github.com/fjosw/pyerrors.git ``` ## Usage From 97434f8a65a31ca381285259fcfa46377e6c75e7 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:29:58 +0100 Subject: [PATCH 37/76] README shortended --- README.md | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 06ef30a4..19e86dd2 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,12 @@ [![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.5+-blue.svg)](https://www.python.org/downloads/) # pyerrors -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: -* automatic differentiation as suggested in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289) (partly based on the [autograd](https://github.com/HIPS/autograd) package) -* the treatment of slow modes in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228) -* multi ensemble analyses -* non-linear fits with y-errors and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289] -* non-linear fits with x- and y-errors and exact linear error propagation based on automatic differentiation -* matrix valued operations and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) -* implementation of the matrix-pencil-method [IEEE Trans. Acoust. 38, 814-824 (1990)](https://ieeexplore.ieee.org/document/56027) for the extraction of energy levels, especially suited for noisy data and excited states +`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: +* **automatic differentiation** as suggested in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289) (partly based on the [autograd](https://github.com/HIPS/autograd) package) +* **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228) +* coherent **error propagation** for data from **different Markov chains** +* **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289] +* **matrix valued operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) There exist similar implementations of gamma method error analysis suites in - [Fortran](https://gitlab.ift.uam-csic.es/alberto/aderrors) From 355a27d95a9563dc65b3fa2393c8d5dcc9a4211b Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:43:42 +0100 Subject: [PATCH 38/76] test_fits added --- tests/test_fits.py | 104 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 tests/test_fits.py diff --git a/tests/test_fits.py b/tests/test_fits.py new file mode 100644 index 00000000..e8441c96 --- /dev/null +++ b/tests/test_fits.py @@ -0,0 +1,104 @@ +import autograd.numpy as np +import os +import random +import string +import copy +import math +import scipy.optimize +from scipy.odr import ODR, Model, Data, RealData +import pyerrors as pe +import pytest + +def test_standard_fit(): + dim = 10 + int(30 * np.random.rand()) + x = np.arange(dim) + y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) + yerr = 0.1 + 0.1 * np.random.rand(dim) + + oy = [] + for i, item in enumerate(x): + oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) + + def f(x, a, b): + return a * np.exp(-b * x) + + popt, pcov = scipy.optimize.curve_fit(f, x, y, sigma=[o.dvalue for o in oy], absolute_sigma=True) + + def func(a, x): + y = a[0] * np.exp(-a[1] * x) + return y + + beta = pe.fits.standard_fit(x, oy, func) + + pe.Obs.e_tag_global = 5 + for i in range(2): + beta[i].gamma_method(e_tag=5, S=1.0) + assert math.isclose(beta[i].value, popt[i], abs_tol=1e-5) + assert math.isclose(pcov[i, i], beta[i].dvalue ** 2, abs_tol=1e-3) + assert math.isclose(pe.covariance(beta[0], beta[1]), pcov[0, 1], abs_tol=1e-3) + pe.Obs.e_tag_global = 0 + + chi2_pyerrors = np.sum(((f(x, *[o.value for o in beta]) - y) / yerr) ** 2) / (len(x) - 2) + chi2_scipy = np.sum(((f(x, *popt) - y) / yerr) ** 2) / (len(x) - 2) + assert math.isclose(chi2_pyerrors, chi2_scipy, abs_tol=1e-10) + + +def test_odr_fit(): + dim = 10 + int(30 * np.random.rand()) + x = np.arange(dim) + np.random.normal(0.0, 0.15, dim) + xerr = 0.1 + 0.1 * np.random.rand(dim) + y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) + yerr = 0.1 + 0.1 * np.random.rand(dim) + + ox = [] + for i, item in enumerate(x): + ox.append(pe.pseudo_Obs(x[i], xerr[i], str(i))) + + oy = [] + for i, item in enumerate(x): + oy.append(pe.pseudo_Obs(y[i], yerr[i], str(i))) + + def f(x, a, b): + return a * np.exp(-b * x) + + def func(a, x): + y = a[0] * np.exp(-a[1] * x) + return y + + data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy]) + model = Model(func) + odr = ODR(data, model, [0,0], partol=np.finfo(np.float).eps) + odr.set_job(fit_type=0, deriv=1) + output = odr.run() + + beta = pe.fits.odr_fit(ox, oy, func) + + pe.Obs.e_tag_global = 5 + for i in range(2): + beta[i].gamma_method(e_tag=5, S=1.0) + assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5) + assert math.isclose(output.cov_beta[i,i], beta[i].dvalue**2, rel_tol=2.5e-1), str(output.cov_beta[i,i]) + ' ' + str(beta[i].dvalue**2) + assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0,1], rel_tol=2.5e-1) + pe.Obs.e_tag_global = 0 + + +def test_odr_derivatives(): + x = [] + y = [] + x_err = 0.01 + y_err = 0.01 + + for n in np.arange(1, 9, 2): + loc_xvalue = n + np.random.normal(0.0, x_err) + x.append(pe.pseudo_Obs(loc_xvalue, x_err, str(n))) + y.append(pe.pseudo_Obs((lambda x: x ** 2 - 1)(loc_xvalue) + + np.random.normal(0.0, y_err), y_err, str(n))) + + def func(a, x): + return a[0] + a[1] * x ** 2 + + fit1 = pe.fits.odr_fit(x, y, func) + + tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20) + assert np.abs(np.max(np.array(list(fit1[1].deltas.values())) + - np.array(list(tfit[1].deltas.values())))) < 10e-8 From 112890c9d212979afe2137b41bd28f88ed3be368 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:45:32 +0100 Subject: [PATCH 39/76] test_fits fixed --- tests/test_fits.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_fits.py b/tests/test_fits.py index e8441c96..e34e1d14 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -67,7 +67,7 @@ def test_odr_fit(): data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy]) model = Model(func) - odr = ODR(data, model, [0,0], partol=np.finfo(np.float).eps) + odr = ODR(data, model, [0,0], partol=np.finfo(np.float64).eps) odr.set_job(fit_type=0, deriv=1) output = odr.run() From 1b63f5dd3b75f4f7e327268bda8591a84237bafd Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 12:56:07 +0100 Subject: [PATCH 40/76] test_roots added --- tests/test_roots.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 tests/test_roots.py diff --git a/tests/test_roots.py b/tests/test_roots.py new file mode 100644 index 00000000..a6ce3d65 --- /dev/null +++ b/tests/test_roots.py @@ -0,0 +1,14 @@ +import numpy as np +import pyerrors as pe +import pytest + +def test_find_root(): + + def root_function(x, d): + return x - d + + value = np.random.normal(0, 100) + my_obs = pe.pseudo_Obs(value, 0.1, 't') + my_root = pe.roots.find_root(my_obs, root_function) + + assert np.isclose(my_root.value, value) From bb9790acd709dbfeeff3656e5b8db768602cf8cb Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 13:05:00 +0100 Subject: [PATCH 41/76] fixed seed added to tests --- tests/test_fits.py | 2 ++ tests/test_linalg.py | 2 ++ tests/test_pyerrors.py | 2 ++ tests/test_roots.py | 7 ++++++- 4 files changed, 12 insertions(+), 1 deletion(-) diff --git a/tests/test_fits.py b/tests/test_fits.py index e34e1d14..470c6a74 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -9,6 +9,8 @@ from scipy.odr import ODR, Model, Data, RealData import pyerrors as pe import pytest +np.random.seed(0) + def test_standard_fit(): dim = 10 + int(30 * np.random.rand()) x = np.arange(dim) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 9d7f6377..e13c80d5 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -11,6 +11,8 @@ from scipy.odr import ODR, Model, Data, RealData import pyerrors as pe import pytest +np.random.seed(0) + def test_matrix_functions(): dim = 3 + int(4 * np.random.rand()) print(dim) diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 393cde8d..fcc00e6d 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -6,6 +6,8 @@ import copy import pyerrors as pe import pytest +np.random.seed(0) + def test_dump(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) diff --git a/tests/test_roots.py b/tests/test_roots.py index a6ce3d65..dce0eb45 100644 --- a/tests/test_roots.py +++ b/tests/test_roots.py @@ -2,7 +2,9 @@ import numpy as np import pyerrors as pe import pytest -def test_find_root(): +np.random.seed(0) + +def test_root_linear(): def root_function(x, d): return x - d @@ -12,3 +14,6 @@ def test_find_root(): my_root = pe.roots.find_root(my_obs, root_function) assert np.isclose(my_root.value, value) + difference = my_obs - my_root + assert all(np.isclose(0.0, difference.deltas['t'])) + From e46746e4caf14a792aaded70ff47068d9f35f19c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 14:13:06 +0100 Subject: [PATCH 42/76] removed unnecessary imports from tests --- tests/test_fits.py | 17 +++++++---------- tests/test_linalg.py | 10 +--------- tests/test_pyerrors.py | 5 +++-- tests/test_roots.py | 2 +- 4 files changed, 12 insertions(+), 22 deletions(-) diff --git a/tests/test_fits.py b/tests/test_fits.py index 470c6a74..461136ce 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -1,16 +1,13 @@ import autograd.numpy as np -import os -import random -import string -import copy import math import scipy.optimize -from scipy.odr import ODR, Model, Data, RealData +from scipy.odr import ODR, Model, RealData import pyerrors as pe import pytest np.random.seed(0) + def test_standard_fit(): dim = 10 + int(30 * np.random.rand()) x = np.arange(dim) @@ -69,7 +66,7 @@ def test_odr_fit(): data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy]) model = Model(func) - odr = ODR(data, model, [0,0], partol=np.finfo(np.float64).eps) + odr = ODR(data, model, [0, 0], partol=np.finfo(np.float64).eps) odr.set_job(fit_type=0, deriv=1) output = odr.run() @@ -79,8 +76,8 @@ def test_odr_fit(): for i in range(2): beta[i].gamma_method(e_tag=5, S=1.0) assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5) - assert math.isclose(output.cov_beta[i,i], beta[i].dvalue**2, rel_tol=2.5e-1), str(output.cov_beta[i,i]) + ' ' + str(beta[i].dvalue**2) - assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0,1], rel_tol=2.5e-1) + assert math.isclose(output.cov_beta[i, i], beta[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(beta[i].dvalue ** 2) + assert math.isclose(pe.covariance(beta[0], beta[1]), output.cov_beta[0, 1], rel_tol=2.5e-1) pe.Obs.e_tag_global = 0 @@ -94,7 +91,7 @@ def test_odr_derivatives(): loc_xvalue = n + np.random.normal(0.0, x_err) x.append(pe.pseudo_Obs(loc_xvalue, x_err, str(n))) y.append(pe.pseudo_Obs((lambda x: x ** 2 - 1)(loc_xvalue) + - np.random.normal(0.0, y_err), y_err, str(n))) + np.random.normal(0.0, y_err), y_err, str(n))) def func(a, x): return a[0] + a[1] * x ** 2 @@ -103,4 +100,4 @@ def test_odr_derivatives(): tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20) assert np.abs(np.max(np.array(list(fit1[1].deltas.values())) - - np.array(list(tfit[1].deltas.values())))) < 10e-8 + - np.array(list(tfit[1].deltas.values())))) < 10e-8 diff --git a/tests/test_linalg.py b/tests/test_linalg.py index e13c80d5..042c1d3a 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1,18 +1,11 @@ -import sys -sys.path.append('..') import autograd.numpy as np -import os -import random import math -import string -import copy -import scipy.optimize -from scipy.odr import ODR, Model, Data, RealData import pyerrors as pe import pytest np.random.seed(0) + def test_matrix_functions(): dim = 3 + int(4 * np.random.rand()) print(dim) @@ -55,4 +48,3 @@ def test_matrix_functions(): tmp[j].gamma_method() assert math.isclose(tmp[j].value, 0.0, abs_tol=1e-9), 'value ' + str(i) + ',' + str(j) assert math.isclose(tmp[j].dvalue, 0.0, abs_tol=1e-9), 'dvalue ' + str(i) + ',' + str(j) - diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index fcc00e6d..0835fb08 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -8,6 +8,7 @@ import pytest np.random.seed(0) + def test_dump(): value = np.random.normal(5, 10) dvalue = np.abs(np.random.normal(0, 1)) @@ -38,8 +39,8 @@ def test_function_overloading(): lambda x: np.sinh(x[0]), lambda x: np.cosh(x[0]), lambda x: np.tanh(x[0])] for i, f in enumerate(fs): - t1 = f([a,b]) - t2 = pe.derived_observable(f, [a,b]) + t1 = f([a, b]) + t2 = pe.derived_observable(f, [a, b]) c = t2 - t1 assert c.value == 0.0, str(i) assert np.all(np.abs(c.deltas['e1']) < 1e-14), str(i) diff --git a/tests/test_roots.py b/tests/test_roots.py index dce0eb45..c1bdc7ea 100644 --- a/tests/test_roots.py +++ b/tests/test_roots.py @@ -4,6 +4,7 @@ import pytest np.random.seed(0) + def test_root_linear(): def root_function(x, d): @@ -16,4 +17,3 @@ def test_root_linear(): assert np.isclose(my_root.value, value) difference = my_obs - my_root assert all(np.isclose(0.0, difference.deltas['t'])) - From 25d250cd53edd2160f5abcefb4eb9fe363811f3f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 15:01:53 +0100 Subject: [PATCH 43/76] test_correlators added, tests extended --- pyerrors/correlators.py | 4 +++- tests/test_correlators.py | 48 +++++++++++++++++++++++++++++++++++++++ tests/test_pyerrors.py | 15 +++++++++--- tests/test_roots.py | 2 +- 4 files changed, 64 insertions(+), 5 deletions(-) create mode 100644 tests/test_correlators.py diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 01a7d61e..ff899dc8 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -242,7 +242,9 @@ class Corr: Parameters ---------- variant -- log: uses the standard effective mass log(C(t) / C(t+1)) - periodic : Solves C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. See, e.g., arXiv:1205.5380 + cosh : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. + sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. + See, e.g., arXiv:1205.5380 guess -- guess for the root finder, only relevant for the root variant """ if self.N != 1: diff --git a/tests/test_correlators.py b/tests/test_correlators.py new file mode 100644 index 00000000..37b59fba --- /dev/null +++ b/tests/test_correlators.py @@ -0,0 +1,48 @@ +import numpy as np +import pyerrors as pe +import pytest + +np.random.seed(0) + + +def test_function_overloading(): + corr_content_a = [] + corr_content_b = [] + for t in range(24): + corr_content_a.append(pe.pseudo_Obs(np.random.normal(1e-10, 1e-8), 1e-4, 't')) + corr_content_b.append(pe.pseudo_Obs(np.random.normal(1e8, 1e10), 1e7, 't')) + + corr_a = pe.correlators.Corr(corr_content_a) + corr_b = pe.correlators.Corr(corr_content_b) + + fs = [lambda x: x[0] + x[1], lambda x: x[1] + x[0], lambda x: x[0] - x[1], lambda x: x[1] - x[0], + lambda x: x[0] * x[1], lambda x: x[1] * x[0], lambda x: x[0] / x[1], lambda x: x[1] / x[0], + lambda x: np.exp(x[0]), lambda x: np.sin(x[0]), lambda x: np.cos(x[0]), lambda x: np.tan(x[0]), + lambda x: np.log(x[0] + 0.1), lambda x: np.sqrt(np.abs(x[0])), + lambda x: np.sinh(x[0]), lambda x: np.cosh(x[0]), lambda x: np.tanh(x[0])] + + for i, f in enumerate(fs): + t1 = f([corr_a, corr_b]) + for o_a, o_b, con in zip(corr_content_a, corr_content_b, t1.content): + t2 = f([o_a, o_b]) + t2.gamma_method() + assert np.isclose(con[0].value, t2.value) + assert np.isclose(con[0].dvalue, t2.dvalue) + assert np.allclose(con[0].deltas['t'], t2.deltas['t']) + +def test_modify_correlator(): + corr_content = [] + for t in range(24): + exponent = np.random.normal(3, 5) + corr_content.append(pe.pseudo_Obs(2 + 10 ** exponent, 10 ** (exponent - 1), 't')) + + corr = pe.correlators.Corr(corr_content) + + with pytest.warns(RuntimeWarning): + corr.symmetric() + with pytest.warns(RuntimeWarning): + corr.anti_symmetric() + corr.roll(np.random.randint(100)) + corr.deriv(symmetric=True) + corr.deriv(symmetric=False) + corr.second_deriv() diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 0835fb08..4c8cf18f 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -35,7 +35,7 @@ def test_function_overloading(): fs = [lambda x: x[0] + x[1], lambda x: x[1] + x[0], lambda x: x[0] - x[1], lambda x: x[1] - x[0], lambda x: x[0] * x[1], lambda x: x[1] * x[0], lambda x: x[0] / x[1], lambda x: x[1] / x[0], lambda x: np.exp(x[0]), lambda x: np.sin(x[0]), lambda x: np.cos(x[0]), lambda x: np.tan(x[0]), - lambda x: np.log(x[0]), lambda x: np.sqrt(x[0]), + lambda x: np.log(x[0]), lambda x: np.sqrt(np.abs(x[0])), lambda x: np.sinh(x[0]), lambda x: np.cosh(x[0]), lambda x: np.tanh(x[0])] for i, f in enumerate(fs): @@ -47,8 +47,17 @@ def test_function_overloading(): def test_overloading_vectorization(): - a = np.random.randint(0, 100, 10) - b = pe.pseudo_Obs(4, 0.8, 'e1') + a = np.random.randint(1, 100, 10) + b = pe.pseudo_Obs(4, 0.8, 't') + + assert [o.value for o in a * b] == [o.value for o in b * a] + assert [o.value for o in a + b] == [o.value for o in b + a] + assert [o.value for o in a - b] == [-1 * o.value for o in b - a] + assert [o.value for o in a / b] == [o.value for o in [p / b for p in a]] + assert [o.value for o in b / a] == [o.value for o in [b / p for p in a]] + + a = np.random.normal(0.0, 1e10, 10) + b = pe.pseudo_Obs(4, 0.8, 't') assert [o.value for o in a * b] == [o.value for o in b * a] assert [o.value for o in a + b] == [o.value for o in b + a] diff --git a/tests/test_roots.py b/tests/test_roots.py index c1bdc7ea..8d3a8191 100644 --- a/tests/test_roots.py +++ b/tests/test_roots.py @@ -16,4 +16,4 @@ def test_root_linear(): assert np.isclose(my_root.value, value) difference = my_obs - my_root - assert all(np.isclose(0.0, difference.deltas['t'])) + assert np.allclose(0.0, difference.deltas['t']) From 4f8345ac86304cdb401ce637ee3a07df5e8c34b4 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 15:27:33 +0100 Subject: [PATCH 44/76] test_m_eff and test_utility added --- .github/workflows/CI.yml | 2 +- pyerrors/correlators.py | 2 +- tests/test_correlators.py | 25 +++++++++++++++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d0608ae5..67343d6e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: true matrix: - python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] + python-version: ["3.5", "3.6", "3.7", "3.8", "3.9", "3.10"] steps: - name: Checkout source diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index ff899dc8..ee949276 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -427,7 +427,7 @@ class Corr: if isinstance(save, str): fig.savefig(save) else: - raise Exception('safe has to be a string.') + raise Exception("Safe has to be a string.") return diff --git a/tests/test_correlators.py b/tests/test_correlators.py index 37b59fba..502f26fa 100644 --- a/tests/test_correlators.py +++ b/tests/test_correlators.py @@ -1,3 +1,4 @@ +import os import numpy as np import pyerrors as pe import pytest @@ -46,3 +47,27 @@ def test_modify_correlator(): corr.deriv(symmetric=True) corr.deriv(symmetric=False) corr.second_deriv() + +def test_m_eff(): + my_corr = pe.correlators.Corr([pe.pseudo_Obs(10, 0.1, 't'), pe.pseudo_Obs(9, 0.05, 't')]) + my_corr.m_eff('log') + my_corr.m_eff('cosh') + +def test_utility(): + corr_content = [] + for t in range(8): + exponent = np.random.normal(3, 5) + corr_content.append(pe.pseudo_Obs(2 + 10 ** exponent, 10 ** (exponent - 1), 't')) + + corr = pe.correlators.Corr(corr_content) + corr.print() + corr.print([2, 4]) + corr.show() + + corr.dump('test_dump') + new_corr = pe.load_object('test_dump.p') + os.remove('test_dump.p') + for o_a, o_b in zip(corr.content, new_corr.content): + assert np.isclose(o_a[0].value, o_b[0].value) + assert np.isclose(o_a[0].dvalue, o_b[0].dvalue) + assert np.allclose(o_a[0].deltas['t'], o_b[0].deltas['t']) From b491096a956ed48ae379b7275c779f3b047cfc6a Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 15:46:59 +0100 Subject: [PATCH 45/76] utils test added --- .github/workflows/CI.yml | 4 ++-- pyerrors/pyerrors.py | 6 +++--- tests/test_pyerrors.py | 14 ++++++++++++++ 3 files changed, 19 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 67343d6e..3ee55571 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: true matrix: - python-version: ["3.5", "3.6", "3.7", "3.8", "3.9", "3.10"] + python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] steps: - name: Checkout source @@ -27,7 +27,7 @@ jobs: - name: Install run: | python -m pip install --upgrade pip - pip install -e . + pip install . pip install pytest pip install pytest-cov diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 15d01f7f..65613ee8 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -419,7 +419,7 @@ class Obs: arr[r] = (self.r_values[r_name] - sub_r_mean) / (self.e_dvalue[e_name] * np.sqrt(e_N / self.shape[r_name] - 1)) plt.hist(arr, rwidth=0.8, bins=len(self.e_content[e_name])) plt.title('Replica distribution' + e_name + ' (mean=0, var=1), Q=' + str(np.around(self.e_Q[e_name], decimals=2))) - plt.show() + plt.draw() def plot_history(self): """Plot derived Monte Carlo history for each ensemble.""" @@ -440,7 +440,7 @@ class Obs: plt.errorbar(x, y, fmt='.', markersize=3) plt.xlim(-0.5, e_N - 0.5) plt.title(e_name) - plt.show() + plt.draw() def plot_piechart(self): """Plot piechart which shows the fractional contribution of each @@ -454,7 +454,7 @@ class Obs: fig1, ax1 = plt.subplots() ax1.pie(sizes, labels=labels, startangle=90, normalize=True) ax1.axis('equal') - plt.show() + plt.draw() return dict(zip(self.e_names, sizes)) diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 4c8cf18f..c06f73d9 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -171,3 +171,17 @@ def test_overloaded_functions(): assert np.max((ad_obs.deltas['t'] - fd_obs.deltas['t']) / ad_obs.deltas['t']) < 1e-8, item.__name__ assert np.abs((ad_obs.value - item(val)) / ad_obs.value) < 1e-10, item.__name__ assert np.abs(ad_obs.dvalue - dval * np.abs(deriv[i](val))) < 1e-6, item.__name__ + +def test_utils(): + my_obs = pe.pseudo_Obs(1.0, 0.5, 't') + my_obs.print(0) + my_obs.print(1) + my_obs.print(2) + assert not my_obs.zero_within_error() + my_obs.plot_tauint() + my_obs.plot_rho() + my_obs.plot_rep_dist() + my_obs.plot_history() + my_obs.plot_piechart() + assert my_obs > (my_obs - 1) + assert my_obs < (my_obs + 1) From 1972f38cb9b3add88bc762d72903ca75c085f90d Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 15:51:12 +0100 Subject: [PATCH 46/76] Minimal python version changed to 3.6 --- .github/workflows/CI.yml | 2 +- README.md | 2 +- setup.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 3ee55571..f1fddc58 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: true matrix: - python-version: ["3.5", "3.6", "3.7", "3.8", "3.9"] + python-version: ["3.6", "3.7", "3.8", "3.9"] steps: - name: Checkout source diff --git a/README.md b/README.md index 19e86dd2..f6bf3a5d 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.5+-blue.svg)](https://www.python.org/downloads/) +[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.6+-blue.svg)](https://www.python.org/downloads/) # pyerrors `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: diff --git a/setup.py b/setup.py index a5b747fa..988c32a2 100644 --- a/setup.py +++ b/setup.py @@ -8,6 +8,6 @@ setup(name='pyerrors', author='Fabian Joswig', author_email='fabian.joswig@ed.ac.uk', packages=find_packages(), - python_requires='>=3.5.0', + python_requires='>=3.6.0', install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib', 'scipy', 'iminuit<2'] ) From 5f1fb285f561098381eccce737822de8489b9bb9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 15 Oct 2021 16:02:17 +0100 Subject: [PATCH 47/76] workaround for m_eff sinh in the middle of the lattice --- pyerrors/correlators.py | 3 +++ setup.py | 2 +- tests/test_correlators.py | 3 ++- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index ee949276..41c61985 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -274,6 +274,9 @@ class Corr: for t in range(self.T - 1): if (self.content[t] is None) or (self.content[t + 1] is None): newcontent.append(None) + # Fill the two timeslices in the middle of the lattice with their predecessors + elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: + newcontent.append(newcontent[-1]) else: newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) if(all([x is None for x in newcontent])): diff --git a/setup.py b/setup.py index 988c32a2..df32d698 100644 --- a/setup.py +++ b/setup.py @@ -9,5 +9,5 @@ setup(name='pyerrors', author_email='fabian.joswig@ed.ac.uk', packages=find_packages(), python_requires='>=3.6.0', - install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib', 'scipy', 'iminuit<2'] + install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit<2'] ) diff --git a/tests/test_correlators.py b/tests/test_correlators.py index 502f26fa..63f1a755 100644 --- a/tests/test_correlators.py +++ b/tests/test_correlators.py @@ -49,9 +49,10 @@ def test_modify_correlator(): corr.second_deriv() def test_m_eff(): - my_corr = pe.correlators.Corr([pe.pseudo_Obs(10, 0.1, 't'), pe.pseudo_Obs(9, 0.05, 't')]) + my_corr = pe.correlators.Corr([pe.pseudo_Obs(10, 0.1, 't'), pe.pseudo_Obs(9, 0.05, 't'), pe.pseudo_Obs(8, 0.1, 't'), pe.pseudo_Obs(7, 0.05, 't')]) my_corr.m_eff('log') my_corr.m_eff('cosh') + my_corr.m_eff('sinh') def test_utility(): corr_content = [] From ea51df3719ffebffecbe2373d52d1324517ca5e9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 11:02:54 +0100 Subject: [PATCH 48/76] __repr__ of Obs changed to more convenient form --- pyerrors/correlators.py | 4 ++-- pyerrors/pyerrors.py | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 41c61985..9868a5c1 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -402,7 +402,7 @@ class Corr: if plateau: if isinstance(plateau, Obs): - ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=plateau.__repr__()[4:-1]) + ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=plateau.__repr__()) ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') else: raise Exception('plateau must be an Obs') @@ -451,7 +451,7 @@ class Corr: else: content_string += str(i + range[0]) for element in sub_corr: - content_string += '\t' + element.__repr__()[4:-1] + content_string += '\t' + element.__repr__() content_string += '\n' return content_string diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 65613ee8..5b0dfb95 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -474,14 +474,17 @@ class Obs: def __repr__(self): if self.dvalue == 0.0: - return 'Obs[' + str(self.value) + ']' + return str(self.value) fexp = np.floor(np.log10(self.dvalue)) if fexp < 0.0: - return 'Obs[{:{form}}({:2.0f})]'.format(self.value, self.dvalue * 10 ** (-fexp + 1), form='.' + str(-int(fexp) + 1) + 'f') + return '{:{form}}({:2.0f})'.format(self.value, self.dvalue * 10 ** (-fexp + 1), form='.' + str(-int(fexp) + 1) + 'f') elif fexp == 0.0: - return 'Obs[{:.1f}({:1.1f})]'.format(self.value, self.dvalue) + return '{:.1f}({:1.1f})'.format(self.value, self.dvalue) else: - return 'Obs[{:.0f}({:2.0f})]'.format(self.value, self.dvalue) + return '{:.0f}({:2.0f})'.format(self.value, self.dvalue) + + def __str__(self): + return 'Obs[' + str(self) + ']' # Overload comparisons def __lt__(self, other): From 62fe1f9e3a7b680417a204843fbaae439942c5a4 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 11:06:11 +0100 Subject: [PATCH 49/76] recursion depth in __str__ fixed --- pyerrors/pyerrors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 5b0dfb95..5febd9ec 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -484,7 +484,7 @@ class Obs: return '{:.0f}({:2.0f})'.format(self.value, self.dvalue) def __str__(self): - return 'Obs[' + str(self) + ']' + return 'Obs[' + self.__repr__() + ']' # Overload comparisons def __lt__(self, other): From be4d9f722e51d44cbe6ac5cf3ea60902f68ccb30 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 11:21:19 +0100 Subject: [PATCH 50/76] repr and str exchanged --- pyerrors/correlators.py | 4 ++-- pyerrors/pyerrors.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 9868a5c1..0e8b74fa 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -402,7 +402,7 @@ class Corr: if plateau: if isinstance(plateau, Obs): - ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=plateau.__repr__()) + ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') else: raise Exception('plateau must be an Obs') @@ -451,7 +451,7 @@ class Corr: else: content_string += str(i + range[0]) for element in sub_corr: - content_string += '\t' + element.__repr__() + content_string += '\t' + str(element) content_string += '\n' return content_string diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 5febd9ec..e4262af6 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -473,6 +473,9 @@ class Obs: pickle.dump(self, fb) def __repr__(self): + return 'Obs[' + str(self) + ']' + + def __str__(self): if self.dvalue == 0.0: return str(self.value) fexp = np.floor(np.log10(self.dvalue)) @@ -483,9 +486,6 @@ class Obs: else: return '{:.0f}({:2.0f})'.format(self.value, self.dvalue) - def __str__(self): - return 'Obs[' + self.__repr__() + ']' - # Overload comparisons def __lt__(self, other): return self.value < other From dc4514e8d033f9bf15f54a8a09ec98c130a8ed13 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 11:41:44 +0100 Subject: [PATCH 51/76] float special method added to pyerrors --- pyerrors/pyerrors.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index e4262af6..20d64fd5 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -472,6 +472,9 @@ class Obs: with open(file_name, 'wb') as fb: pickle.dump(self, fb) + def __float__(self): + return self.value + def __repr__(self): return 'Obs[' + str(self) + ']' From 286008a20e227dc22b55ee81587d0f575fdaa57b Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 12:07:40 +0100 Subject: [PATCH 52/76] CObs class added --- pyerrors/pyerrors.py | 61 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 20d64fd5..a807871f 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -632,6 +632,67 @@ class Obs: return derived_observable(lambda x: anp.sinc(x[0]), [self]) +class CObs: + def __init__(self, real, imag=0.0): + self.real = real + self.imag = imag + + def gamma_method(self, **kwargs): + self.real.gamma_method(**kwargs) + self.imag.gamma_method(**kwargs) + + def __add__(self, other): + if hasattr(other, 'real') and hasattr(other, 'imag'): + return CObs(self.real + other.real, + self.imag + other.imag) + else: + return CObs(self.real + other, self.imag) + + def __radd__(self, y): + return self + y + + def __sub__(self, other): + if hasattr(other, 'real') and hasattr(other, 'imag'): + return CObs(self.real - other.real, self.imag - other.imag) + else: + return CObs(self.real - other, self.imag) + + def __rsub__(self, other): + return -1 * (self - other) + + def __mul__(self, other): + if hasattr(other, 'real') and hasattr(other, 'imag'): + return CObs(self.real * other.real - self.imag * other.imag, + self.imag * other.real + self.real * other.imag) + else: + return CObs(self.real * other, self.imag * other) + + def __rmul__(self, other): + return self * other + + def __truediv__(self, other): + if hasattr(other, 'real') and hasattr(other, 'imag'): + r = other.real ** 2 + other.imag ** 2 + return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r) + else: + return CObs(self.real / other, self.imag / other) + + def __abs__(self): + return np.sqrt(self.real**2 + self.imag**2) + + def __neg__(other): + return -1 * other + + def __eq__(self, other): + return self.real == other.real and self.imag == other.imag + + def __str__(self): + return '(' + str(self.real) + int(self.imag > - np.finfo(np.float64).eps) * '+' + str(self.imag) + 'j)' + + def __repr__(self): + return 'CObs[' + str(self) + ']' + + def derived_observable(func, data, **kwargs): """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation. From db71644e5aca92a2f5f7927e2270e59a2922d643 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 16 Oct 2021 12:18:46 +0100 Subject: [PATCH 53/76] First tests for CObs added --- tests/test_pyerrors.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index c06f73d9..3be59b04 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -185,3 +185,31 @@ def test_utils(): my_obs.plot_piechart() assert my_obs > (my_obs - 1) assert my_obs < (my_obs + 1) + +def test_cobs(): + obs1 = pe.pseudo_Obs(1.0, 0.1, 't') + obs2 = pe.pseudo_Obs(-0.2, 0.03, 't') + + my_cobs = pe.CObs(obs1, obs2) + np.abs(my_cobs) + fs = [[lambda x: x[0] + x[1], lambda x: x[1] + x[0]], + [lambda x: x[0] * x[1], lambda x: x[1] * x[0]]] + for other in [1, 1.1, (1.1-0.2j), pe.CObs(obs1), pe.CObs(obs1, obs2)]: + for funcs in fs: + ta = funcs[0]([my_cobs, other]) + tb = funcs[1]([my_cobs, other]) + diff = ta - tb + assert np.isclose(0.0, float(diff.real)) + assert np.isclose(0.0, float(diff.imag)) + assert np.allclose(0.0, diff.real.deltas['t']) + assert np.allclose(0.0, diff.imag.deltas['t']) + + ta = my_cobs - other + tb = other - my_cobs + diff = ta + tb + assert np.isclose(0.0, float(diff.real)) + assert np.isclose(0.0, float(diff.imag)) + assert np.allclose(0.0, diff.real.deltas['t']) + assert np.allclose(0.0, diff.imag.deltas['t']) + + div = my_cobs / other From aa4c037e397b7877c14f6a4e1673dbfc07f749e6 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 11:34:53 +0100 Subject: [PATCH 54/76] CObs.gamma_method() works now when real or imag is 0 --- pyerrors/pyerrors.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index a807871f..82bf7154 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -638,8 +638,10 @@ class CObs: self.imag = imag def gamma_method(self, **kwargs): - self.real.gamma_method(**kwargs) - self.imag.gamma_method(**kwargs) + if isinstance(self.real, Obs): + self.real.gamma_method(**kwargs) + if isinstance(self.imag, Obs): + self.imag.gamma_method(**kwargs) def __add__(self, other): if hasattr(other, 'real') and hasattr(other, 'imag'): From 573d4992c9dfd0d12063e2cdd879175c161ca8ee Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 12:23:27 +0100 Subject: [PATCH 55/76] Workaround implemented to deal with matrix operations containing floats or integers --- pyerrors/pyerrors.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 82bf7154..64d58a00 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -517,9 +517,10 @@ class Obs: else: if isinstance(y, np.ndarray): return np.array([self * o for o in y]) + elif isinstance(y, complex): + return CObs(self * y.real, self * y.imag) elif y.__class__.__name__ == 'Corr': return NotImplemented - else: return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y]) @@ -729,6 +730,17 @@ def derived_observable(func, data, **kwargs): data = np.asarray(data) raveled_data = data.ravel() + # Workaround for matrix operations containing non Obs data + for i_data in raveled_data: + if isinstance(i_data, Obs): + first_name = i_data.names[0] + first_shape = i_data.shape[first_name] + break + + for i in range(len(raveled_data)): + if isinstance(raveled_data[i], (int, float)): + raveled_data[i] = Obs([raveled_data[i] + np.zeros(first_shape)], [first_name]) + n_obs = len(raveled_data) new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x])) replicas = len(new_names) From 10b228d434ed8ea61c6f3ee51d702849d57f1326 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 12:28:59 +0100 Subject: [PATCH 56/76] New method is_zero --- pyerrors/correlators.py | 2 +- pyerrors/pyerrors.py | 5 ++++- tests/test_pyerrors.py | 2 +- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 0e8b74fa..8ef26ce2 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -148,7 +148,7 @@ class Corr: if self.T % 2 != 0: raise Exception("Can not symmetrize odd T") - if not all([o.zero_within_error() for o in self.content[0]]): + if not all([o.is_zero_within_error() for o in self.content[0]]): warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning) newcontent = [self.content[0]] diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 64d58a00..64ad12c9 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -341,9 +341,12 @@ class Obs: for e_name in self.e_names: print(e_name, ':', self.e_content[e_name]) - def zero_within_error(self): + def is_zero_within_error(self): return np.abs(self.value) <= self.dvalue + def is_zero(self): + np.isclose(0.0, self.value) and all(np.allclose(0.0, delta) for delta in self.deltas.values()) + def plot_tauint(self, save=None): """Plot integrated autocorrelation time for each ensemble.""" if not self.e_names: diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 3be59b04..982dcab3 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -177,7 +177,7 @@ def test_utils(): my_obs.print(0) my_obs.print(1) my_obs.print(2) - assert not my_obs.zero_within_error() + assert not my_obs.is_zero_within_error() my_obs.plot_tauint() my_obs.plot_rho() my_obs.plot_rep_dist() From 15333d2629e0c6285c79a78aed19d7385f597b8c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 12:45:30 +0100 Subject: [PATCH 57/76] New test for matrix inverse with matrix containing floats r integers --- pyerrors/pyerrors.py | 2 +- tests/test_linalg.py | 12 ++++++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 64ad12c9..6039b412 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -345,7 +345,7 @@ class Obs: return np.abs(self.value) <= self.dvalue def is_zero(self): - np.isclose(0.0, self.value) and all(np.allclose(0.0, delta) for delta in self.deltas.values()) + return np.isclose(0.0, self.value) and all(np.allclose(0.0, delta) for delta in self.deltas.values()) def plot_tauint(self, save=None): """Plot integrated autocorrelation time for each ensemble.""" diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 042c1d3a..27653190 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -6,6 +6,18 @@ import pytest np.random.seed(0) +def test_matrix_inverse(): + content = [] + for t in range(9): + exponent = np.random.normal(3, 5) + content.append(pe.pseudo_Obs(2 + 10 ** exponent, 10 ** (exponent - 1), 't')) + + content.append(1.0) # Add 1.0 as a float + matrix = np.diag(content) + inverse_matrix = pe.linalg.mat_mat_op(np.linalg.inv, matrix) + assert all([o.is_zero() for o in np.diag(matrix) * np.diag(inverse_matrix) - 1]) + + def test_matrix_functions(): dim = 3 + int(4 * np.random.rand()) print(dim) From 157fc1058a3788949d73bc212228244e6d4ab6f2 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 13:03:14 +0100 Subject: [PATCH 58/76] mat_mat_op now also works with complex matrices --- pyerrors/linalg.py | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index f8dd6f8b..a8810091 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -3,7 +3,7 @@ import numpy as np import autograd.numpy as anp # Thinly-wrapped numpy -from .pyerrors import derived_observable +from .pyerrors import derived_observable, CObs # This code block is directly taken from the current master branch of autograd and remains @@ -84,9 +84,30 @@ def scalar_mat_op(op, obs, **kwargs): def mat_mat_op(op, obs, **kwargs): """Computes the matrix to matrix operation op to a given matrix of Obs.""" - if kwargs.get('num_grad') is True: - return _num_diff_mat_mat_op(op, obs, **kwargs) - return derived_observable(lambda x, **kwargs: op(x), obs) + # Use real representation to calculate matrix operations for complex matrices + if isinstance(obs.ravel()[0], CObs): + A = np.empty_like(obs) + B = np.empty_like(obs) + for (n, m), entry in np.ndenumerate(obs): + if hasattr(entry, 'real') and hasattr(entry, 'imag'): + A[n, m] = entry.real + B[n, m] = entry.imag + else: + A[n, m] = entry + B[n, m] = 0.0 + big_matrix = np.bmat([[A, -B], [B, A]]) + if kwargs.get('num_grad') is True: + op_big_matrix = _num_diff_mat_mat_op(op, big_matrix, **kwargs) + else: + op_big_matrix = derived_observable(lambda x, **kwargs: op(x), big_matrix) + dim = op_big_matrix.shape[0] + op_A = op_big_matrix[0: dim // 2, 0: dim // 2] + op_B = op_big_matrix[dim // 2:, 0: dim // 2] + return (1 + 0j) * op_A + 1j * op_B + else: + if kwargs.get('num_grad') is True: + return _num_diff_mat_mat_op(op, obs, **kwargs) + return derived_observable(lambda x, **kwargs: op(x), obs) def eigh(obs, **kwargs): From abdff600e68619006312997ce41a2a47f659ae76 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 13:35:27 +0100 Subject: [PATCH 59/76] Added test_complex_matrix_inverse --- tests/test_linalg.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 27653190..af121912 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -18,6 +18,29 @@ def test_matrix_inverse(): assert all([o.is_zero() for o in np.diag(matrix) * np.diag(inverse_matrix) - 1]) +def test_complex_matrix_inverse(): + dimension = 6 + base_matrix = np.empty((dimension, dimension), dtype=object) + matrix = np.empty((dimension, dimension), dtype=complex) + for (n, m), entry in np.ndenumerate(base_matrix): + exponent_real = np.random.normal(3, 5) + exponent_imag = np.random.normal(3, 5) + base_matrix[n, m] = pe.CObs(pe.pseudo_Obs(2 + 10 ** exponent_real, 10 ** (exponent_real - 1), 't'), + pe.pseudo_Obs(2 + 10 ** exponent_imag, 10 ** (exponent_imag - 1), 't')) + + # Construct invertible matrix + obs_matrix = np.identity(dimension) + base_matrix @ base_matrix.T + + for (n, m), entry in np.ndenumerate(obs_matrix): + matrix[n, m] = entry.real.value + 1j * entry.imag.value + + inverse_matrix = np.linalg.inv(matrix) + inverse_obs_matrix = pe.linalg.mat_mat_op(np.linalg.inv, obs_matrix) + for (n, m), entry in np.ndenumerate(inverse_matrix): + assert np.isclose(inverse_matrix[n, m].real, inverse_obs_matrix[n, m].real.value) + assert np.isclose(inverse_matrix[n, m].imag, inverse_obs_matrix[n, m].imag.value) + + def test_matrix_functions(): dim = 3 + int(4 * np.random.rand()) print(dim) From c927ebf0430337b1c5426e589e5473fd67dcd4d6 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 13:44:24 +0100 Subject: [PATCH 60/76] README and CHANGELOG updated --- CHANGELOG.md | 6 ++++++ README.md | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ae85cf5d..240778ea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,12 @@ All notable changes to this project will be documented in this file. +## [2.0.0] - 2021-??-?? +## Added +- `CObs` class added which can handle complex valued observables +- Matrix to matrix operations like the matrix inverse now also work for complex matrices and matrices containing entries that are not `Obs` but `float` or `int` +- `Obs` objects now have methods `is_zero` and `is_zero_within_error` + ## [1.1.0] - 2021-10-11 ### Added - `Corr` class added diff --git a/README.md b/README.md index f6bf3a5d..df93c050 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/ab * **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228) * coherent **error propagation** for data from **different Markov chains** * **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289] -* **matrix valued operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) +* **real and complex matrix operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) There exist similar implementations of gamma method error analysis suites in - [Fortran](https://gitlab.ift.uam-csic.es/alberto/aderrors) From 87ad43583012551f886b8d46705e7637568fb541 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 17 Oct 2021 13:50:34 +0100 Subject: [PATCH 61/76] CHANGELOG extended --- CHANGELOG.md | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 240778ea..8c9b7d41 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,11 +3,18 @@ All notable changes to this project will be documented in this file. ## [2.0.0] - 2021-??-?? -## Added -- `CObs` class added which can handle complex valued observables +### Added +- `CObs` class added which can handle complex valued Markov chain Monte Carlo data and the corresponding error propagation - Matrix to matrix operations like the matrix inverse now also work for complex matrices and matrices containing entries that are not `Obs` but `float` or `int` - `Obs` objects now have methods `is_zero` and `is_zero_within_error` +### Changed +- The default value for Corr.prange is now `None` +- The `input` module was restructured to contain one submodule per data source + +### Deprecated +- The function `plot_corrs` was deprecated as all its functionality is now contained within `Corr.show` + ## [1.1.0] - 2021-10-11 ### Added - `Corr` class added From eacec6bfdc8c8ab52326c925c85a5fcdf2d628c5 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 09:34:10 +0100 Subject: [PATCH 62/76] Obs.__float__ now always returns a python float --- pyerrors/pyerrors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 6039b412..2e14d2e0 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -476,7 +476,7 @@ class Obs: pickle.dump(self, fb) def __float__(self): - return self.value + return float(self.value) def __repr__(self): return 'Obs[' + str(self) + ']' From e73a99409caee1018f786bc1dd7c4f8756099b0f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 10:00:15 +0100 Subject: [PATCH 63/76] speed up for derived_observable when applied to 1d observables --- pyerrors/pyerrors.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 2e14d2e0..e1a5f7d9 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -758,8 +758,10 @@ def derived_observable(func, data, **kwargs): else: if new_shape[name] != tmp: raise Exception('Shapes of ensemble', name, 'do not match.') - - values = np.vectorize(lambda x: x.value)(data) + if data.ndim == 1: + values = np.array([o.value for o in data]) + else: + values = np.vectorize(lambda x: x.value)(data) new_values = func(values, **kwargs) From 4edea2c3a8ae923757e5b1722eb2524690f621cd Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 10:08:48 +0100 Subject: [PATCH 64/76] attribute e_Q removed --- CHANGELOG.md | 1 + pyerrors/misc.py | 1 + pyerrors/pyerrors.py | 16 +--------------- 3 files changed, 3 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c9b7d41..cb748e37 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ All notable changes to this project will be documented in this file. ### Deprecated - The function `plot_corrs` was deprecated as all its functionality is now contained within `Corr.show` +- Obs no longer have an attribute e_Q ## [1.1.0] - 2021-10-11 ### Added diff --git a/pyerrors/misc.py b/pyerrors/misc.py index a4d71a5f..f122fac2 100644 --- a/pyerrors/misc.py +++ b/pyerrors/misc.py @@ -56,6 +56,7 @@ def ks_test(obs=None): else: obs_list = obs + # TODO: Rework to apply to Q-values of all fits in memory Qs = [] for obs_i in obs_list: for ens in obs_i.e_names: diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index e1a5f7d9..4f8cc9bf 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -86,7 +86,6 @@ class Obs: self.e_tauint = {} self.e_dtauint = {} self.e_windowsize = {} - self.e_Q = {} self.e_rho = {} self.e_drho = {} self.e_n_tauint = {} @@ -298,19 +297,6 @@ class Obs: self.e_windowsize[e_name] = n break - if len(self.e_content[e_name]) > 1 and self.e_dvalue[e_name] > np.finfo(np.float64).eps: - e_mean = 0 - for r_name in self.e_content[e_name]: - e_mean += self.shape[r_name] * self.r_values[r_name] - e_mean /= e_N - xi2 = 0 - for r_name in self.e_content[e_name]: - xi2 += self.shape[r_name] * (self.r_values[r_name] - e_mean) ** 2 - xi2 /= self.e_dvalue[e_name] ** 2 * e_N - self.e_Q[e_name] = 1 - scipy.special.gammainc((len(self.e_content[e_name]) - 1.0) / 2.0, xi2 / 2.0) - else: - self.e_Q[e_name] = None - self.dvalue += self.e_dvalue[e_name] ** 2 self.ddvalue += (self.e_dvalue[e_name] * self.e_ddvalue[e_name]) ** 2 @@ -421,7 +407,7 @@ class Obs: for r, r_name in enumerate(self.e_content[e_name]): arr[r] = (self.r_values[r_name] - sub_r_mean) / (self.e_dvalue[e_name] * np.sqrt(e_N / self.shape[r_name] - 1)) plt.hist(arr, rwidth=0.8, bins=len(self.e_content[e_name])) - plt.title('Replica distribution' + e_name + ' (mean=0, var=1), Q=' + str(np.around(self.e_Q[e_name], decimals=2))) + plt.title('Replica distribution' + e_name + ' (mean=0, var=1)') plt.draw() def plot_history(self): From 4f7b5d23ae2dc83c399cbde69caaa61162151964 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 10:28:08 +0100 Subject: [PATCH 65/76] __slots__ added to Obs to reduce the memory footprint of individual objects and slightly speed up the initialisation of an object --- pyerrors/pyerrors.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 4f8cc9bf..55e40d48 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -35,6 +35,10 @@ class Obs: ensemble. N_sigma_global -- Standard value for N_sigma (default 1.0) """ + __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', 'value', 'dvalue', + 'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma', 'e_names', + 'e_content', 'e_dvalue', 'e_ddvalue', 'e_tauint', 'e_dtauint', + 'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint'] e_tag_global = 0 S_global = 2.0 From 8d7a5daafa3b842b2dbcb6d01663333c684644d4 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 10:33:31 +0100 Subject: [PATCH 66/76] redundant import removed, CHANGELOG updated --- CHANGELOG.md | 1 + pyerrors/pyerrors.py | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cb748e37..b2e4244c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ All notable changes to this project will be documented in this file. - `Obs` objects now have methods `is_zero` and `is_zero_within_error` ### Changed +- Additional attributes can no longer be added to existing Obs. This makes it no longer possible to import Obs created with previous versions of pyerrors - The default value for Corr.prange is now `None` - The `input` module was restructured to contain one submodule per data source diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 55e40d48..fa30c0f4 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -8,7 +8,6 @@ import autograd.numpy as anp # Thinly-wrapped numpy from autograd import jacobian import matplotlib.pyplot as plt import numdifftools as nd -import scipy.special class Obs: From 0c83f7a607d1c153148a163d14f62619b31ef993 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 13:07:33 +0100 Subject: [PATCH 67/76] Further comparison operations and test implemented --- pyerrors/pyerrors.py | 12 ++++++++++++ tests/test_pyerrors.py | 4 ++++ 2 files changed, 16 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index fa30c0f4..2da2ae28 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -485,9 +485,21 @@ class Obs: def __lt__(self, other): return self.value < other + def __le__(self, other): + return self.value <= other + def __gt__(self, other): return self.value > other + def __ge__(self, other): + return self.value >= other + + def __eq__(self, other): + return (self - other).is_zero() + + def __ne__(self, other): + return not (self - other).is_zero() + # Overload math operations def __add__(self, y): if isinstance(y, Obs): diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 982dcab3..5d7fa48a 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -26,6 +26,10 @@ def test_comparison(): test_obs2 = pe.pseudo_Obs(value2, 0.1, 't') assert (value1 > value2) == (test_obs1 > test_obs2) assert (value1 < value2) == (test_obs1 < test_obs2) + assert (value1 >= value2) == (test_obs1 >= test_obs2) + assert (value1 <= value2) == (test_obs1 <= test_obs2) + assert test_obs1 == test_obs1 + assert test_obs1 != test_obs2 def test_function_overloading(): From 9917778093f1896e2af0441254fe835a01409e48 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 13:15:20 +0100 Subject: [PATCH 68/76] Additional comparison asserts added to test_pyerrors.py --- tests/test_pyerrors.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 5d7fa48a..bed06530 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -28,8 +28,12 @@ def test_comparison(): assert (value1 < value2) == (test_obs1 < test_obs2) assert (value1 >= value2) == (test_obs1 >= test_obs2) assert (value1 <= value2) == (test_obs1 <= test_obs2) + assert test_obs1 >= test_obs1 + assert test_obs2 <= test_obs2 assert test_obs1 == test_obs1 + assert test_obs2 == test_obs2 assert test_obs1 != test_obs2 + assert test_obs2 != test_obs1 def test_function_overloading(): From f3b1a2bf3459abb1be8b56d600b56e920051b49e Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 13:28:47 +0100 Subject: [PATCH 69/76] tag attribute added to Obs --- pyerrors/pyerrors.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 2da2ae28..a1f18491 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -37,7 +37,8 @@ class Obs: __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', 'value', 'dvalue', 'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma', 'e_names', 'e_content', 'e_dvalue', 'e_ddvalue', 'e_tauint', 'e_dtauint', - 'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint'] + 'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint', + 'tag'] e_tag_global = 0 S_global = 2.0 @@ -94,6 +95,8 @@ class Obs: self.e_n_tauint = {} self.e_n_dtauint = {} + self.tag = None + def gamma_method(self, **kwargs): """Calculate the error and related properties of the Obs. From 040f6dbd51a627ae53796f8b7c1951e93cf11a67 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 13:31:57 +0100 Subject: [PATCH 70/76] Added tag attribute to CObs, added __slots__ to CObs, decreasing memory footprint --- pyerrors/pyerrors.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index a1f18491..0e8bfeb4 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -641,9 +641,11 @@ class Obs: class CObs: + __slots__ = ['real', 'imag', 'tag'] def __init__(self, real, imag=0.0): self.real = real self.imag = imag + self.tag = None def gamma_method(self, **kwargs): if isinstance(self.real, Obs): From f627a52bba608164b9f9b694840a2ab6aab1e54d Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 13:34:02 +0100 Subject: [PATCH 71/76] docstring added to CObs --- pyerrors/pyerrors.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 0e8bfeb4..f5e58494 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -641,7 +641,9 @@ class Obs: class CObs: + """Class for a complex valued observable.""" __slots__ = ['real', 'imag', 'tag'] + def __init__(self, real, imag=0.0): self.real = real self.imag = imag From 09b2525d2915a66772594fc818b41ef9ba937815 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 14:25:41 +0100 Subject: [PATCH 72/76] bias correction removed --- CHANGELOG.md | 15 ++++++++------- pyerrors/pyerrors.py | 10 +--------- 2 files changed, 9 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b2e4244c..259f61ef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,24 +9,25 @@ All notable changes to this project will be documented in this file. - `Obs` objects now have methods `is_zero` and `is_zero_within_error` ### Changed -- Additional attributes can no longer be added to existing Obs. This makes it no longer possible to import Obs created with previous versions of pyerrors -- The default value for Corr.prange is now `None` +- Additional attributes can no longer be added to existing `Obs`. This makes it no longer possible to import `Obs` created with previous versions of pyerrors +- The default value for `Corr.prange` is now `None` - The `input` module was restructured to contain one submodule per data source ### Deprecated - The function `plot_corrs` was deprecated as all its functionality is now contained within `Corr.show` -- Obs no longer have an attribute e_Q +- The kwarg `bias_correction` in `derived_observable` was removed +- Obs no longer have an attribute `e_Q` ## [1.1.0] - 2021-10-11 ### Added - `Corr` class added -- roots module added which can find the roots of a function that depends on Monte Carlo data via pyerrors Obs -- input/hadrons module added which can read hdf5 files written by [Hadrons](https://github.com/aportelli/Hadrons) -- read_rwms can now read reweighting factors in the format used by openQCD-2.0 +- `roots` module added which can find the roots of a function that depends on Monte Carlo data via pyerrors `Obs` +- `input/hadrons` module added which can read hdf5 files written by [Hadrons](https://github.com/aportelli/Hadrons) +- `read_rwms` can now read reweighting factors in the format used by openQCD-2.0 ## [1.0.1] - 2020-11-03 ### Fixed -- Bug in pyerrors.covariance fixed that appeared when working with several +- Bug in `pyerrors.covariance` fixed that appeared when working with several replica of different length. ## [1.0.0] - 2020-10-13 diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index f5e58494..a1ec8910 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -726,9 +726,6 @@ def derived_observable(func, data, **kwargs): man_grad -- manually supply a list or an array which contains the jacobian of func. Use cautiously, supplying the wrong derivative will not be intercepted. - bias_correction -- if True, the bias correction specified in - hep-lat/0306017 eq. (19) is performed, not recommended. - (Only applicable for more than 1 replicum) Notes ----- @@ -831,12 +828,7 @@ def derived_observable(func, data, **kwargs): new_samples.append(new_deltas[name] + new_r_values[name][i_val]) final_result[i_val] = Obs(new_samples, new_names) - - # Bias correction - if replicas > 1 and kwargs.get('bias_correction'): - final_result[i_val].value = (replicas * new_val - final_result[i_val].value) / (replicas - 1) - else: - final_result[i_val].value = new_val + final_result[i_val].value = new_val if multi == 0: final_result = final_result.item() From 526dceef1cd7cfd474727fc4f0b7dcbb18e61e9f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 14:26:34 +0100 Subject: [PATCH 73/76] derived_observable cleaned up --- pyerrors/pyerrors.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index a1ec8910..555c73d8 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -751,7 +751,6 @@ def derived_observable(func, data, **kwargs): n_obs = len(raveled_data) new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x])) - replicas = len(new_names) new_shape = {} for i_data in raveled_data: From 7b0645bda0cd51cb08821c6c48320cdc77e20ca8 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 14:37:09 +0100 Subject: [PATCH 74/76] Complex matrix operations now work with ndarrays instead of numpy matrices --- pyerrors/linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index a8810091..0394af8e 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -95,7 +95,7 @@ def mat_mat_op(op, obs, **kwargs): else: A[n, m] = entry B[n, m] = 0.0 - big_matrix = np.bmat([[A, -B], [B, A]]) + big_matrix = np.block([[A, -B], [B, A]]) if kwargs.get('num_grad') is True: op_big_matrix = _num_diff_mat_mat_op(op, big_matrix, **kwargs) else: From b1e5777586371deaa191071b5e94152faeba9061 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 14:48:16 +0100 Subject: [PATCH 75/76] conjugate method added to CObs class --- pyerrors/pyerrors.py | 5 ++++- tests/test_pyerrors.py | 4 ++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 555c73d8..5d54b9af 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -655,6 +655,9 @@ class CObs: if isinstance(self.imag, Obs): self.imag.gamma_method(**kwargs) + def conjugate(self): + return CObs(self.real, -self.imag) + def __add__(self, other): if hasattr(other, 'real') and hasattr(other, 'imag'): return CObs(self.real + other.real, @@ -701,7 +704,7 @@ class CObs: return self.real == other.real and self.imag == other.imag def __str__(self): - return '(' + str(self.real) + int(self.imag > - np.finfo(np.float64).eps) * '+' + str(self.imag) + 'j)' + return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)' def __repr__(self): return 'CObs[' + str(self) + ']' diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index bed06530..58ec21f3 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -199,6 +199,10 @@ def test_cobs(): obs2 = pe.pseudo_Obs(-0.2, 0.03, 't') my_cobs = pe.CObs(obs1, obs2) + assert not (my_cobs + my_cobs.conjugate()).real.is_zero() + assert (my_cobs + my_cobs.conjugate()).imag.is_zero() + assert (my_cobs - my_cobs.conjugate()).real.is_zero() + assert not (my_cobs - my_cobs.conjugate()).imag.is_zero() np.abs(my_cobs) fs = [[lambda x: x[0] + x[1], lambda x: x[1] + x[0]], [lambda x: x[0] * x[1], lambda x: x[1] * x[0]]] From 877b1a4467f9d98ab5728d899217353b46ecca78 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 18 Oct 2021 17:36:49 +0100 Subject: [PATCH 76/76] Complex matrix multiplication speed up my explicitly coding the derivatives --- pyerrors/pyerrors.py | 12 +++++++++++- tests/test_pyerrors.py | 4 ++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py index 5d54b9af..84a41f17 100644 --- a/pyerrors/pyerrors.py +++ b/pyerrors/pyerrors.py @@ -655,6 +655,9 @@ class CObs: if isinstance(self.imag, Obs): self.imag.gamma_method(**kwargs) + def is_zero(self): + return self.real == 0.0 and self.imag == 0.0 + def conjugate(self): return CObs(self.real, -self.imag) @@ -678,7 +681,14 @@ class CObs: return -1 * (self - other) def __mul__(self, other): - if hasattr(other, 'real') and hasattr(other, 'imag'): + if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]): + return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3], + [self.real, other.real, self.imag, other.imag], + man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]), + derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3], + [self.real, other.real, self.imag, other.imag], + man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value])) + elif hasattr(other, 'real') and hasattr(other, 'imag'): return CObs(self.real * other.real - self.imag * other.imag, self.imag * other.real + self.real * other.imag) else: diff --git a/tests/test_pyerrors.py b/tests/test_pyerrors.py index 58ec21f3..7d304464 100644 --- a/tests/test_pyerrors.py +++ b/tests/test_pyerrors.py @@ -204,6 +204,10 @@ def test_cobs(): assert (my_cobs - my_cobs.conjugate()).real.is_zero() assert not (my_cobs - my_cobs.conjugate()).imag.is_zero() np.abs(my_cobs) + + assert (my_cobs * my_cobs / my_cobs - my_cobs).is_zero() + assert (my_cobs + my_cobs - 2 * my_cobs).is_zero() + fs = [[lambda x: x[0] + x[1], lambda x: x[1] + x[0]], [lambda x: x[0] * x[1], lambda x: x[1] * x[0]]] for other in [1, 1.1, (1.1-0.2j), pe.CObs(obs1), pe.CObs(obs1, obs2)]: