[Fix] Migrate to odrpack because of scipy.odr deprecation

.github/workflows/pytest.yml

@@ -44,7 +44,7 @@ jobs:
 
       - name: Run tests with -Werror
         if: matrix.python-version != '3.14'
-        run: pytest --cov=pyerrors -vv
+        run: pytest --cov=pyerrors -vv -Werror
 
       - name: Run tests without -Werror for python 3.14
         if: matrix.python-version == '3.14'

pyerrors/__init__.py

@@ -394,7 +394,7 @@ def func(a, x):
 For all available options including combined fits to multiple datasets see `pyerrors.fits.least_squares`.
 
 ## Total least squares fits
-`pyerrors` can also fit data with errors on both the dependent and independent variables using the total least squares method also referred to as orthogonal distance regression as implemented in [scipy](https://docs.scipy.org/doc/scipy/reference/odr.html), see `pyerrors.fits.least_squares`. The syntax is identical to the standard least squares case, the only difference being that `x` also has to be a `list` or `numpy.array` of `Obs`.
+`pyerrors` can also fit data with errors on both the dependent and independent variables using the total least squares method also referred to as orthogonal distance regression as implemented in [odrpack](https://pypi.org/project/odrpack/), see `pyerrors.fits.total_least_squares`. The syntax is identical to the standard least squares case, the only difference being that `x` also has to be a `list` or `numpy.array` of `Obs`.
 
 For the full API see `pyerrors.fits` for fits and `pyerrors.roots` for finding roots of functions.
 

pyerrors/fits.py

@@ -7,7 +7,7 @@
 import scipy.stats
 import matplotlib.pyplot as plt
 from matplotlib import gridspec
-from scipy.odr import ODR, Model, RealData
+from odrpack import odr_fit
 import iminuit
 from autograd import jacobian as auto_jacobian
 from autograd import hessian as auto_hessian
@@ -567,7 +567,7 @@ def func(a, x):
 
     Notes
     -----
-    Based on the orthogonal distance regression module of scipy.
+    Based on the odrpack orthogonal distance regression library.
 
     Returns
     -------
@@ -634,33 +634,57 @@ def func(a, x):
         raise Exception('No y errors available, run the gamma method first.')
 
     if 'initial_guess' in kwargs:
-        x0 = kwargs.get('initial_guess')
+        x0 = np.asarray(kwargs.get('initial_guess'), dtype=np.float64)
         if len(x0) != n_parms:
             raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
     else:
-        x0 = [1] * n_parms
-
-    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
-    model = Model(func)
-    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
-    odr.set_job(fit_type=0, deriv=1)
-    out = odr.run()
+        x0 = np.ones(n_parms, dtype=np.float64)
+
+    # odrpack expects f(x, beta), but pyerrors convention is f(beta, x)
+    def wrapped_func(x, beta):
+        return func(beta, x)
+
+    out = odr_fit(
+        wrapped_func,
+        np.asarray(x_f, dtype=np.float64),
+        np.asarray(y_f, dtype=np.float64),
+        beta0=x0,
+        weight_x=1.0 / np.asarray(dx_f, dtype=np.float64) ** 2,
+        weight_y=1.0 / np.asarray(dy_f, dtype=np.float64) ** 2,
+        partol=np.finfo(np.float64).eps,
+        task='explicit-ODR',
+        diff_scheme='central'
+    )
 
     output.residual_variance = out.res_var
 
     output.method = 'ODR'
 
     output.message = out.stopreason
 
-    output.xplus = out.xplus
+    output.xplus = out.xplusd
 
     if not silent:
         print('Method: ODR')
-        print(*out.stopreason)
+        print(out.stopreason)
         print('Residual variance:', output.residual_variance)
 
-    if out.info > 3:
-        raise Exception('The minimization procedure did not converge.')
+    if not out.success:
+        # ODRPACK95 info code structure (see User Guide §4):
+        #   info % 10        -> convergence: 1=sum-of-sq, 2=param, 3=both
+        #   info // 10 % 10  -> 1 = problem not full rank at solution
+        convergence_status = out.info % 10
+        rank_deficient = (out.info // 10 % 10) == 1
+
+        if convergence_status in [1, 2, 3] and rank_deficient:
+            warnings.warn(
+                f"ODR fit is rank deficient (irank={out.irank}, inv_condnum={out.inv_condnum:.2e}). "
+                "This may indicate a vanishing chi-squared (n_obs == n_parms). "
+                "Results may be unreliable.",
+                RuntimeWarning
+            )
+        else:
+            raise Exception('The minimization procedure did not converge.')
 
     m = x_f.size
 
@@ -679,9 +703,9 @@ def odr_chisquare(p):
 
         number_of_x_parameters = int(m / x_f.shape[-1])
 
-        old_jac = jacobian(func)(out.beta, out.xplus)
+        old_jac = jacobian(func)(out.beta, out.xplusd)
         fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
-        fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplus, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0])))
+        fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplusd, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0])))
         new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
 
         A = W @ new_jac
@@ -690,14 +714,14 @@ def odr_chisquare(p):
         if expected_chisquare <= 0.0:
             warnings.warn("Negative expected_chisquare.", RuntimeWarning)
             expected_chisquare = np.abs(expected_chisquare)
-        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
+        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplusd.ravel()))) / expected_chisquare
         if not silent:
             print('chisquare/expected_chisquare:',
                   output.chisquare_by_expected_chisquare)
 
     fitp = out.beta
     try:
-        hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel())))
+        hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplusd.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
 
@@ -706,7 +730,7 @@ def odr_chisquare_compact_x(d):
         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 = hessian(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.xplusd.ravel(), x_f.ravel())))
 
     # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
     try:
@@ -719,7 +743,7 @@ def odr_chisquare_compact_y(d):
         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 = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f)))
+    jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplusd.ravel(), y_f)))
 
     # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
     try:
@@ -733,7 +757,7 @@ def odr_chisquare_compact_y(d):
 
     output.fit_parameters = result
 
-    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
+    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplusd.ravel())))
     output.dof = x.shape[-1] - n_parms
     output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof)
 

setup.py

@@ -25,7 +25,7 @@
       license="MIT",
       packages=find_packages(),
       python_requires='>=3.10.0',
-      install_requires=['numpy>=2.0', 'autograd>=1.7.0', 'numdifftools>=0.9.41', 'matplotlib>=3.9', 'scipy>=1.13', 'iminuit>=2.28', 'h5py>=3.11', 'lxml>=5.0', 'python-rapidjson>=1.20', 'pandas>=2.2'],
+      install_requires=['numpy>=2.0', 'autograd>=1.7.0', 'numdifftools>=0.9.41', 'matplotlib>=3.9', 'scipy>=1.13', 'iminuit>=2.28', 'h5py>=3.11', 'lxml>=5.0', 'python-rapidjson>=1.20', 'pandas>=2.2', 'odrpack>=0.5'],
       extras_require={'test': ['pytest', 'pytest-cov', 'pytest-benchmark', 'hypothesis', 'nbmake', 'flake8']},
       classifiers=[
           'Development Status :: 5 - Production/Stable',

tests/fits_test.py

@@ -1,9 +1,10 @@
+import warnings
 import numpy as np
 import autograd.numpy as anp
 import matplotlib.pyplot as plt
 import math
 import scipy.optimize
-from scipy.odr import ODR, Model, RealData
+from odrpack import odr_fit
 from scipy.linalg import cholesky
 from scipy.stats import norm
 import iminuit
@@ -397,11 +398,21 @@ def func(a, x):
         y = a[0] * anp.exp(-a[1] * x)
         return y
 
-    data = RealData([o.value for o in ox], [o.value for o in oy], sx=[o.dvalue for o in ox], sy=[o.dvalue for o in oy])
-    model = Model(func)
-    odr = ODR(data, model, [0, 0], partol=np.finfo(np.float64).eps)
-    odr.set_job(fit_type=0, deriv=1)
-    output = odr.run()
+    # odrpack expects f(x, beta), but pyerrors convention is f(beta, x)
+    def wrapped_func(x, beta):
+        return func(beta, x)
+
+    output = odr_fit(
+        wrapped_func,
+        np.array([o.value for o in ox]),
+        np.array([o.value for o in oy]),
+        beta0=np.array([0.0, 0.0]),
+        weight_x=1.0 / np.array([o.dvalue for o in ox]) ** 2,
+        weight_y=1.0 / np.array([o.dvalue for o in oy]) ** 2,
+        partol=np.finfo(np.float64).eps,
+        task='explicit-ODR',
+        diff_scheme='central'
+    )
 
     out = pe.total_least_squares(ox, oy, func, expected_chisquare=True)
     beta = out.fit_parameters
@@ -458,6 +469,21 @@ def func(a, x):
     assert((outc.fit_parameters[1] - betac[1]).is_zero())
 
 
+def test_total_least_squares_vanishing_chisquare():
+    """Test that a saturated fit (n_obs == n_parms) works without exception."""
+    def func(a, x):
+        return a[0] + a[1] * x
+
+    x = [pe.pseudo_Obs(1.0, 0.1, 'x0'), pe.pseudo_Obs(2.0, 0.1, 'x1')]
+    y = [pe.pseudo_Obs(1.0, 0.1, 'y0'), pe.pseudo_Obs(2.0, 0.1, 'y1')]
+
+    # Rank-deficient warning may or may not fire depending on platform/solver numerics.
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", message="ODR fit is rank deficient", category=RuntimeWarning)
+        out = pe.total_least_squares(x, y, func, silent=True)
+    assert len(out.fit_parameters) == 2
+
+
 def test_odr_derivatives():
     x = []
     y = []
@@ -502,7 +528,10 @@ def f(a, x):
     assert np.isclose(fitp[1].value, fitp[1].r_values['a'])
     assert np.isclose(fitp[1].value, fitp[1].r_values['b'])
 
-    fitp = pe.fits.total_least_squares(y, y, f)
+    # Rank-deficient warning may or may not fire depending on platform/solver numerics.
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", message="ODR fit is rank deficient", category=RuntimeWarning)
+        fitp = pe.fits.total_least_squares(y, y, f)
 
     assert np.isclose(fitp[0].value, fitp[0].r_values['a'])
     assert np.isclose(fitp[0].value, fitp[0].r_values['b'])
@@ -1431,11 +1460,11 @@ def func(a, x):
     global print_output, beta0
     print_output = 1
     if 'initial_guess' in kwargs:
-        beta0 = kwargs.get('initial_guess')
+        beta0 = np.asarray(kwargs.get('initial_guess'), dtype=np.float64)
         if len(beta0) != n_parms:
             raise Exception('Initial guess does not have the correct length.')
     else:
-        beta0 = np.arange(n_parms)
+        beta0 = np.arange(n_parms, dtype=np.float64)
 
     if len(x) != len(y):
         raise Exception('x and y have to have the same length')
@@ -1463,23 +1492,45 @@ def do_the_fit(obs, **kwargs):
 
         xerr = kwargs.get('xerr')
 
+        # odrpack expects f(x, beta), but pyerrors convention is f(beta, x)
+        def wrapped_func(x, beta):
+            return func(beta, x)
+
         if length == len(obs):
             assert 'x_constants' in kwargs
-            data = RealData(kwargs.get('x_constants'), obs, sy=yerr)
-            fit_type = 2
+            x_data = np.asarray(kwargs.get('x_constants'))
+            y_data = np.asarray(obs)
+            # Ordinary least squares (no x errors)
+            output = odr_fit(
+                wrapped_func,
+                x_data,
+                y_data,
+                beta0=beta0,
+                weight_y=1.0 / np.asarray(yerr) ** 2,
+                partol=np.finfo(np.float64).eps,
+                task='explicit-OLS',
+                diff_scheme='central'
+            )
         elif length == len(obs) // 2:
-            data = RealData(obs[:length], obs[length:], sx=xerr, sy=yerr)
-            fit_type = 0
+            x_data = np.asarray(obs[:length])
+            y_data = np.asarray(obs[length:])
+            # ODR with x errors
+            output = odr_fit(
+                wrapped_func,
+                x_data,
+                y_data,
+                beta0=beta0,
+                weight_x=1.0 / np.asarray(xerr) ** 2,
+                weight_y=1.0 / np.asarray(yerr) ** 2,
+                partol=np.finfo(np.float64).eps,
+                task='explicit-ODR',
+                diff_scheme='central'
+            )
         else:
             raise Exception('x and y do not fit together.')
 
-        model = Model(func)
-
-        odr = ODR(data, model, beta0, partol=np.finfo(np.float64).eps)
-        odr.set_job(fit_type=fit_type, deriv=1)
-        output = odr.run()
         if print_output and not silent:
-            print(*output.stopreason)
+            print(output.stopreason)
             print('chisquare/d.o.f.:', output.res_var)
             print_output = 0
         beta0 = output.beta

tests/obs_test.py

@@ -511,7 +511,7 @@ def test_merge_obs():
     my_obs2 = pe.Obs([np.random.normal(1, .1, 100)], ['t|2'], idl=[range(1, 200, 2)])
     merged = pe.merge_obs([my_obs1, my_obs2])
     diff = merged - (my_obs2 + my_obs1) / 2
-    assert np.isclose(0, diff.value, atol=1e-16)
+    assert np.isclose(0, diff.value, atol=np.finfo(np.float64).eps)
     with pytest.raises(ValueError):
         pe.merge_obs([my_obs1, my_obs1])
     my_covobs = pe.cov_Obs(1.0, 0.003, 'cov')