diff --git a/.github/workflows/binder.yml b/.github/workflows/binder.yml new file mode 100644 index 00000000..460f0fef --- /dev/null +++ b/.github/workflows/binder.yml @@ -0,0 +1,15 @@ +name: binder +on: + push: + branches: + - develop + +jobs: + Create-MyBinderOrg-Cache: + runs-on: ubuntu-latest + steps: + - name: cache binder build on mybinder.org + uses: jupyterhub/repo2docker-action@master + with: + NO_PUSH: true + MYBINDERORG_TAG: ${{ github.event.ref }} # This builds the container on mybinder.org with the branch that was pushed on. diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 4c074906..5e215110 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -13,15 +13,15 @@ jobs: pytest: runs-on: ${{ matrix.os }} strategy: - fail-fast: true + fail-fast: true matrix: os: [ubuntu-latest] python-version: ["3.7", "3.8", "3.9", "3.10"] include: - os: macos-latest - python-version: 3.9 + python-version: 3.9 - os: windows-latest - python-version: 3.9 + python-version: 3.9 steps: - name: Checkout source @@ -40,6 +40,7 @@ jobs: pip install pytest pip install pytest-cov pip install pytest-benchmark + pip freeze - name: Run tests run: pytest --cov=pyerrors -vv diff --git a/.gitignore b/.gitignore index fdaac56f..a8338e8f 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,10 @@ __pycache__ .cache examples/B1k2_pcac_plateau.p examples/Untitled.* +examples/pcac_plateau_test_ensemble.json.gz core.* *.swp htmlcov +build +pyerrors.egg-info +dist diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ac0ee11..0f900f57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,39 @@ All notable changes to this project will be documented in this file. +## [2.1.3] - 2022-06-13 +### Fixed +- Further bugs in connection with correlator objects which have arrays with None entries as content fixed. + +## [2.1.2] - 2022-06-10 +### Fixed +- Bug in `Corr.matrix_symmetric` fixed which appeared when a time slice contained an array with `None` entries. + +## [2.1.1] - 2022-06-06 +### Fixed +- Bug in error propagation of correlated least square fits fixed. +- `Fit_result.gamma_method` can now be called with kwargs + +## [2.1.0] - 2022-05-31 +### Added +- `obs.covariance` now has the option to smooth small eigenvalues of the matrix with the method described in hep-lat/9412087. +- `Corr.prune` was added which can reduce the size of a correlator matrix before solving the GEVP. +- `Corr.show` has two additional optional arguments. `hide_sigma` to hide data points with large errors and `references` to display reference values as horizontal lines. +- I/O routines for ALPHA dobs format added. +- `input.hadrons` functionality extended. + +### Changed +- The standard behavior of the `Corr.GEVP` method has changed. It now returns all eigenvectors of the system instead of only the specified ones as default. The standard way of sorting the eigenvectors was changed to `Eigenvalue`. The argument `sorted_list` was deprecated in favor of `sort`. +- Before performing a correlated fit the routine first runs an uncorrelated one to obtain a better guess for the initial parameters. + +### Fixed +- `obs.covariance` now also gives correct estimators if data defined on non-identical configurations is passed to the function. +- Rounding errors in estimating derivatives of fit parameters with respect to input data from the inverse hessian reduced. This should only make a difference when the magnitude of the errors of different fit parameters vary vastly. +- Bug in json.gz format fixed which did not properly store the replica mean values. Format version bumped to 1.1. +- The GEVP matrix is now symmetrized before solving the system for all sorting options not only the one with fixed `ts`. +- Automatic range estimation improved in `fits.residual_plot`. +- Bugs in `input.bdio` fixed. + ## [2.0.0] - 2022-03-31 ### Added - The possibility to work with Monte Carlo histories which are evenly or unevenly spaced was added. diff --git a/README.md b/README.md index 67d59be0..225d8128 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![pytest](https://github.com/fjosw/pyerrors/actions/workflows/pytest.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/pytest.yml) [![examples](https://github.com/fjosw/pyerrors/actions/workflows/examples.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/examples.yml) [![docs](https://github.com/fjosw/pyerrors/actions/workflows/docs.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/docs.yml) [![](https://img.shields.io/badge/python-3.7+-blue.svg)](https://www.python.org/downloads/) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) +[![flake8](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![pytest](https://github.com/fjosw/pyerrors/actions/workflows/pytest.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/pytest.yml) [![](https://img.shields.io/badge/python-3.7+-blue.svg)](https://www.python.org/downloads/) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/fjosw/pyerrors/HEAD?labpath=examples) # pyerrors `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. @@ -8,11 +8,12 @@ - **Bug reports:** https://github.com/fjosw/pyerrors/issues ## Installation -To install the current `develop` version run +To install the most recent release run +```bash +pip install pyerrors # Fresh install +pip install -U pyerrors # Upgrade +``` +to install the current `develop` version run ```bash pip install git+https://github.com/fjosw/pyerrors.git@develop ``` -to install the most recent release run -```bash -pip install git+https://github.com/fjosw/pyerrors.git@master -``` diff --git a/examples/01_basic_example.ipynb b/examples/01_basic_example.ipynb index 834f9d39..df1d76f2 100644 --- a/examples/01_basic_example.ipynb +++ b/examples/01_basic_example.ipynb @@ -21,6 +21,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import pyerrors as pe" ] @@ -32,7 +33,8 @@ "outputs": [], "source": [ "plt.style.use('./base_style.mplstyle')\n", - "plt.rc('text', usetex=True)" + "usetex = matplotlib.checkdep_usetex(True)\n", + "plt.rc('text', usetex=usetex)" ] }, { diff --git a/examples/02_correlators.ipynb b/examples/02_correlators.ipynb index a647162e..e2d2a09a 100644 --- a/examples/02_correlators.ipynb +++ b/examples/02_correlators.ipynb @@ -8,6 +8,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import pyerrors as pe" ] @@ -20,7 +21,8 @@ "outputs": [], "source": [ "plt.style.use('./base_style.mplstyle')\n", - "plt.rc('text', usetex=True)" + "usetex = matplotlib.checkdep_usetex(True)\n", + "plt.rc('text', usetex=usetex)" ] }, { @@ -480,7 +482,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -494,7 +496,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.8" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/examples/03_pcac_example.ipynb b/examples/03_pcac_example.ipynb index bf0f0b02..835062a8 100644 --- a/examples/03_pcac_example.ipynb +++ b/examples/03_pcac_example.ipynb @@ -7,6 +7,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import pyerrors as pe" ] @@ -18,7 +19,8 @@ "outputs": [], "source": [ "plt.style.use('./base_style.mplstyle')\n", - "plt.rc('text', usetex=True)" + "usetex = matplotlib.checkdep_usetex(True)\n", + "plt.rc('text', usetex=usetex)" ] }, { diff --git a/examples/04_fit_example.ipynb b/examples/04_fit_example.ipynb index 026cca2d..72d62e9e 100644 --- a/examples/04_fit_example.ipynb +++ b/examples/04_fit_example.ipynb @@ -6,9 +6,10 @@ "metadata": {}, "outputs": [], "source": [ - "import pyerrors as pe\n", "import numpy as np\n", - "import matplotlib.pyplot as plt" + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import pyerrors as pe" ] }, { @@ -18,7 +19,8 @@ "outputs": [], "source": [ "plt.style.use('./base_style.mplstyle')\n", - "plt.rc('text', usetex=True)" + "usetex = matplotlib.checkdep_usetex(True)\n", + "plt.rc('text', usetex=usetex)" ] }, { @@ -439,7 +441,7 @@ "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1" }, "kernelspec": { - "display_name": "Python 3.8.10 64-bit", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, diff --git a/examples/05_matrix_operations.ipynb b/examples/05_matrix_operations.ipynb index a49a7ee0..926d1734 100644 --- a/examples/05_matrix_operations.ipynb +++ b/examples/05_matrix_operations.ipynb @@ -393,7 +393,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -407,7 +407,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.8" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/examples/06_gevp.ipynb b/examples/06_gevp.ipynb index e969008c..268346f6 100644 --- a/examples/06_gevp.ipynb +++ b/examples/06_gevp.ipynb @@ -7,6 +7,7 @@ "metadata": {}, "outputs": [], "source": [ + "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import pyerrors as pe" ] @@ -18,7 +19,9 @@ "metadata": {}, "outputs": [], "source": [ - "plt.style.use('./base_style.mplstyle'); plt.rc('text', usetex=True)" + "plt.style.use('./base_style.mplstyle')\n", + "usetex = matplotlib.checkdep_usetex(True)\n", + "plt.rc('text', usetex=usetex)" ] }, { @@ -305,7 +308,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index 6ec89a04..ef257c86 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -8,6 +8,8 @@ It is based on the gamma method [arXiv:hep-lat/0306017](https://arxiv.org/abs/he - non-linear fits with x- and y-errors and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289). - real and complex matrix operations and their error propagation based on automatic differentiation (Matrix inverse, Cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...). +More detailed examples can found in the [GitHub repository](https://github.com/fjosw/pyerrors/tree/develop/examples) [![badge](https://img.shields.io/badge/-try%20it%20out-579ACA.svg?logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFkAAABZCAMAAABi1XidAAAB8lBMVEX///9XmsrmZYH1olJXmsr1olJXmsrmZYH1olJXmsr1olJXmsrmZYH1olL1olJXmsr1olJXmsrmZYH1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olJXmsrmZYH1olL1olL0nFf1olJXmsrmZYH1olJXmsq8dZb1olJXmsrmZYH1olJXmspXmspXmsr1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olLeaIVXmsrmZYH1olL1olL1olJXmsrmZYH1olLna31Xmsr1olJXmsr1olJXmsrmZYH1olLqoVr1olJXmsr1olJXmsrmZYH1olL1olKkfaPobXvviGabgadXmsqThKuofKHmZ4Dobnr1olJXmsr1olJXmspXmsr1olJXmsrfZ4TuhWn1olL1olJXmsqBi7X1olJXmspZmslbmMhbmsdemsVfl8ZgmsNim8Jpk8F0m7R4m7F5nLB6jbh7jbiDirOEibOGnKaMhq+PnaCVg6qWg6qegKaff6WhnpKofKGtnomxeZy3noG6dZi+n3vCcpPDcpPGn3bLb4/Mb47UbIrVa4rYoGjdaIbeaIXhoWHmZYHobXvpcHjqdHXreHLroVrsfG/uhGnuh2bwj2Hxk17yl1vzmljzm1j0nlX1olL3AJXWAAAAbXRSTlMAEBAQHx8gICAuLjAwMDw9PUBAQEpQUFBXV1hgYGBkcHBwcXl8gICAgoiIkJCQlJicnJ2goKCmqK+wsLC4usDAwMjP0NDQ1NbW3Nzg4ODi5+3v8PDw8/T09PX29vb39/f5+fr7+/z8/Pz9/v7+zczCxgAABC5JREFUeAHN1ul3k0UUBvCb1CTVpmpaitAGSLSpSuKCLWpbTKNJFGlcSMAFF63iUmRccNG6gLbuxkXU66JAUef/9LSpmXnyLr3T5AO/rzl5zj137p136BISy44fKJXuGN/d19PUfYeO67Znqtf2KH33Id1psXoFdW30sPZ1sMvs2D060AHqws4FHeJojLZqnw53cmfvg+XR8mC0OEjuxrXEkX5ydeVJLVIlV0e10PXk5k7dYeHu7Cj1j+49uKg7uLU61tGLw1lq27ugQYlclHC4bgv7VQ+TAyj5Zc/UjsPvs1sd5cWryWObtvWT2EPa4rtnWW3JkpjggEpbOsPr7F7EyNewtpBIslA7p43HCsnwooXTEc3UmPmCNn5lrqTJxy6nRmcavGZVt/3Da2pD5NHvsOHJCrdc1G2r3DITpU7yic7w/7Rxnjc0kt5GC4djiv2Sz3Fb2iEZg41/ddsFDoyuYrIkmFehz0HR2thPgQqMyQYb2OtB0WxsZ3BeG3+wpRb1vzl2UYBog8FfGhttFKjtAclnZYrRo9ryG9uG/FZQU4AEg8ZE9LjGMzTmqKXPLnlWVnIlQQTvxJf8ip7VgjZjyVPrjw1te5otM7RmP7xm+sK2Gv9I8Gi++BRbEkR9EBw8zRUcKxwp73xkaLiqQb+kGduJTNHG72zcW9LoJgqQxpP3/Tj//c3yB0tqzaml05/+orHLksVO+95kX7/7qgJvnjlrfr2Ggsyx0eoy9uPzN5SPd86aXggOsEKW2Prz7du3VID3/tzs/sSRs2w7ovVHKtjrX2pd7ZMlTxAYfBAL9jiDwfLkq55Tm7ifhMlTGPyCAs7RFRhn47JnlcB9RM5T97ASuZXIcVNuUDIndpDbdsfrqsOppeXl5Y+XVKdjFCTh+zGaVuj0d9zy05PPK3QzBamxdwtTCrzyg/2Rvf2EstUjordGwa/kx9mSJLr8mLLtCW8HHGJc2R5hS219IiF6PnTusOqcMl57gm0Z8kanKMAQg0qSyuZfn7zItsbGyO9QlnxY0eCuD1XL2ys/MsrQhltE7Ug0uFOzufJFE2PxBo/YAx8XPPdDwWN0MrDRYIZF0mSMKCNHgaIVFoBbNoLJ7tEQDKxGF0kcLQimojCZopv0OkNOyWCCg9XMVAi7ARJzQdM2QUh0gmBozjc3Skg6dSBRqDGYSUOu66Zg+I2fNZs/M3/f/Grl/XnyF1Gw3VKCez0PN5IUfFLqvgUN4C0qNqYs5YhPL+aVZYDE4IpUk57oSFnJm4FyCqqOE0jhY2SMyLFoo56zyo6becOS5UVDdj7Vih0zp+tcMhwRpBeLyqtIjlJKAIZSbI8SGSF3k0pA3mR5tHuwPFoa7N7reoq2bqCsAk1HqCu5uvI1n6JuRXI+S1Mco54YmYTwcn6Aeic+kssXi8XpXC4V3t7/ADuTNKaQJdScAAAAAElFTkSuQmCC)](https://mybinder.org/v2/gh/fjosw/pyerrors/HEAD?labpath=examples). + There exist similar publicly available implementations of gamma method error analysis suites in [Fortran](https://gitlab.ift.uam-csic.es/alberto/aderrors), [Julia](https://gitlab.ift.uam-csic.es/alberto/aderrors.jl) and [Python](https://github.com/mbruno46/pyobs). ## Basic example @@ -443,8 +445,10 @@ Julia I/O routines for the json.gz format, compatible with [ADerrors.jl](https:/ # Citing If you use `pyerrors` for research that leads to a publication please consider citing: - Ulli Wolff, *Monte Carlo errors with less errors*. Comput.Phys.Commun. 156 (2004) 143-153, Comput.Phys.Commun. 176 (2007) 383 (erratum). -- Stefan Schaefer, Rainer Sommer, Francesco Virotta, *Critical slowing down and error analysis in lattice QCD simulations*. Nucl.Phys.B 845 (2011) 93-119. - Alberto Ramos, *Automatic differentiation for error analysis of Monte Carlo data*. Comput.Phys.Commun. 238 (2019) 19-35. +and +- Stefan Schaefer, Rainer Sommer, Francesco Virotta, *Critical slowing down and error analysis in lattice QCD simulations*. Nucl.Phys.B 845 (2011) 93-119. +where applicable. ''' from .obs import * from .correlators import * diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index ab4f6f50..ecb8d525 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -155,7 +155,7 @@ class Corr: raise Exception("Vectors are of wrong shape!") if normalize: vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) - newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] + newcontent = [None if _check_for_none(self, item) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] else: # There are no checks here yet. There are so many possible scenarios, where this can go wrong. @@ -163,7 +163,7 @@ class Corr: for t in range(self.T): vector_l[t], vector_r[t] = vector_l[t] / np.sqrt((vector_l[t] @ vector_l[t])), vector_r[t] / np.sqrt(vector_r[t] @ vector_r[t]) - newcontent = [None if (self.content[t] is None or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)] + newcontent = [None if (_check_for_none(self, self.content[t]) or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)] return Corr(newcontent) def item(self, i, j): @@ -197,6 +197,8 @@ class Corr: def symmetric(self): """ Symmetrize the correlator around x0=0.""" + if self.N != 1: + raise Exception('symmetric cannot be safely applied to multi-dimensional correlators.') if self.T % 2 != 0: raise Exception("Can not symmetrize odd T") @@ -215,6 +217,8 @@ class Corr: def anti_symmetric(self): """Anti-symmetrize the correlator around x0=0.""" + if self.N != 1: + raise Exception('anti_symmetric cannot be safely applied to multi-dimensional correlators.') if self.T % 2 != 0: raise Exception("Can not symmetrize odd T") @@ -236,7 +240,7 @@ class Corr: def matrix_symmetric(self): """Symmetrizes the correlator matrices on every timeslice.""" if self.N > 1: - transposed = [None if (G is None) else G.T for G in self.content] + transposed = [None if _check_for_none(self, G) else G.T for G in self.content] return 0.5 * (Corr(transposed) + self) if self.N == 1: raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") @@ -419,13 +423,15 @@ class Corr: Can either be an Obs which is correlated with all entries of the correlator or a Corr of same length. """ + if self.N != 1: + raise Exception("Only one-dimensional correlators can be safely correlated.") new_content = [] for x0, t_slice in enumerate(self.content): - if t_slice is None: + if _check_for_none(self, t_slice): new_content.append(None) else: if isinstance(partner, Corr): - if partner.content[x0] is None: + if _check_for_none(partner, partner.content[x0]): new_content.append(None) else: new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) @@ -449,9 +455,11 @@ class Corr: the reweighting factor on all configurations in weight.idl and not on the configurations in obs[i].idl. """ + if self.N != 1: + raise Exception("Reweighting only implemented for one-dimensional correlators.") new_content = [] for t_slice in self.content: - if t_slice is None: + if _check_for_none(self, t_slice): new_content.append(None) else: new_content.append(np.array(reweight(weight, t_slice, **kwargs))) @@ -467,6 +475,8 @@ class Corr: partity : int Parity quantum number of the correlator, can be +1 or -1 """ + if self.N != 1: + raise Exception("T_symmetry only implemented for one-dimensional correlators.") if not isinstance(partner, Corr): raise Exception("T partner has to be a Corr object.") if parity not in [+1, -1]: @@ -494,6 +504,8 @@ class Corr: decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, default: symmetric """ + if self.N != 1: + raise Exception("deriv only implemented for one-dimensional correlators.") if variant == "symmetric": newcontent = [] for t in range(1, self.T - 1): @@ -546,6 +558,8 @@ class Corr: decides which definition of the finite differences derivative is used. Available choice: symmetric, improved, default: symmetric """ + if self.N != 1: + raise Exception("second_deriv only implemented for one-dimensional correlators.") if variant == "symmetric": newcontent = [] for t in range(1, self.T - 1): @@ -588,7 +602,7 @@ class Corr: if variant == 'log': newcontent = [] for t in range(self.T - 1): - if (self.content[t] is None) or (self.content[t + 1] is None): + if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): newcontent.append(None) else: newcontent.append(self.content[t] / self.content[t + 1]) @@ -608,7 +622,7 @@ class Corr: newcontent = [] for t in range(self.T - 1): - if (self.content[t] is None) or (self.content[t + 1] is None): + if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): newcontent.append(None) # Fill the two timeslices in the middle of the lattice with their predecessors elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: @@ -623,7 +637,7 @@ class Corr: elif variant == 'arccosh': newcontent = [] for t in range(1, self.T - 1): - if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): + if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0): newcontent.append(None) else: newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) @@ -758,7 +772,7 @@ class Corr: x, y, y_err = self.plottable() if hide_sigma: - hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 + hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 else: hide_from = None ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) @@ -781,7 +795,7 @@ class Corr: corr.gamma_method() x, y, y_err = corr.plottable() if hide_sigma: - hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 + hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 else: hide_from = None plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) @@ -883,12 +897,14 @@ class Corr: else: raise Exception("Unknown datatype " + str(datatype)) - def print(self, range=[0, None]): - print(self.__repr__(range)) + def print(self, print_range=None): + print(self.__repr__(print_range)) + + def __repr__(self, print_range=None): + if print_range is None: + print_range = [0, None] - def __repr__(self, range=[0, None]): content_string = "" - content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here if self.tag is not None: @@ -896,14 +912,14 @@ class Corr: if self.N != 1: return content_string - if range[1]: - range[1] += 1 + if print_range[1]: + print_range[1] += 1 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' - for i, sub_corr in enumerate(self.content[range[0]:range[1]]): + for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): if sub_corr is None: - content_string += str(i + range[0]) + '\n' + content_string += str(i + print_range[0]) + '\n' else: - content_string += str(i + range[0]) + content_string += str(i + print_range[0]) for element in sub_corr: content_string += '\t' + ' ' * int(element >= 0) + str(element) content_string += '\n' @@ -923,7 +939,7 @@ class Corr: raise Exception("Addition of Corrs with different shape") newcontent = [] for t in range(self.T): - if (self.content[t] is None) or (y.content[t] is None): + if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] + y.content[t]) @@ -932,7 +948,7 @@ class Corr: elif isinstance(y, (Obs, int, float, CObs)): newcontent = [] for t in range(self.T): - if (self.content[t] is None): + if _check_for_none(self, self.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] + y) @@ -951,7 +967,7 @@ class Corr: raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") newcontent = [] for t in range(self.T): - if (self.content[t] is None) or (y.content[t] is None): + if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] * y.content[t]) @@ -960,7 +976,7 @@ class Corr: elif isinstance(y, (Obs, int, float, CObs)): newcontent = [] for t in range(self.T): - if (self.content[t] is None): + if _check_for_none(self, self.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] * y) @@ -979,12 +995,12 @@ class Corr: raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") newcontent = [] for t in range(self.T): - if (self.content[t] is None) or (y.content[t] is None): + if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] / y.content[t]) for t in range(self.T): - if newcontent[t] is None: + if _check_for_none(self, newcontent[t]): continue if np.isnan(np.sum(newcontent[t]).value): newcontent[t] = None @@ -1003,7 +1019,7 @@ class Corr: newcontent = [] for t in range(self.T): - if (self.content[t] is None): + if _check_for_none(self, self.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] / y) @@ -1014,7 +1030,7 @@ class Corr: raise Exception('Division by zero will return undefined correlator') newcontent = [] for t in range(self.T): - if (self.content[t] is None): + if _check_for_none(self, self.content[t]): newcontent.append(None) else: newcontent.append(self.content[t] / y) @@ -1028,7 +1044,7 @@ class Corr: raise TypeError('Corr / wrong type') def __neg__(self): - newcontent = [None if (item is None) else -1. * item for item in self.content] + newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] return Corr(newcontent, prange=self.prange) def __sub__(self, y): @@ -1036,31 +1052,31 @@ class Corr: def __pow__(self, y): if isinstance(y, (Obs, int, float, CObs)): - newcontent = [None if (item is None) else item**y for item in self.content] + newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] return Corr(newcontent, prange=self.prange) else: raise TypeError('Type of exponent not supported') def __abs__(self): - newcontent = [None if (item is None) else np.abs(item) for item in self.content] + newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] return Corr(newcontent, prange=self.prange) # The numpy functions: def sqrt(self): - return self**0.5 + return self ** 0.5 def log(self): - newcontent = [None if (item is None) else np.log(item) for item in self.content] + newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] return Corr(newcontent, prange=self.prange) def exp(self): - newcontent = [None if (item is None) else np.exp(item) for item in self.content] + newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] return Corr(newcontent, prange=self.prange) def _apply_func_to_corr(self, func): - newcontent = [None if (item is None) else func(item) for item in self.content] + newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] for t in range(self.T): - if newcontent[t] is None: + if _check_for_none(self, newcontent[t]): continue if np.isnan(np.sum(newcontent[t]).value): newcontent[t] = None @@ -1222,6 +1238,11 @@ def _sort_vectors(vec_set, ts): return sorted_vec_set +def _check_for_none(corr, entry): + """Checks if entry for correlator corr is None""" + return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 + + def _GEVP_solver(Gt, G0): """Helper function for solving the GEVP and sorting the eigenvectors. diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 420d3da0..997c7f6f 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -8,7 +8,6 @@ import scipy.stats import matplotlib.pyplot as plt from matplotlib import gridspec from scipy.odr import ODR, Model, RealData -from scipy.stats import chi2 import iminuit from autograd import jacobian from autograd import elementwise_grad as egrad @@ -34,9 +33,9 @@ class Fit_result(Sequence): def __len__(self): return len(self.fit_parameters) - def gamma_method(self): + def gamma_method(self, **kwargs): """Apply the gamma method to all fit parameters""" - [o.gamma_method() for o in self.fit_parameters] + [o.gamma_method(**kwargs) for o in self.fit_parameters] def __str__(self): my_str = 'Goodness of fit:\n' @@ -95,7 +94,8 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): If true all output to the console is omitted (default False). initial_guess : list can provide an initial guess for the input parameters. Relevant for - non-linear fits with many parameters. + non-linear fits with many parameters. In case of correlated fits the guess is used to perform + an uncorrelated fit which then serves as guess for the correlated fit. method : str, optional can be used to choose an alternative method for the minimization of chisquare. The possible methods are the ones which can be used for scipy.optimize.minimize and @@ -264,7 +264,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): fitp = out.beta try: - hess_inv = np.linalg.pinv(jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel())))) + hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel()))) except TypeError: raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None @@ -275,7 +275,11 @@ def total_least_squares(x, y, func, silent=False, **kwargs): jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) - deriv_x = -hess_inv @ jac_jac_x[:n_parms + m, n_parms + m:] + # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv + try: + deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:]) + except np.linalg.LinAlgError: + raise Exception("Cannot invert hessian matrix.") def odr_chisquare_compact_y(d): model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) @@ -284,7 +288,11 @@ def total_least_squares(x, y, func, silent=False, **kwargs): jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f))) - deriv_y = -hess_inv @ jac_jac_y[:n_parms + m, n_parms + m:] + # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv + try: + deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:]) + except np.linalg.LinAlgError: + raise Exception("Cannot invert hessian matrix.") result = [] for i in range(n_parms): @@ -294,7 +302,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) output.dof = x.shape[-1] - n_parms - output.p_value = 1 - chi2.cdf(output.odr_chisquare, output.dof) + output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof) return output @@ -469,18 +477,17 @@ def _standard_fit(x, y, func, silent=False, **kwargs): if condn > 1 / np.sqrt(np.finfo(float).eps): warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning) chol = np.linalg.cholesky(corr) - chol_inv = np.linalg.inv(chol) - chol_inv = np.dot(chol_inv, covdiag) + chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True) - def chisqfunc(p): + def chisqfunc_corr(p): model = func(p, x) chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2) return chisq - else: - def chisqfunc(p): - model = func(p, x) - chisq = anp.sum(((y_f - model) / dy_f) ** 2) - return chisq + + def chisqfunc(p): + model = func(p, x) + chisq = anp.sum(((y_f - model) / dy_f) ** 2) + return chisq output.method = kwargs.get('method', 'Levenberg-Marquardt') if not silent: @@ -489,29 +496,38 @@ def _standard_fit(x, y, func, silent=False, **kwargs): if output.method != 'Levenberg-Marquardt': if output.method == 'migrad': fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4) # Stopping criterion 0.002 * tol * errordef + if kwargs.get('correlated_fit') is True: + fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=1e-4) # Stopping criterion 0.002 * tol * errordef output.iterations = fit_result.nfev else: fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12) + if kwargs.get('correlated_fit') is True: + fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=1e-12) output.iterations = fit_result.nit chisquare = fit_result.fun else: if kwargs.get('correlated_fit') is True: - def chisqfunc_residuals(p): + def chisqfunc_residuals_corr(p): model = func(p, x) chisq = anp.dot(chol_inv, (y_f - model)) return chisq - else: - def chisqfunc_residuals(p): - model = func(p, x) - chisq = ((y_f - model) / dy_f) - return chisq + def chisqfunc_residuals(p): + model = func(p, x) + chisq = ((y_f - model) / dy_f) + return chisq fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) + if kwargs.get('correlated_fit') is True: + fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) chisquare = np.sum(fit_result.fun ** 2) + if kwargs.get('correlated_fit') is True: + assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14) + else: + assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) output.iterations = fit_result.nfev @@ -542,7 +558,10 @@ def _standard_fit(x, y, func, silent=False, **kwargs): fitp = fit_result.x try: - hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(fitp)) + if kwargs.get('correlated_fit') is True: + hess = jacobian(jacobian(chisqfunc_corr))(fitp) + else: + hess = jacobian(jacobian(chisqfunc))(fitp) except TypeError: raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None @@ -560,7 +579,11 @@ def _standard_fit(x, y, func, silent=False, **kwargs): jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f))) - deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] + # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv + try: + deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:]) + except np.linalg.LinAlgError: + raise Exception("Cannot invert hessian matrix.") result = [] for i in range(n_parms): @@ -568,9 +591,9 @@ def _standard_fit(x, y, func, silent=False, **kwargs): output.fit_parameters = result - output.chisquare = chisqfunc(fit_result.x) + output.chisquare = chisquare output.dof = x.shape[-1] - n_parms - output.p_value = 1 - chi2.cdf(output.chisquare, output.dof) + output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof) if kwargs.get('resplot') is True: residual_plot(x, y, func, result) diff --git a/pyerrors/input/bdio.py b/pyerrors/input/bdio.py index f54dfbcd..2d30f3b2 100644 --- a/pyerrors/input/bdio.py +++ b/pyerrors/input/bdio.py @@ -61,7 +61,7 @@ def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): return_list = [] print('Reading of bdio file started') - while 1 > 0: + while True: bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) @@ -373,7 +373,7 @@ def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs): fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form)) print('Reading of bdio file started') - while 1 > 0: + while True: bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) if ruinfo < 0: @@ -582,7 +582,7 @@ def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs): fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form)) print('Reading of bdio file started') - while 1 > 0: + while True: bdio_seek_record(fbdio) ruinfo = bdio_get_ruinfo(fbdio) if ruinfo < 0: diff --git a/pyerrors/input/dobs.py b/pyerrors/input/dobs.py index 02418a94..65b78fec 100644 --- a/pyerrors/input/dobs.py +++ b/pyerrors/input/dobs.py @@ -648,7 +648,7 @@ def _dobsdict_to_xmlstring_spaces(d, space=' '): return o -def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags={}): +def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags=None): """Generate the string for the export of a list of Obs or structures containing Obs to a .xml.gz file according to the Zeuthen dobs format. @@ -674,6 +674,8 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N Provide alternative enstag for ensembles in the form enstags = {ename: enstag} Otherwise, the ensemble name is used. """ + if enstags is None: + enstags = {} od = {} r_names = [] for o in obsl: @@ -831,7 +833,7 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N return rs -def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags={}, gz=True): +def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=[], who=None, enstags=None, gz=True): """Export a list of Obs or structures containing Obs to a .xml.gz file according to the Zeuthen dobs format. @@ -861,6 +863,8 @@ def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=[], who=No gz : bool If True, the output is a gzipped XML. If False, the output is a XML file. """ + if enstags is None: + enstags = {} dobsstring = create_dobs_string(obsl, name, spec, origin, symbol, who, enstags=enstags) diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py index 31106bda..4645330d 100644 --- a/pyerrors/input/hadrons.py +++ b/pyerrors/input/hadrons.py @@ -55,8 +55,8 @@ def _get_files(path, filestem, idl): return filtered_files, idx -def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None): - """Read hadrons meson hdf5 file and extract the meson labeled 'meson' +def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None, gammas=None): + r'''Read hadrons meson hdf5 file and extract the meson labeled 'meson' Parameters ----------------- @@ -69,13 +69,31 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', idl=None): meson : str label of the meson to be extracted, standard value meson_0 which corresponds to the pseudoscalar pseudoscalar two-point function. + gammas : tuple of strings + Instrad of a meson label one can also provide a tuple of two strings + indicating the gamma matrices at source and sink. + ("Gamma5", "Gamma5") corresponds to the pseudoscalar pseudoscalar + two-point function. The gammas argument dominateds over meson. idl : range If specified only configurations in the given range are read in. - """ + ''' files, idx = _get_files(path, filestem, idl) tree = meson.rsplit('_')[0] + if gammas is not None: + h5file = h5py.File(path + '/' + files[0], "r") + found_meson = None + for key in h5file[tree].keys(): + if gammas[0] == h5file[tree][key].attrs["gamma_snk"][0].decode() and h5file[tree][key].attrs["gamma_src"][0].decode() == gammas[1]: + found_meson = key + break + h5file.close() + if found_meson: + meson = found_meson + else: + raise Exception("Source Sink combination " + str(gammas) + " not found.") + corr_data = [] infos = [] for hd5_file in files: @@ -153,10 +171,6 @@ def read_DistillationContraction_hd5(path, ens_id, diagrams=["direct"], idl=None identifier = tuple(identifier) # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. - check_traj = h5file["DistillationContraction/Metadata"].attrs.get("Traj")[0] - - assert check_traj == n_traj - for diagram in diagrams: real_data = np.zeros(Nt) for x0 in range(Nt): diff --git a/pyerrors/input/json.py b/pyerrors/input/json.py index ac7963cc..69fa2c2d 100644 --- a/pyerrors/input/json.py +++ b/pyerrors/input/json.py @@ -165,7 +165,7 @@ def create_json_string(ol, description='', indent=1): d = {} d['program'] = 'pyerrors %s' % (pyerrorsversion.__version__) - d['version'] = '1.0' + d['version'] = '1.1' d['who'] = getpass.getuser() d['date'] = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S %z') d['host'] = socket.gethostname() + ', ' + platform.platform() @@ -294,6 +294,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): if od: ret = Obs([[ddi[0] + values[0] for ddi in di] for di in od['deltas']], od['names'], idl=od['idl']) + ret._value = values[0] ret.is_merged = od['is_merged'] else: ret = Obs([], [], means=[]) @@ -319,6 +320,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): for i in range(layout): if od: ret.append(Obs([list(di[:, i] + values[i]) for di in od['deltas']], od['names'], idl=od['idl'])) + ret[-1]._value = values[i] ret[-1].is_merged = od['is_merged'] else: ret.append(Obs([], [], means=[])) @@ -346,6 +348,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): for i in range(N): if od: ret.append(Obs([di[:, i] + values[i] for di in od['deltas']], od['names'], idl=od['idl'])) + ret[-1]._value = values[i] ret[-1].is_merged = od['is_merged'] else: ret.append(Obs([], [], means=[])) diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 6d7b6d50..fe75fb1a 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -17,8 +17,6 @@ def read_pbp(path, prefix, **kwargs): list which contains the last config to be read for each replicum """ - extract_nfct = 1 - ls = [] for (dirpath, dirnames, filenames) in os.walk(path): ls.extend(filenames) @@ -78,14 +76,10 @@ def read_pbp(path, prefix, **kwargs): # This block is necessary for openQCD1.6 ms1 files nfct = [] - if extract_nfct == 1: - for i in range(nrw): - t = fp.read(4) - nfct.append(struct.unpack('i', t)[0]) - print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting - else: - for i in range(nrw): - nfct.append(1) + for i in range(nrw): + t = fp.read(4) + nfct.append(struct.unpack('i', t)[0]) + print('nfct: ', nfct) # Hasenbusch factor, 1 for rat reweighting nsrc = [] for i in range(nrw): @@ -93,7 +87,7 @@ def read_pbp(path, prefix, **kwargs): nsrc.append(struct.unpack('i', t)[0]) # body - while 0 < 1: + while True: t = fp.read(4) if len(t) < 4: break diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 05e9fdb7..c76c0596 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -150,7 +150,7 @@ def read_rwms(path, prefix, version='2.0', names=None, **kwargs): print('something is wrong!') configlist.append([]) - while 0 < 1: + while True: t = fp.read(4) if len(t) < 4: break @@ -362,7 +362,7 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwar Ysl = [] configlist.append([]) - while 0 < 1: + while True: t = fp.read(4) if(len(t) < 4): break @@ -657,7 +657,7 @@ def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): nfl = 1 iobs = 8 * nfl # number of flow observables calculated - while 0 < 1: + while True: t = fp.read(4) if(len(t) < 4): break @@ -682,7 +682,7 @@ def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs): t = fp.read(8) eps = struct.unpack('d', t)[0] - while 0 < 1: + while True: t = fp.read(4) if(len(t) < 4): break diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 3114bcd6..d211c541 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1,5 +1,7 @@ import warnings import pickle +from math import gcd +from functools import reduce import numpy as np import autograd.numpy as anp # Thinly-wrapped numpy from autograd import jacobian @@ -981,6 +983,39 @@ def _merge_idx(idl): return sorted(set().union(*idl)) +def _intersection_idx(idl): + """Returns the intersection of all lists in idl as sorted list + + Parameters + ---------- + idl : list + List of lists or ranges. + """ + + def _lcm(*args): + """Returns the lowest common multiple of args. + + From python 3.9 onwards the math library contains an lcm function.""" + return reduce(lambda a, b: a * b // gcd(a, b), args) + + # Use groupby to efficiently check whether all elements of idl are identical + try: + g = groupby(idl) + if next(g, True) and not next(g, False): + return idl[0] + except Exception: + pass + + if np.all([type(idx) is range for idx in idl]): + if len(set([idx[0] for idx in idl])) == 1: + idstart = max([idx.start for idx in idl]) + idstop = min([idx.stop for idx in idl]) + idstep = _lcm(*[idx.step for idx in idl]) + return range(idstart, idstop, idstep) + + return sorted(set.intersection(*[set(o) for o in idl])) + + def _expand_deltas_for_merge(deltas, idx, shape, new_idx): """Expand deltas defined on idx to the list of configs that is defined by new_idx. New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest @@ -1008,6 +1043,34 @@ def _expand_deltas_for_merge(deltas, idx, shape, new_idx): return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) +def _collapse_deltas_for_merge(deltas, idx, shape, new_idx): + """Collapse deltas defined on idx to the list of configs that is defined by new_idx. + If idx and new_idx are of type range, the smallest + common divisor of the step sizes is used as new step size. + + Parameters + ---------- + deltas : list + List of fluctuations + idx : list + List or range of configs on which the deltas are defined. + Has to be a subset of new_idx and has to be sorted in ascending order. + shape : list + Number of configs in idx. + new_idx : list + List of configs that defines the new range, has to be sorted in ascending order. + """ + + if type(idx) is range and type(new_idx) is range: + if idx == new_idx: + return deltas + ret = np.zeros(new_idx[-1] - new_idx[0] + 1) + for i in range(shape): + if idx[i] in new_idx: + ret[idx[i] - new_idx[0]] = deltas[i] + return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) + + def _filter_zeroes(deltas, idx, eps=Obs.filter_eps): """Filter out all configurations with vanishing fluctuation such that they do not contribute to the error estimate anymore. Returns the new deltas and @@ -1389,13 +1452,6 @@ def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): if isinstance(smooth, int): corr = _smooth_eigenvalues(corr, smooth) - errors = [o.dvalue for o in obs] - cov = np.diag(errors) @ corr @ np.diag(errors) - - eigenvalues = np.linalg.eigh(cov)[0] - if not np.all(eigenvalues >= 0): - warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning) - if visualize: plt.matshow(corr, vmin=-1, vmax=1) plt.set_cmap('RdBu') @@ -1404,8 +1460,15 @@ def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): if correlation is True: return corr - else: - return cov + + errors = [o.dvalue for o in obs] + cov = np.diag(errors) @ corr @ np.diag(errors) + + eigenvalues = np.linalg.eigh(cov)[0] + if not np.all(eigenvalues >= 0): + warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning) + + return cov def _smooth_eigenvalues(corr, E): @@ -1429,8 +1492,8 @@ def _covariance_element(obs1, obs2): """Estimates the covariance of two Obs objects, neglecting autocorrelations.""" def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx): - deltas1 = _expand_deltas_for_merge(deltas1, idx1, len(idx1), new_idx) - deltas2 = _expand_deltas_for_merge(deltas2, idx2, len(idx2), new_idx) + deltas1 = _collapse_deltas_for_merge(deltas1, idx1, len(idx1), new_idx) + deltas2 = _collapse_deltas_for_merge(deltas2, idx2, len(idx2), new_idx) return np.sum(deltas1 * deltas2) if set(obs1.names).isdisjoint(set(obs2.names)): @@ -1450,29 +1513,30 @@ def _covariance_element(obs1, obs2): for r_name in obs1.e_content[e_name]: if r_name not in obs2.e_content[e_name]: continue - idl_d[r_name] = _merge_idx([obs1.idl[r_name], obs2.idl[r_name]]) + idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]]) gamma = 0.0 for r_name in obs1.e_content[e_name]: if r_name not in obs2.e_content[e_name]: continue + if len(idl_d[r_name]) == 0: + continue gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name]) if gamma == 0.0: continue gamma_div = 0.0 - e_N = 0 for r_name in obs1.e_content[e_name]: if r_name not in obs2.e_content[e_name]: continue - gamma_div += calc_gamma(np.ones(obs1.shape[r_name]), np.ones(obs2.shape[r_name]), obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name]) - e_N += len(idl_d[r_name]) - gamma /= max(gamma_div, 1.0) + if len(idl_d[r_name]) == 0: + continue + 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])) + gamma /= gamma_div - # Bias correction hep-lat/0306017 eq. (49) - dvalue += (1 + 1 / e_N) * gamma / e_N + dvalue += gamma for e_name in obs1.cov_names: diff --git a/pyerrors/version.py b/pyerrors/version.py index 6cb8ec9d..0d8c75f0 100644 --- a/pyerrors/version.py +++ b/pyerrors/version.py @@ -1 +1 @@ -__version__ = "2.1.0+dev" +__version__ = "2.2.0+dev" diff --git a/setup.py b/setup.py index bc4533ef..0c00aad5 100644 --- a/setup.py +++ b/setup.py @@ -1,16 +1,40 @@ from setuptools import setup, find_packages from pathlib import Path +from distutils.util import convert_path + this_directory = Path(__file__).parent long_description = (this_directory / "README.md").read_text() +version = {} +with open(convert_path('pyerrors/version.py')) as ver_file: + exec(ver_file.read(), version) + setup(name='pyerrors', - version='2.1.0+dev', + version=version['__version__'], description='Error analysis for lattice QCD', long_description=long_description, long_description_content_type='text/markdown', + url="https://github.com/fjosw/pyerrors", + project_urls= { + 'Documentation': 'https://fjosw.github.io/pyerrors/pyerrors.html', + 'Bug Tracker': 'https://github.com/fjosw/pyerrors/issues', + 'Changelog' : 'https://github.com/fjosw/pyerrors/blob/master/CHANGELOG.md' + }, author='Fabian Joswig', author_email='fabian.joswig@ed.ac.uk', + license="MIT", packages=find_packages(), python_requires='>=3.6.0', - install_requires=['numpy>=1.16', 'autograd>=1.4', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit>=2', 'h5py', 'lxml', 'python-rapidjson'] + install_requires=['numpy>=1.16', 'autograd>=1.4', 'numdifftools', 'matplotlib>=3.3', 'scipy>=1', 'iminuit>=2', 'h5py>=3', 'lxml>=4', 'python-rapidjson>=1'], + classifiers=[ + 'Development Status :: 5 - Production/Stable', + 'Intended Audience :: Science/Research', + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python :: 3', + 'Programming Language :: Python :: 3.7', + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', + 'Programming Language :: Python :: 3.10', + 'Topic :: Scientific/Engineering :: Physics' + ], ) diff --git a/tests/correlators_test.py b/tests/correlators_test.py index 97bbebb1..3cfa230b 100644 --- a/tests/correlators_test.py +++ b/tests/correlators_test.py @@ -1,6 +1,7 @@ import os import numpy as np import scipy +import matplotlib.pyplot as plt import pyerrors as pe import pytest @@ -246,6 +247,21 @@ def test_matrix_corr(): corr_mat.Eigenvalue(2, state=0) +def test_corr_none_entries(): + a = pe.pseudo_Obs(1.0, 0.1, 'a') + l = np.asarray([[a, a], [a, a]]) + n = np.asarray([[None, None], [None, None]]) + x = [l, n] + matr = pe.Corr(x) + matr.projected(np.asarray([1.0, 0.0])) + + matr * 2 - 2 * matr + matr * matr + matr ** 2 / matr + + for func in [np.sqrt, np.log, np.exp, np.sin, np.cos, np.tan, np.sinh, np.cosh, np.tanh]: + func(matr) + + def test_GEVP_warnings(): corr_aa = _gen_corr(1) corr_ab = 0.5 * corr_aa @@ -332,6 +348,15 @@ def test_matrix_symmetric(): assert np.all([np.all(o == o.T) for o in sym_corr_mat]) + t_obs = pe.pseudo_Obs(1.0, 0.1, 'test') + o_mat = np.array([[t_obs, t_obs], [t_obs, t_obs]]) + corr1 = pe.Corr([o_mat, None, o_mat]) + corr2 = pe.Corr([o_mat, np.array([[None, None], [None, None]]), o_mat]) + corr3 = pe.Corr([o_mat, np.array([[t_obs, None], [None, t_obs]], dtype=object), o_mat]) + corr1.matrix_symmetric() + corr2.matrix_symmetric() + corr3.matrix_symmetric() + def test_GEVP_solver(): @@ -347,6 +372,17 @@ def test_GEVP_solver(): assert np.allclose(sp_vecs, pe.correlators._GEVP_solver(mat1, mat2), atol=1e-14) +def test_GEVP_none_entries(): + t_obs = pe.pseudo_Obs(1.0, 0.1, 'test') + t_obs2 = pe.pseudo_Obs(0.1, 0.1, 'test') + + o_mat = np.array([[t_obs, t_obs2], [t_obs2, t_obs2]]) + n_arr = np.array([[None, None], [None, None]]) + + corr = pe.Corr([o_mat, o_mat, o_mat, o_mat, o_mat, o_mat, None, o_mat, n_arr, None, o_mat]) + corr.GEVP(t0=2) + + def test_hankel(): corr_content = [] for t in range(8): @@ -405,6 +441,7 @@ def test_spaghetti_plot(): corr.spaghetti_plot(True) corr.spaghetti_plot(False) + plt.close('all') def _gen_corr(val, samples=2000): diff --git a/tests/fits_test.py b/tests/fits_test.py index 374550e4..8dd6677d 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -149,8 +149,10 @@ def test_correlated_fit(): return p[1] * anp.exp(-p[0] * x) fitp = pe.least_squares(x, data, fitf, expected_chisquare=True) + assert np.isclose(fitp.chisquare / fitp.dof, fitp.chisquare_by_dof, atol=1e-14) fitpc = pe.least_squares(x, data, fitf, correlated_fit=True) + assert np.isclose(fitpc.chisquare / fitpc.dof, fitpc.chisquare_by_dof, atol=1e-14) for i in range(2): diff = fitp[i] - fitpc[i] diff.gamma_method() @@ -171,12 +173,35 @@ def test_fit_corr_independent(): y = a[0] * anp.exp(-a[1] * x) return y - out = pe.least_squares(x, oy, func) - out_corr = pe.least_squares(x, oy, func, correlated_fit=True) + for method in ["Levenberg-Marquardt", "migrad", "Nelder-Mead"]: + out = pe.least_squares(x, oy, func, method=method) + out_corr = pe.least_squares(x, oy, func, correlated_fit=True, method=method) - assert np.isclose(out.chisquare, out_corr.chisquare) - assert (out[0] - out_corr[0]).is_zero(atol=1e-5) - assert (out[1] - out_corr[1]).is_zero(atol=1e-5) + assert np.isclose(out.chisquare, out_corr.chisquare) + assert np.isclose(out.dof, out_corr.dof) + assert np.isclose(out.chisquare_by_dof, out_corr.chisquare_by_dof) + assert (out[0] - out_corr[0]).is_zero(atol=1e-5) + assert (out[1] - out_corr[1]).is_zero(atol=1e-5) + + +def test_linear_fit_guesses(): + for err in [10, 0.1, 0.001]: + xvals = [] + yvals = [] + for x in range(1, 8, 2): + xvals.append(x) + yvals.append(pe.pseudo_Obs(x + np.random.normal(0.0, err), err, 'test1') + pe.pseudo_Obs(0, err / 100, 'test2', samples=87)) + lin_func = lambda a, x: a[0] + a[1] * x + with pytest.raises(Exception): + pe.least_squares(xvals, yvals, lin_func) + [o.gamma_method() for o in yvals]; + with pytest.raises(Exception): + pe.least_squares(xvals, yvals, lin_func, initial_guess=[5]) + + bad_guess = pe.least_squares(xvals, yvals, lin_func, initial_guess=[999, 999]) + good_guess = pe.least_squares(xvals, yvals, lin_func, initial_guess=[0, 1]) + assert np.isclose(bad_guess.chisquare, good_guess.chisquare, atol=1e-8) + assert np.all([(go - ba).is_zero(atol=1e-6) for (go, ba) in zip(good_guess, bad_guess)]) def test_total_least_squares(): @@ -218,7 +243,7 @@ def test_total_least_squares(): beta[i].gamma_method(S=1.0) assert math.isclose(beta[i].value, output.beta[i], rel_tol=1e-5) assert math.isclose(output.cov_beta[i, i], beta[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(beta[i].dvalue ** 2) - assert math.isclose(pe.covariance([beta[0], beta[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1) + assert math.isclose(pe.covariance([beta[0], beta[1]])[0, 1], output.cov_beta[0, 1], rel_tol=3.5e-1) out = pe.total_least_squares(ox, oy, func, const_par=[beta[1]]) @@ -241,7 +266,7 @@ def test_total_least_squares(): betac[i].gamma_method(S=1.0) assert math.isclose(betac[i].value, output.beta[i], rel_tol=1e-3) assert math.isclose(output.cov_beta[i, i], betac[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(betac[i].dvalue ** 2) - assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1) + assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=3.5e-1) outc = pe.total_least_squares(oxc, oyc, func, const_par=[betac[1]]) @@ -256,7 +281,7 @@ def test_total_least_squares(): betac[i].gamma_method(S=1.0) assert math.isclose(betac[i].value, output.beta[i], rel_tol=1e-3) assert math.isclose(output.cov_beta[i, i], betac[i].dvalue ** 2, rel_tol=2.5e-1), str(output.cov_beta[i, i]) + ' ' + str(betac[i].dvalue ** 2) - assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=2.5e-1) + assert math.isclose(pe.covariance([betac[0], betac[1]])[0, 1], output.cov_beta[0, 1], rel_tol=3.5e-1) outc = pe.total_least_squares(oxc, oy, func, const_par=[betac[1]]) @@ -376,6 +401,80 @@ def test_error_band(): pe.fits.error_band(x, f, fitp.fit_parameters) +def test_fit_vs_jackknife(): + od = 0.9999999999 + cov1 = np.array([[1, od, od], [od, 1.0, od], [od, od, 1.0]]) + cov1 *= 0.05 + nod = -0.4 + cov2 = np.array([[1, nod, nod], [nod, 1.0, nod], [nod, nod, 1.0]]) + cov2 *= 0.05 + cov3 = np.identity(3) + cov3 *= 0.05 + samples = 500 + + for i, cov in enumerate([cov1, cov2, cov3]): + dat = pe.misc.gen_correlated_data(np.arange(1, 4), cov, 'test', 0.5, samples=samples) + [o.gamma_method(S=0) for o in dat]; + func = lambda a, x: a[0] + a[1] * x + fr = pe.least_squares(np.arange(1, 4), dat, func) + fr.gamma_method(S=0) + + jd = np.array([o.export_jackknife() for o in dat]).T + jfr = [] + for jacks in jd: + + def chisqfunc_residuals(p): + model = func(p, np.arange(1, 4)) + chisq = ((jacks - model) / [o.dvalue for o in dat]) + return chisq + + tf = scipy.optimize.least_squares(chisqfunc_residuals, [0.0, 0.0], method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) + jfr.append(tf.x) + ajfr = np.array(jfr).T + err = np.array([np.sqrt(np.var(ajfr[j][1:], ddof=0) * (samples - 1)) for j in range(2)]) + assert np.allclose(err, [o.dvalue for o in fr], atol=1e-8) + +def test_correlated_fit_vs_jackknife(): + od = 0.999999 + cov1 = np.array([[1, od, od], [od, 1.0, od], [od, od, 1.0]]) + cov1 *= 0.1 + nod = -0.44 + cov2 = np.array([[1, nod, nod], [nod, 1.0, nod], [nod, nod, 1.0]]) + cov2 *= 0.1 + cov3 = np.identity(3) + cov3 *= 0.01 + + samples = 250 + x_val = np.arange(1, 6, 2) + for i, cov in enumerate([cov1, cov2, cov3]): + dat = pe.misc.gen_correlated_data(x_val + x_val ** 2 + np.random.normal(0.0, 0.1, 3), cov, 'test', 0.5, samples=samples) + [o.gamma_method(S=0) for o in dat]; + dat + func = lambda a, x: a[0] * x + a[1] * x ** 2 + fr = pe.least_squares(x_val, dat, func, correlated_fit=True, silent=True) + [o.gamma_method(S=0) for o in fr] + + cov = pe.covariance(dat) + chol = np.linalg.cholesky(cov) + chol_inv = np.linalg.inv(chol) + + jd = np.array([o.export_jackknife() for o in dat]).T + jfr = [] + for jacks in jd: + + def chisqfunc_residuals(p): + model = func(p, x_val) + chisq = np.dot(chol_inv, (jacks - model)) + return chisq + + tf = scipy.optimize.least_squares(chisqfunc_residuals, [0.0, 0.0], method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) + jfr.append(tf.x) + ajfr = np.array(jfr).T + err = np.array([np.sqrt(np.var(ajfr[j][1:], ddof=0) * (samples - 1)) for j in range(2)]) + assert np.allclose(err, [o.dvalue for o in fr], atol=1e-7) + assert np.allclose(ajfr.T[0], [o.value for o in fr], atol=1e-8) + + def test_fit_no_autograd(): dim = 10 x = np.arange(dim) diff --git a/tests/json_io_test.py b/tests/json_io_test.py index 00de3bfc..5474c8ce 100644 --- a/tests/json_io_test.py +++ b/tests/json_io_test.py @@ -354,3 +354,55 @@ def test_dobsio(): if isinstance(ol[i], pe.Obs): for name in ol[i].r_values: assert(np.isclose(ol[i].r_values[name], rl[i].r_values[name])) + + +def test_reconstruct_non_linear_r_obs(tmp_path): + to = pe.Obs([np.random.rand(500), np.random.rand(500), np.random.rand(111)], + ["e|r1", "e|r2", "my_new_ensemble_54^£$|8'[@124435%6^7&()~#"], + idl=[range(1, 501), range(0, 500), range(1, 999, 9)]) + to = np.log(to ** 2) / to + to.dump((tmp_path / "test_equality").as_posix()) + ro = pe.input.json.load_json((tmp_path / "test_equality").as_posix()) + assert assert_equal_Obs(to, ro) + + +def test_reconstruct_non_linear_r_obs_list(tmp_path): + to = pe.Obs([np.random.rand(500), np.random.rand(500), np.random.rand(111)], + ["e|r1", "e|r2", "my_new_ensemble_54^£$|8'[@124435%6^7&()~#"], + idl=[range(1, 501), range(0, 500), range(1, 999, 9)]) + to = np.log(to ** 2) / to + for to_list in [[to, to, to], np.array([to, to, to])]: + pe.input.json.dump_to_json(to_list, (tmp_path / "test_equality_list").as_posix()) + ro_list = pe.input.json.load_json((tmp_path / "test_equality_list").as_posix()) + for oa, ob in zip(to_list, ro_list): + assert assert_equal_Obs(oa, ob) + + +def assert_equal_Obs(to, ro): + for kw in ["N", "cov_names", "covobs", "ddvalue", "dvalue", "e_content", + "e_names", "idl", "mc_names", "names", + "reweighted", "shape", "tag"]: + if not getattr(to, kw) == getattr(ro, kw): + print(kw, "does not match.") + return False + + for kw in ["value"]: + if not np.isclose(getattr(to, kw), getattr(ro, kw), atol=1e-14): + print(kw, "does not match.") + return False + + + for kw in ["r_values", "deltas"]: + for (k, v), (k2, v2) in zip(getattr(to, kw).items(), getattr(ro, kw).items()): + assert k == k2 + if not np.allclose(v, v2, atol=1e-14): + print(kw, "does not match.") + return False + + m_to = getattr(to, "is_merged") + m_ro = getattr(ro, "is_merged") + if not m_to == m_ro: + if not (all(value is False for value in m_ro.values()) and all(value is False for value in m_to.values())): + print("is_merged", "does not match.") + return False + return True diff --git a/tests/obs_test.py b/tests/obs_test.py index 381ecdcb..5b6a8509 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -1,6 +1,7 @@ import autograd.numpy as np import os import copy +import matplotlib.pyplot as plt import pyerrors as pe import pytest @@ -56,6 +57,7 @@ def test_Obs_exceptions(): one.gamma_method() with pytest.raises(Exception): one.plot_piechart() + plt.close('all') def test_dump(): value = np.random.normal(5, 10) @@ -368,6 +370,7 @@ def test_utils(): assert my_obs < (my_obs + 1) float(my_obs) str(my_obs) + plt.close('all') def test_cobs(): @@ -515,6 +518,35 @@ def test_merge_idx(): assert pe.obs._merge_idx([range(500, 6050, 50), range(500, 6250, 250)]) == range(500, 6250, 50) +def test_intersection_idx(): + assert pe.obs._intersection_idx([range(1, 100), range(1, 100), range(1, 100)]) == range(1, 100) + assert pe.obs._intersection_idx([range(1, 100, 10), range(1, 100, 2)]) == range(1, 100, 10) + assert pe.obs._intersection_idx([range(10, 1010, 10), range(10, 1010, 50)]) == range(10, 1010, 50) + assert pe.obs._intersection_idx([range(500, 6050, 50), range(500, 6250, 250)]) == range(500, 6050, 250) + + for ids in [[list(range(1, 80, 3)), list(range(1, 100, 2))], [range(1, 80, 3), range(1, 100, 2), range(1, 100, 7)]]: + assert list(pe.obs._intersection_idx(ids)) == pe.obs._intersection_idx([list(o) for o in ids]) + +def test_merge_intersection(): + for idl_list in [[range(1, 100), range(1, 100), range(1, 100)], + [range(4, 80, 6), range(4, 80, 6)], + [[0, 2, 8, 19, 205], [0, 2, 8, 19, 205]]]: + assert pe.obs._merge_idx(idl_list) == pe.obs._intersection_idx(idl_list) + + +def test_intersection_collapse(): + range1 = range(1, 2000, 2) + range2 = range(2, 2001, 8) + + obs1 = pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + obs_merge = obs1 + pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2]) + + intersection = pe.obs._intersection_idx([o.idl["ens"] for o in [obs1, obs_merge]]) + coll = pe.obs._collapse_deltas_for_merge(obs_merge.deltas["ens"], obs_merge.idl["ens"], len(obs_merge.idl["ens"]), range1) + + assert np.all(coll == obs1.deltas["ens"]) + + def test_irregular_error_propagation(): obs_list = [pe.Obs([np.random.rand(100)], ['t']), pe.Obs([np.random.rand(50)], ['t'], idl=[range(1, 100, 2)]), @@ -619,6 +651,26 @@ def test_covariance_is_variance(): assert np.isclose(test_obs.dvalue ** 2, pe.covariance([test_obs, test_obs])[0, 1]) +def test_covariance_vs_numpy(): + N = 1078 + data1 = np.random.normal(2.5, 0.2, N) + data2 = np.random.normal(0.5, 0.08, N) + data3 = np.random.normal(-178, 5, N) + uncorr = np.row_stack([data1, data2, data3]) + corr = np.random.multivariate_normal([0.0, 17, -0.0487], [[1.0, 0.6, -0.22], [0.6, 0.8, 0.01], [-0.22, 0.01, 1.9]], N).T + + for X in [uncorr, corr]: + obs1 = pe.Obs([X[0]], ["ens1"]) + obs2 = pe.Obs([X[1]], ["ens1"]) + obs3 = pe.Obs([X[2]], ["ens1"]) + obs1.gamma_method(S=0.0) + obs2.gamma_method(S=0.0) + obs3.gamma_method(S=0.0) + pe_cov = pe.covariance([obs1, obs2, obs3]) + np_cov = np.cov(X) / N + assert np.allclose(pe_cov, np_cov, atol=1e-14) + + def test_covariance_symmetry(): value1 = np.random.normal(5, 10) dvalue1 = np.abs(np.random.normal(0, 1)) @@ -687,6 +739,23 @@ def test_covariance_factorizing(): assert np.isclose(pe.covariance([mt0, tt[1]])[0, 1], -pe.covariance(tt)[0, 1]) +def test_covariance_smooth_eigenvalues(): + for c_coeff in range(0, 14, 2): + length = 14 + sm = 5 + t_fac = 1.5 + tt = pe.misc.gen_correlated_data(np.zeros(length), 1 - 0.1 ** c_coeff * np.ones((length, length)) + 0.1 ** c_coeff * np.diag(np.ones(length)), 'test', tau=0.5 + t_fac * np.random.rand(length), samples=200) + [o.gamma_method() for o in tt] + full_corr = pe.covariance(tt, correlation=True) + cov = pe.covariance(tt, smooth=sm, correlation=True) + + full_evals = np.linalg.eigh(full_corr)[0] + sm_length = np.where(full_evals < np.mean(full_evals[:-sm]))[0][-1] + + evals = np.linalg.eigh(cov)[0] + assert np.all(np.isclose(evals[:sm_length], evals[0], atol=1e-8)) + + def test_covariance_alternation(): length = 12 t_fac = 2.5 @@ -729,6 +798,86 @@ def test_covariance_idl(): pe.covariance([obs1, obs2]) +def test_correlation_intersection_of_idls(): + range1 = range(1, 2000, 2) + range2 = range(2, 2001, 2) + + obs1 = pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + obs2_a = 0.4 * pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + 0.6 * obs1 + obs1.gamma_method() + obs2_a.gamma_method() + + cov1 = pe.covariance([obs1, obs2_a]) + corr1 = pe.covariance([obs1, obs2_a], correlation=True) + + obs2_b = obs2_a + pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2]) + obs2_b.gamma_method() + + cov2 = pe.covariance([obs1, obs2_b]) + corr2 = pe.covariance([obs1, obs2_b], correlation=True) + + assert np.isclose(corr1[0, 1], corr2[0, 1], atol=1e-14) + assert cov1[0, 1] > cov2[0, 1] + + obs2_c = pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2]) + obs2_c.gamma_method() + assert np.isclose(0, pe.covariance([obs1, obs2_c])[0, 1], atol=1e-14) + + +def test_covariance_non_identical_objects(): + obs1 = pe.Obs([np.random.normal(1.0, 0.1, 1000), np.random.normal(1.0, 0.1, 1000), np.random.normal(1.0, 0.1, 732)], ["ens|r1", "ens|r2", "ens2"]) + obs1.gamma_method() + obs2 = obs1 + 1e-18 + obs2.gamma_method() + assert obs1 == obs2 + assert obs1 is not obs2 + assert np.allclose(np.ones((2, 2)), pe.covariance([obs1, obs2], correlation=True), atol=1e-14) + + +def test_covariance_additional_non_overlapping_data(): + range1 = range(1, 20, 2) + + data2 = np.random.normal(0.0, 0.1, len(range1)) + + obs1 = pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + obs2_a = pe.Obs([data2], ["ens"], idl=[range1]) + obs1.gamma_method() + obs2_a.gamma_method() + + corr1 = pe.covariance([obs1, obs2_a], correlation=True) + + added_data = np.random.normal(0.0, 0.1, len(range1)) + added_data -= np.mean(added_data) - np.mean(data2) + data2_extended = np.ravel([data2, added_data], 'F') + + obs2_b = pe.Obs([data2_extended], ["ens"]) + obs2_b.gamma_method() + + corr2 = pe.covariance([obs1, obs2_b], correlation=True) + + assert np.isclose(corr1[0, 1], corr2[0, 1], atol=1e-14) + + +def test_coavariance_reorder_non_overlapping_data(): + range1 = range(1, 20, 2) + range2 = range(1, 41, 2) + + obs1 = pe.Obs([np.random.normal(1.0, 0.1, len(range1))], ["ens"], idl=[range1]) + obs2_b = pe.Obs([np.random.normal(1.0, 0.1, len(range2))], ["ens"], idl=[range2]) + obs1.gamma_method() + obs2_b.gamma_method() + + corr1 = pe.covariance([obs1, obs2_b], correlation=True) + + deltas = list(obs2_b.deltas['ens'][:len(range1)]) + sorted(obs2_b.deltas['ens'][len(range1):]) + obs2_a = pe.Obs([obs2_b.value + np.array(deltas)], ["ens"], idl=[range2]) + obs2_a.gamma_method() + + corr2 = pe.covariance([obs1, obs2_a], correlation=True) + + assert np.isclose(corr1[0, 1], corr2[0, 1], atol=1e-14) + + def test_empty_obs(): o = pe.Obs([np.random.rand(100)], ['test']) q = o + pe.Obs([], [], means=[])