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