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