Compute the `xt::median()` over any axes

[DavisVaughan]

I think xt::median() should have the same kind of interface as a reducer. It should be able to be computed over any axes, not just 1 axis.

[ghost]

I need this functionality. I may contribute a version implementing it similar to NumPy where it reduces the median axis down to a size of 1. @tdegeus any implementation hints on naming or namespace?

[tdegeus]

The question about naming, namespace, etc. That's easy. You can just copy from average and make sure that the API is the same as NumPy.

What is more challenging could be (efficient) implementation. And the NumPy docs are vague about it, and I even doubt that their definition is correct.

I my mind, taking median over a sequence of axes is the same as recursively calling median on one axis (each time reducing the dimension by one). In this case the output values (not the output shape) would change if the order the input axes changes. However, what I think that NumPy is doing is flatting according the input axes and then taking the median once. Here is an example:

import numpy as np

A = np.random.random((2, 3, 4))

print(np.median(A, (0, 1)))
print(np.median(A, (1, 0)))
print(np.median(A.T.reshape(4, -1), 1))

print(np.median(np.median(A, 1), 0))
print(np.median(np.median(A, 0), 0))

which prints

[0.37402117 0.20483371 0.51644919 0.75283429]
[0.37402117 0.20483371 0.51644919 0.75283429]
[0.37402117 0.20483371 0.51644919 0.75283429]
[0.4531213  0.40388562 0.50749118 0.75283429]
[0.47376498 0.32346399 0.5106028  0.68236615]

When in doubt I think that our policy is to copy NumPy. But I could imagine having two versions.