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

Merge branch 'develop' into documentation

Browse source
This commit is contained in:
fjosw 2023-02-02 15:09:30 +00:00
commit 822ce887c1
8 changed files with 597 additions and 24 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
examples/01_basic_example.ipynb
View file

@ -33,7 +33,8 @@
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"usetex = matplotlib.checkdep_usetex(True)\n",
"import shutil\n",
"usetex = shutil.which('latex') != ''\n",
"plt.rc('text', usetex=usetex)"
]
},

3
examples/02_correlators.ipynb
View file

@ -21,7 +21,8 @@
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"usetex = matplotlib.checkdep_usetex(True)\n",
"import shutil\n",
"usetex = shutil.which('latex') != ''\n",
"plt.rc('text', usetex=usetex)"
]
},

3
examples/03_pcac_example.ipynb
View file

@ -19,7 +19,8 @@
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"usetex = matplotlib.checkdep_usetex(True)\n",
"import shutil\n",
"usetex = shutil.which('latex') != ''\n",
"plt.rc('text', usetex=usetex)"
]
},

3
examples/04_fit_example.ipynb
View file

@ -19,7 +19,8 @@
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"usetex = matplotlib.checkdep_usetex(True)\n",
"import shutil\n",
"usetex = shutil.which('latex') != ''\n",
"plt.rc('text', usetex=usetex)"
]
},

3
examples/06_gevp.ipynb
View file

@ -20,7 +20,8 @@
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"usetex = matplotlib.checkdep_usetex(True)\n",
"import shutil\n",
"usetex = shutil.which('latex') != ''\n",
"plt.rc('text', usetex=usetex)"
]
},

134
examples/08_combined_fit_example.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

211
pyerrors/fits.py
View file

@ -72,9 +72,12 @@ class Fit_result(Sequence):
def least_squares(x, y, func, priors=None, silent=False, **kwargs):
r'''Performs a non-linear fit to y = func(x).
```
Parameters
----------
For an uncombined fit:
x : list
list of floats.
y : list
@ -96,9 +99,32 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
(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.
OR For a combined fit:
x : dict
dict of lists.
y : dict
dict of lists of Obs.
funcs : dict
dict of objects
fit functions have to be of the form (here a[0] is the common fit parameter)
```python
import autograd.numpy as anp
funcs = {"a": func_a,
"b": func_b}
def func_a(a, x):
return a[1] * anp.exp(-a[0] * x)
def func_b(a, x):
return a[2] * anp.exp(-a[0] * x)
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
will not work.
priors : list, optional
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
@ -114,6 +140,11 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
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.
tol: float, optional
can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence
to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly
invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values
The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max)
correlated_fit : bool
If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`.
@ -137,6 +168,11 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
'''
if priors is not None:
return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
elif (type(x) == dict and type(y) == dict and type(func) == dict):
return _combined_fit(x, y, func, silent=silent, **kwargs)
elif (type(x) == dict or type(y) == dict or type(func) == dict):
raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.")
else:
return _standard_fit(x, y, func, silent=silent, **kwargs)
@ -474,7 +510,6 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs):
def _standard_fit(x, y, func, silent=False, **kwargs):
output = Fit_result()
output.fit_function = func
@ -668,6 +703,180 @@ def _standard_fit(x, y, func, silent=False, **kwargs):
return output
def _combined_fit(x, y, func, silent=False, **kwargs):
if kwargs.get('correlated_fit') is True:
raise Exception("Correlated fit has not been implemented yet")
output = Fit_result()
output.fit_function = func
if kwargs.get('num_grad') is True:
jacobian = num_jacobian
hessian = num_hessian
else:
jacobian = auto_jacobian
hessian = auto_hessian
key_ls = sorted(list(x.keys()))
if sorted(list(y.keys())) != key_ls:
raise Exception('x and y dictionaries do not contain the same keys.')
if sorted(list(func.keys())) != key_ls:
raise Exception('x and func dictionaries do not contain the same keys.')
x_all = np.concatenate([np.array(x[key]) for key in key_ls])
y_all = np.concatenate([np.array(y[key]) for key in key_ls])
y_f = [o.value for o in y_all]
dy_f = [o.dvalue for o in y_all]
if len(x_all.shape) > 2:
raise Exception('Unknown format for x values')
# number of fit parameters
n_parms_ls = []
for key in key_ls:
if not callable(func[key]):
raise TypeError('func (key=' + key + ') is not a function.')
if len(x[key]) != len(y[key]):
raise Exception('x and y input (key=' + key + ') do not have the same length')
for i in range(100):
try:
func[key](np.arange(i), x_all.T[0])
except TypeError:
continue
except IndexError:
continue
else:
break
else:
raise RuntimeError("Fit function (key=" + key + ") is not valid.")
n_parms = i
n_parms_ls.append(n_parms)
n_parms = max(n_parms_ls)
if not silent:
print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
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: %d vs. %d' % (len(x0), n_parms))
else:
x0 = [0.1] * n_parms
def chisqfunc(p):
func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
model = anp.array([func_list[i](p, x_all[i]) for i in range(len(x_all))])
chisq = anp.sum(((y_f - model) / dy_f) ** 2)
return chisq
output.method = kwargs.get('method', 'Levenberg-Marquardt')
if not silent:
print('Method:', output.method)
if output.method != 'Levenberg-Marquardt':
if output.method == 'migrad':
tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic
if 'tol' in kwargs:
tolerance = kwargs.get('tol')
fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
output.iterations = fit_result.nfev
else:
tolerance = 1e-12
if 'tol' in kwargs:
tolerance = kwargs.get('tol')
fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=tolerance)
output.iterations = fit_result.nit
chisquare = fit_result.fun
else:
def chisqfunc_residuals(p):
model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls])
chisq = ((y_f - model) / dy_f)
return chisq
if 'tol' in kwargs:
print('tol cannot be set for Levenberg-Marquardt')
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)
assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
output.iterations = fit_result.nfev
output.message = fit_result.message
if not fit_result.success:
raise Exception('The minimization procedure did not converge.')
if x_all.shape[-1] - n_parms > 0:
output.chisquare = chisquare
output.dof = x_all.shape[-1] - n_parms
output.chisquare_by_dof = output.chisquare / output.dof
output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
else:
output.chisquare_by_dof = float('nan')
if not silent:
print(fit_result.message)
print('chisquare/d.o.f.:', output.chisquare_by_dof)
print('fit parameters', fit_result.x)
def chisqfunc_compact(d):
func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
return chisq
def prepare_hat_matrix():
hat_vector = []
for key in key_ls:
x_array = np.asarray(x[key])
if (len(x_array) != 0):
hat_vector.append(jacobian(func[key])(fit_result.x, x_array))
hat_vector = [item for sublist in hat_vector for item in sublist]
return hat_vector
fitp = fit_result.x
if np.any(np.asarray(dy_f) <= 0.0):
raise Exception('No y errors available, run the gamma method first.')
try:
hess = hessian(chisqfunc)(fitp)
except TypeError:
raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
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
try:
deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
except np.linalg.LinAlgError:
raise Exception("Cannot invert hessian matrix.")
if kwargs.get('expected_chisquare') is True:
if kwargs.get('correlated_fit') is not True:
W = np.diag(1 / np.asarray(dy_f))
cov = covariance(y_all)
hat_vector = prepare_hat_matrix()
A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)'
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)
output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
if not silent:
print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
result = []
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])))
output.fit_parameters = result
return output
def fit_lin(x, y, **kwargs):
"""Performs a linear fit to y = n + m * x and returns two Obs n, m.

261
tests/fits_test.py
View file

@ -108,24 +108,6 @@ def test_prior_fit_num_grad():
auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False, piors=y[:2])
def test_least_squares_num_grad():
x = []
y = []
for i in range(2, 5):
x.append(i * 0.01)
y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens"))
num = pe.fits.least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True)
auto = pe.fits.least_squares(x, y, lambda a, x: anp.exp(a[0] * x) + a[1], num_grad=False)
assert(num[0] == auto[0])
assert(num[1] == auto[1])
assert(num[0] == auto[0])
assert(num[1] == auto[1])
def test_total_least_squares_num_grad():
x = []
y = []
@ -608,6 +590,249 @@ def test_ks_test():
pe.fits.ks_test(fit_res)
def test_combined_fit_list_v_array():
res = []
for y_test in [{'a': [pe.Obs([np.random.normal(i, 0.5, 1000)], ['ensemble1']) for i in range(1, 7)]},
{'a': np.array([pe.Obs([np.random.normal(i, 0.5, 1000)], ['ensemble1']) for i in range(1, 7)])}]:
for x_test in [{'a': [0, 1, 2, 3, 4, 5]}, {'a': np.arange(6)}]:
for key in y_test.keys():
[item.gamma_method() for item in y_test[key]]
def func_a(a, x):
return a[1] * x + a[0]
funcs_test = {"a": func_a}
res.append(pe.fits.least_squares(x_test, y_test, funcs_test))
assert (res[0][0] - res[1][0]).is_zero(atol=1e-8)
assert (res[0][1] - res[1][1]).is_zero(atol=1e-8)
def test_combined_fit_vs_standard_fit():
x_const = {'a':[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 'b':np.arange(10, 20)}
y_const = {'a':[pe.Obs([np.random.normal(1, val, 1000)], ['ensemble1'])
for val in [0.25, 0.3, 0.01, 0.2, 0.5, 1.3, 0.26, 0.4, 0.1, 1.0]],
'b':[pe.Obs([np.random.normal(1, val, 1000)], ['ensemble1'])
for val in [0.5, 1.12, 0.26, 0.25, 0.3, 0.01, 0.2, 1.0, 0.38, 0.1]]}
for key in y_const.keys():
[item.gamma_method() for item in y_const[key]]
y_const_ls = np.concatenate([np.array(o) for o in y_const.values()])
x_const_ls = np.arange(0, 20)
def func_const(a,x):
return 0 * x + a[0]
funcs_const = {"a": func_const,"b": func_const}
for method_kw in ['Levenberg-Marquardt', 'migrad', 'Powell', 'Nelder-Mead']:
res = []
res.append(pe.fits.least_squares(x_const, y_const, funcs_const, method = method_kw, expected_chisquare=True))
res.append(pe.fits.least_squares(x_const_ls, y_const_ls, func_const, method = method_kw, expected_chisquare=True))
[item.gamma_method for item in res]
assert np.isclose(0.0, (res[0].chisquare_by_dof - res[1].chisquare_by_dof), 1e-14, 1e-8)
assert np.isclose(0.0, (res[0].chisquare_by_expected_chisquare - res[1].chisquare_by_expected_chisquare), 1e-14, 1e-8)
assert np.isclose(0.0, (res[0].p_value - res[1].p_value), 1e-14, 1e-8)
assert (res[0][0] - res[1][0]).is_zero(atol=1e-8)
def test_combined_fit_no_autograd():
def func_exp1(x):
return 0.3*np.exp(0.5*x)
def func_exp2(x):
return 0.3*np.exp(0.8*x)
xvals_b = np.arange(0,6)
xvals_a = np.arange(0,8)
def func_a(a,x):
return a[0]*np.exp(a[1]*x)
def func_b(a,x):
return a[0]*np.exp(a[2]*x)
funcs = {'a':func_a, 'b':func_b}
xs = {'a':xvals_a, 'b':xvals_b}
ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)],
'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]}
for key in funcs.keys():
[item.gamma_method() for item in ys[key]]
with pytest.raises(Exception):
pe.least_squares(xs, ys, funcs)
pe.least_squares(xs, ys, funcs, num_grad=True)
def test_combined_fit_invalid_fit_functions():
def func1(a, x):
return a[0] + a[1] * x + a[2] * anp.sinh(x) + a[199]
def func2(a, x, y):
return a[0] + a[1] * x
def func3(x):
return x
def func_valid(a,x):
return a[0] + a[1] * x
xvals =[]
yvals =[]
err = 0.1
for x in range(1, 8, 2):
xvals.append(x)
yvals.append(pe.pseudo_Obs(x + np.random.normal(0.0, err), err, 'test1') + pe.pseudo_Obs(0, err / 100, 'test2', samples=87))
[o.gamma_method() for o in yvals]
for func in [func1, func2, func3]:
with pytest.raises(Exception):
pe.least_squares({'a':xvals}, {'a':yvals}, {'a':func})
with pytest.raises(Exception):
pe.least_squares({'a':xvals, 'b':xvals}, {'a':yvals, 'b':yvals}, {'a':func, 'b':func_valid})
with pytest.raises(Exception):
pe.least_squares({'a':xvals, 'b':xvals}, {'a':yvals, 'b':yvals}, {'a':func_valid, 'b':func})
def test_combined_fit_invalid_input():
xvals =[]
yvals =[]
err = 0.1
def func_valid(a,x):
return a[0] + a[1] * x
for x in range(1, 8, 2):
xvals.append(x)
yvals.append(pe.pseudo_Obs(x + np.random.normal(0.0, err), err, 'test1') + pe.pseudo_Obs(0, err / 100, 'test2', samples=87))
with pytest.raises(Exception):
pe.least_squares({'a':xvals}, {'b':yvals}, {'a':func_valid})
def test_combined_fit_no_autograd():
def func_exp1(x):
return 0.3*np.exp(0.5*x)
def func_exp2(x):
return 0.3*np.exp(0.8*x)
xvals_b = np.arange(0,6)
xvals_a = np.arange(0,8)
def func_a(a,x):
return a[0]*np.exp(a[1]*x)
def func_b(a,x):
return a[0]*np.exp(a[2]*x)
funcs = {'a':func_a, 'b':func_b}
xs = {'a':xvals_a, 'b':xvals_b}
ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)],
'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]}
for key in funcs.keys():
[item.gamma_method() for item in ys[key]]
with pytest.raises(Exception):
pe.least_squares(xs, ys, funcs)
pe.least_squares(xs, ys, funcs, num_grad=True)
def test_combined_fit_num_grad():
def func_exp1(x):
return 0.3*np.exp(0.5*x)
def func_exp2(x):
return 0.3*np.exp(0.8*x)
xvals_b = np.arange(0,6)
xvals_a = np.arange(0,8)
def func_num_a(a,x):
return a[0]*np.exp(a[1]*x)
def func_num_b(a,x):
return a[0]*np.exp(a[2]*x)
def func_auto_a(a,x):
return a[0]*anp.exp(a[1]*x)
def func_auto_b(a,x):
return a[0]*anp.exp(a[2]*x)
funcs_num = {'a':func_num_a, 'b':func_num_b}
funcs_auto = {'a':func_auto_a, 'b':func_auto_b}
xs = {'a':xvals_a, 'b':xvals_b}
ys = {'a':[pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)],
'b':[pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]}
for key in funcs_num.keys():
[item.gamma_method() for item in ys[key]]
num = pe.fits.least_squares(xs, ys, funcs_num, num_grad=True)
auto = pe.fits.least_squares(xs, ys, funcs_auto, num_grad=False)
assert(num[0] == auto[0])
assert(num[1] == auto[1])
def test_combined_fit_dictkeys_no_order():
def func_exp1(x):
return 0.3*np.exp(0.5*x)
def func_exp2(x):
return 0.3*np.exp(0.8*x)
xvals_b = np.arange(0,6)
xvals_a = np.arange(0,8)
def func_num_a(a,x):
return a[0]*np.exp(a[1]*x)
def func_num_b(a,x):
return a[0]*np.exp(a[2]*x)
def func_auto_a(a,x):
return a[0]*anp.exp(a[1]*x)
def func_auto_b(a,x):
return a[0]*anp.exp(a[2]*x)
funcs = {'a':func_auto_a, 'b':func_auto_b}
funcs_no_order = {'b':func_auto_b, 'a':func_auto_a}
xs = {'a':xvals_a, 'b':xvals_b}
xs_no_order = {'b':xvals_b, 'a':xvals_a}
yobs_a = [pe.Obs([np.random.normal(item, item*1.5, 1000)],['ensemble1']) for item in func_exp1(xvals_a)]
yobs_b = [pe.Obs([np.random.normal(item, item*1.4, 1000)],['ensemble1']) for item in func_exp2(xvals_b)]
ys = {'a': yobs_a, 'b': yobs_b}
ys_no_order = {'b': yobs_b, 'a': yobs_a}
for key in funcs.keys():
[item.gamma_method() for item in ys[key]]
[item.gamma_method() for item in ys_no_order[key]]
for method_kw in ['Levenberg-Marquardt', 'migrad', 'Powell', 'Nelder-Mead']:
order = pe.fits.least_squares(xs, ys, funcs,method = method_kw)
no_order_func = pe.fits.least_squares(xs, ys, funcs_no_order,method = method_kw)
no_order_x = pe.fits.least_squares(xs_no_order, ys, funcs,method = method_kw)
no_order_y = pe.fits.least_squares(xs, ys_no_order, funcs,method = method_kw)
no_order_func_x = pe.fits.least_squares(xs_no_order, ys, funcs_no_order,method = method_kw)
no_order_func_y = pe.fits.least_squares(xs, ys_no_order, funcs_no_order,method = method_kw)
no_order_x_y = pe.fits.least_squares(xs_no_order, ys_no_order, funcs,method = method_kw)
assert(no_order_func[0] == order[0])
assert(no_order_func[1] == order[1])
assert(no_order_x[0] == order[0])
assert(no_order_x[1] == order[1])
assert(no_order_y[0] == order[0])
assert(no_order_y[1] == order[1])
assert(no_order_func_x[0] == order[0])
assert(no_order_func_x[1] == order[1])
assert(no_order_func_y[0] == order[0])
assert(no_order_func_y[1] == order[1])
assert(no_order_x_y[0] == order[0])
assert(no_order_x_y[1] == order[1])
def fit_general(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.