Repository: scipy/scipy
File: scipy/integrate/_quadrature.py (Python)

Code before:
    N = y.shape[axis]
    last_dx = dx
    if x is not None:
        if x.ndim == 1:
            shapex = [1] * nd
            shapex[axis] = x.shape[0]
            x = xp.reshape(x, tuple(shapex))
        elif x.ndim != y.ndim:
            raise ValueError("If given, shape of x must be 1-D or the same as y.")
        if x.shape[axis] != N:
            raise ValueError("If given, length of x along axis must be the same as y.")

    if N % 2 == 0:
        val = 0.0
        result = 0.0
        slice_all = (slice(None),) * nd

        if N == 2:
            # need at least 3 points in integration axis to form parabolic
            # segment. If there are two points then any of 'avg', 'first',
            # 'last' should give the same result.
            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            if x is not None:
                last_dx = x[slice1] - x[slice2]
            val = val + 0.5 * last_dx * (y[slice1] + y[slice2])
        else:
            # use Simpson's rule on first intervals
            result = _basic_simpson(y, 0, N-3, x, dx, axis, xp=xp)

            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            slice3 = tupleset(slice_all, axis, -3)

            h = xp.stack([dx, dx])
            if x is not None:
                # grab the last two spacings from the appropriate axis
                hm2 = tupleset(slice_all, axis, slice(-2, -1, 1))
                hm1 = tupleset(slice_all, axis, slice(-1, None, 1))

                diffs = xp.diff(x, axis=axis)
                h = [xp.squeeze(diffs[hm2], axis=axis),
                     xp.squeeze(diffs[hm1], axis=axis)]

            # This is the correction for the last interval according to
            # Cartwright.
            # However, I used the equations given at
            # https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
            # A footnote on Wikipedia says:
            # Cartwright 2017, Equation 8. The equation in Cartwright is
            # calculating the first interval whereas the equations in the

Review comment (suggestion):
I suppose this is fine, since `val = 0.0` above.
```suggestion
            val = 0.5 * last_dx * (y[slice1] + y[slice2])
```

Code after:

    """
    xp = array_namespace(y, x, dx)
    y, x, dx = xp_promote(y, x, dx, xp=xp)
    nd = y.ndim
    N = y.shape[axis]
    last_dx = dx
    if x is not None:
        if x.ndim == 1:
            shapex = [1] * nd
            shapex[axis] = x.shape[0]
            x = xp.reshape(x, tuple(shapex))
        elif x.ndim != y.ndim:
            raise ValueError("If given, shape of x must be 1-D or the same as y.")
        if x.shape[axis] != N:
            raise ValueError("If given, length of x along axis must be the same as y.")

    if N % 2 == 0:
        val = 0.0
        result = 0.0
        slice_all = (slice(None),) * nd

        if N == 2:
            # need at least 3 points in integration axis to form parabolic
            # segment. If there are two points then any of 'avg', 'first',
            # 'last' should give the same result.
            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            if x is not None:
                last_dx = x[slice1] - x[slice2]
            val = 0.5 * last_dx * (y[slice1] + y[slice2])
        else:
            # use Simpson's rule on first intervals
            result = _basic_simpson(y, 0, N-3, x, dx, axis, xp=xp)

            slice1 = tupleset(slice_all, axis, -1)
            slice2 = tupleset(slice_all, axis, -2)
            slice3 = tupleset(slice_all, axis, -3)

            h = xp.stack([dx, dx])
            if x is not None:
                # grab the last two spacings from the appropriate axis
                hm2 = tupleset(slice_all, axis, slice(-2, -1, 1))
                hm1 = tupleset(slice_all, axis, slice(-1, None, 1))

                diffs = xp.diff(x, axis=axis)
                h = [xp.squeeze(diffs[hm2], axis=axis),
                     xp.squeeze(diffs[hm1], axis=axis)]

            # This is the correction for the last interval according to
            # Cartwright.