pyerrors.misc

  1import pickle
  2import numpy as np
  3from .obs import Obs
  4
  5
  6def dump_object(obj, name, **kwargs):
  7    """Dump object into pickle file.
  8
  9    Parameters
 10    ----------
 11    obj : object
 12        object to be saved in the pickle file
 13    name : str
 14        name of the file
 15    path : str
 16        specifies a custom path for the file (default '.')
 17
 18    Returns
 19    -------
 20    None
 21    """
 22    if 'path' in kwargs:
 23        file_name = kwargs.get('path') + '/' + name + '.p'
 24    else:
 25        file_name = name + '.p'
 26    with open(file_name, 'wb') as fb:
 27        pickle.dump(obj, fb)
 28
 29
 30def load_object(path):
 31    """Load object from pickle file.
 32
 33    Parameters
 34    ----------
 35    path : str
 36        path to the file
 37
 38    Returns
 39    -------
 40    object : Obs
 41        Loaded Object
 42    """
 43    with open(path, 'rb') as file:
 44        return pickle.load(file)
 45
 46
 47def pseudo_Obs(value, dvalue, name, samples=1000):
 48    """Generate an Obs object with given value, dvalue and name for test purposes
 49
 50    Parameters
 51    ----------
 52    value : float
 53        central value of the Obs to be generated.
 54    dvalue : float
 55        error of the Obs to be generated.
 56    name : str
 57        name of the ensemble for which the Obs is to be generated.
 58    samples: int
 59        number of samples for the Obs (default 1000).
 60
 61    Returns
 62    -------
 63    res : Obs
 64        Generated Observable
 65    """
 66    if dvalue <= 0.0:
 67        return Obs([np.zeros(samples) + value], [name])
 68    else:
 69        for _ in range(100):
 70            deltas = [np.random.normal(0.0, dvalue * np.sqrt(samples), samples)]
 71            deltas -= np.mean(deltas)
 72            deltas *= dvalue / np.sqrt((np.var(deltas) / samples)) / np.sqrt(1 + 3 / samples)
 73            deltas += value
 74            res = Obs(deltas, [name])
 75            res.gamma_method(S=2, tau_exp=0)
 76            if abs(res.dvalue - dvalue) < 1e-10 * dvalue:
 77                break
 78
 79        res._value = float(value)
 80
 81        return res
 82
 83
 84def gen_correlated_data(means, cov, name, tau=0.5, samples=1000):
 85    """ Generate observables with given covariance and autocorrelation times.
 86
 87    Parameters
 88    ----------
 89    means : list
 90        list containing the mean value of each observable.
 91    cov : numpy.ndarray
 92        covariance matrix for the data to be generated.
 93    name : str
 94        ensemble name for the data to be geneated.
 95    tau : float or list
 96        can either be a real number or a list with an entry for
 97        every dataset.
 98    samples : int
 99        number of samples to be generated for each observable.
100
101    Returns
102    -------
103    corr_obs : list[Obs]
104        Generated observable list
105    """
106
107    assert len(means) == cov.shape[-1]
108    tau = np.asarray(tau)
109    if np.min(tau) < 0.5:
110        raise Exception('All integrated autocorrelations have to be >= 0.5.')
111
112    a = (2 * tau - 1) / (2 * tau + 1)
113    rand = np.random.multivariate_normal(np.zeros_like(means), cov * samples, samples)
114
115    # Normalize samples such that sample variance matches input
116    norm = np.array([np.var(o, ddof=1) / samples for o in rand.T])
117    rand = rand @ np.diag(np.sqrt(np.diag(cov))) @ np.diag(1 / np.sqrt(norm))
118
119    data = [rand[0]]
120    for i in range(1, samples):
121        data.append(np.sqrt(1 - a ** 2) * rand[i] + a * data[-1])
122    corr_data = np.array(data) - np.mean(data, axis=0) + means
123    return [Obs([dat], [name]) for dat in corr_data.T]
124
125
126def _assert_equal_properties(ol, otype=Obs):
127    otype = type(ol[0])
128    for o in ol[1:]:
129        if not isinstance(o, otype):
130            raise Exception("Wrong data type in list.")
131        for attr in ["reweighted", "e_content", "idl"]:
132            if hasattr(ol[0], attr):
133                if not getattr(ol[0], attr) == getattr(o, attr):
134                    raise Exception(f"All Obs in list have to have the same state '{attr}'.")