From ee2944d5b02bcae67b91c12add5f6073f17d2071 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 1 Mar 2023 16:12:31 +0000 Subject: [PATCH 1/3] refactor: chisqfunc rewritten as sum over residuals. --- pyerrors/fits.py | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 7f57c5df..2f635b4a 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -580,6 +580,17 @@ def _combined_fit(x, y, func, silent=False, **kwargs): else: x0 = [0.1] * n_parms + if kwargs.get('correlated_fit') is True: + def chisqfunc_residuals_corr(p): + model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) + chisq = anp.dot(chol_inv, (y_f - model)) + return chisq + + def chisqfunc_residuals(p): + model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) + chisq = ((y_f - model) / dy_f) + return chisq + if kwargs.get('correlated_fit') is True: corr = covariance(y_all, correlation=True, **kwargs) covdiag = np.diag(1 / np.asarray(dy_f)) @@ -592,15 +603,10 @@ def _combined_fit(x, y, func, silent=False, **kwargs): chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True) def chisqfunc_corr(p): - model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) - chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2) - return chisq + return anp.sum(chisqfunc_residuals_corr(p) ** 2) def chisqfunc(p): - func_list = np.concatenate([[funcd[k]] * len(xd[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 + return anp.sum(chisqfunc_residuals(p) ** 2) output.method = kwargs.get('method', 'Levenberg-Marquardt') if not silent: @@ -627,17 +633,6 @@ def _combined_fit(x, y, func, silent=False, **kwargs): chisquare = fit_result.fun else: - if kwargs.get('correlated_fit') is True: - def chisqfunc_residuals_corr(p): - model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) - chisq = anp.dot(chol_inv, (y_f - model)) - return chisq - - def chisqfunc_residuals(p): - model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) 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') From a140b2ab39d1914a2dd0fcc0b69bbef3d8361fe5 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 1 Mar 2023 16:26:37 +0000 Subject: [PATCH 2/3] refactor: removed redundant formulations of the chisquare function in least_squares. --- pyerrors/fits.py | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 2f635b4a..9aa9e8bf 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -580,17 +580,6 @@ def _combined_fit(x, y, func, silent=False, **kwargs): else: x0 = [0.1] * n_parms - if kwargs.get('correlated_fit') is True: - def chisqfunc_residuals_corr(p): - model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) - chisq = anp.dot(chol_inv, (y_f - model)) - return chisq - - def chisqfunc_residuals(p): - model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls]) - chisq = ((y_f - model) / dy_f) - return chisq - if kwargs.get('correlated_fit') is True: corr = covariance(y_all, correlation=True, **kwargs) covdiag = np.diag(1 / np.asarray(dy_f)) @@ -602,9 +591,24 @@ def _combined_fit(x, y, func, silent=False, **kwargs): chol = np.linalg.cholesky(corr) chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True) + def general_chisqfunc_corr(p, ivars): + model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls]) + return anp.dot(chol_inv, (ivars - model)) + + def general_chisqfunc(p, ivars): + model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls]) + return ((ivars - model) / dy_f) + + if kwargs.get('correlated_fit') is True: + def chisqfunc_residuals_corr(p): + return general_chisqfunc_corr(p, y_f) + def chisqfunc_corr(p): return anp.sum(chisqfunc_residuals_corr(p) ** 2) + def chisqfunc_residuals(p): + return general_chisqfunc(p, y_f) + def chisqfunc(p): return anp.sum(chisqfunc_residuals(p) ** 2) @@ -700,16 +704,10 @@ def _combined_fit(x, y, func, silent=False, **kwargs): if kwargs.get('correlated_fit') is True: def chisqfunc_compact(d): - func_list = np.concatenate([[funcd[k]] * len(xd[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(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2) - return chisq + return anp.sum(general_chisqfunc_corr(d[:n_parms], d[n_parms:]) ** 2) else: def chisqfunc_compact(d): - func_list = np.concatenate([[funcd[k]] * len(xd[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 + return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms:]) ** 2) jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f))) From 2ac38515b6ec3a913ddb1410dc66554841435923 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 2 Mar 2023 18:39:57 +0000 Subject: [PATCH 3/3] refactor: logic in least square fits simplified. Co-authored-by: Simon Kuberski --- pyerrors/fits.py | 67 ++++++++++++++++++++++-------------------------- 1 file changed, 30 insertions(+), 37 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 9aa9e8bf..dbbebc61 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -580,6 +580,13 @@ def _combined_fit(x, y, func, silent=False, **kwargs): else: x0 = [0.1] * n_parms + def general_chisqfunc_uncorr(p, ivars): + model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls]) + return ((ivars - model) / dy_f) + + def chisqfunc_uncorr(p): + return anp.sum(general_chisqfunc_uncorr(p, y_f) ** 2) + if kwargs.get('correlated_fit') is True: corr = covariance(y_all, correlation=True, **kwargs) covdiag = np.diag(1 / np.asarray(dy_f)) @@ -591,26 +598,15 @@ def _combined_fit(x, y, func, silent=False, **kwargs): chol = np.linalg.cholesky(corr) chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True) - def general_chisqfunc_corr(p, ivars): + def general_chisqfunc(p, ivars): model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls]) return anp.dot(chol_inv, (ivars - model)) - def general_chisqfunc(p, ivars): - model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls]) - return ((ivars - model) / dy_f) - - if kwargs.get('correlated_fit') is True: - def chisqfunc_residuals_corr(p): - return general_chisqfunc_corr(p, y_f) - - def chisqfunc_corr(p): - return anp.sum(chisqfunc_residuals_corr(p) ** 2) - - def chisqfunc_residuals(p): - return general_chisqfunc(p, y_f) - - def chisqfunc(p): - return anp.sum(chisqfunc_residuals(p) ** 2) + def chisqfunc(p): + return anp.sum(general_chisqfunc(p, y_f) ** 2) + else: + general_chisqfunc = general_chisqfunc_uncorr + chisqfunc = chisqfunc_uncorr output.method = kwargs.get('method', 'Levenberg-Marquardt') if not silent: @@ -621,17 +617,17 @@ def _combined_fit(x, y, func, silent=False, **kwargs): 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 + fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef if kwargs.get('correlated_fit') is True: - fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=tolerance) + fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance) 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) + fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance) if kwargs.get('correlated_fit') is True: - fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=tolerance) + fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance) output.iterations = fit_result.nit chisquare = fit_result.fun @@ -640,15 +636,19 @@ def _combined_fit(x, y, func, silent=False, **kwargs): 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) + def chisqfunc_residuals_uncorr(p): + return general_chisqfunc_uncorr(p, y_f) + + fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) if kwargs.get('correlated_fit') is True: - fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) + + def chisqfunc_residuals(p): + return general_chisqfunc(p, y_f) + + fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) chisquare = np.sum(fit_result.fun ** 2) - if kwargs.get('correlated_fit') is True: - assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14) - else: - assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) + assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) output.iterations = fit_result.nfev @@ -695,19 +695,12 @@ def _combined_fit(x, y, func, silent=False, **kwargs): raise Exception('No y errors available, run the gamma method first.') try: - if kwargs.get('correlated_fit') is True: - hess = hessian(chisqfunc_corr)(fitp) - else: - hess = hessian(chisqfunc)(fitp) + 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 - if kwargs.get('correlated_fit') is True: - def chisqfunc_compact(d): - return anp.sum(general_chisqfunc_corr(d[:n_parms], d[n_parms:]) ** 2) - else: - def chisqfunc_compact(d): - return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms:]) ** 2) + def chisqfunc_compact(d): + return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms:]) ** 2) jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))