diff --git a/examples/my_db.sqlite b/examples/my_db.sqlite new file mode 100644 index 00000000..880c4275 Binary files /dev/null and b/examples/my_db.sqlite differ diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index 1949e68a..a6fefe8d 100644 --- a/pyerrors/__init__.py +++ b/pyerrors/__init__.py @@ -457,7 +457,6 @@ from .obs import * from .correlators import * from .fits import * from .misc import * -from . import combined_fits from . import dirac from . import input from . import linalg diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 180cc440..722f9bc0 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -70,13 +70,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 @@ -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 will not work. - + OR For a combined fit: - + 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. - + x : dict dict of lists. y : dict @@ -117,16 +116,16 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): 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 @@ -160,10 +159,10 @@ 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): + + elif (type(x) == dict and type(y) == dict and type(func) == dict): return _combined_fit(x, y, func, silent=silent, **kwargs) - + else: return _standard_fit(x, y, func, silent=silent, **kwargs) @@ -688,39 +687,40 @@ def _standard_fit(x, y, func, silent=False, **kwargs): 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: 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 - + x_all = [] y_all = [] for key in x.keys(): - x_all+=x[key] - y_all+=y[key] - + x_all += x[key] + y_all += y[key] + x_all = np.asarray(x_all) - + if len(x_all.shape) > 2: raise Exception('Unknown format for x values') - + # number of fit parameters n_parms_ls = [] for key in func.keys(): 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]): - 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): try: func[key](np.arange(i), x_all.T[0]) @@ -731,41 +731,41 @@ def _combined_fit(x,y,func,silent=False,**kwargs): else: break else: - raise RuntimeError("Fit function (key="+ key + ") is not valid.") + 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 - + output.method = kwargs.get('method', 'migrad') if not silent: print('Method:', output.method) - + def chisqfunc(p): chisq = 0.0 for key in func.keys(): 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_f = [o.value 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 - chisq += anp.sum((y_f - model)@ C_inv @(y_f - model)) + 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)) return chisq - + if output.method == 'migrad': tolerance = 1e-4 if 'tol' in kwargs: 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 else: tolerance = 1e-12 @@ -776,23 +776,23 @@ def _combined_fit(x,y,func,silent=False,**kwargs): chisquare = fit_result.fun 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 = chisqfunc(fit_result.x) 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) 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) - + print('chisquare/d.o.f.:', output.chisquare_by_dof) + print('fit parameters', fit_result.x) + def chisqfunc_compact(d): chisq = 0.0 list_tmp = [] @@ -800,38 +800,37 @@ def _combined_fit(x,y,func,silent=False,**kwargs): c2 = 0 for key in func.keys(): x_array = np.asarray(x[key]) - c2+=len(x_array) - model = anp.array(func[key](d[:n_parms],x_array)) + c2 += len(x_array) + model = anp.array(func[key](d[:n_parms], x_array)) y_obs = y[key] - y_f = [o.value 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 - 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) + 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))) + c1 += len(x_array) chisq = anp.sum(list_tmp) return chisq - + def prepare_hat_matrix(): hat_vector = [] for key in func.keys(): 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 = [item for sublist in hat_vector for item in sublist] return hat_vector - + 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] - + 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 @@ -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:]) 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)' + 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 = 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 + output.fit_parameters = result + + return output def fit_lin(x, y, **kwargs):