pyerrors.obs

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