Merge branch 'develop' of github.com:fjosw/pyerrors into develop

.github/workflows/docs.yml

@@ -0,0 +1,34 @@
+name: docs
+
+on:
+  push:
+    branches:
+      - develop
+
+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

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
 ```

README.md

@@ -1,17 +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/) [![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.
-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
+- **Documentation:** https://fjosw.github.io/pyerrors/pyerrors.html
-- [Fortran](https://gitlab.ift.uam-csic.es/alberto/aderrors)
+- **Examples**: https://github.com/fjosw/pyerrors/tree/develop/examples
-- [Julia](https://gitlab.ift.uam-csic.es/alberto/aderrors.jl)
+- **Contributing:** https://github.com/fjosw/pyerrors/blob/develop/CONTRIBUTING.md
-- [Python 3](https://github.com/mbruno46/pyobs)
+- **Bug reports:** https://github.com/fjosw/pyerrors/issues
 
 ## Installation
 To install the most recent release of `pyerrors` run
@@ -23,31 +17,8 @@ to install the current `develop` version run
 pip install git+https://github.com/fjosw/pyerrors.git@develop
 ```
 
-## Usage
+## Other implementations
-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.
+There exist similar implementations of gamma method error analysis suites in
-```python
+- [Fortran](https://gitlab.ift.uam-csic.es/alberto/aderrors)
-import pyerrors as pe
+- [Julia](https://gitlab.ift.uam-csic.es/alberto/aderrors.jl)
-
+- [Python](https://github.com/mbruno46/pyobs)
-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)
-
-## License
-[MIT](https://choosealicense.com/licenses/mit/)

pyerrors/__init__.py

@@ -1,4 +1,78 @@
-from .pyerrors import *
+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
+
+```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
+
+# 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 .fits import *
 from . import dirac

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
@@ -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
         """

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):
@@ -57,9 +57,9 @@ 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).
 
-    Arguments:
+    Parameters
     ----------
     x : list
         list of floats.
@@ -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
@@ -87,173 +92,35 @@ 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).
-
+    initial_guess : list
-
+        can provide an initial guess for the input parameters. Relevant for
-    Keyword arguments
-    -----------------
-    initial_guess -- 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.
+    method : str
-              The possible methods are the ones which can be used for scipy.optimize.minimize and
+        can be used to choose an alternative method for the minimization of chisquare.
-              migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
+        The possible methods are the ones which can be used for scipy.optimize.minimize and
-              Reliable alternatives are migrad, Powell and Nelder-Mead.
+        migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
-    resplot -- If true, a plot which displays fit, data and residuals is generated (default False).
+        Reliable alternatives are migrad, Powell and Nelder-Mead.
-    qqplot -- If true, a quantile-quantile plot of the fit result is generated (default False).
+    resplot : bool
-    expected_chisquare -- If true prints the expected chisquare which is
+        If true, a plot which displays fit, data and residuals is generated (default False).
-                          corrected by effects caused by correlated input data.
+    qqplot : bool
-                          This can take a while as the full correlation matrix
+        If true, a quantile-quantile plot of the fit result is generated (default False).
-                          has to be calculated (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)
     else:
         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
+    ----------
     x : list
         list of Obs, or a tuple of lists of Obs
     y : list
@@ -261,31 +128,35 @@ 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.
+    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).
 
-    Keyword arguments
+    Based on the orthogonal distance regression module of scipy
-    -----------------
+    '''
-    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).
-    """
 
     output = Fit_result()
 
@@ -539,6 +410,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.
 

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):

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
 

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):

pyerrors/input/openQCD.py

@@ -6,22 +6,23 @@ 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
 
 
 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
+    Parameters
-    -----------------
+    ----------
-    version -- version of openQCD, default 2.0
+    version : str
-
+        version of openQCD, default 2.0
-    Keyword arguments
+    r_start : list
-    -----------------
+        list which contains the first config to be read for each replicum
-    r_start -- list which contains the first config to be read for each replicum
+    r_stop : list
-    r_stop -- list which contains the last config to be read for each replicum
+        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
+    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
+    path : str
-    prefix -- Ensemble prefix
+        Path to .ms.dat files
-    dtr_read -- Determines how many trajectories should be skipped when reading the ms.dat files.
+    prefix : str
-                Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
+        Ensemble prefix
-    xmin -- First timeslice where the boundary effects have sufficiently decayed.
+    dtr_read : int
-    spatial_extent -- spatial extent of the lattice, required for normalization.
+        Determines how many trajectories should be skipped when reading the ms.dat files.
-    fit_range -- Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5)
+        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
-
+    xmin : int
-    Keyword arguments
+        First timeslice where the boundary effects have sufficiently decayed.
-    -----------------
+    spatial_extent : int
-    r_start -- list which contains the first config to be read for each replicum.
+        spatial extent of the lattice, required for normalization.
-    r_stop -- list which contains the last config to be read for each replicum.
+    fit_range : int
-    plaquette -- If true extract the plaquette estimate of t0 instead.
+        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 = []

pyerrors/input/sfcf.py

@@ -5,14 +5,14 @@ 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):
     """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.
     """

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
@@ -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.

pyerrors/misc.py

@@ -2,14 +2,14 @@
 # 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):
     """ 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.

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))

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:
+        Parameters
-        inv_prop -- Inverted 12x12 quark propagator
+        ----------
-        fermion -- Fermion type for which the tree-level propagator is used
+        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()

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
+        Parameters
-        -----------------
+        ----------
-        path -- specifies a custom path for the file (default '.')
+        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
+    names : list
-    deltas -- Dict lists of fluctuations
+        List of names
-    idx    -- Dict of lists or ranges of configs on which the deltas are defined.
+    deltas : dict
-               Has to be a subset of new_idx.
+        Dict lists of fluctuations
-
+    idx : dict
-    Optional parameters
+        Dict of lists or ranges of configs on which the deltas are defined.
-    ----------
+        Has to be a subset of new_idx.
-    eps    -- Prefactor that enters the filter criterion.
+    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
+    deltas : list
-    idx_old -- List or range of configs on which the deltas are defined
+        List of fluctuations
-    idx_new -- List of configs for which we want to extract the deltas.
+    idx_old : list
-               Has to be a subset of idx_old.
+        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
+    Parameters
-    -----------------
+    ----------
-    correlation -- if true the correlation instead of the covariance is
+    correlation : bool
-                   returned (default False)
+        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
+    Parameters
-    -----------------
+    ----------
-    path -- specifies a custom path for the file (default '.')
+    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]

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):

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']
      )