lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-03-15 14:50:25 +01:00
Code Issues Projects Releases Packages Wiki Activity

fix: flak8 & pytest

Browse source
This commit is contained in:
ppetrak 2022-12-19 14:03:45 +01:00
parent 500c5234cf
commit 3d6ec7b397
3 changed files with 54 additions and 59 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

BIN
examples/my_db.sqlite Normal file
View file

Binary file not shown.

1
pyerrors/__init__.py
View file

@ -457,7 +457,6 @@ from .obs import *
from .correlators import * from .correlators import *
from .fits import * from .fits import *
from .misc import * from .misc import *
from . import combined_fits
from . import dirac from . import dirac
from . import input from . import input
from . import linalg from . import linalg

112
pyerrors/fits.py
View file

@ -70,13 +70,12 @@ class Fit_result(Sequence):
def least_squares(x, y, func, priors=None, silent=False, **kwargs): def least_squares(x, y, func, priors=None, silent=False, **kwargs):
r'''Performs a non-linear fit to y = func(x). r'''Performs a non-linear fit to y = func(x).
``` ```
Parameters Parameters
---------- ----------
For an uncombined fit: For an uncombined fit:
x : list x : list
list of floats. list of floats.
y : list y : list
@ -100,12 +99,12 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
``` ```
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
will not work. will not work.
OR For a combined fit: OR For a combined fit:
Do not need to use ordered dictionaries: python version >= 3.7: Dictionary order is guaranteed to be insertion order. Do not need to use ordered dictionaries: python version >= 3.7: Dictionary order is guaranteed to be insertion order.
(https://docs.python.org/3/library/stdtypes.html#dict-views) Ensures that x, y and func values are mapped correctly. (https://docs.python.org/3/library/stdtypes.html#dict-views) Ensures that x, y and func values are mapped correctly.
x : dict x : dict
dict of lists. dict of lists.
y : dict y : dict
@ -117,16 +116,16 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
import autograd.numpy as anp import autograd.numpy as anp
funcs = {"a": func_a, funcs = {"a": func_a,
"b": func_b} "b": func_b}
def func_a(a, x): def func_a(a, x):
return a[1] * anp.exp(-a[0] * x) return a[1] * anp.exp(-a[0] * x)
def func_b(a, x): def func_b(a, x):
return a[2] * anp.exp(-a[0] * x) return a[2] * anp.exp(-a[0] * x)
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
will not work. will not work.
priors : list, optional priors : list, optional
priors has to be a list with an entry for every parameter in the fit. The entries can either be priors has to be a list with an entry for every parameter in the fit. The entries can either be
Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like
@ -160,10 +159,10 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
''' '''
if priors is not None: if priors is not None:
return _prior_fit(x, y, func, priors, silent=silent, **kwargs) return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
elif (type(x)==dict and type(y)==dict and type(func)==dict): elif (type(x) == dict and type(y) == dict and type(func) == dict):
return _combined_fit(x, y, func, silent=silent, **kwargs) return _combined_fit(x, y, func, silent=silent, **kwargs)
else: else:
return _standard_fit(x, y, func, silent=silent, **kwargs) return _standard_fit(x, y, func, silent=silent, **kwargs)
@ -688,39 +687,40 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
return output return output
def _combined_fit(x,y,func,silent=False,**kwargs):
def _combined_fit(x, y, func, silent=False, **kwargs):
if kwargs.get('correlated_fit') is True: if kwargs.get('correlated_fit') is True:
raise Exception("Correlated fit has not been implemented yet") raise Exception("Correlated fit has not been implemented yet")
output = Fit_result() output = Fit_result()
output.fit_function = func output.fit_function = func
if kwargs.get('num_grad') is True: if kwargs.get('num_grad') is True:
jacobian = num_jacobian jacobian = num_jacobian
hessian = num_hessian hessian = num_hessian
else: else:
jacobian = auto_jacobian jacobian = auto_jacobian
hessian = auto_hessian hessian = auto_hessian
x_all = [] x_all = []
y_all = [] y_all = []
for key in x.keys(): for key in x.keys():
x_all+=x[key] x_all += x[key]
y_all+=y[key] y_all += y[key]
x_all = np.asarray(x_all) x_all = np.asarray(x_all)
if len(x_all.shape) > 2: if len(x_all.shape) > 2:
raise Exception('Unknown format for x values') raise Exception('Unknown format for x values')
# number of fit parameters # number of fit parameters
n_parms_ls = [] n_parms_ls = []
for key in func.keys(): for key in func.keys():
if not callable(func[key]): if not callable(func[key]):
raise TypeError('func (key='+ key + ') is not a function.') raise TypeError('func (key=' + key + ') is not a function.')
if len(x[key]) != len(y[key]): if len(x[key]) != len(y[key]):
raise Exception('x and y input (key='+ key + ') do not have the same length') raise Exception('x and y input (key=' + key + ') do not have the same length')
for i in range(42): for i in range(42):
try: try:
func[key](np.arange(i), x_all.T[0]) func[key](np.arange(i), x_all.T[0])
@ -731,41 +731,41 @@ def _combined_fit(x,y,func,silent=False,**kwargs):
else: else:
break break
else: else:
raise RuntimeError("Fit function (key="+ key + ") is not valid.") raise RuntimeError("Fit function (key=" + key + ") is not valid.")
n_parms = i n_parms = i
n_parms_ls.append(n_parms) n_parms_ls.append(n_parms)
n_parms = max(n_parms_ls) n_parms = max(n_parms_ls)
if not silent: if not silent:
print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
if 'initial_guess' in kwargs: if 'initial_guess' in kwargs:
x0 = kwargs.get('initial_guess') x0 = kwargs.get('initial_guess')
if len(x0) != n_parms: if len(x0) != n_parms:
raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
else: else:
x0 = [0.1] * n_parms x0 = [0.1] * n_parms
output.method = kwargs.get('method', 'migrad') output.method = kwargs.get('method', 'migrad')
if not silent: if not silent:
print('Method:', output.method) print('Method:', output.method)
def chisqfunc(p): def chisqfunc(p):
chisq = 0.0 chisq = 0.0
for key in func.keys(): for key in func.keys():
x_array = np.asarray(x[key]) x_array = np.asarray(x[key])
model = anp.array(func[key](p,x_array)) model = anp.array(func[key](p, x_array))
y_obs = y[key] y_obs = y[key]
y_f = [o.value for o in y_obs] y_f = [o.value for o in y_obs]
dy_f = [o.dvalue for o in y_obs] dy_f = [o.dvalue for o in y_obs]
C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f C_inv = np.diag(np.diag(np.ones((len(x_array), len(x_array))))) / dy_f / dy_f
chisq += anp.sum((y_f - model)@ C_inv @(y_f - model)) chisq += anp.sum((y_f - model) @ C_inv @ (y_f - model))
return chisq return chisq
if output.method == 'migrad': if output.method == 'migrad':
tolerance = 1e-4 tolerance = 1e-4
if 'tol' in kwargs: if 'tol' in kwargs:
tolerance = kwargs.get('tol') tolerance = kwargs.get('tol')
fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
output.iterations = fit_result.nfev output.iterations = fit_result.nfev
else: else:
tolerance = 1e-12 tolerance = 1e-12
@ -776,23 +776,23 @@ def _combined_fit(x,y,func,silent=False,**kwargs):
chisquare = fit_result.fun chisquare = fit_result.fun
output.message = fit_result.message output.message = fit_result.message
if not fit_result.success: if not fit_result.success:
raise Exception('The minimization procedure did not converge.') raise Exception('The minimization procedure did not converge.')
if x_all.shape[-1] - n_parms > 0: if x_all.shape[-1] - n_parms > 0:
output.chisquare = chisqfunc(fit_result.x) output.chisquare = chisqfunc(fit_result.x)
output.dof = x_all.shape[-1] - n_parms output.dof = x_all.shape[-1] - n_parms
output.chisquare_by_dof = output.chisquare/output.dof output.chisquare_by_dof = output.chisquare / output.dof
output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof) output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
else: else:
output.chisquare_by_dof = float('nan') output.chisquare_by_dof = float('nan')
if not silent: if not silent:
print(fit_result.message) print(fit_result.message)
print('chisquare/d.o.f.:', output.chisquare_by_dof ) print('chisquare/d.o.f.:', output.chisquare_by_dof)
print('fit parameters',fit_result.x) print('fit parameters', fit_result.x)
def chisqfunc_compact(d): def chisqfunc_compact(d):
chisq = 0.0 chisq = 0.0
list_tmp = [] list_tmp = []
@ -800,38 +800,37 @@ def _combined_fit(x,y,func,silent=False,**kwargs):
c2 = 0 c2 = 0
for key in func.keys(): for key in func.keys():
x_array = np.asarray(x[key]) x_array = np.asarray(x[key])
c2+=len(x_array) c2 += len(x_array)
model = anp.array(func[key](d[:n_parms],x_array)) model = anp.array(func[key](d[:n_parms], x_array))
y_obs = y[key] y_obs = y[key]
y_f = [o.value for o in y_obs]
dy_f = [o.dvalue for o in y_obs] dy_f = [o.dvalue for o in y_obs]
C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f C_inv = np.diag(np.diag(np.ones((len(x_array), len(x_array))))) / dy_f / dy_f
list_tmp.append(anp.sum((d[n_parms+c1:n_parms+c2]- model)@ C_inv @(d[n_parms+c1:n_parms+c2]- model))) list_tmp.append(anp.sum((d[n_parms + c1:n_parms + c2] - model) @ C_inv @ (d[n_parms + c1:n_parms + c2] - model)))
c1+=len(x_array) c1 += len(x_array)
chisq = anp.sum(list_tmp) chisq = anp.sum(list_tmp)
return chisq return chisq
def prepare_hat_matrix(): def prepare_hat_matrix():
hat_vector = [] hat_vector = []
for key in func.keys(): for key in func.keys():
x_array = np.asarray(x[key]) x_array = np.asarray(x[key])
if (len(x_array)!= 0): if (len(x_array) != 0):
hat_vector.append(jacobian(func[key])(fit_result.x, x_array)) hat_vector.append(jacobian(func[key])(fit_result.x, x_array))
hat_vector = [item for sublist in hat_vector for item in sublist] hat_vector = [item for sublist in hat_vector for item in sublist]
return hat_vector return hat_vector
fitp = fit_result.x fitp = fit_result.x
y_f = [o.value for o in y_all] y_f = [o.value for o in y_all]
dy_f = [o.dvalue for o in y_all] dy_f = [o.dvalue for o in y_all]
if np.any(np.asarray(dy_f) <= 0.0): if np.any(np.asarray(dy_f) <= 0.0):
raise Exception('No y errors available, run the gamma method first.') raise Exception('No y errors available, run the gamma method first.')
try: try:
hess = hessian(chisqfunc)(fitp) hess = hessian(chisqfunc)(fitp)
except TypeError: except TypeError:
raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f))) jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
# Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
@ -839,29 +838,26 @@ def _combined_fit(x,y,func,silent=False,**kwargs):
deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:]) deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
except np.linalg.LinAlgError: except np.linalg.LinAlgError:
raise Exception("Cannot invert hessian matrix.") raise Exception("Cannot invert hessian matrix.")
if kwargs.get('expected_chisquare') is True: if kwargs.get('expected_chisquare') is True:
if kwargs.get('correlated_fit') is not True: if kwargs.get('correlated_fit') is not True:
W = np.diag(1 / np.asarray(dy_f)) W = np.diag(1 / np.asarray(dy_f))
cov = covariance(y_all) cov = covariance(y_all)
hat_vector = prepare_hat_matrix() hat_vector = prepare_hat_matrix()
A = W @ hat_vector #hat_vector = 'jacobian(func)(fit_result.x, x)' A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)'
P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W) expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
output.chisquare_by_expected_chisquare = chisquare / expected_chisquare output.chisquare_by_expected_chisquare = chisquare / expected_chisquare
if not silent: if not silent:
print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare) print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
result = [] result = []
for i in range(n_parms): for i in range(n_parms):
result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all), man_grad=list(deriv_y[i]))) result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all), man_grad=list(deriv_y[i])))
output.fit_parameters = result
return output
output.fit_parameters = result
return output
def fit_lin(x, y, **kwargs): def fit_lin(x, y, **kwargs):