Repository: scipy/scipy
Pull request: BUG: cluster: Avoid OOB write when distances are NaN in `centroid`
If a distance matrix with negative distances is passed, to centroid or median linkage, this can result in NaN distances written back to the distance matrix. If this happens, none of those distances will be considered less than inf, so find_min_distance() will return -1, which will eventually be used to index into an array. This causes an OOB write and a crash. Forbid this and throw an exception if this ever happens.

<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
http://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue

Closes #21011

#### What does this implement/fix?

Example test case:

```
from scipy.cluster.hierarchy import centroid

y = [0, 0, -1]
centroid(y)
```

Output:

```
double free or corruption (out)
Aborted (core dumped)
```
<!--Please explain your changes.-->

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
BUG: Double free or corruption when passing invalid arguments to some functions in the `scipy.cluster.hierarchy` module
### Describe your issue.

I was working on Pynguin, an automated unit test generation tool for Python, and the tool found that if we pass invalid arguments to some functions in the `scipy.cluster.hierarchy` module, Scipy crashes with a "double free or corruption" error. It's not a particularly serious bug but it might be a good idea to fix it so that it is more user-friendly in the event of a user error.

### Reproducing Code Example

```python
#######################################################################
# Example 1
#######################################################################

from scipy.cluster.hierarchy import centroid

y = [-23, 17, 48, -76, -84, -88, -23, -74, 12, 35]
centroid(y)

#######################################################################
# Example 2
#######################################################################

from scipy.cluster.hierarchy import median

y = [-23, 17, 48, -76, -84, -88, -23, -74, 12, 35]
median(y)
```


### Error message

```shell
double free or corruption (out)
Aborted (core dumped)
```


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

```shell
1.13.1 2.0.0 sys.version_info(major=3, minor=10, micro=14, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.27
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.27
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-mnl4e8vy/overlay/lib/python3.10/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp310-cp310/bin/python
  version: '3.10'
```

Hints from issue discussion:
I can confirm that this fails on current main.

I minimized the centroid test case to the following:

```
from scipy.cluster.hierarchy import centroid

y = [0, 0, -1]
centroid(y)
```

Here is what I have figured out about this failure:

If _centroid recomputes the distance, and the expression 

```
cdef double _centroid(double d_xi, double d_yi, double d_xy,
                      int size_x, int size_y, int size_i) noexcept:
    return sqrt((((size_x * d_xi * d_xi) + (size_y * d_yi * d_yi)) -
                 (size_x * size_y * d_xy * d_xy) / (size_x + size_y)) /
                (size_x + size_y))
```

returns NaN, that will be saved to the distance matrix. (I believe this can only happen if the triangle inequality is not satisfied for the given distance matrix.)

If find_min_distance() fails to find any node closer than inf, then it will return a node index of -1.

https://github.com/scipy/scipy/blob/v1.13.1/scipy/cluster/_hierarchy.pyx#L757-L772

Since NaN is not smaller than inf, then -1 will be returned. That will be saved into neighbors as the index of the nearest neighbor.

Then, in the following expression, y = 0 and x = -1 is passed into the following expression:

```
            if dist == D[condensed_index(n, x, y)]:
```

This performs an out of bounds memory read. There are also out of bounds memory writes later on.

It seems like the following hardening steps could be taken to fix this:

* In `find_min_distance()` if -1 would be returned as the neighbor index, (i.e. no neighbors are close by) then raise an exception instead.
* If `new_dist` returns NaN, then raise an exception.
@BergLucas

PS: Really cool project, and I see a lot of potential in the automatic generation of test-cases. It occurs to me that you could probably find more issues if you turned on cython bounds checking. For example, in this case, the array D[] is being written to out-of-bounds, but I bet there are other cases where there's an out-of-bounds read, or an out-of-bounds write where the corruption is not noticed.

Commit message: BUG: Avoid OOB write when distances are NaN

If a distance matrix with negative distances is passed, to centroid or
median linkage, this can result in NaN distances written back to
the distance matrix. If this happens, none of those distances will be
considered less than inf, so find_min_distance() will return -1, which
will eventually be used to index into an array. This causes an OOB write
and a crash. Forbid this and throw an exception if this ever happens.

Patch to review:
diff --git a/scipy/cluster/_hierarchy.pyx b/scipy/cluster/_hierarchy.pyx
index 511e47fa0e6b..d801236888dd 100644
--- a/scipy/cluster/_hierarchy.pyx
+++ b/scipy/cluster/_hierarchy.pyx
@@ -769,6 +769,12 @@ cdef Pair find_min_dist(int n, double[:] D, int[:] size, int x):
             current_min = dist
             y = i
 
+    if y == -1:
+        raise ValueError(
+            "find_min_dist cannot find any neighbors closer than inf away. "
+            "Check that distances contain no negative/infinite/NaN entries. "
+        )
+
     return Pair(y, current_min)
 
 
diff --git a/scipy/cluster/tests/test_hierarchy.py b/scipy/cluster/tests/test_hierarchy.py
index 0b6ca2da121f..019bc0ec50fc 100644
--- a/scipy/cluster/tests/test_hierarchy.py
+++ b/scipy/cluster/tests/test_hierarchy.py
@@ -1260,3 +1260,12 @@ def test_Heap(xp):
     pair = heap.get_min()
     assert_equal(pair['key'], 1)
     assert_equal(pair['value'], 10)
+
+
+@skip_xp_backends(np_only=True, reasons=['scipy specific corner case'])
+def test_centroid_neg_distance(xp):
+    # gh-21011
+    values = xp.asarray([0, 0, -1])
+    with pytest.raises(ValueError):
+        # This is just checking that this doesn't crash
+        linkage(values, method='centroid')


Code review comments:
diff_hunk: @@ -1260,3 +1260,12 @@ def test_Heap(xp):
     pair = heap.get_min()
     assert_equal(pair['key'], 1)
     assert_equal(pair['value'], 10)
+
+
+@skip_xp_backends(np_only=True, reasons=['scipy specific corner case'])
path: scipy/cluster/tests/test_hierarchy.py
text: @user2:
This function supports cpu backends so I think we should do this?
```suggestion
@user1(cpu_only=True)
```

@author:
Thanks, applied in dc8d1b524fea8234ff6609ac748ec8317de4e9c1.