Allow to choose median function for mad_std and median_absolute_deviation by saimn · Pull Request #5835 · astropy/astropy

saimn · 2017-02-23T10:27:07Z

The idea here is to allow the user to choose which function must be used for the median computation for mad_std and median_absolute_deviation (similarly to sigma_clip), which would be useful for arrays containing NaNs. Using np.nanmedian can be much faster than having to use masked arrays (see example below).
Also, np.nanmedian allows to use a multi-dimensional axis, as well as np.median, but np.ma.median cannot do this (though it's possible to bypass this limitation with a reshape).

Example on a data cube:

In [3]: cube.shape
Out[3]: (3681, 337, 336)

Masked array version:

In [4]: %time std = mad_std(mcube.reshape(mcube.shape[0], -1), axis=1)
CPU times: user 1min 56s, sys: 22.4 s, total: 2min 19s
Wall time: 2min 19s

With this PR, np.nanmedian and a multi-dimensional axis:

In [5]: %time std2 = mad_std(cube, axis=(1, 2), func=np.nanmedian)        
CPU times: user 14 s, sys: 2.41 s, total: 16.4 s                             
Wall time: 16.4 s
                                                                       
In [6]: np.allclose(std, std2)                                                              
Out[6]: True

With bottleneck (does not support multi-dimensional axis):

In [9]: %time std3 = mad_std(cube.reshape(mcube.shape[0], -1), axis=1, func=bn.nanmedian)
CPU times: user 12.2 s, sys: 1.82 s, total: 14.1 s
Wall time: 14.1 s

The PR is not polished yet, it needs tests, docstring, changelog etc., but first I wanted to open the discussion (@larrybradley , anyone else ?)

saimn · 2017-02-23T10:28:12Z

Oh, I could have skipped the CI, if someone wants to kill the builds ...

larrybradley · 2017-02-23T12:54:26Z

This looks great to me, please continue! Thanks, @saimn.

MSeifert04 · 2017-02-23T12:56:00Z

This has some overlap with #5232.

I like the possibility of providing a function (🎉). But it seems weird that any providided function is overridden if the input is a masked array.

larrybradley · 2017-02-23T13:04:06Z

Ah, yes, thanks for reminding me about #5232.

For a masked array input, we could turn the masked values to NaNs and use np.nanmedian.

EDIT: I'm not sure if we can give the user a choice for masked arrays.

saimn · 2017-02-23T15:36:02Z

Oh indeed it's quite the same as #5232, but with a different approach. I prefer a simple approach like here where the user has the responsibility to choose the right function (nanmedian is not available before Numpy 1.9, but one could use bottleneck instead), instead of a function that tries to handle everything as in #5232 (and as such has to perform more operations, like np.isnan calls, or potentially re-masking the NaNs in a masked array). @keflavich - What's your opinion ?

For the masked array case, currently it overrides the provided function, which is I think the best choice here.
One could fill the array with NaNs and use np.nanmedian, but np.nanmedian is Numpy 1.9+, and we need to give back a masked array to the user so it needs some extra steps.

saimn · 2017-02-23T15:46:25Z

Reading a bit more #5232, in particular #5232 (comment) and the benchmark part (#5232 (comment) and below). For me it strengthens the idea that it is difficult to make good default choices, and the best is to give the user to possibility to choose :)

keflavich · 2017-02-23T15:56:50Z

@saimn I lean the other way - I don't think it's reasonable to expect users to understand all the intricacies of interactions between different versions of numpy and different array types. Personally, if I ever had to figure out which variant of median to use, I'd go back to the source code in #5232.

crawfordsm · 2017-03-03T16:12:22Z

I definitely like the possibility of the user having the option to chose the function, but I see the argument in that the average user who will just want the result doesn't want to go through all the work to figure out which to use.

So I think the best solution would be a hybrid of the two. Something where the user could specify the function they want to use (and then it would be use that function regardless of what else the user did), but the default would not be an numpy.median but an astropy.stats function that would select the correct median to call depending on the version and the type of array. This then moves the decision tree into one function where it can be maintained and also re-used if necessary.

I have to say I don't like how this one changes the function after the user has specified it depending on the array type as that seems non-intuitive from the users perspective to me and something that is done without any warning to the user.

saimn · 2017-03-07T22:05:35Z

I have changed the default function to None which allows to respect the user's choice even for masked arrays.
Having a astropy.stats.median function which provides a good default choice could indeed be a solution, as long as there is the possibility to override this. Maybe a solution based on #5232 can be added later ?

crawfordsm · 2017-03-10T19:57:14Z

@keflavich are you happy with this merging and then in another pull request adding an astropy.stats.median that chooses the appropriate median?

keflavich · 2017-03-10T20:38:00Z

Yes, I think this approach is fine for now and we can hopefully enhance it later.

saimn · 2017-03-10T22:19:57Z

Ok great ! Added changelog and a test for the multi-dimensional axis case, it should be ready now.

crawfordsm · 2017-03-13T22:07:52Z

@saimn -- it looks like your test for multi-dimensional axis is failing in older versions of numpy.

saimn · 2017-03-13T22:26:30Z

@crawfordsm - Thanks for the notice, it was indeed introduced in Numpy 1.9 for np.median. I will skip the test for older Numpy versions.

saimn · 2017-03-14T09:27:42Z

Failure is due to a timeout in the documentation build:

Warning: Embedding the documentation hyperlinks requires Internet access.
Please check your network connection.
Unable to continue embedding `numpy` links due to a URL Error:
(TimeoutError(110, 'Connection timed out'),)

saimn · 2017-03-28T20:28:43Z

ping @crawfordsm - I think this one is ready (modulo the doc build failure, which you may want to restart ?)

MSeifert04 · 2017-03-28T20:30:59Z

I restarted it. That was the time the numpy documentation was down, right? 😄

MSeifert04 · 2017-03-28T20:35:08Z

    # broadcast the median array before subtraction
    if axis is not None:
-        a_median = np.expand_dims(a_median, axis=axis)
+        if isinstance(axis, (tuple, list)):


Are lists possible? I thought it was only tuple or int. Additional thought: Could this be replaced with a try / except TypeError?

Ah fun, lists are valid for np.median (docstring speaks of a sequence of int) but not for np.mean. I usually prefer to use try/except for real errors, here using isinstance is more explicit imo.
And thanks for restarting 😉

Yeah, that sequence of int is the problem. What if I used axis=np.array([1, 2]) as input? arrays are valid axis inputs for median.

Yes sure, it's valid for this case. And after more testing it seems safe to catch TypeError, so if it's your pleasure to use a Numpy array here let's go for it.

Except that np.expand_dims([0], axis=np.array([1,2])) raises a ValueError, so two exceptions to catch.

MSeifert04 · 2017-03-28T21:40:31Z

+        try:
+            a_median = np.expand_dims(a_median, axis=axis)
+        except (ValueError, TypeError):
+            for ax in sorted(list(axis)):


Does sorted work correctly even for negative axis?

MSeifert04 · 2017-03-28T21:42:50Z

-        a_median = np.expand_dims(a_median, axis=axis)
+        try:
+            a_median = np.expand_dims(a_median, axis=axis)
+        except (ValueError, TypeError):


I thought about:

try: len(axis) except TypeError: axis = [axis] for ax in ....

Replaced with astropy.utils.isiterable which is quite similar.

saimn · 2017-03-29T10:57:57Z

Nice, 2 unrelated failures this time:

a timeout error in test_skip_remote_data_astropy
an error with numpy-dev in TestUfuncCoverage.test_coverage related to the addition of a new heaviside ufunc (ENH: Add the 'heaviside' ufunc. numpy/numpy#8795) cc @mhvk

…tion

saimn · 2017-03-30T15:28:53Z

Rebased to get the heaviside fix, now Travis is green (thanks @mhvk !)

crawfordsm · 2017-03-30T15:44:52Z

Merging! Thanks @saimn !

larrybradley added Enhancement stats labels Feb 23, 2017

larrybradley added this to the v2.0.0 milestone Feb 23, 2017

MSeifert04 reviewed Mar 28, 2017

View reviewed changes

saimn mentioned this pull request Mar 29, 2017

Support heaviside ufunc. #5920

Merged

saimn and others added 5 commits March 30, 2017 14:20

Allow to choose median function for mad_std and median_absolute_devia…

ba258f4

…tion

Don't force the use of np.ma.median for masked arrays

612d16b

Add a test for NaNs and axis

4b442f1

Add changelog entry [skip ci]

24b92b0

Add a test for the multi-dimensional axis case

e9dbf75

saimn added 3 commits March 30, 2017 14:20

Skip multi-dimensional test for Numpy < 1.9

8c802c7

Allow more options for axis, as np.median is more flexible

986d94b

Use utils.isiterable

2a553d7

saimn force-pushed the mad-choose-func branch from 6227d06 to 2a553d7 Compare March 30, 2017 12:22

crawfordsm merged commit d9c8238 into astropy:master Mar 30, 2017

saimn deleted the mad-choose-func branch March 30, 2017 20:25

saimn mentioned this pull request Feb 22, 2018

sigma_clip fails when given a tuple for the axis parameter #7227

Closed

Conversation

saimn commented Feb 23, 2017

Uh oh!

saimn commented Feb 23, 2017

Uh oh!

larrybradley commented Feb 23, 2017

Uh oh!

MSeifert04 commented Feb 23, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

larrybradley commented Feb 23, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

saimn commented Feb 23, 2017

Uh oh!

saimn commented Feb 23, 2017

Uh oh!

keflavich commented Feb 23, 2017

Uh oh!

crawfordsm commented Mar 3, 2017

Uh oh!

saimn commented Mar 7, 2017

Uh oh!

crawfordsm commented Mar 10, 2017

Uh oh!

keflavich commented Mar 10, 2017

Uh oh!

saimn commented Mar 10, 2017

Uh oh!

crawfordsm commented Mar 13, 2017

Uh oh!

saimn commented Mar 13, 2017

Uh oh!

saimn commented Mar 14, 2017

Uh oh!

saimn commented Mar 28, 2017

Uh oh!

MSeifert04 commented Mar 28, 2017

Uh oh!

MSeifert04 Mar 28, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

saimn Mar 28, 2017

Choose a reason for hiding this comment

Uh oh!

MSeifert04 Mar 28, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

saimn Mar 28, 2017

Choose a reason for hiding this comment

Uh oh!

saimn Mar 28, 2017

Choose a reason for hiding this comment

Uh oh!

MSeifert04 Mar 28, 2017

Choose a reason for hiding this comment

Uh oh!

saimn Mar 29, 2017

Choose a reason for hiding this comment

Uh oh!

MSeifert04 Mar 28, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

saimn Mar 29, 2017

Choose a reason for hiding this comment

Uh oh!

saimn commented Mar 29, 2017

Uh oh!

saimn commented Mar 30, 2017

Uh oh!

crawfordsm commented Mar 30, 2017

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

MSeifert04 commented Feb 23, 2017 •

edited

Loading

larrybradley commented Feb 23, 2017 •

edited

Loading

MSeifert04 Mar 28, 2017 •

edited

Loading

MSeifert04 Mar 28, 2017 •

edited

Loading

MSeifert04 Mar 28, 2017 •

edited

Loading