np.average strips units on sum of weights when weights and input have different shapes #19054
Description
Description
The np.average function returns incorrect output for the optionally-returned sum of weights when the input array and weights do not have the same shape (e.g. you want the weighted sum of every row in a 2D array, but the Nth element of each row has the same weight, so the weights are a 1D array).
Expected behavior
The sum of the weights returned by np.average should have the same units as the weights passed as input.
How to Reproduce
Here's a minimal example:
import numpy as np
from astropy import units as u
a = np.arange(9).reshape((3, 3)) * u.m
w = np.arange(3) * u.s
avg, wsum = np.average(a, weights=w, returned=True, axis=-1)
assert np.all(wsum == np.sum(w)) # fails! wsum should have units, but is a bare array
Versions
platform
--------
platform.platform() = 'Linux-5.14.0-570.42.2.el9_6.x86_64-x86_64-with-glibc2.34'
platform.version() = '#1 SMP PREEMPT_DYNAMIC Sun Sep 14 13:59:34 UTC 2025'
platform.python_version() = '3.13.9'
packages
--------
astropy 7.2.0
numpy 2.3.3
scipy 1.16.2
matplotlib 3.10.5
pandas 2.3.1
pyerfa 2.0.1.5