add ignore_nan option to median_absolute_deviation by keflavich · Pull Request #5232 · astropy/astropy

keflavich · 2016-08-11T08:33:02Z

Simple PR: use np.nanmedian instead of np.median if ignore_nan is set. I would favor it defaulting to True, but to avoid breaking backward compatibility I set it to False. Feedback on that option would be helpful.

bsipocz · 2016-08-11T11:14:01Z

        to compute the MAD of the flattened array.
+    ignore_nan : bool
+        Ignore NaN values (treat them as if they are not in the array) when
+        computing the median.  This will use `np.nanmedian`, which requires


this causes sphinx failure, please change to numpy.nanmedian

bsipocz · 2016-08-11T11:18:03Z

GH shows all green, but this just an annoying bug to show all green after [ci skip] is used (they call it feature tough).

@keflavich - Could you please fix the failures? np.nanmedian was added in 1.10, so the test should be skipped for earlier versions. The sphinx failure is also a trivial one.

bsipocz · 2016-08-11T11:39:44Z

 from numpy.testing import assert_equal
 from numpy.testing.utils import assert_allclose

+NUMPY_LT_1P9 = [int(x) for x in np.__version__.split('.')[:2]] < [1, 9]


Most of these are already available from astropy.utils.compat

pllim · 2016-08-11T13:32:38Z

@larrybradley , do you wish to have a look as well?

keflavich · 2016-08-11T13:35:05Z

The one failure looks like an issue in numpy 1.9: apparently, np.median([1,2,3,4,5,np.nan]) was not nan in that version!

Would anyone like to confirm that interpretation of the error? I'll set the minimum np version to 1.10 as @bsipocz recommended earlier.

astrofrog · 2016-08-11T13:40:27Z

@keflavich - try np.median([1,2,np.nan,4,5]) with numpy 1.9. I think what happened back then was that if there were nans, sorting didn't work?

keflavich · 2016-08-11T13:53:36Z

@astrofrog: agreed that there's a sorting error, but it seems to just put the nan at the end or beginning, not sure which:

>       assert np.isnan(funcs.mad_std([1,2,np.nan,4,5]))
E       assert <ufunc 'isnan'>(2.9652044370112041)
E        +  where <ufunc 'isnan'> = np.isnan
E        +  and   2.9652044370112041 = <function mad_std at 0x10bf2ecf8>([1, 2, nan, 4, 5])
E        +    where <function mad_std at 0x10bf2ecf8> = funcs.mad_std

EDIT: so, I think we just want to keep the version bumped to 1.10 for the test

astrofrog · 2016-08-11T14:10:49Z

@keflavich - I think this needs to work for all Numpys we support - so just add a workaround function for Numpy < 1.10 which does the nanmedian manually

astrofrog · 2016-08-11T14:12:42Z

    if isinstance(a, np.ma.MaskedArray):
        func = np.ma.median
+    elif ignore_nan:
+        func = np.nanmedian


Something like this could do the trick:

try: func = np.nanmedian except AttributeError: func = lambda x: np.median(x[~np.isnan(x)])

I think this would fail if an axis argument is given because boolean indexing returns a 1d array

MSeifert04 · 2016-08-11T14:16:25Z

The function already supports masked arrays so one could also mask the nans:

# mask nan/inf
if ignore_invalid:
    a = np.ma.masked_invalid(a)

# or just ignore nans
# if ignore_nan:
#    a = np.ma.masked_where(np.isnan(a), a)

This would work for masked and unmasked arrays alike and work across different numpy versions. It's however a bit slower to use np.ma.median instead of np.nanmedian and may return a masked value if no not-nan(inf) was present.

larrybradley · 2016-08-11T14:30:54Z

I suggest the same as @MSeifert04. If no invalid values are present we could ensure the output is not masked (with axis=None).

keflavich · 2016-08-11T14:34:06Z

@MSeifert04 I've implemented your suggestion in f887005, which also led me to tackle another corner case (if ignore_nan is specified and nans aren't masked out already). Is this reasonable?

@larrybradley do you want an extra check at the end that the return is not a masked array?

larrybradley · 2016-08-11T14:42:02Z

@keflavich Currently, there is a bug in np.ma.median in that it will return a single-valued masked array instead of a scalar (numpy/numpy#7835). This will be fixed in numpy 1.11.2, but I suggest we implement a workaround to make the behavior consistent across all numpy versions (e.g. see https://github.com/astropy/astropy/blob/master/astropy/stats/sigma_clipping.py#L312).

keflavich · 2016-08-11T14:55:40Z

@larrybradley would you like me to add that same workaround to this PR, then?

larrybradley · 2016-08-11T15:02:46Z

@keflavich Yes. For example, I would expect median_absolute_deviation(np.arange(10)) to always return a scalar (here, 2.5), not a masked array (perhaps there is even a test for that).

MSeifert04 · 2016-08-11T15:32:05Z

@keflavich
I do have to correct my statement about np.ma.median being slower than np.nanmedian. It's a bit more complicated than that. In case axis=None nanmedian is faster but with a specified axis masked_median is faster. (see https://gist.github.com/MSeifert04/1551e80e1aa85b48c9191e727c38f9a9)

But I thought it might be better always to use np.ma.median if the user specifies ignore_nan (btw masked_invalid also masks infinities). That has the advantage that the return type doesn't depend on the numpy version being used. Currently numpy 1.9 users get a masked-array while numpy 1.10 users get a normal array if he specifies an axis.

More general if we want to return scalars instead of returning masked-arrays wouldn't it be better to ship some compat/monkeypatch-function instead of adding the workaround to several functions?

larrybradley · 2016-08-11T16:31:21Z

More general if we want to return scalars instead of returning masked-arrays wouldn't it be better to ship some compat/monkeypatch-function instead of adding the workaround to several functions?

Yes, I should have suggested that before. For photutils, I use a private function to smooth out the inconsistencies:
https://github.com/astropy/photutils/blob/master/photutils/background/core.py#L35

crawfordsm · 2016-08-11T16:57:11Z

Yes, if the test isn't already there, I would suggest adding the test with this addition that median_absolute_deviation(np.arange(10)) should always return a scalar and not a masked array.

MSeifert04 · 2016-08-12T09:36:48Z



-def median_absolute_deviation(a, axis=None):
+def median_absolute_deviation(data, axis=None, ignore_nan=False):


If you're going to rename the parameter #5214 would be ideal here 😄

(sorry for double posting this)

Sorry, I should have mentioned that the PR is not yet reviewed nor merged.

If you are going to use stuff from #5214, then the tag from this PR needs to be changed from "affect release" to "affect dev", right?

juliantaylor · 2016-08-17T15:14:04Z

curious, np.nanmedian internally switches to ma.median if the reduced axis is too small. The threshold is 400 elements and used to work well, now the threshold is around 1200 elements. Something in the other path must have regressed in performance.

crawfordsm · 2016-10-17T13:51:46Z

@keflavich Let me know if you still want to update this and have it reviewed. Looks like some tests failed and now a conflict with some other code updates.

see astropy#5232 (review)

keflavich · 2017-04-26T16:42:09Z

This is ready for further review. I can't see any remaining problems with sphinx, these warnings all seem unassociated:

WARNING: The sphinx_gallery extension is not installed, so the gallery will not be built.  You will probably see additional warnings about undefined references due to this.
WARNING: FITSFixedWarning: LONPOLE2= 180.000000000 /lonpole
WARNING: FITSFixedWarning: LATPOLE2= 0.00000000000 /latpole
/Users/adam/repos/astropy/docs/index.rst:55: WARNING: toctree contains reference to nonexisting document 'generated/examples/index'
/Users/adam/repos/astropy/docs/visualization/wcsaxes/overlays.rst:163: WARNING: Exception occurred in plotting overlays-6
/Users/adam/repos/astropy/docs/visualization/wcsaxes/overlays.rst:174: WARNING: Exception occurred in plotting overlays-7
/Users/adam/repos/astropy/docs/coordinates/frames.rst:318: WARNING: undefined label: sphx_glr_generated_examples_coordinates_plot_sgr-coordinate-frame.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/coordinates/index.rst:129: WARNING: undefined label: sphx_glr_generated_examples_coordinates_plot_obs-planning.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/coordinates/index.rst:292: WARNING: undefined label: sphx_glr_generated_examples_coordinates_plot_obs-planning.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/appendix/faq.rst:283: WARNING: undefined label: sphx_glr_generated_examples_io_skip_create-large-fits.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/appendix/faq.rst:307: WARNING: undefined label: sphx_glr_generated_examples_io_create-mef.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/index.rst:273: WARNING: undefined label: sphx_glr_generated_examples_io_modify-fits-header.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/index.rst:351: WARNING: undefined label: sphx_glr_generated_examples_io_plot_fits-image.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/index.rst:466: WARNING: undefined label: sphx_glr_generated_examples_io_fits-tables.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/io/fits/index.rst:644: WARNING: undefined label: sphx_glr_generated_examples_io_create-mef.py (if the link has no caption the label must precede a section header)
/Users/adam/repos/astropy/docs/whatsnew/1.0.rst:68: WARNING: undefined label: sphx_glr_generated_examples_coordinates_plot_obs-planning.py (if the link has no caption the label must precede a section header)

see astropy#5232 (review)

keflavich · 2017-05-08T20:56:51Z

Ready for review again!

mirca · 2017-05-08T21:01:09Z

+
+    if not NUMPY_LT_1_10:
+        assert np.isnan(funcs.mad_std([1,2,3,4,5,np.nan]))
+    assert_allclose(funcs.mad_std([1,2,3,4,5,np.nan], ignore_nan=True),


nitpick on those two lines: add blank space after commas :)

saimn

Apart from the minor inline comments, I think the main point is what we discussed in #5835 (with @crawfordsm): it would be useful to extract the logic to choose which median to use in a astropy.stats.median function, and it would simplify the code.
Also something should be done is a func is provided and ignore_nan is True. Just issue a warning or raise an exception ? This could also be stated more clearly in the docstring.

saimn · 2017-05-09T10:04:05Z

+    that will use ``np.ma.median`` with nans masked out or ``np.nanmedian``
+    instead of ``np.median`` when computing the median. [#5232]
+
+- ``astropy.sphinx``


Residual from rebase ?

saimn · 2017-05-09T10:10:52Z

+            is_masked = True
            func = np.ma.median
+            if ignore_nan:
+                data = np.ma.masked_where(np.isnan(data), data)


You could use np.ma.masked_invalid instead (it is already used below).

saimn · 2017-05-09T10:20:25Z

 NUMPY_LT_1_10_4 = not minversion('numpy', '1.10.4')
 NUMPY_LT_1_11 = not minversion('numpy', '1.11.0')
 NUMPY_LT_1_12 = not minversion('numpy', '1.12')
+NUMPY_GE_1_13 = minversion('numpy', '1.13dev')


This is not used, can be removed ?

mhvk

@keflavich - looks good except that I think now that we dropped numpy 1.7 and 1.8, this can be quite a bit simplified (since nanmedian is in 1.9, contrary to the earlier comment).

But wow, what a hassle!

mhvk · 2017-05-09T12:59:24Z

        # for normal arrays should not be done (summary: np.ma.median always
        # returns an masked array even if the result should be scalar). (#4658)
-        if isinstance(a, np.ma.MaskedArray):
+        if isinstance(data, np.ma.MaskedArray):


If func is explicitly given, this whole stanza should be ignored.

Also, maybe use np.ma.isMaskedArray here?

mhvk · 2017-05-09T13:04:13Z

+                data = np.ma.masked_where(np.isnan(data), data)
+        elif ignore_nan:
+            is_masked = False
+            if NUMPY_LT_1_10:


I'm confused - help on np.nanmedian states this has been present since numpy 1.9, and I just checked that this is indeed the case (contrary to what was stated in #5232 (comment)).. Since we just dropped support for 1.7 and 1.8, it would seem this means we can get rid of this stanza!

mhvk · 2017-05-09T13:07:08Z

+    data = np.asanyarray(data)
+    data_median = func(data, axis=axis)

    # broadcast the median array before subtraction


np.nanmedian has a keepdims argument... (it was only added in 1.10 for np.ma.median). But since any user function can be passed in, it is probably best to keep this (maybe add a comment, though).

mhvk · 2017-05-09T13:10:09Z

@@ -835,7 +885,8 @@ def mad_std(data, axis=None, func=None):
    """

    # NOTE: 1. / scipy.stats.norm.ppf(0.75) = 1.482602218505602


This is unrelated to this PR, but it would be good to mention this factor in the docstring; if I understand correctly, it essentially is the exact number for the following:

1./np.median(np.abs(np.random.normal(size=100000))) # 1.4832115241968997

@mhvk I did describe this factor in the docstring:
http://docs.astropy.org/en/stable/api/astropy.stats.mad_std.html#astropy.stats.mad_std

The factor is the normal inverse cumulative distribution function evaluated at P=3/4. It's exactly 1. / scipy.stats.norm.ppf(0.75).

@larrybradley - sorry, don't know how I could have missed this.

No worries!

mhvk · 2017-05-09T13:12:01Z

+        data[5,5] = np.nan
+        rslt = funcs.mad_std(data, ignore_nan=True)
+        assert np.isscalar(rslt)
+        with catch_warnings() as warns:


Are any warnings still given with ignore_nan=True?

I think there should be but I am admittedly fuzzy-headed at the moment.

mhvk · 2017-05-09T13:15:42Z


-__all__ = ['NUMPY_LT_1_9_1', 'NUMPY_LT_1_10', 'NUMPY_LT_1_10_4',
-           'NUMPY_LT_1_11', 'NUMPY_LT_1_12']
+__all__ = ['NUMPY_LT_1_9_1', 'NUMPY_LT_1_10',


Just leave the line as is? Then this file does not need to be touched at all.

Probably right; I think this was from a failed merge during rebasing.

see astropy#5232 (review)

larrybradley · 2017-05-09T15:03:09Z

+from astropy.utils.compat import NUMPY_LT_1_10, NUMPY_LT_1_12
+from astropy.utils.decorators import deprecated_renamed_argument
+
+


Should just be one blank line.

larrybradley · 2017-05-09T16:18:07Z

    a multi-dimensional ``axis``. [#5835]

- ``astropy.sphinx``
-


@keflavich You inadvertently deleted these lines instead of the ones below (line 551).

larrybradley · 2017-05-09T16:37:18Z

+from warnings import warn
+
+from astropy.utils.compat import NUMPY_LT_1_10
+from astropy.utils.decorators import deprecated_renamed_argument


I would write these as

from ..utils.compat import NUMPY_LT_1_10 from ..utils.decorators import deprecated_renamed_argument

and then remove the follow two blank lines.

keflavich · 2017-05-09T20:21:03Z

I... should probably squash all this if we ever decide to merge it. I think the commit log is now longer than the changes?

crawfordsm · 2017-05-10T10:11:08Z

I'm happy to merge this after any additional comments. Thanks @keflavich for taking all the time to push this through.

@bsipocz

make mad_std nan-compatible mad_std too add changelog entry [ci skip] np->numpy in docstrings & check numpy version (using 1.9, because officially it was added to numpy then according to their docs...) use astropy.utils.compat (thanks @bsipocz) 1_9, not 1P9 bump to 1.10 to avoid np 1.9 weirdness attempt to work around older numpy versions using ma.median add tests for scalarness and use np.ma.median except when axis=None. Also, change 'a' to 'data'. 'a' is a bad variable name and led to a few mistakes add another test of madstd with nan & axes add a warning about NaNs in arrays for np<1.10. Remove a test for those versions. update docstrings & changelog. Add deprecation of old argument correct the code to match the comment change order of bool checking to do expensive one last Test appropriate numpy versions, check warnings, and use asked array comparison switch to catch_warnigns for numpy 1.13, treatment of NaNs in masked arrays changed add GE_13 to __ALL__ xfail the masked array test to avoid supporting incorrect behavior see astropy#5232 (review) udpate the special cases as per @mwcraig's helpful chart! [ci skip] fix the tests: the behavior of the function when given masked arrays vs non-masked arrays has changed: we're now forcing the correct behavior for earlier versions implement @juliantaylor's suggetions test should no longer use np.ma.allclose (though I don't know why ma wouldn't just work....) this commit will fail: I've fixed the underlying issue but kept the tests flexible using the wrong data type... tests for array type double backticks a second instance of bacticks had a case wrong: one of the warn cases should *not* return a NaN for np<1.11 (which is why we're catching a warning there) numpy 1.10 np.ma.median([1,2,nan,4,5]) returns nan by default no more old numpy supports. Mixed feelings about this fix an import cleanup to address the numpy <=1.8 deprecation more cleanup fix changelog add a note about keepdims whitespace fix fix imports trailing whitespace?!!?!?!

@bsipocz

(this is a squash of 38 commits with messages below:) make mad_std nan-compatible mad_std too add changelog entry np->numpy in docstrings & check numpy version (using 1.9, because officially it was added to numpy then according to their docs...) use astropy.utils.compat (thanks @bsipocz) 1_9, not 1P9 bump to 1.10 to avoid np 1.9 weirdness attempt to work around older numpy versions using ma.median add tests for scalarness and use np.ma.median except when axis=None. Also, change 'a' to 'data'. 'a' is a bad variable name and led to a few mistakes add another test of madstd with nan & axes add a warning about NaNs in arrays for np<1.10. Remove a test for those versions. update docstrings & changelog. Add deprecation of old argument correct the code to match the comment change order of bool checking to do expensive one last Test appropriate numpy versions, check warnings, and use asked array comparison switch to catch_warnigns for numpy 1.13, treatment of NaNs in masked arrays changed add GE_13 to __ALL__ xfail the masked array test to avoid supporting incorrect behavior see astropy#5232 (review) udpate the special cases as per @mwcraig's helpful chart! fix the tests: the behavior of the function when given masked arrays vs non-masked arrays has changed: we're now forcing the correct behavior for earlier versions implement @juliantaylor's suggetions test should no longer use np.ma.allclose (though I don't know why ma wouldn't just work....) this commit will fail: I've fixed the underlying issue but kept the tests flexible using the wrong data type... tests for array type double backticks a second instance of bacticks had a case wrong: one of the warn cases should *not* return a NaN for np<1.11 (which is why we're catching a warning there) numpy 1.10 np.ma.median([1,2,nan,4,5]) returns nan by default no more old numpy supports. Mixed feelings about this fix an import cleanup to address the numpy <=1.8 deprecation more cleanup fix changelog add a note about keepdims whitespace fix fix imports trailing whitespace?!!?!?!

keflavich · 2017-05-15T18:19:13Z

Rebased and squashed.

bsipocz · 2017-05-16T10:51:33Z

@crawfordsm - Could you clarify, is this ready to be merged or we're still waiting for others to comment?

crawfordsm · 2017-05-16T12:15:21Z

I was primarily just waiting on the rebase and squashing. Merged.

Thanks @keflavich !

bsipocz reviewed Aug 11, 2016
View reviewed changes

pllim added stats Affects-release labels Aug 11, 2016

pllim added this to the v1.3.0 milestone Aug 11, 2016

astrofrog reviewed Aug 11, 2016
View reviewed changes

MSeifert04 reviewed Aug 12, 2016
View reviewed changes

MSeifert04 mentioned this pull request Feb 23, 2017

Allow to choose median function for mad_std and median_absolute_deviation #5835

Merged

keflavich force-pushed the madstdnan branch from 6643478 to 86314c3 Compare April 25, 2017 23:51

keflavich added a commit to keflavich/astropy that referenced this pull request Apr 25, 2017

xfail the masked array test to avoid supporting incorrect behavior

aa22a8b

see astropy#5232 (review)

keflavich force-pushed the madstdnan branch from 12fe88a to 3c8f699 Compare May 8, 2017 20:56

keflavich added a commit to keflavich/astropy that referenced this pull request May 8, 2017

xfail the masked array test to avoid supporting incorrect behavior

8392f54

see astropy#5232 (review)

mirca reviewed May 8, 2017

View reviewed changes

bsipocz requested review from larrybradley, mhvk and saimn May 9, 2017 09:57

saimn reviewed May 9, 2017

View reviewed changes

mhvk reviewed May 9, 2017

View reviewed changes

keflavich force-pushed the madstdnan branch from 1f7b578 to 55c43fb Compare May 9, 2017 15:12

keflavich added a commit to keflavich/astropy that referenced this pull request May 9, 2017

xfail the masked array test to avoid supporting incorrect behavior

3d5221d

see astropy#5232 (review)

larrybradley reviewed May 9, 2017

View reviewed changes

keflavich force-pushed the madstdnan branch from 62e512c to aacf5f2 Compare May 15, 2017 18:17

keflavich force-pushed the madstdnan branch from aacf5f2 to e96efa7 Compare May 15, 2017 18:19

crawfordsm merged commit 1149f0b into astropy:master May 16, 2017

keflavich deleted the madstdnan branch May 16, 2017 12:37

pllim mentioned this pull request Mar 8, 2018

Remove xpass cases and turn xpass into failures #7267

Merged

pllim mentioned this pull request Jun 20, 2019

Remove deprecated functionality for 4.0 release #8888

Closed

16 tasks done



		def median_absolute_deviation(a, axis=None):
		def median_absolute_deviation(data, axis=None, ignore_nan=False):

		@@ -835,7 +885,8 @@ def mad_std(data, axis=None, func=None):
		"""

		# NOTE: 1. / scipy.stats.norm.ppf(0.75) = 1.482602218505602

		from astropy.utils.compat import NUMPY_LT_1_10, NUMPY_LT_1_12
		from astropy.utils.decorators import deprecated_renamed_argument

Conversation

keflavich commented Aug 11, 2016

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bsipocz commented Aug 11, 2016

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pllim commented Aug 11, 2016

Uh oh!

keflavich commented Aug 11, 2016

Uh oh!

astrofrog commented Aug 11, 2016

Uh oh!

keflavich commented Aug 11, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

astrofrog commented Aug 11, 2016

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MSeifert04 commented Aug 11, 2016 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

larrybradley commented Aug 11, 2016

Uh oh!

keflavich commented Aug 11, 2016

Uh oh!

larrybradley commented Aug 11, 2016

Uh oh!

keflavich commented Aug 11, 2016

Uh oh!

larrybradley commented Aug 11, 2016

Uh oh!

MSeifert04 commented Aug 11, 2016

Uh oh!

larrybradley commented Aug 11, 2016

Uh oh!

crawfordsm commented Aug 11, 2016

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MSeifert04 Aug 13, 2016 • 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!

Choose a reason for hiding this comment

Uh oh!

juliantaylor commented Aug 17, 2016

Uh oh!

crawfordsm commented Oct 17, 2016

Uh oh!

keflavich commented Apr 26, 2017

Uh oh!

keflavich commented May 8, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

saimn left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

keflavich commented Aug 11, 2016 •

edited

Loading

MSeifert04 commented Aug 11, 2016 •

edited

Loading

MSeifert04 Aug 13, 2016 •

edited

Loading

keflavich commented May 15, 2017 •

edited

Loading