From 9805ee2b4b4297d523c1f0890279bad802a23e54 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 11:24:56 +0000 Subject: [PATCH 01/12] pyerrors.py renamed to obs.py, h5py added as dependency --- pyerrors/__init__.py | 2 +- pyerrors/correlators.py | 2 +- pyerrors/fits.py | 2 +- pyerrors/input/bdio.py | 2 +- pyerrors/input/hadrons.py | 2 +- pyerrors/input/misc.py | 2 +- pyerrors/input/openQCD.py | 2 +- pyerrors/input/sfcf.py | 2 +- pyerrors/linalg.py | 2 +- pyerrors/misc.py | 2 +- pyerrors/mpm.py | 4 +--- pyerrors/{pyerrors.py => obs.py} | 0 pyerrors/roots.py | 2 +- setup.py | 2 +- 14 files changed, 13 insertions(+), 15 deletions(-) rename pyerrors/{pyerrors.py => obs.py} (100%) diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index 855606da..db7b6859 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -1,4 +1,4 @@ -from .pyerrors import * +from .obs import * from .correlators import * from .fits import * from . import dirac diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 1d11bedb..8459f1f5 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -3,7 +3,7 @@ import numpy as np import autograd.numpy as anp import matplotlib.pyplot as plt import scipy.linalg -from .pyerrors import Obs, dump_object, reweight, correlate +from .obs import Obs, dump_object, reweight, correlate from .fits import least_squares from .linalg import eigh, inv, cholesky from .roots import find_root diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 06d8abf4..82a83f3a 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -11,7 +11,7 @@ from scipy.odr import ODR, Model, RealData import iminuit from autograd import jacobian from autograd import elementwise_grad as egrad -from .pyerrors import Obs, derived_observable, covariance, pseudo_Obs +from .obs import Obs, derived_observable, covariance, pseudo_Obs class Fit_result(Sequence): diff --git a/pyerrors/input/bdio.py b/pyerrors/input/bdio.py index 4c4eec03..3e29924b 100644 --- a/pyerrors/input/bdio.py +++ b/pyerrors/input/bdio.py @@ -4,7 +4,7 @@ import ctypes import hashlib import autograd.numpy as np # Thinly-wrapped numpy -from ..pyerrors import Obs +from ..obs import Obs def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs): diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py index ac1bb78a..c184ad8f 100644 --- a/pyerrors/input/hadrons.py +++ b/pyerrors/input/hadrons.py @@ -4,7 +4,7 @@ import os import h5py import numpy as np -from ..pyerrors import Obs, CObs +from ..obs import Obs, CObs from ..correlators import Corr from ..npr import Npr_matrix diff --git a/pyerrors/input/misc.py b/pyerrors/input/misc.py index 0cb10f74..b1328241 100644 --- a/pyerrors/input/misc.py +++ b/pyerrors/input/misc.py @@ -6,7 +6,7 @@ import fnmatch import re import struct import numpy as np # Thinly-wrapped numpy -from ..pyerrors import Obs +from ..obs import Obs def read_pbp(path, prefix, **kwargs): diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 9aaf30be..1710a25e 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -6,7 +6,7 @@ import fnmatch import re import struct import numpy as np # Thinly-wrapped numpy -from ..pyerrors import Obs +from ..obs import Obs from ..fits import fit_lin diff --git a/pyerrors/input/sfcf.py b/pyerrors/input/sfcf.py index c6309cfe..3f079b79 100644 --- a/pyerrors/input/sfcf.py +++ b/pyerrors/input/sfcf.py @@ -5,7 +5,7 @@ import os import fnmatch import re import numpy as np # Thinly-wrapped numpy -from ..pyerrors import Obs +from ..obs import Obs def read_sfcf(path, prefix, name, **kwargs): diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index 0f08d72b..6f6d757a 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -4,7 +4,7 @@ import numpy as np from autograd import jacobian import autograd.numpy as anp # Thinly-wrapped numpy -from .pyerrors import derived_observable, CObs, Obs +from .obs import derived_observable, CObs, Obs from functools import partial from autograd.extend import defvjp diff --git a/pyerrors/misc.py b/pyerrors/misc.py index 920518e5..04f29dc2 100644 --- a/pyerrors/misc.py +++ b/pyerrors/misc.py @@ -2,7 +2,7 @@ # coding: utf-8 import numpy as np -from .pyerrors import Obs +from .obs import Obs def gen_correlated_data(means, cov, name, tau=0.5, samples=1000): diff --git a/pyerrors/mpm.py b/pyerrors/mpm.py index c835b1e2..00ace703 100644 --- a/pyerrors/mpm.py +++ b/pyerrors/mpm.py @@ -3,7 +3,7 @@ import numpy as np import scipy.linalg -from .pyerrors import Obs +from .obs import Obs from .linalg import svd, eig, pinv @@ -105,8 +105,6 @@ def matrix_pencil_method_old(data, p, noise_level=None, verbose=1, **kwargs): # Moore–Penrose pseudoinverse pinv_y1 = pinv(y1) - # Note: Automatic differentiation of eig is implemented in the git of autograd - # but not yet released to PyPi (1.3). The code is currently part of pyerrors e = eig((pinv_y1 @ y2), **kwargs) energy_levels = -np.log(np.abs(e)) return sorted(energy_levels, key=lambda x: abs(x.value)) diff --git a/pyerrors/pyerrors.py b/pyerrors/obs.py similarity index 100% rename from pyerrors/pyerrors.py rename to pyerrors/obs.py diff --git a/pyerrors/roots.py b/pyerrors/roots.py index 0c7c1566..792b2e63 100644 --- a/pyerrors/roots.py +++ b/pyerrors/roots.py @@ -3,7 +3,7 @@ import scipy.optimize from autograd import jacobian -from .pyerrors import derived_observable, pseudo_Obs +from .obs import derived_observable, pseudo_Obs def find_root(d, func, guess=1.0, **kwargs): diff --git a/setup.py b/setup.py index 20e4803b..d632fdcb 100644 --- a/setup.py +++ b/setup.py @@ -9,5 +9,5 @@ setup(name='pyerrors', author_email='fabian.joswig@ed.ac.uk', packages=find_packages(), python_requires='>=3.6.0', - install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit<2'] + install_requires=['numpy>=1.16', 'autograd>=1.2', 'numdifftools', 'matplotlib>=3.3', 'scipy', 'iminuit<2', 'h5py'] ) From 70240c71aa7b39735245b62a05e8db6b56d2c4c1 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 11:30:41 +0000 Subject: [PATCH 02/12] CONTRIBUTING updated --- CONTRIBUTING.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 299e5891..a6705807 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,6 +1,6 @@ # Development ### Setup -If you want to contribute to `pyerrors` please fork `pyerrors` on Github, clone the current `develop` branch +If you want to contribute to `pyerrors` please [fork](https://docs.github.com/en/get-started/quickstart/fork-a-repo) `pyerrors` on Github, clone the current `develop` branch ``` git clone http://github.com/my_username/pyerrors.git --branch develop ``` From bb7935b870448d0bebf50ff4abccad7e07d694e6 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 17:12:48 +0000 Subject: [PATCH 03/12] Documentation extended --- pyerrors/__init__.py | 59 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index db7b6859..dc02fd04 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -1,3 +1,62 @@ +r''' +# What is pyerrors? +`pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. + +## Getting started + +```python +import numpy as np +import pyerrors as pe + +my_obs = pe.Obs([samples], ['ensemble_name']) +my_new_obs = 2 * np.log(my_obs) / my_obs +my_new_obs.gamma_method() +my_new_obs.details() +print(my_new_obs) +``` +# The `Obs` class +`pyerrors.obs.Obs` +```python +import pyerrors as pe + +my_obs = pe.Obs([samples], ['ensemble_name']) +``` + +## Multiple ensembles/replica + +## Irregular Monte Carlo chains + +# Error propagation +Automatic differentiation, cite Alberto, + +numpy overloaded +```python +import numpy as np +import pyerrors as pe + +my_obs = pe.Obs([samples], ['ensemble_name']) +my_new_obs = 2 * np.log(my_obs) / my_obs +my_new_obs.gamma_method() +my_new_obs.details() +``` + +# Error estimation +`pyerrors.obs.Obs.gamma_method` + +$\delta_i\delta_j$ + +## Exponential tails + +## Covariance + +# Optimization / fits / roots + +# Complex observables + +# Matrix operations + +# Input +''' from .obs import * from .correlators import * from .fits import * From 2c004d303687ccc2aefe74b3436547b1c32c78e0 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 20:50:55 +0000 Subject: [PATCH 04/12] docs workflow added --- .github/workflows/docs.yml | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 .github/workflows/docs.yml diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 00000000..7bc2bc13 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,34 @@ +name: docs + +on: + push: + branches: + - feature/docs + +jobs: + docs: + runs-on: ubuntu-latest + + steps: + - name: Set up Python environment + uses: actions/setup-python@v2 + with: + python-version: "3.8" + - uses: actions/checkout@v2 + - name: Updated documentation + run: | + git config --global user.email "${{ github.actor }}@users.noreply.github.com" + git config --global user.name "${{ github.actor }}" + git fetch origin documentation + git checkout documentation + git pull + git merge --allow-unrelated-histories -X theirs develop + python -m pip install --upgrade pip + pip install wheel + pip install . + pip install pdoc + echo $(ls -l docs) + pdoc --docformat numpy --math -o ./docs ./pyerrors + echo $(ls -l docs) + git add docs + if [ -n "$(git diff --cached --exit-code)" ]; then git commit -am "Documentation updated"; git push; fi From fa6b5029db8786493a54c1989801c5db9dcfa6f9 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 20:52:17 +0000 Subject: [PATCH 05/12] docs workflow corrected --- .github/workflows/docs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 7bc2bc13..d4782abe 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -3,7 +3,7 @@ name: docs on: push: branches: - - feature/docs + - develop jobs: docs: From c41ab8d299561b729f2307e0ab54fc8636f4d18b Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 21:03:08 +0000 Subject: [PATCH 06/12] README updated --- README.md | 42 +++++++----------------------------------- pyerrors/__init__.py | 6 ++++++ 2 files changed, 13 insertions(+), 35 deletions(-) diff --git a/README.md b/README.md index 6e53ea35..abf327a3 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,10 @@ [![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.6+-blue.svg)](https://www.python.org/downloads/) # pyerrors `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. -It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017). Some of its features are: -* **automatic differentiation** as suggested in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289) (partly based on the [autograd](https://github.com/HIPS/autograd) package) -* **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228) -* coherent **error propagation** for data from **different Markov chains** -* **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289] -* **real and complex matrix operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) -There exist similar 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) -- [Python 3](https://github.com/mbruno46/pyobs) +- **Documentation:** https://fjosw.github.io/pyerrors/pyerrors.html +- **Examples**: https://github.com/fjosw/pyerrors/tree/develop/examples +- **Contributing:** https://github.com/fjosw/pyerrors/blob/develop/CONTRIBUTING.md ## Installation To install the most recent release of `pyerrors` run @@ -23,31 +16,10 @@ to install the current `develop` version run pip install git+https://github.com/fjosw/pyerrors.git@develop ``` -## Usage -The basic objects of a pyerrors analysis are instances of the class `Obs`. They can be initialized with an array of Monte Carlo data (e.g. `samples1`) and a name for the given ensemble (e.g. `'ensemble1'`). The `gamma_method` can then be used to compute the statistical error, taking into account autocorrelations. The `print` method outputs a human readable result. -```python -import pyerrors as pe - -obs1 = pe.Obs([samples1], ['ensemble1']) -obs1.gamma_method() -obs1.print() -``` -Often one is interested in secondary observables which can be arbitrary functions of primary observables. `pyerrors` overloads most basic math operations and `numpy` functions such that the user can work with `Obs` objects as if they were floats -```python -import numpy as np - -obs3 = 12.0 / obs1 ** 2 - np.exp(-1.0 / obs2) -obs3.gamma_method() -obs3.print() -``` - -More detailed examples can be found in the `examples` folder: - -* [01_basic_example](examples/01_basic_example.ipynb) -* [02_correlators](examples/02_correlators.ipynb) -* [03_pcac_example](examples/03_pcac_example.ipynb) -* [04_fit_example](examples/04_fit_example.ipynb) -* [05_matrix_operations](examples/05_matrix_operations.ipynb) +There exist similar 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) +- [Python 3](https://github.com/mbruno46/pyobs) ## License [MIT](https://choosealicense.com/licenses/mit/) diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index dc02fd04..b8f3d149 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -1,6 +1,12 @@ r''' # What is pyerrors? `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. +It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017). Some of its features are: +- **automatic differentiation** as suggested in [arXiv:1809.01289](https://arxiv.org/abs/1809.01289) (partly based on the [autograd](https://github.com/HIPS/autograd) package) +- **treatment of slow modes** in the simulation as suggested in [arXiv:1009.5228](https://arxiv.org/abs/1009.5228) +- coherent **error propagation** for data from **different Markov chains** +- **non-linear fits with x- and y-errors** and exact linear error propagation based on automatic differentiation as introduced in [arXiv:1809.01289] +- **real and complex matrix operations** and their error propagation based on automatic differentiation (cholesky decomposition, calculation of eigenvalues and eigenvectors, singular value decomposition...) ## Getting started From 11d2f788cb43c8e06cebda78a17d416c109a7d51 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 21:08:40 +0000 Subject: [PATCH 07/12] Links added to documentation --- pyerrors/__init__.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index b8f3d149..0c39e998 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -55,13 +55,22 @@ $\delta_i\delta_j$ ## Covariance +# Correlators +`pyerrors.correlators.Corr` + # Optimization / fits / roots +`pyerrors.fits` +`pyerrors.roots` + # Complex observables +`pyerrors.obs.CObs` # Matrix operations +`pyerrors.linalg` # Input +`pyerrors.input` ''' from .obs import * from .correlators import * From 9d0135c5e48fd0c7150d613358dc300a62bfa7e4 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 21:17:06 +0000 Subject: [PATCH 08/12] README updated --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index abf327a3..4974c8ec 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,11 @@ -[![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.6+-blue.svg)](https://www.python.org/downloads/) +[![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) [![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.6+-blue.svg)](https://www.python.org/downloads/) # pyerrors `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. - **Documentation:** https://fjosw.github.io/pyerrors/pyerrors.html - **Examples**: https://github.com/fjosw/pyerrors/tree/develop/examples - **Contributing:** https://github.com/fjosw/pyerrors/blob/develop/CONTRIBUTING.md +- **Bug reports:** https://github.com/fjosw/pyerrors/issues ## Installation To install the most recent release of `pyerrors` run @@ -16,10 +17,11 @@ to install the current `develop` version run pip install git+https://github.com/fjosw/pyerrors.git@develop ``` +## Other implementations There exist similar 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) -- [Python 3](https://github.com/mbruno46/pyobs) +- [Python](https://github.com/mbruno46/pyobs) ## License [MIT](https://choosealicense.com/licenses/mit/) From a23a97aed16777744f8a7e8dbc406715b72325dc Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 21:18:36 +0000 Subject: [PATCH 09/12] License badge added --- README.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/README.md b/README.md index 4974c8ec..1b34e861 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) [![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.6+-blue.svg)](https://www.python.org/downloads/) +[![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) [![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.6+-blue.svg)](https://www.python.org/downloads/) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) # pyerrors `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. @@ -22,6 +22,3 @@ There exist similar 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) - [Python](https://github.com/mbruno46/pyobs) - -## License -[MIT](https://choosealicense.com/licenses/mit/) From effccb1cc8ca89698f9030e32d064153f95311ed Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sun, 7 Nov 2021 21:44:22 +0000 Subject: [PATCH 10/12] docstrings updated --- pyerrors/correlators.py | 27 +++--- pyerrors/fits.py | 52 ++++++------ pyerrors/input/openQCD.py | 51 +++++++----- pyerrors/input/sfcf.py | 15 ++-- pyerrors/linalg.py | 3 - pyerrors/misc.py | 4 +- pyerrors/npr.py | 9 +- pyerrors/obs.py | 171 ++++++++++++++++++++------------------ 8 files changed, 173 insertions(+), 159 deletions(-) diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py index 8459f1f5..6a2108f4 100644 --- a/pyerrors/correlators.py +++ b/pyerrors/correlators.py @@ -226,8 +226,8 @@ class Corr: def roll(self, dt): """Periodically shift the correlator by dt timeslices - Attributes: - ----------- + Parameters + ---------- dt : int number of timeslices """ @@ -264,9 +264,6 @@ class Corr: weight : Obs Reweighting factor. An Observable that has to be defined on a superset of the configurations in obs[i].idl for all i. - - Keyword arguments - ----------------- 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 @@ -283,8 +280,8 @@ class Corr: def T_symmetry(self, partner, parity=+1): """Return the time symmetry average of the correlator and its partner - Attributes: - ----------- + Parameters + ---------- partner : Corr Time symmetry partner of the Corr partity : int @@ -309,8 +306,8 @@ class Corr: def deriv(self, symmetric=True): """Return the first derivative of the correlator with respect to x0. - Attributes: - ----------- + Parameters + ---------- symmetric : bool decides whether symmertic of simple finite differences are used. Default: True """ @@ -414,8 +411,8 @@ class Corr: def fit(self, function, fitrange=None, silent=False, **kwargs): """Fits function to the data - Attributes: - ----------- + Parameters + ---------- function : obj function to fit to the data. See fits.least_squares for details. fitrange : list @@ -447,8 +444,8 @@ class Corr: def plateau(self, plateau_range=None, method="fit"): """ Extract a plateu value from a Corr object - Attributes: - ----------- + Parameters + ---------- plateau_range : list list with two entries, indicating the first and the last timeslice of the plateau region. @@ -578,8 +575,8 @@ class Corr: def dump(self, filename): """Dumps the Corr into a pickel file - Attributes: - ----------- + Parameters + ---------- filename : str Name of the file """ diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 82a83f3a..44cf4acf 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -59,7 +59,7 @@ class Fit_result(Sequence): def least_squares(x, y, func, priors=None, silent=False, **kwargs): """Performs a non-linear fit to y = func(x). - Arguments: + Parameters ---------- x : list list of floats. @@ -87,22 +87,23 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): enough. silent : bool, optional If true all output to the console is omitted (default False). - - - Keyword arguments - ----------------- - initial_guess -- can provide an initial guess for the input parameters. Relevant for + initial_guess : list + can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters. - method -- 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 - migrad of iminuit. If no method is specified, Levenberg-Marquard is used. - Reliable alternatives are migrad, Powell and Nelder-Mead. - resplot -- If true, a plot which displays fit, data and residuals is generated (default False). - qqplot -- If true, a quantile-quantile plot of the fit result is generated (default False). - expected_chisquare -- If true prints the expected chisquare which is - corrected by effects caused by correlated input data. - This can take a while as the full correlation matrix - has to be calculated (default False). + method : str + 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 + migrad of iminuit. If no method is specified, Levenberg-Marquard is used. + Reliable alternatives are migrad, Powell and Nelder-Mead. + resplot : bool + If true, a plot which displays fit, data and residuals is generated (default False). + qqplot : bool + If true, a quantile-quantile plot of the fit result is generated (default False). + expected_chisquare : bool + If true prints the expected chisquare which is + corrected by effects caused by correlated input data. + This can take a while as the full correlation matrix + has to be calculated (default False). """ if priors is not None: return _prior_fit(x, y, func, priors, silent=silent, **kwargs) @@ -254,6 +255,8 @@ def odr_fit(x, y, func, silent=False, **kwargs): def total_least_squares(x, y, func, silent=False, **kwargs): """Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. + Parameters + ---------- x : list list of Obs, or a tuple of lists of Obs y : list @@ -276,15 +279,14 @@ def total_least_squares(x, y, func, silent=False, **kwargs): silent : bool, optional If true all output to the console is omitted (default False). Based on the orthogonal distance regression module of scipy - - Keyword arguments - ----------------- - initial_guess -- can provide an initial guess for the input parameters. Relevant for non-linear - fits with many parameters. - expected_chisquare -- If true prints the expected chisquare which is - corrected by effects caused by correlated input data. - This can take a while as the full correlation matrix - has to be calculated (default False). + initial_guess : list + can provide an initial guess for the input parameters. Relevant for non-linear + fits with many parameters. + expected_chisquare : bool + If true prints the expected chisquare which is + corrected by effects caused by correlated input data. + This can take a while as the full correlation matrix + has to be calculated (default False). """ output = Fit_result() diff --git a/pyerrors/input/openQCD.py b/pyerrors/input/openQCD.py index 1710a25e..5e1c8d49 100644 --- a/pyerrors/input/openQCD.py +++ b/pyerrors/input/openQCD.py @@ -13,15 +13,16 @@ from ..fits import fit_lin def read_rwms(path, prefix, version='2.0', names=None, **kwargs): """Read rwms format from given folder structure. Returns a list of length nrw - Attributes - ----------------- - version -- version of openQCD, default 2.0 - - Keyword arguments - ----------------- - r_start -- list which contains the first config to be read for each replicum - r_stop -- list which contains the last config to be read for each replicum - postfix -- postfix of the file to read, e.g. '.ms1' for openQCD-files + Parameters + ---------- + version : str + version of openQCD, default 2.0 + r_start : list + list which contains the first config to be read for each replicum + r_stop : list + list which contains the last config to be read for each replicum + postfix : str + postfix of the file to read, e.g. '.ms1' for openQCD-files """ known_oqcd_versions = ['1.4', '1.6', '2.0'] if not (version in known_oqcd_versions): @@ -174,19 +175,25 @@ def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwar Parameters ---------- - path -- Path to .ms.dat files - prefix -- Ensemble prefix - dtr_read -- Determines how many trajectories should be skipped when reading the ms.dat files. - Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. - xmin -- First timeslice where the boundary effects have sufficiently decayed. - spatial_extent -- spatial extent of the lattice, required for normalization. - fit_range -- Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5) - - Keyword arguments - ----------------- - r_start -- list which contains the first config to be read for each replicum. - r_stop -- list which contains the last config to be read for each replicum. - plaquette -- If true extract the plaquette estimate of t0 instead. + path : str + Path to .ms.dat files + prefix : str + Ensemble prefix + dtr_read : int + Determines how many trajectories should be skipped when reading the ms.dat files. + Corresponds to dtr_cnfg / dtr_ms in the openQCD input file. + xmin : int + First timeslice where the boundary effects have sufficiently decayed. + spatial_extent : int + spatial extent of the lattice, required for normalization. + fit_range : int + Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5) + r_start : list + list which contains the first config to be read for each replicum. + r_stop: list + list which contains the last config to be read for each replicum. + plaquette : bool + If true extract the plaquette estimate of t0 instead. """ ls = [] diff --git a/pyerrors/input/sfcf.py b/pyerrors/input/sfcf.py index 3f079b79..e48bdd16 100644 --- a/pyerrors/input/sfcf.py +++ b/pyerrors/input/sfcf.py @@ -11,8 +11,8 @@ from ..obs import Obs def read_sfcf(path, prefix, name, **kwargs): """Read sfcf C format from given folder structure. - Keyword arguments - ----------------- + Parameters + ---------- im -- if True, read imaginary instead of real part of the correlation function. single -- if True, read a boundary-to-boundary correlation function with a single value b2b -- if True, read a time-dependent boundary-to-boundary correlation function @@ -113,15 +113,12 @@ def read_sfcf(path, prefix, name, **kwargs): def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwargs): """Read sfcf c format from given folder structure. - Arguments - ----------------- + Parameters + ---------- quarks -- Label of the quarks used in the sfcf input file noffset -- Offset of the source (only relevant when wavefunctions are used) wf -- ID of wave function wf2 -- ID of the second wavefunction (only relevant for boundary-to-boundary correlation functions) - - Keyword arguments - ----------------- im -- if True, read imaginary instead of real part of the correlation function. b2b -- if True, read a time-dependent boundary-to-boundary correlation function names -- Alternative labeling for replicas/ensembles. Has to have the appropriate length @@ -236,8 +233,8 @@ def read_sfcf_c(path, prefix, name, quarks='.*', noffset=0, wf=0, wf2=0, **kwarg def read_qtop(path, prefix, **kwargs): """Read qtop format from given folder structure. - Keyword arguments - ----------------- + Parameters + ---------- target -- specifies the topological sector to be reweighted to (default 0) full -- if true read the charge instead of the reweighting factor. """ diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index 6f6d757a..bc265efa 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -20,9 +20,6 @@ def derived_array(func, data, **kwargs): automatic differentiation to work, all numpy functions have to have the autograd wrapper (use 'import autograd.numpy as anp'). data -- list of Obs, e.g. [obs1, obs2, obs3]. - - Keyword arguments - ----------------- man_grad -- manually supply a list or an array which contains the jacobian of func. Use cautiously, supplying the wrong derivative will not be intercepted. diff --git a/pyerrors/misc.py b/pyerrors/misc.py index 04f29dc2..7fd0af58 100644 --- a/pyerrors/misc.py +++ b/pyerrors/misc.py @@ -8,8 +8,8 @@ from .obs import Obs def gen_correlated_data(means, cov, name, tau=0.5, samples=1000): """ Generate observables with given covariance and autocorrelation times. - Arguments - ----------------- + Parameters + ---------- means -- list containing the mean value of each observable. cov -- covariance matrix for the data to be geneated. name -- ensemble name for the data to be geneated. diff --git a/pyerrors/npr.py b/pyerrors/npr.py index b1c32c8a..013e0d25 100644 --- a/pyerrors/npr.py +++ b/pyerrors/npr.py @@ -73,9 +73,12 @@ def inv_propagator(prop): def Zq(inv_prop, fermion='Wilson'): """ Calculates the quark field renormalization constant Zq - Attributes: - inv_prop -- Inverted 12x12 quark propagator - fermion -- Fermion type for which the tree-level propagator is used + Parameters + ---------- + inv_prop : array + Inverted 12x12 quark propagator + fermion : str + Fermion type for which the tree-level propagator is used in the calculation of Zq. Default Wilson. """ _check_geometry() diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 225424f7..794e7287 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -53,7 +53,7 @@ class Obs: def __init__(self, samples, names, idl=None, means=None, **kwargs): """ Initialize Obs object. - Attributes + Parameters ---------- samples : list list of numpy arrays containing the Monte Carlo samples @@ -150,57 +150,11 @@ class Obs: res[e_name].append(e_name) return res - def expand_deltas(self, deltas, idx, shape): - """Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0. - If idx is of type range, the deltas are not changed - - Parameters - ---------- - deltas -- List of fluctuations - idx -- List or range of configs on which the deltas are defined. - shape -- Number of configs in idx. - """ - if type(idx) is range: - return deltas - else: - ret = np.zeros(idx[-1] - idx[0] + 1) - for i in range(shape): - ret[idx[i] - idx[0]] = deltas[i] - return ret - - def calc_gamma(self, deltas, idx, shape, w_max, fft): - """Calculate Gamma_{AA} from the deltas, which are defined on idx. - idx is assumed to be a contiguous range (possibly with a stepsize != 1) - - Parameters - ---------- - deltas -- List of fluctuations - idx -- List or range of configs on which the deltas are defined. - shape -- Number of configs in idx. - w_max -- Upper bound for the summation window - fft -- boolean, which determines whether the fft algorithm is used for - the computation of the autocorrelation function - """ - gamma = np.zeros(w_max) - deltas = self.expand_deltas(deltas, idx, shape) - new_shape = len(deltas) - if fft: - max_gamma = min(new_shape, w_max) - # The padding for the fft has to be even - padding = new_shape + max_gamma + (new_shape + max_gamma) % 2 - gamma[:max_gamma] += np.fft.irfft(np.abs(np.fft.rfft(deltas, padding)) ** 2)[:max_gamma] - else: - for n in range(w_max): - if new_shape - n >= 0: - gamma[n] += deltas[0:new_shape - n].dot(deltas[n:new_shape]) - - return gamma - def gamma_method(self, **kwargs): """Calculate the error and related properties of the Obs. - Keyword arguments - ----------------- + Parameters + ---------- S : float specifies a custom value for the parameter S (default 2.0), can be a float or an array of floats for different ensembles @@ -378,6 +332,52 @@ class Obs: self.ddvalue = np.sqrt(self.ddvalue) / self.dvalue return + def expand_deltas(self, deltas, idx, shape): + """Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0. + If idx is of type range, the deltas are not changed + + Parameters + ---------- + deltas -- List of fluctuations + idx -- List or range of configs on which the deltas are defined. + shape -- Number of configs in idx. + """ + if type(idx) is range: + return deltas + else: + ret = np.zeros(idx[-1] - idx[0] + 1) + for i in range(shape): + ret[idx[i] - idx[0]] = deltas[i] + return ret + + def calc_gamma(self, deltas, idx, shape, w_max, fft): + """Calculate Gamma_{AA} from the deltas, which are defined on idx. + idx is assumed to be a contiguous range (possibly with a stepsize != 1) + + Parameters + ---------- + deltas -- List of fluctuations + idx -- List or range of configs on which the deltas are defined. + shape -- Number of configs in idx. + w_max -- Upper bound for the summation window + fft -- boolean, which determines whether the fft algorithm is used for + the computation of the autocorrelation function + """ + gamma = np.zeros(w_max) + deltas = self.expand_deltas(deltas, idx, shape) + new_shape = len(deltas) + if fft: + max_gamma = min(new_shape, w_max) + # The padding for the fft has to be even + padding = new_shape + max_gamma + (new_shape + max_gamma) % 2 + gamma[:max_gamma] += np.fft.irfft(np.abs(np.fft.rfft(deltas, padding)) ** 2)[:max_gamma] + else: + for n in range(w_max): + if new_shape - n >= 0: + gamma[n] += deltas[0:new_shape - n].dot(deltas[n:new_shape]) + + return gamma + def print(self, level=1): warnings.warn("Method 'print' renamed to 'details'", DeprecationWarning) self.details(level > 1) @@ -539,9 +539,10 @@ class Obs: def dump(self, name, **kwargs): """Dump the Obs to a pickle file 'name'. - Keyword arguments - ----------------- - path -- specifies a custom path for the file (default '.') + Parameters + ---------- + path : str + specifies a custom path for the file (default '.') """ if 'path' in kwargs: file_name = kwargs.get('path') + '/' + name + '.p' @@ -836,7 +837,8 @@ def merge_idx(idl): Parameters ---------- - idl -- List of lists or ranges. + idl : list + List of lists or ranges. """ # Use groupby to efficiently check whether all elements of idl are identical @@ -893,14 +895,15 @@ def filter_zeroes(names, deltas, idl, eps=Obs.filter_eps): Parameters ---------- - names -- List of names - deltas -- Dict lists of fluctuations - idx -- Dict of lists or ranges of configs on which the deltas are defined. - Has to be a subset of new_idx. - - Optional parameters - ---------- - eps -- Prefactor that enters the filter criterion. + names : list + List of names + deltas : dict + Dict lists of fluctuations + idx : dict + Dict of lists or ranges of configs on which the deltas are defined. + Has to be a subset of new_idx. + eps : float + Prefactor that enters the filter criterion. """ new_names = [] new_deltas = {} @@ -931,9 +934,6 @@ def derived_observable(func, data, **kwargs): the autograd wrapper (use 'import autograd.numpy as anp'). data : list list of Obs, e.g. [obs1, obs2, obs3]. - - Keyword arguments - ----------------- num_grad : bool if True, numerical derivatives are used instead of autograd (default False). To control the numerical differentiation the @@ -1072,10 +1072,13 @@ def reduce_deltas(deltas, idx_old, idx_new): Parameters ---------- - deltas -- List of fluctuations - idx_old -- List or range of configs on which the deltas are defined - idx_new -- List of configs for which we want to extract the deltas. - Has to be a subset of idx_old. + deltas : list + List of fluctuations + idx_old : list + List or range of configs on which the deltas are defined + idx_new : list + List of configs for which we want to extract the deltas. + Has to be a subset of idx_old. """ if not len(deltas) == len(idx_old): raise Exception('Lenght of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old))) @@ -1109,9 +1112,6 @@ def reweight(weight, obs, **kwargs): configurations in obs[i].idl for all i. obs : list list of Obs, e.g. [obs1, obs2, obs3]. - - Keyword arguments - ----------------- 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 @@ -1146,8 +1146,8 @@ def reweight(weight, obs, **kwargs): def correlate(obs_a, obs_b): """Correlate two observables. - Attributes: - ----------- + Parameters + ---------- obs_a : Obs First observable obs_b : Obs @@ -1193,10 +1193,11 @@ def covariance(obs1, obs2, correlation=False, **kwargs): is constrained to the maximum value in order to make sure that covariance matrices are positive semidefinite. - Keyword arguments - ----------------- - correlation -- if true the correlation instead of the covariance is - returned (default False) + Parameters + ---------- + correlation : bool + if true the correlation instead of the covariance is + returned (default False) """ for name in sorted(set(obs1.names + obs2.names)): @@ -1450,9 +1451,14 @@ def pseudo_Obs(value, dvalue, name, samples=1000): def dump_object(obj, name, **kwargs): """Dump object into pickle file. - Keyword arguments - ----------------- - path -- specifies a custom path for the file (default '.') + Parameters + ---------- + obj : object + object to be saved in the pickle file + name : str + name of the file + path : str + specifies a custom path for the file (default '.') """ if 'path' in kwargs: file_name = kwargs.get('path') + '/' + name + '.p' @@ -1471,6 +1477,11 @@ def load_object(path): def merge_obs(list_of_obs): """Combine all observables in list_of_obs into one new observable + Parameters + ---------- + list_of_obs : list + list of the Obs object to be combined + It is not possible to combine obs which are based on the same replicum """ replist = [item for obs in list_of_obs for item in obs.names] From 2a594dfb4b34203f7b80880edd6d41975e2860ac Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 8 Nov 2021 09:50:51 +0000 Subject: [PATCH 11/12] total_least_squares docstring updated --- pyerrors/fits.py | 293 ++++++++++++++++++++++++----------------------- 1 file changed, 149 insertions(+), 144 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 44cf4acf..4787841c 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -111,149 +111,8 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): return _standard_fit(x, y, func, silent=silent, **kwargs) -def standard_fit(x, y, func, silent=False, **kwargs): - warnings.warn("standard_fit renamed to least_squares", DeprecationWarning) - return least_squares(x, y, func, silent=silent, **kwargs) - - -def _standard_fit(x, y, func, silent=False, **kwargs): - - output = Fit_result() - - output.fit_function = func - - x = np.asarray(x) - - if x.shape[-1] != len(y): - raise Exception('x and y input have to have the same length') - - if len(x.shape) > 2: - raise Exception('Unkown format for x values') - - if not callable(func): - raise TypeError('func has to be a function.') - - for i in range(25): - try: - func(np.arange(i), x.T[0]) - except: - pass - else: - break - - n_parms = i - - if not silent: - print('Fit with', n_parms, 'parameters') - - y_f = [o.value for o in y] - dy_f = [o.dvalue for o in y] - - if np.any(np.asarray(dy_f) <= 0.0): - raise Exception('No y errors available, run the gamma method first.') - - if 'initial_guess' in kwargs: - x0 = kwargs.get('initial_guess') - if len(x0) != n_parms: - raise Exception('Initial guess does not have the correct length.') - else: - x0 = [0.1] * n_parms - - def chisqfunc(p): - model = func(p, x) - chisq = anp.sum(((y_f - model) / dy_f) ** 2) - return chisq - - if 'method' in kwargs: - output.method = kwargs.get('method') - if not silent: - print('Method:', kwargs.get('method')) - if kwargs.get('method') == 'migrad': - fit_result = iminuit.minimize(chisqfunc, x0) - fit_result = iminuit.minimize(chisqfunc, fit_result.x) - else: - fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method')) - fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=1e-12) - - chisquare = fit_result.fun - - output.iterations = fit_result.nit - else: - output.method = 'Levenberg-Marquardt' - if not silent: - print('Method: Levenberg-Marquardt') - - 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) - - chisquare = np.sum(fit_result.fun ** 2) - - output.iterations = fit_result.nfev - - if not fit_result.success: - raise Exception('The minimization procedure did not converge.') - - if x.shape[-1] - n_parms > 0: - output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms) - else: - output.chisquare_by_dof = float('nan') - - output.message = fit_result.message - if not silent: - print(fit_result.message) - print('chisquare/d.o.f.:', output.chisquare_by_dof) - - if kwargs.get('expected_chisquare') is True: - W = np.diag(1 / np.asarray(dy_f)) - cov = covariance_matrix(y) - A = W @ jacobian(func)(fit_result.x, x) - P_phi = A @ np.linalg.inv(A.T @ A) @ A.T - expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W) - output.chisquare_by_expected_chisquare = chisquare / expected_chisquare - if not silent: - print('chisquare/expected_chisquare:', - output.chisquare_by_expected_chisquare) - - hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(fit_result.x)) - - def chisqfunc_compact(d): - model = func(d[:n_parms], x) - chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2) - return chisq - - jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fit_result.x, y_f))) - - deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] - - result = [] - for i in range(n_parms): - result.append(derived_observable(lambda x, **kwargs: x[0], [pseudo_Obs(fit_result.x[i], 0.0, y[0].names[0], y[0].shape[y[0].names[0]])] + list(y), man_grad=[0] + list(deriv[i]))) - - output.fit_parameters = result - - output.chisquare = chisqfunc(fit_result.x) - output.dof = x.shape[-1] - n_parms - - if kwargs.get('resplot') is True: - residual_plot(x, y, func, result) - - if kwargs.get('qqplot') is True: - qqplot(x, y, func, result) - - return output - - -def odr_fit(x, y, func, silent=False, **kwargs): - warnings.warn("odr_fit renamed to total_least_squares", DeprecationWarning) - return total_least_squares(x, y, func, silent=silent, **kwargs) - - def total_least_squares(x, y, func, silent=False, **kwargs): - """Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. + r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. Parameters ---------- @@ -264,21 +123,24 @@ def total_least_squares(x, y, func, silent=False, **kwargs): func : object func has to be of the form + ```python def func(a, x): y = a[0] + a[1] * x + a[2] * anp.sinh(x) return y + ``` For multiple x values func can be of the form + ```python def func(a, x): (x1, x2) = x return a[0] * x1 ** 2 + a[1] * x2 + ``` It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work. silent : bool, optional If true all output to the console is omitted (default False). - Based on the orthogonal distance regression module of scipy initial_guess : list can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters. @@ -287,7 +149,9 @@ def total_least_squares(x, y, func, silent=False, **kwargs): corrected by effects caused by correlated input data. This can take a while as the full correlation matrix has to be calculated (default False). - """ + + Based on the orthogonal distance regression module of scipy + ''' output = Fit_result() @@ -541,6 +405,147 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): return output +def standard_fit(x, y, func, silent=False, **kwargs): + warnings.warn("standard_fit renamed to least_squares", DeprecationWarning) + return least_squares(x, y, func, silent=silent, **kwargs) + + +def _standard_fit(x, y, func, silent=False, **kwargs): + + output = Fit_result() + + output.fit_function = func + + x = np.asarray(x) + + if x.shape[-1] != len(y): + raise Exception('x and y input have to have the same length') + + if len(x.shape) > 2: + raise Exception('Unkown format for x values') + + if not callable(func): + raise TypeError('func has to be a function.') + + for i in range(25): + try: + func(np.arange(i), x.T[0]) + except: + pass + else: + break + + n_parms = i + + if not silent: + print('Fit with', n_parms, 'parameters') + + y_f = [o.value for o in y] + dy_f = [o.dvalue for o in y] + + if np.any(np.asarray(dy_f) <= 0.0): + raise Exception('No y errors available, run the gamma method first.') + + if 'initial_guess' in kwargs: + x0 = kwargs.get('initial_guess') + if len(x0) != n_parms: + raise Exception('Initial guess does not have the correct length.') + else: + x0 = [0.1] * n_parms + + def chisqfunc(p): + model = func(p, x) + chisq = anp.sum(((y_f - model) / dy_f) ** 2) + return chisq + + if 'method' in kwargs: + output.method = kwargs.get('method') + if not silent: + print('Method:', kwargs.get('method')) + if kwargs.get('method') == 'migrad': + fit_result = iminuit.minimize(chisqfunc, x0) + fit_result = iminuit.minimize(chisqfunc, fit_result.x) + else: + fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method')) + fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=1e-12) + + chisquare = fit_result.fun + + output.iterations = fit_result.nit + else: + output.method = 'Levenberg-Marquardt' + if not silent: + print('Method: Levenberg-Marquardt') + + 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) + + chisquare = np.sum(fit_result.fun ** 2) + + output.iterations = fit_result.nfev + + if not fit_result.success: + raise Exception('The minimization procedure did not converge.') + + if x.shape[-1] - n_parms > 0: + output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms) + else: + output.chisquare_by_dof = float('nan') + + output.message = fit_result.message + if not silent: + print(fit_result.message) + print('chisquare/d.o.f.:', output.chisquare_by_dof) + + if kwargs.get('expected_chisquare') is True: + W = np.diag(1 / np.asarray(dy_f)) + cov = covariance_matrix(y) + A = W @ jacobian(func)(fit_result.x, x) + P_phi = A @ np.linalg.inv(A.T @ A) @ A.T + expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W) + output.chisquare_by_expected_chisquare = chisquare / expected_chisquare + if not silent: + print('chisquare/expected_chisquare:', + output.chisquare_by_expected_chisquare) + + hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(fit_result.x)) + + def chisqfunc_compact(d): + model = func(d[:n_parms], x) + chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2) + return chisq + + jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fit_result.x, y_f))) + + deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] + + result = [] + for i in range(n_parms): + result.append(derived_observable(lambda x, **kwargs: x[0], [pseudo_Obs(fit_result.x[i], 0.0, y[0].names[0], y[0].shape[y[0].names[0]])] + list(y), man_grad=[0] + list(deriv[i]))) + + output.fit_parameters = result + + output.chisquare = chisqfunc(fit_result.x) + output.dof = x.shape[-1] - n_parms + + if kwargs.get('resplot') is True: + residual_plot(x, y, func, result) + + if kwargs.get('qqplot') is True: + qqplot(x, y, func, result) + + return output + + +def odr_fit(x, y, func, silent=False, **kwargs): + warnings.warn("odr_fit renamed to total_least_squares", DeprecationWarning) + return total_least_squares(x, y, func, silent=silent, **kwargs) + + def fit_lin(x, y, **kwargs): """Performs a linear fit to y = n + m * x and returns two Obs n, m. From bc819828bc193419fd7c552561e903dda03538f3 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 8 Nov 2021 10:01:26 +0000 Subject: [PATCH 12/12] least_squares docstring updated --- pyerrors/fits.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 4787841c..686ed81f 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -57,7 +57,7 @@ class Fit_result(Sequence): def least_squares(x, y, func, priors=None, silent=False, **kwargs): - """Performs a non-linear fit to y = func(x). + r'''Performs a non-linear fit to y = func(x). Parameters ---------- @@ -68,14 +68,19 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): func : object fit function, has to be of the form + ```python def func(a, x): - return a[0] + a[1] * x + a[2] * anp.sinh(x) + y = a[0] + a[1] * x + a[2] * anp.sinh(x) + return y + ``` For multiple x values func can be of the form + ```python def func(a, x): (x1, x2) = x return a[0] * x1 ** 2 + a[1] * x2 + ``` It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work @@ -104,7 +109,7 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): corrected by effects caused by correlated input data. This can take a while as the full correlation matrix has to be calculated (default False). - """ + ''' if priors is not None: return _prior_fit(x, y, func, priors, silent=silent, **kwargs) else: