BUG: scipy.spatial.Delaunay, scipy.interpolate.LinearNDInterpolator segfault on invalid shapes

[lpsinger]

Describe your issue.

scipy.spatial.Delaunay and scipy.interpolate.LinearNDInterpolator segfault on inputs of invalid shapes.

Reproducing Code Example

>>> import numpy as np
>>> from scipy.spatial import Delaunay
>>> Delaunay(np.random.randn(4, 10, 3))

or

>>> import numpy as np
>>> from scipy.interpolate import LinearNDInterpolator
>>> LinearNDInterpolator(np.random.randn(4, 10, 3), np.random.randn(10))

Error message

Segmentation fault: 11

SciPy/NumPy/Python version and system information

1.13.0 1.26.4 sys.version_info(major=3, minor=11, micro=9, 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= SANDYBRIDGE MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  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= SANDYBRIDGE MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 13.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 13.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 11.3.0
  pythran:
    include directory: ../../../../../../private/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/pip-build-env-mo36y8b_/overlay/lib/python3.11/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: darwin
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: darwin
Python Information:
  path: /private/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/cibw-run-1iia41ej/cp311-macosx_x86_64/build/venv/bin/python
  version: '3.11'

[lucascolley]

I see Qhull and segfault - @andfoy related to gh-20619, maybe?

[andfoy]

This issue is not related to #20619, since it is not related to concurrency. However, I can reproduce it, which implies that some other kind of error is occurring, I can further investigate this one.

[tylerjereddy]

As the cross-linked PR indicates, we just need to error out for the spatial cases since those are invalid inputs. I didn't touch the interpolate case, though it may or may not benefit from my patch if it passes through the same code path (didn't check).

[AtsushiSakai]

In the interpolate case, the LinearNDInterpolator is using the qhull.Delaunay in its initializer.

scipy/scipy/interpolate/interpnd.pyx

Lines 329 to 330 in 4a537b7

	def _calculate_triangulation(self, points):
	self.tri = qhull.Delaunay(points)

In the latest code, I have confirmed that using the code reported in this issue does not cause a segmentation fault, thanks to the validation introduced in PR #20633:

In [1]: import numpy as np

In [2]: from scipy.interpolate import LinearNDInterpolator

In [3]: LinearNDInterpolator(np.random.randn(4, 10, 3), np.random.randn(10))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 LinearNDInterpolator(np.random.randn(4, 10, 3), np.random.randn(10))

File interpnd.pyx:326, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__()

File interpnd.pyx:92, in scipy.interpolate.interpnd.NDInterpolatorBase.__init__()

File interpnd.pyx:330, in scipy.interpolate.interpnd.LinearNDInterpolator._calculate_triangulation()

File _qhull.pyx:1808, in scipy.spatial._qhull.Delaunay.__init__()

ValueError: Input points array must have 2 dimensions.

So, I think we can close this issue.