Repository: scipy/scipy
Issue category: Bug Report

Issue:
BUG: integrate.trapezoid: broadcasting failure after #21524
### Describe your issue.

`scipy.integrate.trapezoid` fails when 2D arrays of identical shape are passed as `x` and `y` arguments.

This failure started in recent nightlies, probably from https://github.com/scipy/scipy/pull/21524.

The same example works with scipy 1.14.1 and earlier.

@Armavica @lucascolley 

### Reproducing Code Example

```python
In [1]: import numpy as np, scipy

In [2]: x = np.random.randn(3, 10)

In [3]: scipy.integrate.trapezoid(x, x)
```


### Error message

```shell
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 scipy.integrate.trapezoid(x, x)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/integrate/_quadrature.py:142, in trapezoid(y, x, dx, axis)
    140     shape = [1] * y.ndim
    141     shape[axis] = y.shape[axis]
--> 142     x = xp.reshape(x, shape)
    143     d = x[tuple(slice1)] - x[tuple(slice2)]
    144 try:

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/_lib/array_api_compat/_internal.py:28, in get_xp.<locals>.inner.<locals>.wrapped_f(*args, **kwargs)
     26 @wraps(f)
     27 def wrapped_f(*args, **kwargs):
---> 28     return f(*args, xp=xp, **kwargs)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/_lib/array_api_compat/common/_aliases.py:382, in reshape(x, shape, xp, copy, **kwargs)
    380     y.shape = shape
    381     return y
--> 382 return xp.reshape(x, shape, **kwargs)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/numpy/core/fromnumeric.py:285, in reshape(a, newshape, order)
    200 @array_function_dispatch(_reshape_dispatcher)
    201 def reshape(a, newshape, order='C'):
    202     """
    203     Gives a new shape to an array without changing its data.
    204
   (...)
    283            [5, 6]])
    284     """
--> 285     return _wrapfunc(a, 'reshape', newshape, order=order)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/numpy/core/fromnumeric.py:59, in _wrapfunc(obj, method, *args, **kwds)
     56     return _wrapit(obj, method, *args, **kwds)
     58 try:
---> 59     return bound(*args, **kwds)
     60 except TypeError:
     61     # A TypeError occurs if the object does have such a method in its
     62     # class, but its signature is not identical to that of NumPy's. This
   (...)
     66     # Call _wrapit from within the except clause to ensure a potential
     67     # exception has a traceback chain.
     68     return _wrapit(obj, method, *args, **kwds)

ValueError: cannot reshape array of size 30 into shape (1,10)
```


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

```shell
1.15.0.dev0+git20241118.0ea90cf 1.26.4 sys.version_info(major=3, minor=12, micro=3, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  lapack:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.13.6
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 15.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 15.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 13.3.0
  pythran:
    include directory: ../../../../../../private/var/folders/g6/rgtlsw6n123b0gt5483s5_cm0000gn/T/pip-build-env-zkcj2f79/overlay/lib/python3.12/site-packages/pythran
    version: 0.17.0
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /private/var/folders/g6/rgtlsw6n123b0gt5483s5_cm0000gn/T/cibw-run-7dkhp5o6/cp312-macosx_arm64/build/venv/bin/python
  version: '3.12'
```



Hints:
Can reproduce.

```
import numpy as np
from scipy.integrate import trapezoid
import scipy
print(scipy.version.version)

x = np.linspace(0, 29, 30).reshape(3, 10)
y = np.linspace(0, 29, 30).reshape(3, 10)

print(trapezoid(y, x=x, axis=0))
print(trapezoid(y, x=x, axis=1))
```
@lucascolley @Armavica, it's not clear what [these lines do](https://github.com/scipy/scipy/blob/main/scipy/integrate/_quadrature.py#L139):

```
        # reshape to correct shape
        shape = [1] * y.ndim
        shape[axis] = y.shape[axis]
        x = xp.reshape(x, shape)
```

why is a reshape necessary?
the pre-existing code included a similar reshape, but of `d` and only when `x.ndim == 1`. I think the new reshape of `x` is serving the same purpose, as `d` is then constructed from the reshaped `x` rather than via `np.diff`. At a guess, @Armavica has accurately replicated `np.diff` for `x.ndim == 1`, but in higher dimensions, where the pre-existing code did not reshape, the new reshape is causing trouble.

`xp.diff` is coming in the next version of the standard, so at the worst we can just revert and wait for that to become available in array-api-compat. Maybe the fix here is easy, though.

Patch:
diff --git a/scipy/integrate/_quadrature.py b/scipy/integrate/_quadrature.py
index 4f423da63895..65faaaa79b85 100644
--- a/scipy/integrate/_quadrature.py
+++ b/scipy/integrate/_quadrature.py
@@ -136,11 +136,16 @@ def trapezoid(y, x=None, dx=1.0, axis=-1):
         d = dx
     else:
         x = _asarray(x, xp=xp, subok=True)
-        # reshape to correct shape
-        shape = [1] * y.ndim
-        shape[axis] = y.shape[axis]
-        x = xp.reshape(x, shape)
-        d = x[tuple(slice1)] - x[tuple(slice2)]
+        if x.ndim == 1:
+            d = x[1:] - x[:-1]
+            # make d broadcastable to y
+            slice3 = [None] * nd
+            slice3[axis] = slice(None)
+            d = d[tuple(slice3)]
+        else:
+            # if x is n-D it should be broadcastable to y
+            x = xp.broadcast_to(x, y.shape)
+            d = x[tuple(slice1)] - x[tuple(slice2)]
     try:
         ret = xp.sum(
             d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0,


Test patch:
diff --git a/scipy/integrate/tests/test_quadrature.py b/scipy/integrate/tests/test_quadrature.py
index e49ded996a2e..ec554d3b0bc8 100644
--- a/scipy/integrate/tests/test_quadrature.py
+++ b/scipy/integrate/tests/test_quadrature.py
@@ -309,14 +309,6 @@ def test_ndim(self, xp):
         r = trapezoid(q, x=z[None, None,:], axis=2)
         xp_assert_close(r, qz)
 
-        # n-d `x` but not the same as `y`
-        r = trapezoid(q, x=xp.reshape(x[:, None, None], (3, 1)), axis=0)
-        xp_assert_close(r, qx)
-        r = trapezoid(q, x=xp.reshape(y[None,:, None], (8, 1)), axis=1)
-        xp_assert_close(r, qy)
-        r = trapezoid(q, x=xp.reshape(z[None, None,:], (13, 1)), axis=2)
-        xp_assert_close(r, qz)
-
         # 1-d `x`
         r = trapezoid(q, x=x, axis=0)
         xp_assert_close(r, qx)
@@ -325,6 +317,33 @@ def test_ndim(self, xp):
         r = trapezoid(q, x=z, axis=2)
         xp_assert_close(r, qz)
 
+    @skip_xp_backends('jax.numpy',
+                      reasons=["JAX arrays do not support item assignment"])
+    @pytest.mark.usefixtures("skip_xp_backends")
+    def test_gh21908(self, xp):
+        # extended testing for n-dim arrays
+        x = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
+        y = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
+
+        out0 = xp.linspace(200, 380, 10)
+        xp_assert_close(trapezoid(y, x=x, axis=0), out0)
+        xp_assert_close(trapezoid(y, x=xp.asarray([0, 10., 20.]), axis=0), out0)
+        # x needs to be broadcastable against y
+        xp_assert_close(
+            trapezoid(y, x=xp.asarray([0, 10., 20.])[:, None], axis=0),
+            out0
+        )
+        with pytest.raises(Exception):
+            # x is not broadcastable against y
+            trapezoid(y, x=xp.asarray([0, 10., 20.])[None, :], axis=0)
+
+        out1 = xp.asarray([ 40.5, 130.5, 220.5])
+        xp_assert_close(trapezoid(y, x=x, axis=1), out1)
+        xp_assert_close(
+            trapezoid(y, x=xp.linspace(0, 9, 10), axis=1),
+            out1
+        )
+
     @skip_xp_invalid_arg
     def test_masked(self, xp):
         # Testing that masked arrays behave as if the function is 0 where