pyerrors.obs

   1import warnings
   2import hashlib
   3import pickle
   4from math import gcd
   5from functools import reduce
   6import numpy as np
   7import autograd.numpy as anp  # Thinly-wrapped numpy
   8from autograd import jacobian
   9import matplotlib.pyplot as plt
  10from scipy.stats import skew, skewtest, kurtosis, kurtosistest
  11import numdifftools as nd
  12from itertools import groupby
  13from .covobs import Covobs
  14
  15# Improve print output of numpy.ndarrays containing Obs objects.
  16np.set_printoptions(formatter={'object': lambda x: str(x)})
  17
  18
  19class Obs:
  20    """Class for a general observable.
  21
  22    Instances of Obs are the basic objects of a pyerrors error analysis.
  23    They are initialized with a list which contains arrays of samples for
  24    different ensembles/replica and another list of same length which contains
  25    the names of the ensembles/replica. Mathematical operations can be
  26    performed on instances. The result is another instance of Obs. The error of
  27    an instance can be computed with the gamma_method. Also contains additional
  28    methods for output and visualization of the error calculation.
  29
  30    Attributes
  31    ----------
  32    S_global : float
  33        Standard value for S (default 2.0)
  34    S_dict : dict
  35        Dictionary for S values. If an entry for a given ensemble
  36        exists this overwrites the standard value for that ensemble.
  37    tau_exp_global : float
  38        Standard value for tau_exp (default 0.0)
  39    tau_exp_dict : dict
  40        Dictionary for tau_exp values. If an entry for a given ensemble exists
  41        this overwrites the standard value for that ensemble.
  42    N_sigma_global : float
  43        Standard value for N_sigma (default 1.0)
  44    N_sigma_dict : dict
  45        Dictionary for N_sigma values. If an entry for a given ensemble exists
  46        this overwrites the standard value for that ensemble.
  47    """
  48    __slots__ = ['names', 'shape', 'r_values', 'deltas', 'N', '_value', '_dvalue',
  49                 'ddvalue', 'reweighted', 'S', 'tau_exp', 'N_sigma',
  50                 'e_dvalue', 'e_ddvalue', 'e_tauint', 'e_dtauint',
  51                 'e_windowsize', 'e_rho', 'e_drho', 'e_n_tauint', 'e_n_dtauint',
  52                 'idl', 'tag', '_covobs', '__dict__']
  53
  54    S_global = 2.0
  55    S_dict = {}
  56    tau_exp_global = 0.0
  57    tau_exp_dict = {}
  58    N_sigma_global = 1.0
  59    N_sigma_dict = {}
  60
  61    def __init__(self, samples, names, idl=None, **kwargs):
  62        """ Initialize Obs object.
  63
  64        Parameters
  65        ----------
  66        samples : list
  67            list of numpy arrays containing the Monte Carlo samples
  68        names : list
  69            list of strings labeling the individual samples
  70        idl : list, optional
  71            list of ranges or lists on which the samples are defined
  72        """
  73
  74        if kwargs.get("means") is None and len(samples):
  75            if len(samples) != len(names):
  76                raise Exception('Length of samples and names incompatible.')
  77            if idl is not None:
  78                if len(idl) != len(names):
  79                    raise Exception('Length of idl incompatible with samples and names.')
  80            name_length = len(names)
  81            if name_length > 1:
  82                if name_length != len(set(names)):
  83                    raise Exception('names are not unique.')
  84                if not all(isinstance(x, str) for x in names):
  85                    raise TypeError('All names have to be strings.')
  86            else:
  87                if not isinstance(names[0], str):
  88                    raise TypeError('All names have to be strings.')
  89            if min(len(x) for x in samples) <= 4:
  90                raise Exception('Samples have to have at least 5 entries.')
  91
  92        self.names = sorted(names)
  93        self.shape = {}
  94        self.r_values = {}
  95        self.deltas = {}
  96        self._covobs = {}
  97
  98        self._value = 0
  99        self.N = 0
 100        self.idl = {}
 101        if idl is not None:
 102            for name, idx in sorted(zip(names, idl)):
 103                if isinstance(idx, range):
 104                    self.idl[name] = idx
 105                elif isinstance(idx, (list, np.ndarray)):
 106                    dc = np.unique(np.diff(idx))
 107                    if np.any(dc < 0):
 108                        raise Exception("Unsorted idx for idl[%s]" % (name))
 109                    if len(dc) == 1:
 110                        self.idl[name] = range(idx[0], idx[-1] + dc[0], dc[0])
 111                    else:
 112                        self.idl[name] = list(idx)
 113                else:
 114                    raise Exception('incompatible type for idl[%s].' % (name))
 115        else:
 116            for name, sample in sorted(zip(names, samples)):
 117                self.idl[name] = range(1, len(sample) + 1)
 118
 119        if kwargs.get("means") is not None:
 120            for name, sample, mean in sorted(zip(names, samples, kwargs.get("means"))):
 121                self.shape[name] = len(self.idl[name])
 122                self.N += self.shape[name]
 123                self.r_values[name] = mean
 124                self.deltas[name] = sample
 125        else:
 126            for name, sample in sorted(zip(names, samples)):
 127                self.shape[name] = len(self.idl[name])
 128                self.N += self.shape[name]
 129                if len(sample) != self.shape[name]:
 130                    raise Exception('Incompatible samples and idx for %s: %d vs. %d' % (name, len(sample), self.shape[name]))
 131                self.r_values[name] = np.mean(sample)
 132                self.deltas[name] = sample - self.r_values[name]
 133                self._value += self.shape[name] * self.r_values[name]
 134            self._value /= self.N
 135
 136        self._dvalue = 0.0
 137        self.ddvalue = 0.0
 138        self.reweighted = False
 139
 140        self.tag = None
 141
 142    @property
 143    def value(self):
 144        return self._value
 145
 146    @property
 147    def dvalue(self):
 148        return self._dvalue
 149
 150    @property
 151    def e_names(self):
 152        return sorted(set([o.split('|')[0] for o in self.names]))
 153
 154    @property
 155    def cov_names(self):
 156        return sorted(set([o for o in self.covobs.keys()]))
 157
 158    @property
 159    def mc_names(self):
 160        return sorted(set([o.split('|')[0] for o in self.names if o not in self.cov_names]))
 161
 162    @property
 163    def e_content(self):
 164        res = {}
 165        for e, e_name in enumerate(self.e_names):
 166            res[e_name] = sorted(filter(lambda x: x.startswith(e_name + '|'), self.names))
 167            if e_name in self.names:
 168                res[e_name].append(e_name)
 169        return res
 170
 171    @property
 172    def covobs(self):
 173        return self._covobs
 174
 175    def gamma_method(self, **kwargs):
 176        """Estimate the error and related properties of the Obs.
 177
 178        Parameters
 179        ----------
 180        S : float
 181            specifies a custom value for the parameter S (default 2.0).
 182            If set to 0 it is assumed that the data exhibits no
 183            autocorrelation. In this case the error estimates coincides
 184            with the sample standard error.
 185        tau_exp : float
 186            positive value triggers the critical slowing down analysis
 187            (default 0.0).
 188        N_sigma : float
 189            number of standard deviations from zero until the tail is
 190            attached to the autocorrelation function (default 1).
 191        fft : bool
 192            determines whether the fft algorithm is used for the computation
 193            of the autocorrelation function (default True)
 194        """
 195
 196        e_content = self.e_content
 197        self.e_dvalue = {}
 198        self.e_ddvalue = {}
 199        self.e_tauint = {}
 200        self.e_dtauint = {}
 201        self.e_windowsize = {}
 202        self.e_n_tauint = {}
 203        self.e_n_dtauint = {}
 204        e_gamma = {}
 205        self.e_rho = {}
 206        self.e_drho = {}
 207        self._dvalue = 0
 208        self.ddvalue = 0
 209
 210        self.S = {}
 211        self.tau_exp = {}
 212        self.N_sigma = {}
 213
 214        if kwargs.get('fft') is False:
 215            fft = False
 216        else:
 217            fft = True
 218
 219        def _parse_kwarg(kwarg_name):
 220            if kwarg_name in kwargs:
 221                tmp = kwargs.get(kwarg_name)
 222                if isinstance(tmp, (int, float)):
 223                    if tmp < 0:
 224                        raise Exception(kwarg_name + ' has to be larger or equal to 0.')
 225                    for e, e_name in enumerate(self.e_names):
 226                        getattr(self, kwarg_name)[e_name] = tmp
 227                else:
 228                    raise TypeError(kwarg_name + ' is not in proper format.')
 229            else:
 230                for e, e_name in enumerate(self.e_names):
 231                    if e_name in getattr(Obs, kwarg_name + '_dict'):
 232                        getattr(self, kwarg_name)[e_name] = getattr(Obs, kwarg_name + '_dict')[e_name]
 233                    else:
 234                        getattr(self, kwarg_name)[e_name] = getattr(Obs, kwarg_name + '_global')
 235
 236        _parse_kwarg('S')
 237        _parse_kwarg('tau_exp')
 238        _parse_kwarg('N_sigma')
 239
 240        for e, e_name in enumerate(self.mc_names):
 241            r_length = []
 242            for r_name in e_content[e_name]:
 243                if isinstance(self.idl[r_name], range):
 244                    r_length.append(len(self.idl[r_name]))
 245                else:
 246                    r_length.append((self.idl[r_name][-1] - self.idl[r_name][0] + 1))
 247
 248            e_N = np.sum([self.shape[r_name] for r_name in e_content[e_name]])
 249            w_max = max(r_length) // 2
 250            e_gamma[e_name] = np.zeros(w_max)
 251            self.e_rho[e_name] = np.zeros(w_max)
 252            self.e_drho[e_name] = np.zeros(w_max)
 253
 254            for r_name in e_content[e_name]:
 255                e_gamma[e_name] += self._calc_gamma(self.deltas[r_name], self.idl[r_name], self.shape[r_name], w_max, fft)
 256
 257            gamma_div = np.zeros(w_max)
 258            for r_name in e_content[e_name]:
 259                gamma_div += self._calc_gamma(np.ones((self.shape[r_name])), self.idl[r_name], self.shape[r_name], w_max, fft)
 260            gamma_div[gamma_div < 1] = 1.0
 261            e_gamma[e_name] /= gamma_div[:w_max]
 262
 263            if np.abs(e_gamma[e_name][0]) < 10 * np.finfo(float).tiny:  # Prevent division by zero
 264                self.e_tauint[e_name] = 0.5
 265                self.e_dtauint[e_name] = 0.0
 266                self.e_dvalue[e_name] = 0.0
 267                self.e_ddvalue[e_name] = 0.0
 268                self.e_windowsize[e_name] = 0
 269                continue
 270
 271            gaps = []
 272            for r_name in e_content[e_name]:
 273                if isinstance(self.idl[r_name], range):
 274                    gaps.append(1)
 275                else:
 276                    gaps.append(np.min(np.diff(self.idl[r_name])))
 277
 278            if not np.all([gi == gaps[0] for gi in gaps]):
 279                raise Exception(f"Replica for ensemble {e_name} are not equally spaced.", gaps)
 280            else:
 281                gapsize = gaps[0]
 282
 283            self.e_rho[e_name] = e_gamma[e_name][:w_max] / e_gamma[e_name][0]
 284            self.e_n_tauint[e_name] = np.cumsum(np.concatenate(([0.5], self.e_rho[e_name][1:])))
 285            # Make sure no entry of tauint is smaller than 0.5
 286            self.e_n_tauint[e_name][self.e_n_tauint[e_name] <= 0.5] = 0.5 + np.finfo(np.float64).eps
 287            # hep-lat/0306017 eq. (42)
 288            self.e_n_dtauint[e_name] = self.e_n_tauint[e_name] * 2 * np.sqrt(np.abs(np.arange(w_max) / gapsize + 0.5 - self.e_n_tauint[e_name]) / e_N)
 289            self.e_n_dtauint[e_name][0] = 0.0
 290
 291            def _compute_drho(i):
 292                tmp = self.e_rho[e_name][i + 1:w_max] + np.concatenate([self.e_rho[e_name][i - 1::-1], self.e_rho[e_name][1:w_max - 2 * i]]) - 2 * self.e_rho[e_name][i] * self.e_rho[e_name][1:w_max - i]
 293                self.e_drho[e_name][i] = np.sqrt(np.sum(tmp ** 2) / e_N)
 294
 295            if self.tau_exp[e_name] > 0:
 296                _compute_drho(gapsize)
 297                texp = self.tau_exp[e_name]
 298                # Critical slowing down analysis
 299                if w_max // 2 <= 1:
 300                    raise Exception("Need at least 8 samples for tau_exp error analysis")
 301                for n in range(gapsize, w_max // 2, gapsize):
 302                    _compute_drho(n + gapsize)
 303                    if (self.e_rho[e_name][n] - self.N_sigma[e_name] * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:
 304                        # Bias correction hep-lat/0306017 eq. (49) included
 305                        self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 1) / e_N) / (1 + 1 / e_N) + texp * np.abs(self.e_rho[e_name][n + 1])  # The absolute makes sure, that the tail contribution is always positive
 306                        self.e_dtauint[e_name] = np.sqrt(self.e_n_dtauint[e_name][n] ** 2 + texp ** 2 * self.e_drho[e_name][n + 1] ** 2)
 307                        # Error of tau_exp neglected so far, missing term: self.e_rho[e_name][n + 1] ** 2 * d_tau_exp ** 2
 308                        self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)
 309                        self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n / gapsize + 0.5) / e_N)
 310                        self.e_windowsize[e_name] = n
 311                        break
 312            else:
 313                if self.S[e_name] == 0.0:
 314                    self.e_tauint[e_name] = 0.5
 315                    self.e_dtauint[e_name] = 0.0
 316                    self.e_dvalue[e_name] = np.sqrt(e_gamma[e_name][0] / (e_N - 1))
 317                    self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt(0.5 / e_N)
 318                    self.e_windowsize[e_name] = 0
 319                else:
 320                    # Standard automatic windowing procedure
 321                    tau = self.S[e_name] / np.log((2 * self.e_n_tauint[e_name][gapsize::gapsize] + 1) / (2 * self.e_n_tauint[e_name][gapsize::gapsize] - 1))
 322                    g_w = np.exp(- np.arange(1, len(tau) + 1) / tau) - tau / np.sqrt(np.arange(1, len(tau) + 1) * e_N)
 323                    for n in range(1, w_max):
 324                        if g_w[n - 1] < 0 or n >= w_max - 1:
 325                            _compute_drho(gapsize * n)
 326                            n *= gapsize
 327                            self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 1) / e_N) / (1 + 1 / e_N)  # Bias correction hep-lat/0306017 eq. (49)
 328                            self.e_dtauint[e_name] = self.e_n_dtauint[e_name][n]
 329                            self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)
 330                            self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n / gapsize + 0.5) / e_N)
 331                            self.e_windowsize[e_name] = n
 332                            break
 333
 334            self._dvalue += self.e_dvalue[e_name] ** 2
 335            self.ddvalue += (self.e_dvalue[e_name] * self.e_ddvalue[e_name]) ** 2
 336
 337        for e_name in self.cov_names:
 338            self.e_dvalue[e_name] = np.sqrt(self.covobs[e_name].errsq())
 339            self.e_ddvalue[e_name] = 0
 340            self._dvalue += self.e_dvalue[e_name]**2
 341
 342        self._dvalue = np.sqrt(self._dvalue)
 343        if self._dvalue == 0.0:
 344            self.ddvalue = 0.0
 345        else:
 346            self.ddvalue = np.sqrt(self.ddvalue) / self._dvalue
 347        return
 348
 349    gm = gamma_method
 350
 351    def _calc_gamma(self, deltas, idx, shape, w_max, fft):
 352        """Calculate Gamma_{AA} from the deltas, which are defined on idx.
 353           idx is assumed to be a contiguous range (possibly with a stepsize != 1)
 354
 355        Parameters
 356        ----------
 357        deltas : list
 358            List of fluctuations
 359        idx : list
 360            List or range of configurations on which the deltas are defined.
 361        shape : int
 362            Number of configurations in idx.
 363        w_max : int
 364            Upper bound for the summation window.
 365        fft : bool
 366            determines whether the fft algorithm is used for the computation
 367            of the autocorrelation function.
 368        """
 369        gamma = np.zeros(w_max)
 370        deltas = _expand_deltas(deltas, idx, shape)
 371        new_shape = len(deltas)
 372        if fft:
 373            max_gamma = min(new_shape, w_max)
 374            # The padding for the fft has to be even
 375            padding = new_shape + max_gamma + (new_shape + max_gamma) % 2
 376            gamma[:max_gamma] += np.fft.irfft(np.abs(np.fft.rfft(deltas, padding)) ** 2)[:max_gamma]
 377        else:
 378            for n in range(w_max):
 379                if new_shape - n >= 0:
 380                    gamma[n] += deltas[0:new_shape - n].dot(deltas[n:new_shape])
 381
 382        return gamma
 383
 384    def details(self, ens_content=True):
 385        """Output detailed properties of the Obs.
 386
 387        Parameters
 388        ----------
 389        ens_content : bool
 390            print details about the ensembles and replica if true.
 391        """
 392        if self.tag is not None:
 393            print("Description:", self.tag)
 394        if not hasattr(self, 'e_dvalue'):
 395            print('Result\t %3.8e' % (self.value))
 396        else:
 397            if self.value == 0.0:
 398                percentage = np.nan
 399            else:
 400                percentage = np.abs(self._dvalue / self.value) * 100
 401            print('Result\t %3.8e +/- %3.8e +/- %3.8e (%3.3f%%)' % (self.value, self._dvalue, self.ddvalue, percentage))
 402            if len(self.e_names) > 1:
 403                print(' Ensemble errors:')
 404            e_content = self.e_content
 405            for e_name in self.mc_names:
 406                if isinstance(self.idl[e_content[e_name][0]], range):
 407                    gap = self.idl[e_content[e_name][0]].step
 408                else:
 409                    gap = np.min(np.diff(self.idl[e_content[e_name][0]]))
 410
 411                if len(self.e_names) > 1:
 412                    print('', e_name, '\t %3.6e +/- %3.6e' % (self.e_dvalue[e_name], self.e_ddvalue[e_name]))
 413                tau_string = " \N{GREEK SMALL LETTER TAU}_int\t " + _format_uncertainty(self.e_tauint[e_name], self.e_dtauint[e_name])
 414                tau_string += f" in units of {gap} config"
 415                if gap > 1:
 416                    tau_string += "s"
 417                if self.tau_exp[e_name] > 0:
 418                    tau_string = f"{tau_string: <45}" + '\t(\N{GREEK SMALL LETTER TAU}_exp=%3.2f, N_\N{GREEK SMALL LETTER SIGMA}=%1.0i)' % (self.tau_exp[e_name], self.N_sigma[e_name])
 419                else:
 420                    tau_string = f"{tau_string: <45}" + '\t(S=%3.2f)' % (self.S[e_name])
 421                print(tau_string)
 422            for e_name in self.cov_names:
 423                print('', e_name, '\t %3.8e' % (self.e_dvalue[e_name]))
 424        if ens_content is True:
 425            if len(self.e_names) == 1:
 426                print(self.N, 'samples in', len(self.e_names), 'ensemble:')
 427            else:
 428                print(self.N, 'samples in', len(self.e_names), 'ensembles:')
 429            my_string_list = []
 430            for key, value in sorted(self.e_content.items()):
 431                if key not in self.covobs:
 432                    my_string = '  ' + "\u00B7 Ensemble '" + key + "' "
 433                    if len(value) == 1:
 434                        my_string += f': {self.shape[value[0]]} configurations'
 435                        if isinstance(self.idl[value[0]], range):
 436                            my_string += f' (from {self.idl[value[0]].start} to {self.idl[value[0]][-1]}' + int(self.idl[value[0]].step != 1) * f' in steps of {self.idl[value[0]].step}' + ')'
 437                        else:
 438                            my_string += f' (irregular range from {self.idl[value[0]][0]} to {self.idl[value[0]][-1]})'
 439                    else:
 440                        sublist = []
 441                        for v in value:
 442                            my_substring = '    ' + "\u00B7 Replicum '" + v[len(key) + 1:] + "' "
 443                            my_substring += f': {self.shape[v]} configurations'
 444                            if isinstance(self.idl[v], range):
 445                                my_substring += f' (from {self.idl[v].start} to {self.idl[v][-1]}' + int(self.idl[v].step != 1) * f' in steps of {self.idl[v].step}' + ')'
 446                            else:
 447                                my_substring += f' (irregular range from {self.idl[v][0]} to {self.idl[v][-1]})'
 448                            sublist.append(my_substring)
 449
 450                        my_string += '\n' + '\n'.join(sublist)
 451                else:
 452                    my_string = '  ' + "\u00B7 Covobs   '" + key + "' "
 453                my_string_list.append(my_string)
 454            print('\n'.join(my_string_list))
 455
 456    def reweight(self, weight):
 457        """Reweight the obs with given rewighting factors.
 458
 459        Parameters
 460        ----------
 461        weight : Obs
 462            Reweighting factor. An Observable that has to be defined on a superset of the
 463            configurations in obs[i].idl for all i.
 464        all_configs : bool
 465            if True, the reweighted observables are normalized by the average of
 466            the reweighting factor on all configurations in weight.idl and not
 467            on the configurations in obs[i].idl. Default False.
 468        """
 469        return reweight(weight, [self])[0]
 470
 471    def is_zero_within_error(self, sigma=1):
 472        """Checks whether the observable is zero within 'sigma' standard errors.
 473
 474        Parameters
 475        ----------
 476        sigma : int
 477            Number of standard errors used for the check.
 478
 479        Works only properly when the gamma method was run.
 480        """
 481        return self.is_zero() or np.abs(self.value) <= sigma * self._dvalue
 482
 483    def is_zero(self, atol=1e-10):
 484        """Checks whether the observable is zero within a given tolerance.
 485
 486        Parameters
 487        ----------
 488        atol : float
 489            Absolute tolerance (for details see numpy documentation).
 490        """
 491        return np.isclose(0.0, self.value, 1e-14, atol) and all(np.allclose(0.0, delta, 1e-14, atol) for delta in self.deltas.values()) and all(np.allclose(0.0, delta.errsq(), 1e-14, atol) for delta in self.covobs.values())
 492
 493    def plot_tauint(self, save=None):
 494        """Plot integrated autocorrelation time for each ensemble.
 495
 496        Parameters
 497        ----------
 498        save : str
 499            saves the figure to a file named 'save' if.
 500        """
 501        if not hasattr(self, 'e_dvalue'):
 502            raise Exception('Run the gamma method first.')
 503
 504        for e, e_name in enumerate(self.mc_names):
 505            fig = plt.figure()
 506            plt.xlabel(r'$W$')
 507            plt.ylabel(r'$\tau_\mathrm{int}$')
 508            length = int(len(self.e_n_tauint[e_name]))
 509            if self.tau_exp[e_name] > 0:
 510                base = self.e_n_tauint[e_name][self.e_windowsize[e_name]]
 511                x_help = np.arange(2 * self.tau_exp[e_name])
 512                y_help = (x_help + 1) * np.abs(self.e_rho[e_name][self.e_windowsize[e_name] + 1]) * (1 - x_help / (2 * (2 * self.tau_exp[e_name] - 1))) + base
 513                x_arr = np.arange(self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name])
 514                plt.plot(x_arr, y_help, 'C' + str(e), linewidth=1, ls='--', marker=',')
 515                plt.errorbar([self.e_windowsize[e_name] + 2 * self.tau_exp[e_name]], [self.e_tauint[e_name]],
 516                             yerr=[self.e_dtauint[e_name]], fmt='C' + str(e), linewidth=1, capsize=2, marker='o', mfc=plt.rcParams['axes.facecolor'])
 517                xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
 518                label = e_name + r', $\tau_\mathrm{exp}$=' + str(np.around(self.tau_exp[e_name], decimals=2))
 519            else:
 520                label = e_name + ', S=' + str(np.around(self.S[e_name], decimals=2))
 521                xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
 522
 523            plt.errorbar(np.arange(length)[:int(xmax) + 1], self.e_n_tauint[e_name][:int(xmax) + 1], yerr=self.e_n_dtauint[e_name][:int(xmax) + 1], linewidth=1, capsize=2, label=label)
 524            plt.axvline(x=self.e_windowsize[e_name], color='C' + str(e), alpha=0.5, marker=',', ls='--')
 525            plt.legend()
 526            plt.xlim(-0.5, xmax)
 527            ylim = plt.ylim()
 528            plt.ylim(bottom=0.0, top=max(1.0, ylim[1]))
 529            plt.draw()
 530            if save:
 531                fig.savefig(save + "_" + str(e))
 532
 533    def plot_rho(self, save=None):
 534        """Plot normalized autocorrelation function time for each ensemble.
 535
 536        Parameters
 537        ----------
 538        save : str
 539            saves the figure to a file named 'save' if.
 540        """
 541        if not hasattr(self, 'e_dvalue'):
 542            raise Exception('Run the gamma method first.')
 543        for e, e_name in enumerate(self.mc_names):
 544            fig = plt.figure()
 545            plt.xlabel('W')
 546            plt.ylabel('rho')
 547            length = int(len(self.e_drho[e_name]))
 548            plt.errorbar(np.arange(length), self.e_rho[e_name][:length], yerr=self.e_drho[e_name][:], linewidth=1, capsize=2)
 549            plt.axvline(x=self.e_windowsize[e_name], color='r', alpha=0.25, ls='--', marker=',')
 550            if self.tau_exp[e_name] > 0:
 551                plt.plot([self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name]],
 552                         [self.e_rho[e_name][self.e_windowsize[e_name] + 1], 0], 'k-', lw=1)
 553                xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
 554                plt.title('Rho ' + e_name + r', tau\_exp=' + str(np.around(self.tau_exp[e_name], decimals=2)))
 555            else:
 556                xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
 557                plt.title('Rho ' + e_name + ', S=' + str(np.around(self.S[e_name], decimals=2)))
 558            plt.plot([-0.5, xmax], [0, 0], 'k--', lw=1)
 559            plt.xlim(-0.5, xmax)
 560            plt.draw()
 561            if save:
 562                fig.savefig(save + "_" + str(e))
 563
 564    def plot_rep_dist(self):
 565        """Plot replica distribution for each ensemble with more than one replicum."""
 566        if not hasattr(self, 'e_dvalue'):
 567            raise Exception('Run the gamma method first.')
 568        for e, e_name in enumerate(self.mc_names):
 569            if len(self.e_content[e_name]) == 1:
 570                print('No replica distribution for a single replicum (', e_name, ')')
 571                continue
 572            r_length = []
 573            sub_r_mean = 0
 574            for r, r_name in enumerate(self.e_content[e_name]):
 575                r_length.append(len(self.deltas[r_name]))
 576                sub_r_mean += self.shape[r_name] * self.r_values[r_name]
 577            e_N = np.sum(r_length)
 578            sub_r_mean /= e_N
 579            arr = np.zeros(len(self.e_content[e_name]))
 580            for r, r_name in enumerate(self.e_content[e_name]):
 581                arr[r] = (self.r_values[r_name] - sub_r_mean) / (self.e_dvalue[e_name] * np.sqrt(e_N / self.shape[r_name] - 1))
 582            plt.hist(arr, rwidth=0.8, bins=len(self.e_content[e_name]))
 583            plt.title('Replica distribution' + e_name + ' (mean=0, var=1)')
 584            plt.draw()
 585
 586    def plot_history(self, expand=True):
 587        """Plot derived Monte Carlo history for each ensemble
 588
 589        Parameters
 590        ----------
 591        expand : bool
 592            show expanded history for irregular Monte Carlo chains (default: True).
 593        """
 594        for e, e_name in enumerate(self.mc_names):
 595            plt.figure()
 596            r_length = []
 597            tmp = []
 598            tmp_expanded = []
 599            for r, r_name in enumerate(self.e_content[e_name]):
 600                tmp.append(self.deltas[r_name] + self.r_values[r_name])
 601                if expand:
 602                    tmp_expanded.append(_expand_deltas(self.deltas[r_name], list(self.idl[r_name]), self.shape[r_name]) + self.r_values[r_name])
 603                    r_length.append(len(tmp_expanded[-1]))
 604                else:
 605                    r_length.append(len(tmp[-1]))
 606            e_N = np.sum(r_length)
 607            x = np.arange(e_N)
 608            y_test = np.concatenate(tmp, axis=0)
 609            if expand:
 610                y = np.concatenate(tmp_expanded, axis=0)
 611            else:
 612                y = y_test
 613            plt.errorbar(x, y, fmt='.', markersize=3)
 614            plt.xlim(-0.5, e_N - 0.5)
 615            plt.title(e_name + f'\nskew: {skew(y_test):.3f} (p={skewtest(y_test).pvalue:.3f}), kurtosis: {kurtosis(y_test):.3f} (p={kurtosistest(y_test).pvalue:.3f})')
 616            plt.draw()
 617
 618    def plot_piechart(self, save=None):
 619        """Plot piechart which shows the fractional contribution of each
 620        ensemble to the error and returns a dictionary containing the fractions.
 621
 622        Parameters
 623        ----------
 624        save : str
 625            saves the figure to a file named 'save' if.
 626        """
 627        if not hasattr(self, 'e_dvalue'):
 628            raise Exception('Run the gamma method first.')
 629        if np.isclose(0.0, self._dvalue, atol=1e-15):
 630            raise Exception('Error is 0.0')
 631        labels = self.e_names
 632        sizes = [self.e_dvalue[name] ** 2 for name in labels] / self._dvalue ** 2
 633        fig1, ax1 = plt.subplots()
 634        ax1.pie(sizes, labels=labels, startangle=90, normalize=True)
 635        ax1.axis('equal')
 636        plt.draw()
 637        if save:
 638            fig1.savefig(save)
 639
 640        return dict(zip(self.e_names, sizes))
 641
 642    def dump(self, filename, datatype="json.gz", description="", **kwargs):
 643        """Dump the Obs to a file 'name' of chosen format.
 644
 645        Parameters
 646        ----------
 647        filename : str
 648            name of the file to be saved.
 649        datatype : str
 650            Format of the exported file. Supported formats include
 651            "json.gz" and "pickle"
 652        description : str
 653            Description for output file, only relevant for json.gz format.
 654        path : str
 655            specifies a custom path for the file (default '.')
 656        """
 657        if 'path' in kwargs:
 658            file_name = kwargs.get('path') + '/' + filename
 659        else:
 660            file_name = filename
 661
 662        if datatype == "json.gz":
 663            from .input.json import dump_to_json
 664            dump_to_json([self], file_name, description=description)
 665        elif datatype == "pickle":
 666            with open(file_name + '.p', 'wb') as fb:
 667                pickle.dump(self, fb)
 668        else:
 669            raise Exception("Unknown datatype " + str(datatype))
 670
 671    def export_jackknife(self):
 672        """Export jackknife samples from the Obs
 673
 674        Returns
 675        -------
 676        numpy.ndarray
 677            Returns a numpy array of length N + 1 where N is the number of samples
 678            for the given ensemble and replicum. The zeroth entry of the array contains
 679            the mean value of the Obs, entries 1 to N contain the N jackknife samples
 680            derived from the Obs. The current implementation only works for observables
 681            defined on exactly one ensemble and replicum. The derived jackknife samples
 682            should agree with samples from a full jackknife analysis up to O(1/N).
 683        """
 684
 685        if len(self.names) != 1:
 686            raise Exception("'export_jackknife' is only implemented for Obs defined on one ensemble and replicum.")
 687
 688        name = self.names[0]
 689        full_data = self.deltas[name] + self.r_values[name]
 690        n = full_data.size
 691        mean = self.value
 692        tmp_jacks = np.zeros(n + 1)
 693        tmp_jacks[0] = mean
 694        tmp_jacks[1:] = (n * mean - full_data) / (n - 1)
 695        return tmp_jacks
 696
 697    def __float__(self):
 698        return float(self.value)
 699
 700    def __repr__(self):
 701        return 'Obs[' + str(self) + ']'
 702
 703    def __str__(self):
 704        return _format_uncertainty(self.value, self._dvalue)
 705
 706    def __hash__(self):
 707        hash_tuple = (np.array([self.value]).astype(np.float32).data.tobytes(),)
 708        hash_tuple += tuple([o.astype(np.float32).data.tobytes() for o in self.deltas.values()])
 709        hash_tuple += tuple([np.array([o.errsq()]).astype(np.float32).data.tobytes() for o in self.covobs.values()])
 710        hash_tuple += tuple([o.encode() for o in self.names])
 711        m = hashlib.md5()
 712        [m.update(o) for o in hash_tuple]
 713        return int(m.hexdigest(), 16) & 0xFFFFFFFF
 714
 715    # Overload comparisons
 716    def __lt__(self, other):
 717        return self.value < other
 718
 719    def __le__(self, other):
 720        return self.value <= other
 721
 722    def __gt__(self, other):
 723        return self.value > other
 724
 725    def __ge__(self, other):
 726        return self.value >= other
 727
 728    def __eq__(self, other):
 729        return (self - other).is_zero()
 730
 731    def __ne__(self, other):
 732        return not (self - other).is_zero()
 733
 734    # Overload math operations
 735    def __add__(self, y):
 736        if isinstance(y, Obs):
 737            return derived_observable(lambda x, **kwargs: x[0] + x[1], [self, y], man_grad=[1, 1])
 738        else:
 739            if isinstance(y, np.ndarray):
 740                return np.array([self + o for o in y])
 741            elif y.__class__.__name__ in ['Corr', 'CObs']:
 742                return NotImplemented
 743            else:
 744                return derived_observable(lambda x, **kwargs: x[0] + y, [self], man_grad=[1])
 745
 746    def __radd__(self, y):
 747        return self + y
 748
 749    def __mul__(self, y):
 750        if isinstance(y, Obs):
 751            return derived_observable(lambda x, **kwargs: x[0] * x[1], [self, y], man_grad=[y.value, self.value])
 752        else:
 753            if isinstance(y, np.ndarray):
 754                return np.array([self * o for o in y])
 755            elif isinstance(y, complex):
 756                return CObs(self * y.real, self * y.imag)
 757            elif y.__class__.__name__ in ['Corr', 'CObs']:
 758                return NotImplemented
 759            else:
 760                return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y])
 761
 762    def __rmul__(self, y):
 763        return self * y
 764
 765    def __sub__(self, y):
 766        if isinstance(y, Obs):
 767            return derived_observable(lambda x, **kwargs: x[0] - x[1], [self, y], man_grad=[1, -1])
 768        else:
 769            if isinstance(y, np.ndarray):
 770                return np.array([self - o for o in y])
 771            elif y.__class__.__name__ in ['Corr', 'CObs']:
 772                return NotImplemented
 773            else:
 774                return derived_observable(lambda x, **kwargs: x[0] - y, [self], man_grad=[1])
 775
 776    def __rsub__(self, y):
 777        return -1 * (self - y)
 778
 779    def __pos__(self):
 780        return self
 781
 782    def __neg__(self):
 783        return -1 * self
 784
 785    def __truediv__(self, y):
 786        if isinstance(y, Obs):
 787            return derived_observable(lambda x, **kwargs: x[0] / x[1], [self, y], man_grad=[1 / y.value, - self.value / y.value ** 2])
 788        else:
 789            if isinstance(y, np.ndarray):
 790                return np.array([self / o for o in y])
 791            elif y.__class__.__name__ in ['Corr', 'CObs']:
 792                return NotImplemented
 793            else:
 794                return derived_observable(lambda x, **kwargs: x[0] / y, [self], man_grad=[1 / y])
 795
 796    def __rtruediv__(self, y):
 797        if isinstance(y, Obs):
 798            return derived_observable(lambda x, **kwargs: x[0] / x[1], [y, self], man_grad=[1 / self.value, - y.value / self.value ** 2])
 799        else:
 800            if isinstance(y, np.ndarray):
 801                return np.array([o / self for o in y])
 802            elif y.__class__.__name__ in ['Corr', 'CObs']:
 803                return NotImplemented
 804            else:
 805                return derived_observable(lambda x, **kwargs: y / x[0], [self], man_grad=[-y / self.value ** 2])
 806
 807    def __pow__(self, y):
 808        if isinstance(y, Obs):
 809            return derived_observable(lambda x: x[0] ** x[1], [self, y])
 810        else:
 811            return derived_observable(lambda x: x[0] ** y, [self])
 812
 813    def __rpow__(self, y):
 814        if isinstance(y, Obs):
 815            return derived_observable(lambda x: x[0] ** x[1], [y, self])
 816        else:
 817            return derived_observable(lambda x: y ** x[0], [self])
 818
 819    def __abs__(self):
 820        return derived_observable(lambda x: anp.abs(x[0]), [self])
 821
 822    # Overload numpy functions
 823    def sqrt(self):
 824        return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
 825
 826    def log(self):
 827        return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
 828
 829    def exp(self):
 830        return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
 831
 832    def sin(self):
 833        return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
 834
 835    def cos(self):
 836        return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
 837
 838    def tan(self):
 839        return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
 840
 841    def arcsin(self):
 842        return derived_observable(lambda x: anp.arcsin(x[0]), [self])
 843
 844    def arccos(self):
 845        return derived_observable(lambda x: anp.arccos(x[0]), [self])
 846
 847    def arctan(self):
 848        return derived_observable(lambda x: anp.arctan(x[0]), [self])
 849
 850    def sinh(self):
 851        return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
 852
 853    def cosh(self):
 854        return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
 855
 856    def tanh(self):
 857        return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
 858
 859    def arcsinh(self):
 860        return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
 861
 862    def arccosh(self):
 863        return derived_observable(lambda x: anp.arccosh(x[0]), [self])
 864
 865    def arctanh(self):
 866        return derived_observable(lambda x: anp.arctanh(x[0]), [self])
 867
 868
 869class CObs:
 870    """Class for a complex valued observable."""
 871    __slots__ = ['_real', '_imag', 'tag']
 872
 873    def __init__(self, real, imag=0.0):
 874        self._real = real
 875        self._imag = imag
 876        self.tag = None
 877
 878    @property
 879    def real(self):
 880        return self._real
 881
 882    @property
 883    def imag(self):
 884        return self._imag
 885
 886    def gamma_method(self, **kwargs):
 887        """Executes the gamma_method for the real and the imaginary part."""
 888        if isinstance(self.real, Obs):
 889            self.real.gamma_method(**kwargs)
 890        if isinstance(self.imag, Obs):
 891            self.imag.gamma_method(**kwargs)
 892
 893    def is_zero(self):
 894        """Checks whether both real and imaginary part are zero within machine precision."""
 895        return self.real == 0.0 and self.imag == 0.0
 896
 897    def conjugate(self):
 898        return CObs(self.real, -self.imag)
 899
 900    def __add__(self, other):
 901        if isinstance(other, np.ndarray):
 902            return other + self
 903        elif hasattr(other, 'real') and hasattr(other, 'imag'):
 904            return CObs(self.real + other.real,
 905                        self.imag + other.imag)
 906        else:
 907            return CObs(self.real + other, self.imag)
 908
 909    def __radd__(self, y):
 910        return self + y
 911
 912    def __sub__(self, other):
 913        if isinstance(other, np.ndarray):
 914            return -1 * (other - self)
 915        elif hasattr(other, 'real') and hasattr(other, 'imag'):
 916            return CObs(self.real - other.real, self.imag - other.imag)
 917        else:
 918            return CObs(self.real - other, self.imag)
 919
 920    def __rsub__(self, other):
 921        return -1 * (self - other)
 922
 923    def __mul__(self, other):
 924        if isinstance(other, np.ndarray):
 925            return other * self
 926        elif hasattr(other, 'real') and hasattr(other, 'imag'):
 927            if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
 928                return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
 929                                               [self.real, other.real, self.imag, other.imag],
 930                                               man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
 931                            derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
 932                                               [self.real, other.real, self.imag, other.imag],
 933                                               man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
 934            elif getattr(other, 'imag', 0) != 0:
 935                return CObs(self.real * other.real - self.imag * other.imag,
 936                            self.imag * other.real + self.real * other.imag)
 937            else:
 938                return CObs(self.real * other.real, self.imag * other.real)
 939        else:
 940            return CObs(self.real * other, self.imag * other)
 941
 942    def __rmul__(self, other):
 943        return self * other
 944
 945    def __truediv__(self, other):
 946        if isinstance(other, np.ndarray):
 947            return 1 / (other / self)
 948        elif hasattr(other, 'real') and hasattr(other, 'imag'):
 949            r = other.real ** 2 + other.imag ** 2
 950            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
 951        else:
 952            return CObs(self.real / other, self.imag / other)
 953
 954    def __rtruediv__(self, other):
 955        r = self.real ** 2 + self.imag ** 2
 956        if hasattr(other, 'real') and hasattr(other, 'imag'):
 957            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
 958        else:
 959            return CObs(self.real * other / r, -self.imag * other / r)
 960
 961    def __abs__(self):
 962        return np.sqrt(self.real**2 + self.imag**2)
 963
 964    def __pos__(self):
 965        return self
 966
 967    def __neg__(self):
 968        return -1 * self
 969
 970    def __eq__(self, other):
 971        return self.real == other.real and self.imag == other.imag
 972
 973    def __str__(self):
 974        return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
 975
 976    def __repr__(self):
 977        return 'CObs[' + str(self) + ']'
 978
 979
 980def _format_uncertainty(value, dvalue):
 981    """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)"""
 982    if dvalue == 0.0:
 983        return str(value)
 984    fexp = np.floor(np.log10(dvalue))
 985    if fexp < 0.0:
 986        return '{:{form}}({:2.0f})'.format(value, dvalue * 10 ** (-fexp + 1), form='.' + str(-int(fexp) + 1) + 'f')
 987    elif fexp == 0.0:
 988        return '{:.1f}({:1.1f})'.format(value, dvalue)
 989    else:
 990        return '{:.0f}({:2.0f})'.format(value, dvalue)
 991
 992
 993def _expand_deltas(deltas, idx, shape):
 994    """Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0.
 995       If idx is of type range, the deltas are not changed
 996
 997    Parameters
 998    ----------
 999    deltas : list
1000        List of fluctuations
1001    idx : list
1002        List or range of configs on which the deltas are defined, has to be sorted in ascending order.
1003    shape : int
1004        Number of configs in idx.
1005    """
1006    if isinstance(idx, range):
1007        return deltas
1008    else:
1009        ret = np.zeros(idx[-1] - idx[0] + 1)
1010        for i in range(shape):
1011            ret[idx[i] - idx[0]] = deltas[i]
1012        return ret
1013
1014
1015def _merge_idx(idl):
1016    """Returns the union of all lists in idl as sorted list
1017
1018    Parameters
1019    ----------
1020    idl : list
1021        List of lists or ranges.
1022    """
1023
1024    # Use groupby to efficiently check whether all elements of idl are identical
1025    try:
1026        g = groupby(idl)
1027        if next(g, True) and not next(g, False):
1028            return idl[0]
1029    except Exception:
1030        pass
1031
1032    if np.all([type(idx) is range for idx in idl]):
1033        if len(set([idx[0] for idx in idl])) == 1:
1034            idstart = min([idx.start for idx in idl])
1035            idstop = max([idx.stop for idx in idl])
1036            idstep = min([idx.step for idx in idl])
1037            return range(idstart, idstop, idstep)
1038
1039    return sorted(set().union(*idl))
1040
1041
1042def _intersection_idx(idl):
1043    """Returns the intersection of all lists in idl as sorted list
1044
1045    Parameters
1046    ----------
1047    idl : list
1048        List of lists or ranges.
1049    """
1050
1051    def _lcm(*args):
1052        """Returns the lowest common multiple of args.
1053
1054        From python 3.9 onwards the math library contains an lcm function."""
1055        return reduce(lambda a, b: a * b // gcd(a, b), args)
1056
1057    # Use groupby to efficiently check whether all elements of idl are identical
1058    try:
1059        g = groupby(idl)
1060        if next(g, True) and not next(g, False):
1061            return idl[0]
1062    except Exception:
1063        pass
1064
1065    if np.all([type(idx) is range for idx in idl]):
1066        if len(set([idx[0] for idx in idl])) == 1:
1067            idstart = max([idx.start for idx in idl])
1068            idstop = min([idx.stop for idx in idl])
1069            idstep = _lcm(*[idx.step for idx in idl])
1070            return range(idstart, idstop, idstep)
1071
1072    return sorted(set.intersection(*[set(o) for o in idl]))
1073
1074
1075def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
1076    """Expand deltas defined on idx to the list of configs that is defined by new_idx.
1077       New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
1078       common divisor of the step sizes is used as new step size.
1079
1080    Parameters
1081    ----------
1082    deltas : list
1083        List of fluctuations
1084    idx : list
1085        List or range of configs on which the deltas are defined.
1086        Has to be a subset of new_idx and has to be sorted in ascending order.
1087    shape : list
1088        Number of configs in idx.
1089    new_idx : list
1090        List of configs that defines the new range, has to be sorted in ascending order.
1091    """
1092
1093    if type(idx) is range and type(new_idx) is range:
1094        if idx == new_idx:
1095            return deltas
1096    ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
1097    for i in range(shape):
1098        ret[idx[i] - new_idx[0]] = deltas[i]
1099    return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
1100
1101
1102def derived_observable(func, data, array_mode=False, **kwargs):
1103    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
1104
1105    Parameters
1106    ----------
1107    func : object
1108        arbitrary function of the form func(data, **kwargs). For the
1109        automatic differentiation to work, all numpy functions have to have
1110        the autograd wrapper (use 'import autograd.numpy as anp').
1111    data : list
1112        list of Obs, e.g. [obs1, obs2, obs3].
1113    num_grad : bool
1114        if True, numerical derivatives are used instead of autograd
1115        (default False). To control the numerical differentiation the
1116        kwargs of numdifftools.step_generators.MaxStepGenerator
1117        can be used.
1118    man_grad : list
1119        manually supply a list or an array which contains the jacobian
1120        of func. Use cautiously, supplying the wrong derivative will
1121        not be intercepted.
1122
1123    Notes
1124    -----
1125    For simple mathematical operations it can be practical to use anonymous
1126    functions. For the ratio of two observables one can e.g. use
1127
1128    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
1129    """
1130
1131    data = np.asarray(data)
1132    raveled_data = data.ravel()
1133
1134    # Workaround for matrix operations containing non Obs data
1135    if not all(isinstance(x, Obs) for x in raveled_data):
1136        for i in range(len(raveled_data)):
1137            if isinstance(raveled_data[i], (int, float)):
1138                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
1139
1140    allcov = {}
1141    for o in raveled_data:
1142        for name in o.cov_names:
1143            if name in allcov:
1144                if not np.allclose(allcov[name], o.covobs[name].cov):
1145                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
1146            else:
1147                allcov[name] = o.covobs[name].cov
1148
1149    n_obs = len(raveled_data)
1150    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
1151    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
1152    new_sample_names = sorted(set(new_names) - set(new_cov_names))
1153
1154    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
1155
1156    if data.ndim == 1:
1157        values = np.array([o.value for o in data])
1158    else:
1159        values = np.vectorize(lambda x: x.value)(data)
1160
1161    new_values = func(values, **kwargs)
1162
1163    multi = int(isinstance(new_values, np.ndarray))
1164
1165    new_r_values = {}
1166    new_idl_d = {}
1167    for name in new_sample_names:
1168        idl = []
1169        tmp_values = np.zeros(n_obs)
1170        for i, item in enumerate(raveled_data):
1171            tmp_values[i] = item.r_values.get(name, item.value)
1172            tmp_idl = item.idl.get(name)
1173            if tmp_idl is not None:
1174                idl.append(tmp_idl)
1175        if multi > 0:
1176            tmp_values = np.array(tmp_values).reshape(data.shape)
1177        new_r_values[name] = func(tmp_values, **kwargs)
1178        new_idl_d[name] = _merge_idx(idl)
1179
1180    if 'man_grad' in kwargs:
1181        deriv = np.asarray(kwargs.get('man_grad'))
1182        if new_values.shape + data.shape != deriv.shape:
1183            raise Exception('Manual derivative does not have correct shape.')
1184    elif kwargs.get('num_grad') is True:
1185        if multi > 0:
1186            raise Exception('Multi mode currently not supported for numerical derivative')
1187        options = {
1188            'base_step': 0.1,
1189            'step_ratio': 2.5}
1190        for key in options.keys():
1191            kwarg = kwargs.get(key)
1192            if kwarg is not None:
1193                options[key] = kwarg
1194        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
1195        if tmp_df.size == 1:
1196            deriv = np.array([tmp_df.real])
1197        else:
1198            deriv = tmp_df.real
1199    else:
1200        deriv = jacobian(func)(values, **kwargs)
1201
1202    final_result = np.zeros(new_values.shape, dtype=object)
1203
1204    if array_mode is True:
1205
1206        class _Zero_grad():
1207            def __init__(self, N):
1208                self.grad = np.zeros((N, 1))
1209
1210        new_covobs_lengths = dict(set([y for x in [[(n, o.covobs[n].N) for n in o.cov_names] for o in raveled_data] for y in x]))
1211        d_extracted = {}
1212        g_extracted = {}
1213        for name in new_sample_names:
1214            d_extracted[name] = []
1215            ens_length = len(new_idl_d[name])
1216            for i_dat, dat in enumerate(data):
1217                d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
1218        for name in new_cov_names:
1219            g_extracted[name] = []
1220            zero_grad = _Zero_grad(new_covobs_lengths[name])
1221            for i_dat, dat in enumerate(data):
1222                g_extracted[name].append(np.array([o.covobs.get(name, zero_grad).grad for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (new_covobs_lengths[name], 1)))
1223
1224    for i_val, new_val in np.ndenumerate(new_values):
1225        new_deltas = {}
1226        new_grad = {}
1227        if array_mode is True:
1228            for name in new_sample_names:
1229                ens_length = d_extracted[name][0].shape[-1]
1230                new_deltas[name] = np.zeros(ens_length)
1231                for i_dat, dat in enumerate(d_extracted[name]):
1232                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
1233            for name in new_cov_names:
1234                new_grad[name] = 0
1235                for i_dat, dat in enumerate(g_extracted[name]):
1236                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
1237        else:
1238            for j_obs, obs in np.ndenumerate(data):
1239                for name in obs.names:
1240                    if name in obs.cov_names:
1241                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
1242                    else:
1243                        new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name])
1244
1245        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
1246
1247        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
1248            raise Exception('The same name has been used for deltas and covobs!')
1249        new_samples = []
1250        new_means = []
1251        new_idl = []
1252        new_names_obs = []
1253        for name in new_names:
1254            if name not in new_covobs:
1255                new_samples.append(new_deltas[name])
1256                new_idl.append(new_idl_d[name])
1257                new_means.append(new_r_values[name][i_val])
1258                new_names_obs.append(name)
1259        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
1260        for name in new_covobs:
1261            final_result[i_val].names.append(name)
1262        final_result[i_val]._covobs = new_covobs
1263        final_result[i_val]._value = new_val
1264        final_result[i_val].reweighted = reweighted
1265
1266    if multi == 0:
1267        final_result = final_result.item()
1268
1269    return final_result
1270
1271
1272def _reduce_deltas(deltas, idx_old, idx_new):
1273    """Extract deltas defined on idx_old on all configs of idx_new.
1274
1275    Assumes, that idx_old and idx_new are correctly defined idl, i.e., they
1276    are ordered in an ascending order.
1277
1278    Parameters
1279    ----------
1280    deltas : list
1281        List of fluctuations
1282    idx_old : list
1283        List or range of configs on which the deltas are defined
1284    idx_new : list
1285        List of configs for which we want to extract the deltas.
1286        Has to be a subset of idx_old.
1287    """
1288    if not len(deltas) == len(idx_old):
1289        raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
1290    if type(idx_old) is range and type(idx_new) is range:
1291        if idx_old == idx_new:
1292            return deltas
1293    # Use groupby to efficiently check whether all elements of idx_old and idx_new are identical
1294    try:
1295        g = groupby([idx_old, idx_new])
1296        if next(g, True) and not next(g, False):
1297            return deltas
1298    except Exception:
1299        pass
1300    indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1]
1301    if len(indices) < len(idx_new):
1302        raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old')
1303    return np.array(deltas)[indices]
1304
1305
1306def reweight(weight, obs, **kwargs):
1307    """Reweight a list of observables.
1308
1309    Parameters
1310    ----------
1311    weight : Obs
1312        Reweighting factor. An Observable that has to be defined on a superset of the
1313        configurations in obs[i].idl for all i.
1314    obs : list
1315        list of Obs, e.g. [obs1, obs2, obs3].
1316    all_configs : bool
1317        if True, the reweighted observables are normalized by the average of
1318        the reweighting factor on all configurations in weight.idl and not
1319        on the configurations in obs[i].idl. Default False.
1320    """
1321    result = []
1322    for i in range(len(obs)):
1323        if len(obs[i].cov_names):
1324            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
1325        if not set(obs[i].names).issubset(weight.names):
1326            raise Exception('Error: Ensembles do not fit')
1327        for name in obs[i].names:
1328            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
1329                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
1330        new_samples = []
1331        w_deltas = {}
1332        for name in sorted(obs[i].names):
1333            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
1334            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
1335        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
1336
1337        if kwargs.get('all_configs'):
1338            new_weight = weight
1339        else:
1340            new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
1341
1342        result.append(tmp_obs / new_weight)
1343        result[-1].reweighted = True
1344
1345    return result
1346
1347
1348def correlate(obs_a, obs_b):
1349    """Correlate two observables.
1350
1351    Parameters
1352    ----------
1353    obs_a : Obs
1354        First observable
1355    obs_b : Obs
1356        Second observable
1357
1358    Notes
1359    -----
1360    Keep in mind to only correlate primary observables which have not been reweighted
1361    yet. The reweighting has to be applied after correlating the observables.
1362    Currently only works if ensembles are identical (this is not strictly necessary).
1363    """
1364
1365    if sorted(obs_a.names) != sorted(obs_b.names):
1366        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
1367    if len(obs_a.cov_names) or len(obs_b.cov_names):
1368        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
1369    for name in obs_a.names:
1370        if obs_a.shape[name] != obs_b.shape[name]:
1371            raise Exception('Shapes of ensemble', name, 'do not fit')
1372        if obs_a.idl[name] != obs_b.idl[name]:
1373            raise Exception('idl of ensemble', name, 'do not fit')
1374
1375    if obs_a.reweighted is True:
1376        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
1377    if obs_b.reweighted is True:
1378        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
1379
1380    new_samples = []
1381    new_idl = []
1382    for name in sorted(obs_a.names):
1383        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
1384        new_idl.append(obs_a.idl[name])
1385
1386    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
1387    o.reweighted = obs_a.reweighted or obs_b.reweighted
1388    return o
1389
1390
1391def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
1392    r'''Calculates the error covariance matrix of a set of observables.
1393
1394    WARNING: This function should be used with care, especially for observables with support on multiple
1395             ensembles with differing autocorrelations. See the notes below for details.
1396
1397    The gamma method has to be applied first to all observables.
1398
1399    Parameters
1400    ----------
1401    obs : list or numpy.ndarray
1402        List or one dimensional array of Obs
1403    visualize : bool
1404        If True plots the corresponding normalized correlation matrix (default False).
1405    correlation : bool
1406        If True the correlation matrix instead of the error covariance matrix is returned (default False).
1407    smooth : None or int
1408        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
1409        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
1410        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
1411        small ones.
1412
1413    Notes
1414    -----
1415    The error covariance is defined such that it agrees with the squared standard error for two identical observables
1416    $$\operatorname{cov}(a,a)=\sum_{s=1}^N\delta_a^s\delta_a^s/N^2=\Gamma_{aa}(0)/N=\operatorname{var}(a)/N=\sigma_a^2$$
1417    in the absence of autocorrelation.
1418    The error covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
1419    $$\sum_{i,j}v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v\in\mathbb{R}^M$, while such an identity does not hold for larger windows/lags.
1420    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
1421    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
1422    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
1423    '''
1424
1425    length = len(obs)
1426
1427    max_samples = np.max([o.N for o in obs])
1428    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
1429        warnings.warn(f"The dimension of the covariance matrix ({length}) is larger or equal to the number of samples ({max_samples}). This will result in a rank deficient matrix.", RuntimeWarning)
1430
1431    cov = np.zeros((length, length))
1432    for i in range(length):
1433        for j in range(i, length):
1434            cov[i, j] = _covariance_element(obs[i], obs[j])
1435    cov = cov + cov.T - np.diag(np.diag(cov))
1436
1437    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
1438
1439    if isinstance(smooth, int):
1440        corr = _smooth_eigenvalues(corr, smooth)
1441
1442    if visualize:
1443        plt.matshow(corr, vmin=-1, vmax=1)
1444        plt.set_cmap('RdBu')
1445        plt.colorbar()
1446        plt.draw()
1447
1448    if correlation is True:
1449        return corr
1450
1451    errors = [o.dvalue for o in obs]
1452    cov = np.diag(errors) @ corr @ np.diag(errors)
1453
1454    eigenvalues = np.linalg.eigh(cov)[0]
1455    if not np.all(eigenvalues >= 0):
1456        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
1457
1458    return cov
1459
1460
1461def _smooth_eigenvalues(corr, E):
1462    """Eigenvalue smoothing as described in hep-lat/9412087
1463
1464    corr : np.ndarray
1465        correlation matrix
1466    E : integer
1467        Number of eigenvalues to be left substantially unchanged
1468    """
1469    if not (2 < E < corr.shape[0] - 1):
1470        raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).")
1471    vals, vec = np.linalg.eigh(corr)
1472    lambda_min = np.mean(vals[:-E])
1473    vals[vals < lambda_min] = lambda_min
1474    vals /= np.mean(vals)
1475    return vec @ np.diag(vals) @ vec.T
1476
1477
1478def _covariance_element(obs1, obs2):
1479    """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
1480
1481    def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
1482        deltas1 = _reduce_deltas(deltas1, idx1, new_idx)
1483        deltas2 = _reduce_deltas(deltas2, idx2, new_idx)
1484        return np.sum(deltas1 * deltas2)
1485
1486    if set(obs1.names).isdisjoint(set(obs2.names)):
1487        return 0.0
1488
1489    if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
1490        raise Exception('The gamma method has to be applied to both Obs first.')
1491
1492    dvalue = 0.0
1493
1494    for e_name in obs1.mc_names:
1495
1496        if e_name not in obs2.mc_names:
1497            continue
1498
1499        idl_d = {}
1500        for r_name in obs1.e_content[e_name]:
1501            if r_name not in obs2.e_content[e_name]:
1502                continue
1503            idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
1504
1505        gamma = 0.0
1506
1507        for r_name in obs1.e_content[e_name]:
1508            if r_name not in obs2.e_content[e_name]:
1509                continue
1510            if len(idl_d[r_name]) == 0:
1511                continue
1512            gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
1513
1514        if gamma == 0.0:
1515            continue
1516
1517        gamma_div = 0.0
1518        for r_name in obs1.e_content[e_name]:
1519            if r_name not in obs2.e_content[e_name]:
1520                continue
1521            if len(idl_d[r_name]) == 0:
1522                continue
1523            gamma_div += np.sqrt(calc_gamma(obs1.deltas[r_name], obs1.deltas[r_name], obs1.idl[r_name], obs1.idl[r_name], idl_d[r_name]) * calc_gamma(obs2.deltas[r_name], obs2.deltas[r_name], obs2.idl[r_name], obs2.idl[r_name], idl_d[r_name]))
1524        gamma /= gamma_div
1525
1526        dvalue += gamma
1527
1528    for e_name in obs1.cov_names:
1529
1530        if e_name not in obs2.cov_names:
1531            continue
1532
1533        dvalue += float(np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)))
1534
1535    return dvalue
1536
1537
1538def import_jackknife(jacks, name, idl=None):
1539    """Imports jackknife samples and returns an Obs
1540
1541    Parameters
1542    ----------
1543    jacks : numpy.ndarray
1544        numpy array containing the mean value as zeroth entry and
1545        the N jackknife samples as first to Nth entry.
1546    name : str
1547        name of the ensemble the samples are defined on.
1548    """
1549    length = len(jacks) - 1
1550    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
1551    samples = jacks[1:] @ prj
1552    mean = np.mean(samples)
1553    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
1554    new_obs._value = jacks[0]
1555    return new_obs
1556
1557
1558def merge_obs(list_of_obs):
1559    """Combine all observables in list_of_obs into one new observable
1560
1561    Parameters
1562    ----------
1563    list_of_obs : list
1564        list of the Obs object to be combined
1565
1566    Notes
1567    -----
1568    It is not possible to combine obs which are based on the same replicum
1569    """
1570    replist = [item for obs in list_of_obs for item in obs.names]
1571    if (len(replist) == len(set(replist))) is False:
1572        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
1573    if any([len(o.cov_names) for o in list_of_obs]):
1574        raise Exception('Not possible to merge data that contains covobs!')
1575    new_dict = {}
1576    idl_dict = {}
1577    for o in list_of_obs:
1578        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
1579                        for key in set(o.deltas) | set(o.r_values)})
1580        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
1581
1582    names = sorted(new_dict.keys())
1583    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
1584    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
1585    return o
1586
1587
1588def cov_Obs(means, cov, name, grad=None):
1589    """Create an Obs based on mean(s) and a covariance matrix
1590
1591    Parameters
1592    ----------
1593    mean : list of floats or float
1594        N mean value(s) of the new Obs
1595    cov : list or array
1596        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
1597    name : str
1598        identifier for the covariance matrix
1599    grad : list or array
1600        Gradient of the Covobs wrt. the means belonging to cov.
1601    """
1602
1603    def covobs_to_obs(co):
1604        """Make an Obs out of a Covobs
1605
1606        Parameters
1607        ----------
1608        co : Covobs
1609            Covobs to be embedded into the Obs
1610        """
1611        o = Obs([], [], means=[])
1612        o._value = co.value
1613        o.names.append(co.name)
1614        o._covobs[co.name] = co
1615        o._dvalue = np.sqrt(co.errsq())
1616        return o
1617
1618    ol = []
1619    if isinstance(means, (float, int)):
1620        means = [means]
1621
1622    for i in range(len(means)):
1623        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
1624    if ol[0].covobs[name].N != len(means):
1625        raise Exception('You have to provide %d mean values!' % (ol[0].N))
1626    if len(ol) == 1:
1627        return ol[0]
1628    return ol