pyerrors.input.misc

  1import fnmatch
  2import os
  3import re
  4import struct
  5import warnings
  6
  7import matplotlib.pyplot as plt
  8import numpy as np  # Thinly-wrapped numpy
  9from matplotlib import gridspec
 10
 11from ..fits import fit_lin
 12from ..obs import Obs
 13
 14
 15def fit_t0(t2E_dict, fit_range, plot_fit=False, observable='t0'):
 16    """Compute the root of (flow-based) data based on a dictionary that contains
 17    the necessary information in key-value pairs a la (flow time: observable at flow time).
 18
 19    It is assumed that the data is monotonically increasing and passes zero from below.
 20    No exception is thrown if this is not the case (several roots, no monotonic increase).
 21    An exception is thrown if no root can be found in the data.
 22
 23    A linear fit in the vicinity of the root is performed to exctract the root from the
 24    two fit parameters.
 25
 26    Parameters
 27    ----------
 28    t2E_dict : dict
 29        Dictionary with pairs of (flow time: observable at flow time) where the flow times
 30        are of type float and the observables of type Obs.
 31    fit_range : int
 32        Number of data points left and right of the zero
 33        crossing to be included in the linear fit.
 34    plot_fit : bool
 35        If true, the fit for the extraction of t0 is shown together with the data. (Default: False)
 36    observable: str
 37        Keyword to identify the observable to print the correct ylabel (if plot_fit is True)
 38        for the observables 't0' and 'w0'. No y label is printed otherwise. (Default: 't0')
 39
 40    Returns
 41    -------
 42    root : Obs
 43        The root of the data series.
 44    """
 45
 46    zero_crossing = np.argmax(np.array(
 47        [o.value for o in t2E_dict.values()]) > 0.0)
 48
 49    if zero_crossing == 0:
 50        raise Exception('Desired flow time not in data')
 51
 52    x = list(t2E_dict.keys())[zero_crossing - fit_range:
 53                              zero_crossing + fit_range]
 54    y = list(t2E_dict.values())[zero_crossing - fit_range:
 55                                zero_crossing + fit_range]
 56    [o.gamma_method() for o in y]
 57
 58    if len(x) < 2 * fit_range:
 59        warnings.warn(f'Fit range smaller than expected! Fitting from {x[0]:1.2e} to {x[-1]:1.2e}', stacklevel=2)
 60
 61    fit_result = fit_lin(x, y)
 62
 63    if plot_fit is True:
 64        plt.figure()
 65        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
 66        ax0 = plt.subplot(gs[0])
 67        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
 68        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
 69        [o.gamma_method() for o in ymore]
 70        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
 71        xplot = np.linspace(np.min(x), np.max(x))
 72        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
 73        [yi.gamma_method() for yi in yplot]
 74        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
 75        retval = (-fit_result[0] / fit_result[1])
 76        retval.gamma_method()
 77        ylim = ax0.get_ylim()
 78        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
 79        ax0.set_ylim(ylim)
 80        if observable == 't0':
 81            ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
 82        elif observable == 'w0':
 83            ax0.set_ylabel(r'$t d(t^2 \langle E(t) \rangle)/dt - 0.3 $')
 84        xlim = ax0.get_xlim()
 85
 86        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
 87        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
 88        ax1 = plt.subplot(gs[1])
 89        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
 90        ax1.tick_params(direction='out')
 91        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
 92        ax1.axhline(y=0.0, ls='--', color='k')
 93        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
 94        ax1.set_xlim(xlim)
 95        ax1.set_ylabel('Residuals')
 96        ax1.set_xlabel(r'$t/a^2$')
 97
 98        plt.draw()
 99    return -fit_result[0] / fit_result[1]
100
101
102def read_pbp(path, prefix, **kwargs):
103    """Read pbp format from given folder structure.
104
105    Parameters
106    ----------
107    r_start : list
108        list which contains the first config to be read for each replicum
109    r_stop : list
110        list which contains the last config to be read for each replicum
111
112    Returns
113    -------
114    result : list[Obs]
115        list of observables read
116    """
117
118    ls = []
119    for (_dirpath, _dirnames, filenames) in os.walk(path):
120        ls.extend(filenames)
121        break
122
123    if not ls:
124        raise Exception('Error, directory not found')
125
126    # Exclude files with different names
127    for exc in ls:
128        if not fnmatch.fnmatch(exc, prefix + '*.dat'):
129            ls = list(set(ls) - set([exc]))
130    if len(ls) > 1:
131        ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
132    replica = len(ls)
133
134    if 'r_start' in kwargs:
135        r_start = kwargs.get('r_start')
136        if len(r_start) != replica:
137            raise Exception('r_start does not match number of replicas')
138        # Adjust Configuration numbering to python index
139        r_start = [o - 1 if o else None for o in r_start]
140    else:
141        r_start = [None] * replica
142
143    if 'r_stop' in kwargs:
144        r_stop = kwargs.get('r_stop')
145        if len(r_stop) != replica:
146            raise Exception('r_stop does not match number of replicas')
147    else:
148        r_stop = [None] * replica
149
150    print(r'Read <bar{psi}\psi> from', prefix[:-1], ',', replica, 'replica', end='')
151
152    print_err = 0
153    if 'print_err' in kwargs:
154        print_err = 1
155        print()
156
157    deltas = []
158
159    for rep in range(replica):
160        tmp_array = []
161        with open(path + '/' + ls[rep], 'rb') as fp:
162
163            t = fp.read(4)  # number of reweighting factors
164            if rep == 0:
165                nrw = struct.unpack('i', t)[0]
166                for _ in range(nrw):
167                    deltas.append([])
168            else:
169                if nrw != struct.unpack('i', t)[0]:
170                    raise Exception('Error: different number of factors for replicum', rep)
171
172            for _ in range(nrw):
173                tmp_array.append([])
174
175            # This block is necessary for openQCD1.6 ms1 files
176            nfct = []
177            for _ in range(nrw):
178                t = fp.read(4)
179                nfct.append(struct.unpack('i', t)[0])
180            print('nfct: ', nfct)  # Hasenbusch factor, 1 for rat reweighting
181
182            nsrc = []
183            for _ in range(nrw):
184                t = fp.read(4)
185                nsrc.append(struct.unpack('i', t)[0])
186
187            # body
188            while True:
189                t = fp.read(4)
190                if len(t) < 4:
191                    break
192                if print_err:
193                    config_no = struct.unpack('i', t)
194                for i in range(nrw):
195                    tmp_nfct = 1.0
196                    for j in range(nfct[i]):
197                        t = fp.read(8 * nsrc[i])
198                        t = fp.read(8 * nsrc[i])
199                        tmp_rw = struct.unpack('d' * nsrc[i], t)
200                        tmp_nfct *= np.mean(np.asarray(tmp_rw))
201                        if print_err:
202                            print(config_no, i, j, np.mean(np.asarray(tmp_rw)), np.std(np.asarray(tmp_rw)))
203                            print('Sources:', np.asarray(tmp_rw))
204                            print('Partial factor:', tmp_nfct)
205                    tmp_array[i].append(tmp_nfct)
206
207            for k in range(nrw):
208                deltas[k].append(tmp_array[k][r_start[rep]:r_stop[rep]])
209
210    rep_names = []
211    for entry in ls:
212        truncated_entry = entry.split('.')[0]
213        idx = truncated_entry.index('r')
214        rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
215    print(',', nrw, r'<bar{psi}\psi> with', nsrc, 'sources')
216    result = []
217    for t in range(nrw):
218        result.append(Obs(deltas[t], rep_names))
219
220    return result