pyerrors.misc

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