pyerrors.misc

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