Merge pull request #5203 from njsmith/master

BUG: copy inherited masks in MaskedArray.__array_finalize__

4 files changed, 62 insertions, 6 deletions

diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index a8c63dff1..35195338f 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -1037,7 +1037,8 @@ def gradient(f, *varargs, **kwargs):
             out[slice1] = (3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0
 
         # divide by step size
-        outvals.append(out / dx[axis])
+        out /= dx[axis]
+        outvals.append(out)
 
         # reset the slice object in this dimension to ":"
         slice1[axis] = slice(None)
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 063281feb..eb0b04057 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -526,8 +526,18 @@ class TestGradient(TestCase):
 
     def test_masked(self):
         # Make sure that gradient supports subclasses like masked arrays
-        x = np.ma.array([[1, 1], [3, 4]])
-        assert_equal(type(gradient(x)[0]), type(x))
+        x = np.ma.array([[1, 1], [3, 4]],
+                        mask=[[False, False], [False, False]])
+        out = gradient(x)[0]
+        assert_equal(type(out), type(x))
+        # And make sure that the output and input don't have aliased mask
+        # arrays
+        assert_(x.mask is not out.mask)
+        # Also check that edge_order=2 doesn't alter the original mask
+        x2 = np.ma.arange(5)
+        x2[2] = np.ma.masked
+        np.gradient(x2, edge_order=2)
+        assert_array_equal(x2.mask, [False, False, True, False, False])
 
     def test_datetime64(self):
         # Make sure gradient() can handle special types like datetime64
diff --git a/numpy/ma/core.py b/numpy/ma/core.py
index e0d81419b..501177bad 100644
--- a/numpy/ma/core.py
+++ b/numpy/ma/core.py
@@ -2783,12 +2783,50 @@ class MaskedArray(ndarray):
         """
         # Get main attributes .........
         self._update_from(obj)
+        # We have to decide how to initialize self.mask, based on
+        # obj.mask. This is very difficult.  There might be some
+        # correspondence between the elements in the array we are being
+        # created from (= obj) and us. Or... there might not. This method can
+        # be called in all kinds of places for all kinds of reasons -- could
+        # be empty_like, could be slicing, could be a ufunc, could be a view,
+        # ... The numpy subclassing interface simply doesn't give us any way
+        # to know, which means that at best this method will be based on
+        # guesswork and heuristics. To make things worse, there isn't even any
+        # clear consensus about what the desired behavior is. For instance,
+        # most users think that np.empty_like(marr) -- which goes via this
+        # method -- should return a masked array with an empty mask (see
+        # gh-3404 and linked discussions), but others disagree, and they have
+        # existing code which depends on empty_like returning an array that
+        # matches the input mask.
+        #
+        # Historically our algorithm was: if the template object mask had the
+        # same *number of elements* as us, then we used *it's mask object
+        # itself* as our mask, so that writes to us would also write to the
+        # original array. This is horribly broken in multiple ways.
+        #
+        # Now what we do instead is, if the template object mask has the same
+        # number of elements as us, and we do not have the same base pointer
+        # as the template object (b/c views like arr[...] should keep the same
+        # mask), then we make a copy of the template object mask and use
+        # that. This is also horribly broken but somewhat less so. Maybe.
         if isinstance(obj, ndarray):
-            odtype = obj.dtype
-            if odtype.names:
-                _mask = getattr(obj, '_mask', make_mask_none(obj.shape, odtype))
+            # XX: This looks like a bug -- shouldn't it check self.dtype
+            # instead?
+            if obj.dtype.names:
+                _mask = getattr(obj, '_mask',
+                                make_mask_none(obj.shape, obj.dtype))
             else:
                 _mask = getattr(obj, '_mask', nomask)
+            # If self and obj point to exactly the same data, then probably
+            # self is a simple view of obj (e.g., self = obj[...]), so they
+            # should share the same mask. (This isn't 100% reliable, e.g. self
+            # could be the first row of obj, or have strange strides, but as a
+            # heuristic it's not bad.) In all other cases, we make a copy of
+            # the mask, so that future modifications to 'self' do not end up
+            # side-effecting 'obj' as well.
+            if (obj.__array_interface__["data"][0]
+                    != self.__array_interface__["data"][0]):
+                _mask = _mask.copy()
         else:
             _mask = nomask
         self._mask = _mask
diff --git a/numpy/ma/tests/test_core.py b/numpy/ma/tests/test_core.py
index 34951875d..4ac3465aa 100644
--- a/numpy/ma/tests/test_core.py
+++ b/numpy/ma/tests/test_core.py
@@ -2226,6 +2226,13 @@ class TestMaskedArrayMethods(TestCase):
         assert_equal(b.shape, a.shape)
         assert_equal(b.fill_value, a.fill_value)
 
+        # check empty_like mask handling
+        a = masked_array([1, 2, 3], mask=[False, True, False])
+        b = empty_like(a)
+        assert_(not np.may_share_memory(a.mask, b.mask))
+        b = a.view(masked_array)
+        assert_(np.may_share_memory(a.mask, b.mask))
+
     def test_put(self):
         # Tests put.
         d = arange(5)