Repository: scipy/scipy
File: scipy/stats/_stats_py.py (Python)

Code before:
    >>> gmean([1, 4, 7], weights=[3, 1, 3])
    2.80668351922014

    """
    if not isinstance(p, int | float):
        raise ValueError("Power mean only defined for exponent of type int or "
                         "float.")
    if p == 0:
        return gmean(a, axis=axis, dtype=dtype, weights=weights)
    elif math.isinf(p):
        message = "Power mean only implemented for finite `p`"
        raise NotImplementedError(message)

    xp = array_namespace(a, weights)
    dtype = (xp_result_type(a, p, weights, force_floating=True, xp=xp)
             if dtype is None else dtype)
    a = xp.asarray(a, dtype=dtype)

    if weights is not None:
        weights = xp.asarray(weights, dtype=dtype)

    # Linearized approximation for small p to avoid numerical issues; see gh-23407
    p_threshold = 2e-6
    if xp.abs(p) <= p_threshold:
        return _linearized_pmean(a, p, axis=axis, weights=weights, xp=xp)

    negative_mask = a < 0
    if xp.any(negative_mask):
        # `where` avoids having to be careful about dtypes and will work with
        # JAX. This is the exceptional case, so it's OK to be a little slower.
        # Won't work for array_api_strict for now, but see data-apis/array-api#807
        a = xp.where(negative_mask, np.nan, a)
        message = ("The power mean is only defined if all elements are "
                   "non-negative; otherwise, the result is NaN.")
        warnings.warn(message, RuntimeWarning, stacklevel=2)

    with np.errstate(divide='ignore', invalid='ignore'):
        return _xp_mean(a**float(p), axis=axis, weights=weights)**(1/p)


ModeResult = namedtuple('ModeResult', ('mode', 'count'))


def _mode_result(mode, count):
    # When a slice is empty, `_axis_nan_policy` automatically produces
    # NaN for `mode` and `count`. This is a reasonable convention for `mode`,
    # but `count` should not be NaN; it should be zero.
    xp = array_namespace(mode, count)
    i = xp.isnan(count)
    if i.shape == ():
        count = xp.asarray(0, dtype=count.dtype)[()] if i else count

Review comment (bug):
This function doesn't support negative `p`. Please move this below the `negative_mask` block, right above the `np.errstate` block. 

With that done, you can remove the `xp.abs(p)`, so it will be back to the code I suggested. That should fix the `torch` failures we're seing in the [GPU workflow](https://github.com/scipy/scipy/actions/runs/20763285674/job/59623407948?pr=24287).

Code after:
    >>> gmean([1, 4, 7], weights=[3, 1, 3])
    2.80668351922014

    """
    if not isinstance(p, int | float):
        raise ValueError("Power mean only defined for exponent of type int or "
                         "float.")
    if p == 0:
        return gmean(a, axis=axis, dtype=dtype, weights=weights)
    elif math.isinf(p):
        message = "Power mean only implemented for finite `p`"
        raise NotImplementedError(message)

    xp = array_namespace(a, weights)
    dtype = (xp_result_type(a, p, weights, force_floating=True, xp=xp)
             if dtype is None else dtype)
    a = xp.asarray(a, dtype=dtype)

    if weights is not None:
        weights = xp.asarray(weights, dtype=dtype)

    negative_mask = a < 0
    if xp.any(negative_mask):
        # `where` avoids having to be careful about dtypes and will work with
        # JAX. This is the exceptional case, so it's OK to be a little slower.
        # Won't work for array_api_strict for now, but see data-apis/array-api#807
        a = xp.where(negative_mask, np.nan, a)
        message = ("The power mean is only defined if all elements are "
                   "non-negative; otherwise, the result is NaN.")
        warnings.warn(message, RuntimeWarning, stacklevel=2)

    # Linearized approximation for small p to avoid numerical issues; see gh-23407
    p_threshold = 2e-6
    if abs(p) <= p_threshold:
        return _linearized_pmean(a, p, axis=axis, weights=weights, xp=xp)

    with np.errstate(divide='ignore', invalid='ignore'):
        return _xp_mean(a**float(p), axis=axis, weights=weights)**(1/p)


ModeResult = namedtuple('ModeResult', ('mode', 'count'))


def _mode_result(mode, count):
    # When a slice is empty, `_axis_nan_policy` automatically produces
    # NaN for `mode` and `count`. This is a reasonable convention for `mode`,
    # but `count` should not be NaN; it should be zero.
    xp = array_namespace(mode, count)
    i = xp.isnan(count)
    if i.shape == ():
        count = xp.asarray(0, dtype=count.dtype)[()] if i else count