pyerrors.input.misc

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