From ee4149e4983a747bfe02fc255b95db1bd2b25cc7 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 5 Oct 2022 17:44:38 +0100 Subject: [PATCH 01/10] feat: least_squares fit error propagation can now also be performed via numerical derivatives. --- pyerrors/fits.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index f97665c8..fd49135a 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -9,8 +9,9 @@ import matplotlib.pyplot as plt from matplotlib import gridspec from scipy.odr import ODR, Model, RealData import iminuit -from autograd import jacobian +from autograd import jacobian as auto_jacobian from autograd import elementwise_grad as egrad +from numdifftools import Jacobian as num_jacobian from .obs import Obs, derived_observable, covariance, cov_Obs @@ -114,6 +115,8 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): If True, a plot which displays fit, data and residuals is generated (default False). qqplot : bool If True, a quantile-quantile plot of the fit result is generated (default False). + num_grad : bool + Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). ''' if priors is not None: return _prior_fit(x, y, func, priors, silent=silent, **kwargs) @@ -174,6 +177,8 @@ def total_least_squares(x, y, func, silent=False, **kwargs): x_shape = x.shape + jacobian = auto_jacobian + if not callable(func): raise TypeError('func has to be a function.') @@ -318,6 +323,11 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): x = np.asarray(x) + if kwargs.get('num_grad') is True: + jacobian = num_jacobian + else: + jacobian = auto_jacobian + if not callable(func): raise TypeError('func has to be a function.') @@ -441,6 +451,11 @@ def _standard_fit(x, y, func, silent=False, **kwargs): x = np.asarray(x) + if kwargs.get('num_grad') is True: + jacobian = num_jacobian + else: + jacobian = auto_jacobian + if x.shape[-1] != len(y): raise Exception('x and y input have to have the same length') From 99e130d33c1a3f5e1235726c65ebec81ee829f1b Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 5 Oct 2022 17:54:25 +0100 Subject: [PATCH 02/10] fix: index of num diff jacobian in least squares fit corrected. --- pyerrors/fits.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index fd49135a..fc99f7fd 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -416,7 +416,10 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): if not m.fmin.is_valid: raise Exception('The minimization procedure did not converge.') - hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(params)) + hess = jacobian(jacobian(chisqfunc))(params) + if kwargs.get('num_grad') is True: + hess = hess[0] + hess_inv = np.linalg.pinv(hess) def chisqfunc_compact(d): model = func(d[:n_parms], x) @@ -424,6 +427,8 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): return chisq jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f))) + if kwargs.get('num_grad') is True: + jac_jac = jac_jac[0] deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] @@ -591,6 +596,8 @@ def _standard_fit(x, y, func, silent=False, **kwargs): hess = jacobian(jacobian(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('num_grad') is True: + hess = hess[0] if kwargs.get('correlated_fit') is True: def chisqfunc_compact(d): @@ -605,6 +612,8 @@ def _standard_fit(x, y, func, silent=False, **kwargs): return chisq jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f))) + if kwargs.get('num_grad') is True: + jac_jac = jac_jac[0] # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv try: From f22614f9993cbdaaadf5f410bf92e5e2e57d2fbe Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Wed, 5 Oct 2022 18:41:28 +0100 Subject: [PATCH 03/10] tests: tests for num diff least square fits added. --- tests/fits_test.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/tests/fits_test.py b/tests/fits_test.py index 013239f5..9c43e375 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -83,8 +83,22 @@ def test_least_squares(): assert math.isclose(pcov[i, i], betac[i].dvalue ** 2, abs_tol=1e-3) +def test_fit_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: a[0] * np.exp(x) + a[1], num_grad=True) + auto = pe.fits.least_squares(x, y, lambda a, x: a[0] * anp.exp(x) + a[1], num_grad=False) + + assert(num[0] == auto[0]) + assert(num[1] == auto[1]) + + def test_alternative_solvers(): - dim = 192 + dim = 92 x = np.arange(dim) y = 2 * np.exp(-0.06 * x) + np.random.normal(0.0, 0.15, dim) yerr = 0.1 + 0.1 * np.random.rand(dim) @@ -158,7 +172,7 @@ def test_correlated_fit(): def test_fit_corr_independent(): - dim = 50 + dim = 30 x = np.arange(dim) y = 0.84 * np.exp(-0.12 * x) + np.random.normal(0.0, 0.1, dim) yerr = [0.1] * dim @@ -470,7 +484,7 @@ def test_correlated_fit_vs_jackknife(): def test_fit_no_autograd(): - dim = 10 + dim = 3 x = np.arange(dim) y = 2 * np.exp(-0.08 * x) + np.random.normal(0.0, 0.15, dim) yerr = 0.1 + 0.1 * np.random.rand(dim) @@ -486,6 +500,8 @@ def test_fit_no_autograd(): with pytest.raises(Exception): pe.least_squares(x, oy, func) + pe.least_squares(x, oy, func, num_grad=True) + with pytest.raises(Exception): pe.total_least_squares(oy, oy, func) From b84ea1cc3ba4902f10657f01445495805844630f Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 6 Oct 2022 10:44:06 +0100 Subject: [PATCH 04/10] feat: double jacobian in standard fit replaced by hessian This greatly improves performance for numerical derivatives and helps with readability. --- pyerrors/fits.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index fc99f7fd..0ba7db55 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -10,8 +10,10 @@ from matplotlib import gridspec from scipy.odr import ODR, Model, RealData import iminuit from autograd import jacobian as auto_jacobian +from autograd import hessian as auto_hessian from autograd import elementwise_grad as egrad from numdifftools import Jacobian as num_jacobian +from numdifftools import Hessian as num_hessian from .obs import Obs, derived_observable, covariance, cov_Obs @@ -458,8 +460,10 @@ def _standard_fit(x, y, func, silent=False, **kwargs): if kwargs.get('num_grad') is True: jacobian = num_jacobian + hessian = num_hessian else: jacobian = auto_jacobian + hessian = auto_hessian if x.shape[-1] != len(y): raise Exception('x and y input have to have the same length') @@ -591,13 +595,11 @@ def _standard_fit(x, y, func, silent=False, **kwargs): fitp = fit_result.x try: if kwargs.get('correlated_fit') is True: - hess = jacobian(jacobian(chisqfunc_corr))(fitp) + hess = hessian(chisqfunc_corr)(fitp) else: - hess = jacobian(jacobian(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('num_grad') is True: - hess = hess[0] if kwargs.get('correlated_fit') is True: def chisqfunc_compact(d): @@ -611,9 +613,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs): chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2) return chisq - jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f))) - if kwargs.get('num_grad') is True: - jac_jac = jac_jac[0] + jac_jac = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f))) # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv try: From 8cf15b651f0ed2056a8da69bee1e9a68bd499571 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 6 Oct 2022 17:55:04 +0100 Subject: [PATCH 05/10] tests: num_grad fit test improved. --- tests/fits_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/fits_test.py b/tests/fits_test.py index 9c43e375..698acd26 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -90,8 +90,8 @@ def test_fit_num_grad(): 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: a[0] * np.exp(x) + a[1], num_grad=True) - auto = pe.fits.least_squares(x, y, lambda a, x: a[0] * anp.exp(x) + a[1], num_grad=False) + 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]) From 58b8383c1ad5c2a4be6cc497872497735d4cb190 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 6 Oct 2022 18:07:19 +0100 Subject: [PATCH 06/10] feat: hessian added to prior fit and odr fit routines. --- pyerrors/fits.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 0ba7db55..d4087caf 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -165,6 +165,8 @@ def total_least_squares(x, y, func, silent=False, **kwargs): corrected by effects caused by correlated input data. This can take a while as the full correlation matrix has to be calculated (default False). + num_grad : bool + Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). Notes ----- @@ -179,7 +181,12 @@ def total_least_squares(x, y, func, silent=False, **kwargs): x_shape = x.shape - jacobian = auto_jacobian + if kwargs.get('num_grad') is True: + jacobian = num_jacobian + hessian = num_hessian + else: + jacobian = auto_jacobian + hessian = auto_hessian if not callable(func): raise TypeError('func has to be a function.') @@ -275,7 +282,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): fitp = out.beta try: - hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel()))) + hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel()))) except TypeError: raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None @@ -284,7 +291,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) return chisq - jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) + jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv try: @@ -297,7 +304,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) return chisq - jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f))) + jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f))) # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv try: @@ -326,9 +333,9 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): x = np.asarray(x) if kwargs.get('num_grad') is True: - jacobian = num_jacobian + hessian = num_hessian else: - jacobian = auto_jacobian + hessian = auto_hessian if not callable(func): raise TypeError('func has to be a function.') @@ -418,7 +425,7 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): if not m.fmin.is_valid: raise Exception('The minimization procedure did not converge.') - hess = jacobian(jacobian(chisqfunc))(params) + hess = hessian(chisqfunc)(params) if kwargs.get('num_grad') is True: hess = hess[0] hess_inv = np.linalg.pinv(hess) @@ -428,7 +435,7 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): chisq = anp.sum(((d[n_parms: n_parms + len(x)] - model) / dy_f) ** 2) + anp.sum(((d[n_parms + len(x):] - d[:n_parms]) / dp_f) ** 2) return chisq - jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f))) + jac_jac = hessian(chisqfunc_compact)(np.concatenate((params, y_f, p_f))) if kwargs.get('num_grad') is True: jac_jac = jac_jac[0] From e4d18fc8c8af09e40c510a30a024a4d011566ae3 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 6 Oct 2022 18:09:37 +0100 Subject: [PATCH 07/10] tests: num_grad test for total least squares added. --- tests/fits_test.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/tests/fits_test.py b/tests/fits_test.py index 698acd26..4018abb8 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -83,7 +83,7 @@ def test_least_squares(): assert math.isclose(pcov[i, i], betac[i].dvalue ** 2, abs_tol=1e-3) -def test_fit_num_grad(): +def test_least_squares_num_grad(): x = [] y = [] for i in range(2, 5): @@ -97,6 +97,20 @@ def test_fit_num_grad(): assert(num[1] == auto[1]) +def test_total_least_squares_num_grad(): + x = [] + y = [] + for i in range(2, 5): + x.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens")) + y.append(pe.pseudo_Obs(i * 0.01, 0.0001, "ens")) + + num = pe.fits.total_least_squares(x, y, lambda a, x: np.exp(a[0] * x) + a[1], num_grad=True) + auto = pe.fits.total_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]) + + def test_alternative_solvers(): dim = 92 x = np.arange(dim) From ca4f83f5d100b095cec30aafc2a1fe61b7618782 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 7 Oct 2022 17:57:02 +0100 Subject: [PATCH 08/10] fix: bug in num_grad version of prior fit fixed. --- pyerrors/fits.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index d4087caf..79c11d39 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -426,8 +426,6 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): raise Exception('The minimization procedure did not converge.') hess = hessian(chisqfunc)(params) - if kwargs.get('num_grad') is True: - hess = hess[0] hess_inv = np.linalg.pinv(hess) def chisqfunc_compact(d): From 89125c91e2713862ae392b435df94226fcbe13a6 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 7 Oct 2022 18:02:14 +0100 Subject: [PATCH 09/10] fix: another bug in num_grad version of prior fit fixed. --- pyerrors/fits.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 79c11d39..fdc11feb 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -434,8 +434,6 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): return chisq jac_jac = hessian(chisqfunc_compact)(np.concatenate((params, y_f, p_f))) - if kwargs.get('num_grad') is True: - jac_jac = jac_jac[0] deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] From f5c9a30266ad6c2a5d835aab00e0025194cca227 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Fri, 7 Oct 2022 18:02:27 +0100 Subject: [PATCH 10/10] tests: tests for num_grad version of prior fits added. --- tests/fits_test.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/fits_test.py b/tests/fits_test.py index 4018abb8..374f590a 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -97,6 +97,35 @@ def test_least_squares_num_grad(): assert(num[1] == auto[1]) +def test_prior_fit_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, priors=y[:2]) + 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 = []