BUG: Cannot call np.median on datetime array

[guilhermeleobas]

Describe the issue:

Not sure if this is by design or actually a bug but one cannot call np.median on datetime arrays.

This was also reported on StackOverflow some years ago.

Reproduce the code example:

import numpy as np
a = np.array([np.datetime64('2005-02-25'), np.datetime64(1, 'Y')])
np.median(a)

Error message:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<__array_function__ internals>", line 180, in median
  File "/Users/guilhermeleobas/micromamba/envs/numba/lib/python3.9/site-packages/numpy/lib/function_base.py", line 3793, in median
    r, k = _ureduce(a, func=_median, axis=axis, out=out,
  File "/Users/guilhermeleobas/micromamba/envs/numba/lib/python3.9/site-packages/numpy/lib/function_base.py", line 3702, in _ureduce
    r = func(a, **kwargs)
  File "/Users/guilhermeleobas/micromamba/envs/numba/lib/python3.9/site-packages/numpy/lib/function_base.py", line 3847, in _median
    rout = mean(part[indexer], axis=axis, out=out)
  File "<__array_function__ internals>", line 180, in mean
  File "/Users/guilhermeleobas/micromamba/envs/numba/lib/python3.9/site-packages/numpy/core/fromnumeric.py", line 3474, in mean
    return _methods._mean(a, axis=axis, dtype=dtype,
  File "/Users/guilhermeleobas/micromamba/envs/numba/lib/python3.9/site-packages/numpy/core/_methods.py", line 179, in _mean
    ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)
numpy.core._exceptions._UFuncBinaryResolutionError: ufunc 'add' cannot use operands with types dtype('<M8[D]') and dtype('<M8[D]')

Python and NumPy Versions:

1.22.4
3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:35:41)
[Clang 16.0.6 ]

Also tried on Numpy 1.24, 1.25 and 1.26.

Runtime Environment:

No response

Context for the issue:

No response

[mwanink]

I could take a look at this if nobody else has yet!

[seberg]

This may need a bit of thought/discussion (i.e. may be a bit slow). Since similar things happen in mean also. You cannot implement mean as sum()/count or the interpolation here using addition for datetimes.
OTOH, if you can rephrase the computation to support it without breaking/slowing down the (fragile) logic for other types, I think that would definitely be OK.

(Although, maybe nobody will mind a small special case for datetimes.)

[kasmith11]

@seberg is this still an open discussion? If so, do you feel as if this would be a good first contribution?