Repository: scipy/scipy
Pull request: BUG: interpolate.Akima1DInterpolator: Identical Rearrangement Fixes
The subset of changes proposed in #22888 that is strictly an identical numerical stability improvement. The tests for the fix to #22882 are marked as skipped for now.

#### Reference issue
closes #22880 

#### What does this implement/fix?
As requested in the comments of  #22888, I have split the pull request into a smaller one that contains only the identical numerical improvement that reduce the potential for overflow errors, while the other changes that either slightly break backwards compatibility or require an additional keyword argument (or both) remain under discussion. 

I marked the 3 tests for #22882 as skipped. Two of the tests would pass with only the addition of the explicit handling of the m2=m3 case, for the third to pass either `break_mult` would have to be set to exactly zero in:

` ind = np.nonzero(f12 > break_mult * np.max(f12, initial=-np.inf))`

or the `np.max(f12, initial=-np.inf))` would need to be set to a replaced with a similar quantity that depends only on nearby points, as discussed further in #22888.



Problem statement:
BUG: interpolate: Akima1DInterpolator raises ValueError for valid inputs
### Describe your issue.

It is possible to get a `ValueError: `dydx` must contain only finite values.` when calling Akima1DInterpolator with finite-valued x and ys, for example if the input contains a sharp feature at a region where the value of y is also large. This can happen even if (y[i+1]-y[i])/(x[i+1]-x[i]) is finite everywhere and well within the limits of 64 bit floating point precision.  The error message itself is also confusing, because Akima1DInterpolator does not take dydx as an argument at all. 

This appears to occur because of a floating point overflow in line 539/540 of `_cubic.py`, where the multiplication of f1 and m in the numerator can overflow a 64 bit float even if the denominator f12 is also large and would therefore cancel the overflow. It is straightforward to rearrange this equation to avoid the potential for overflow (e.g. (f1/f12)*m1+(f2/f12)*m2 could be one possibility for an alternate arrangement). 

I attached a simple example I created to replicate the behavior. GSL's akima spline implementation performs as expected for this test case, so it should be a case that a similar spline implementation could be reasonably expected to handle without a floating point overflow.



### Reproducing Code Example

```python
import numpy as np
from scipy.interpolate import Akima1DInterpolator

xs = np.arange(1,10)
ys = 1.e160*np.heaviside(xs-4,0.5)

m = (ys[1:]-ys[:-1])/(xs[1:]-xs[:-1])

assert np.all(np.isfinite(m))
assert np.all(m < 1.e-148 * np.finfo(np.float64).max)

akima_spline = Akima1DInterpolator(xs, ys, extrapolate=False,method='makima')
```

### Error message

```shell
conda/envs/Akima/lib/python3.12/site-packages/scipy/interpolate/_cubic.py:539: RuntimeWarning: overflow encountered in multiply
  t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
conda/envs/Akima/lib/python3.12/site-packages/scipy/interpolate/_cubic.py:540: RuntimeWarning: overflow encountered in multiply
  f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
Traceback (most recent call last):
  File "scipy_spline_bug_demo.py", line 12, in <module>
    akima_spline = Akima1DInterpolator(xs, ys, extrapolate=False,method='akima')
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "conda/envs/Akima/lib/python3.12/site-packages/scipy/interpolate/_cubic.py", line 542, in __init__
    super().__init__(x, y, t, axis=0, extrapolate=extrapolate)
  File "conda/envs/Akima/lib/python3.12/site-packages/scipy/interpolate/_cubic.py", line 143, in __init__
    x, dx, y, axis, dydx = prepare_input(x, y, axis, dydx)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "conda/envs/Akima/lib/python3.12/site-packages/scipy/interpolate/_cubic.py", line 59, in prepare_input
    raise ValueError("`dydx` must contain only finite values.")
ValueError: `dydx` must contain only finite values.
```

### SciPy/NumPy/Python version and system information

```shell
{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "18.1.8",
      "commands": "arm64-apple-darwin20.0.0-clang",
      "args": "-ftree-vectorize, -fPIC, -fstack-protector-strong, -O2, -pipe, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix, -D_FORTIFY_SOURCE=2, -isystem, /Users/mcdigman/conda/envs/Akima/include, -mmacosx-version-min=11.0, -mmacosx-version-min=11.0",
      "linker args": "-Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs, -Wl,-rpath,/Users/mcdigman/conda/envs/Akima/lib, -L/Users/mcdigman/conda/envs/Akima/lib, -ftree-vectorize, -fPIC, -fstack-protector-strong, -O2, -pipe, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix, -D_FORTIFY_SOURCE=2, -isystem, /Users/mcdigman/conda/envs/Akima/include, -mmacosx-version-min=11.0, -mmacosx-version-min=11.0"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.12",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "18.1.8",
      "commands": "arm64-apple-darwin20.0.0-clang++",
      "args": "-ftree-vectorize, -fPIC, -fstack-protector-strong, -O2, -pipe, -stdlib=libc++, -fvisibility-inlines-hidden, -fmessage-length=0, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix, -D_FORTIFY_SOURCE=2, -isystem, /Users/mcdigman/conda/envs/Akima/include, -mmacosx-version-min=11.0, -mmacosx-version-min=11.0",
      "linker args": "-Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs, -Wl,-rpath,/Users/mcdigman/conda/envs/Akima/lib, -L/Users/mcdigman/conda/envs/Akima/lib, -ftree-vectorize, -fPIC, -fstack-protector-strong, -O2, -pipe, -stdlib=libc++, -fvisibility-inlines-hidden, -fmessage-length=0, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix, -D_FORTIFY_SOURCE=2, -isystem, /Users/mcdigman/conda/envs/Akima/include, -mmacosx-version-min=11.0, -mmacosx-version-min=11.0"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld64",
      "version": "13.2.0",
      "commands": "/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/_build_env/bin/arm64-apple-darwin20.0.0-gfortran",
      "args": "-march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix",
      "linker args": "-Wl,-headerpad_max_install_names, -Wl,-dead_strip_dylibs, -Wl,-rpath,/Users/mcdigman/conda/envs/Akima/lib, -L/Users/mcdigman/conda/envs/Akima/lib, -march=armv8.3-a, -ftree-vectorize, -fPIC, -fno-stack-protector, -O2, -pipe, -isystem, /Users/mcdigman/conda/envs/Akima/include, -fdebug-prefix-map=/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/work=/usr/local/src/conda/scipy-split-1.15.2, -fdebug-prefix-map=/Users/mcdigman/conda/envs/Akima=/usr/local/src/conda-prefix"
    },
    "pythran": {
      "version": "0.17.0",
      "include directory": "../../_build_env/venv/lib/python3.12/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "arm64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "darwin"
    },
    "cross-compiled": true
  },
  "Build Dependencies": {
    "blas": {
      "name": "blas",
      "found": true,
      "version": "3.9.0",
      "detection method": "pkgconfig",
      "include directory": "/Users/mcdigman/conda/envs/Akima/include",
      "lib directory": "/Users/mcdigman/conda/envs/Akima/lib",
      "openblas configuration": "unknown",
      "pc file directory": "/Users/mcdigman/conda/envs/Akima/lib/pkgconfig"
    },
    "lapack": {
      "name": "lapack",
      "found": true,
      "version": "3.9.0",
      "detection method": "pkgconfig",
      "include directory": "/Users/mcdigman/conda/envs/Akima/include",
      "lib directory": "/Users/mcdigman/conda/envs/Akima/lib",
      "openblas configuration": "unknown",
      "pc file directory": "/Users/mcdigman/conda/envs/Akima/lib/pkgconfig"
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.13.6",
      "detection method": "pkgconfig",
      "include directory": "/Users/mcdigman/conda/envs/Akima/include"
    }
  },
  "Python Information": {
    "path": "/Users/runner/miniforge3/conda-bld/scipy-split_1739790747129/_build_env/bin/python3.12",
    "version": "3.12"
  }
}
```

Hints from issue discussion:
If a simple rearrangement helps numerical stability, great, would you be interested in submitting a PR with a fix?
@ev-br I can, although in terms of numerical stability the issue I describe in the related issue #22882 (which I found by inspection while trying to diagnose this error) is a more significant problem and can be fixed at the same time. As I describe in my most recent comment there, I think the best fix for that issue would slightly change the behavior in some edge cases. I haven't contributed to scipy before though, and I don't know what scipy's standards (or unit testing) are for something that (slightly) changes the behavior vs purely improving the numerical stability, or if it is a good idea to have a keyword argument as I mention in the other comment. If it is the best way to do it, I can just submit a pull request that fixes both issues the way I would do it, and let others decide if it is good.

Commit message: Update scipy/interpolate/tests/test_interpolate.py

[skip ci]

Patch to review:
diff --git a/scipy/interpolate/_cubic.py b/scipy/interpolate/_cubic.py
index ead4e5ebb7d8..feefcbc6301d 100644
--- a/scipy/interpolate/_cubic.py
+++ b/scipy/interpolate/_cubic.py
@@ -538,13 +538,24 @@ def __init__(self, x, y, axis=0, *, method: Literal["akima", "makima"]="akima",
             else:
                 f1 = dm[2:]
                 f2 = dm[:-2]
+
+            # makima is more numerically stable for small f12,
+            # so a finite cutoff should not improve any behavior
+            # however, akima has a qualitative discontinuity near f12=0
+            # a finite cutoff moves it, but cannot remove it.
+            # the cutoff break_mult could be made a keyword argument
+            # method='akima' also benefits from a check for m2=m3
+            break_mult = 1.e-9
+
             f12 = f1 + f2
+
             # These are the mask of where the slope at breakpoint is defined:
-            ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
+            ind = np.nonzero(f12 > break_mult * np.max(f12, initial=-np.inf))
             x_ind, y_ind = ind[0], ind[1:]
             # Set the slope at breakpoint
-            t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
-                    f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
+            t[ind] = m[(x_ind + 1,) + y_ind] + (
+                    (f2[ind] / f12[ind]) * (
+                        m[(x_ind + 2,) + y_ind] - m[(x_ind + 1,) + y_ind]))
 
         super().__init__(x, y, t, axis=0, extrapolate=extrapolate)
         self.axis = axis
diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py
index 7e058487ce9a..00edf647bf3f 100644
--- a/scipy/interpolate/tests/test_interpolate.py
+++ b/scipy/interpolate/tests/test_interpolate.py
@@ -994,6 +994,21 @@ def test_extrapolate_attr(self):
         xp_assert_close(y_ext, ak_true(x_ext), atol=1e-15)
 
 
+    def test_no_overflow(self):
+        # check a large jump does not cause a float overflow
+        x = np.arange(1, 10)
+        y = 1.e6*np.sqrt(np.finfo(float).max)*np.heaviside(x-4, 0.5)
+
+        ak1 = Akima1DInterpolator(x, y, method='makima')
+        ak2 = Akima1DInterpolator(x, y, method='akima')
+
+        y_eval1 = ak1(x)
+        y_eval2 = ak2(x)
+
+        assert np.isfinite(y_eval1).all()
+        assert np.isfinite(y_eval2).all()
+
+
 @pytest.mark.parametrize("method", [Akima1DInterpolator, PchipInterpolator])
 def test_complex(method):
     # Complex-valued data deprecated


Code review comments:
diff_hunk: @@ -538,13 +538,24 @@ def __init__(self, x, y, axis=0, *, method: Literal["akima", "makima"]="akima",
             else:
                 f1 = dm[2:]
                 f2 = dm[:-2]
+
+            # makima is more numerically stable for small f12,
+            # so a finite cutoff should not improve any behavior
+            # however, akima has a qualitative discontinuity near f12=0
+            # a finite cutoff moves it, but cannot remove it.
+            # the cutoff break_mult could be made a keyword argument
+            # method='akima' also benefits from a check for m2=m3
+            break_mult = 1.e-9
+
             f12 = f1 + f2
+
             # These are the mask of where the slope at breakpoint is defined:
-            ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
+            ind = np.nonzero(f12 > break_mult * np.max(f12, initial=-np.inf))
             x_ind, y_ind = ind[0], ind[1:]
             # Set the slope at breakpoint
-            t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
-                    f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
+            t[ind] = m[(x_ind + 1,) + y_ind] + (
+                    (f2[ind] / f12[ind]) * (
+                        m[(x_ind + 2,) + y_ind] - m[(x_ind + 1,) + y_ind]))
path: scipy/interpolate/_cubic.py
text: @user1:
```suggestion
            t[ind] = m[(x_ind + 1,) + y_ind] + (
                (f2[ind] / f12[ind])
                * (m[(x_ind + 2,) + y_ind] - m[(x_ind + 1,) + y_ind])
            )
```