Repository: dask/dask

Issue:
nan vs mask
I have a use case where I want to perform some statistical processing over a large data set. So, naturally I want to leverage all the awesome goodness that's baked into dask ... however, my data has an MDI.

I'm glad to see that there is nan support already in dask.array, and as a strategy that helps me move forward somewhat - although things are a wee bit clunky in the case of integral data types.

Normally, I'd turn to masked arrays, but that's clearly not an option currently within dask. I've seen peppered requests over time from the community for masked support in dask, but it seems to be a big ask, otherwise I guess it would've been already implemented.

Would it be possible to have some guidance/clarity from the dask devs as to whether such a feature is realistically and technically viable? and just as importantly, is it really desireable from your point of view?

It would be most helpful to understand if there is a genuine appetite to have masked array support within dask, or whether the effort/issues involved simply outweigh the benefit, rendering such a feature request not feasible/worthwhile.

Many thanks in advance

Discussion:
mrocklin: @shoyer @pelson @niallrobinson @rabernat @pwolfram this may interest you or others within your respective organizations
mrocklin: The special-casing of `np.ma` in dask here is somewhat concerning.  This seems fine if masked arrays are a priority (as those cc'ed above may be able to support) but puts us down an odd path where we might have a lot of highly branching code to support a few in-memory arrays.  Are there ways around this through dispatching and such?
rabernat: Looks very cool! However, it will be of limited use to xarray, since xarray does not use masked arrays internally. (It uses NaN to represent masked elements.)
jcrist: This PR is falling out of sync with master. Before I take the time to fix the merge conflicts, is this of actual use to anyone? It would be good to hear from @shoyer, @pelson, @niallrobinson, @pwolfram, or @bjlittle on whether this would be useful for any of the work y'all do.

FWIW I think the maintenance costs here are fairly minimal, but it'd be good to have a real world use case before actually merging.
dkillick: @jcrist Iris is interested in this work! Thanks muchly for the effort you've put into it -- hopefully @bjlittle or @pelson will be able to respond with something more concrete soon :+1: 
mrocklin: How does Iris handle masked arrays today?

Also @njsmith, do you have a sense for the planned longevity of the `numpy.ma` module?
bjlittle: @mrocklin At the moment, `iris` will either use `numpy.ma` for the concrete case, or `biggus` to handle masked arrays in a lazy way. But change is afoot ...

We're in the process of replacing `biggus` with `dask`, and opting for use `numpy.nan` to fudge support for lazy, mask like behaviour ... with pimple on the princess is the lazy, masked integral/bool case, where we need to change `dtype` to use `numpy.nan`. To be honest, if `numpy.nan` was `dtype` agnostic, then that would make life a heck of a lot easier, but instead to support that corner case, we need to deal with it explicitly, which isn't particularly pleasant.
mrocklin: How would life change if dask.array supported masked arrays?  Would this have near-term positive impact on Iris and the Met office?
bjlittle: IMO we're committed to cutting the next release of `iris` with the current `numpy.nan` approach and `dask` as a replacement for `biggus`. If we don't do that, well then, that's news to me. So we've totally bought into the benefits of `dask`, and commited to it as a means to underpin our needs for deferred loading and lazy evaluation.

Naturally, from our next release, users of `iris` will be exposed to `dask`. If they choose to use `dask` natively (without `iris` getting in the way) then they may also be exposed to the `numpy.nan` approach to masking ... and perhaps the lazy, masked integral/bool case stings them, forcing them to deal with that explicitly. I really can't quantify the cost of that to them or others in the community.

So, for me, the near-term positive impact for `dask.array` supporting masked arrays is that it provides a minimal overhead to us in supporting the lazy, masked integral/bool case, which has been the biggest bane of our endeavor to replace `biggus` with `dask` in `iris`, and we would also be able to pass-off all masked array handling naturally to `dask` and not jump through the `numpy.nan` hoop .... an issue that others i.e. `xarray` at least, also have to address. So, I guess the question is, should all users of `dask` accept `numpy.nan` as a given approach to dealing with masking, or can the problem be addressed in one place, but not a the cost of compromising the design of `dask`?
jakirkham: Wanted to add something here. From past experience there are many functions operating on NumPy's masked arrays that have bugs, strange behavior, or are incomplete somehow. Feel free to take a look at NumPy's issue tracker for examples. A relevant search of these is linked below.

So taking on Masked Array support means handling these sorts of cases somehow or at a bare minimum directing feedback to upstream. This isn't an argument against it. Just trying to make sure you are aware of these problems. Admittedly people that are already using `numpy.ma` are probably aware of these shortcomings.

ref: https://github.com/numpy/numpy/issues?utf8=%E2%9C%93&q=is%3Aissue%20is%3Aopen%20label%3A%22component%3A%20numpy.ma%22%20 

> ...do you have a sense for the planned longevity of the numpy.ma module?

While I can't speak for Nathaniel, my understanding is Matplotlib makes use of Masked Arrays throughout. So NumPy couldn't drop Masked Arrays without breaking Matplotlib in a pretty big way, which seems like incentive enough to leave it alone. This in spite of occasional rumblings on the NumPy issue tracker to the contrary.
shoyer: NumPy tries *very* hard to avoid breaking backwards compatibility. You can
count on masked arrays being around for the indefinite future, despite all
of their issues. We haven't even gotten rid of np.matrix, even though all
the NumPy developers agree that it shouldn't be used for new code.

On Thu, Jun 1, 2017 at 8:21 AM, jakirkham <notifications@github.com> wrote:

> Wanted to add something here. From past experience there are many
> functions operating on NumPy's masked arrays that have bugs, strange
> behavior, or are incomplete somehow. Feel free to take a look at NumPy's
> issue tracker for examples. A relevant search of these is linked below.
>
> So taking on Masked Array support means handling these sorts of cases
> somehow or at a bare minimum directing feedback to upstream. This isn't an
> argument against it. Just trying to make sure you are aware of these
> problems. Admittedly people that are already using numpy.ma are probably
> aware of these shortcomings.
>
> ref: https://github.com/numpy/numpy/issues?utf8=%E2%9C%93&q=
> is%3Aissue%20is%3Aopen%20label%3A%22component%3A%20numpy.ma%22%20
>
> ...do you have a sense for the planned longevity of the numpy.ma module?
>
> While I can't speak for Nathaniel, my understanding is Matplotlib makes
> use of Masked Arrays throughout. So NumPy couldn't drop Masked Arrays
> without breaking Matplotlib in a pretty big way, which seems like incentive
> enough to leave it alone. This in spite of occasional rumblings on the
> NumPy issue tracker to the contrary.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/dask/dask/pull/2301#issuecomment-305527352>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/ABKS1qvIJDwQWcIcKWfQ9mVukjvQzyvZks5r_tcMgaJpZM4NRcMo>
> .
>

mrocklin: Perhaps a more optimistic question is "is there likely to be a masked array alternative in the near future?"
njsmith: I would very much like to see proper NA support in numpy. That definitely doesn't replace all use cases for `np.ma` (in particular the cases where people want to toggle elements back and forth between masked/unmasked without losing data), but it sounds like it's what `iris` really wants. I won't try to make any concrete estimate for when that will happen, though. In general, I'd strongly recommend against writing new code using `np.ma` just because the design is kind of inherently broken in several ways and will never be fixed. From the maintenance point of view we generally treat it like a museum piece that needs preservation, though some people do submit the occasional bugfix. As `__numpy_ufunc__` and any future follow-ups land then it will become more possible to write a better version of what `np.ma` does – perhaps someone will take up the challenge. In the mean time, @shoyer is right that it's not going away – at most we might eventually have a better replacement, and might then eventually deprecate it, and it's just barely possible that some time after that we might eventually split it out into a new package and make people update their dependencies. But the code isn't going to disappear.
bjlittle: I'm in agreement with @shoyer and @njsmith 

Our whole `biggus` to `dask` migration and dealing with masking explicitly through `numpy.nan` has really highlighted some not very pleasant areas within `numpy` masked arrays. There are certainly a few dragons living there.

I'm pretty much in awe of @jcrist's efforts to get this PR up :+1:. Totally awesome :beers: 

That said, this PR has forced the discussion on whether `dask` *should* support `numpy` masked arrays. Being impartial, I think it's a decision that shouldn't be taken lightly by the `dask` core devs. I see it as a pretty big, long term commitment - caveat emptor.

jcrist: The changes here were mostly:
- Remove the `package_of` function in favor of a method registry (still think this is a good idea)
- Generalize a few small bits
- Some masked array specific changes that I'm not sure how/if to generalize
- Tests

Of these, I'm only really in favor of adding the first separate from masked-array support. The rest would clutter up the code a bit without a specific reason behind them. I'm fine with letting this idle (even closed) until masked arrays are requested - shouldn't take long to bring back up to date.
marqh: > Wait on Iris and XArray to tell us that they would prefer to use masked arrays if available

Hello @mrocklin et al

I'd like to share my perspective, as an Iris developer.

I think that the adoption of Dask within Iris is a really positive and exciting step for us.  The work to integrate dask into Iris has been fairly intensive, and working around the lack of masked arrays has formed a significant quantity of this work.
Some of the work around code has been accepted as pragmatic, but does not feel ideal to me.
There are concerns from our user community about some of the side effects of this implementation for some use cases.

I think that in the short term, we will be working with dask and patching around the lack of masked array support.  However, I feel that this represents a degree of technical debt for our implementation.

So, in the near future, I'm interested in dask supporting numpy's masked array implementation, as it feels like the widely used implementation and many libraries using numpy rely on this implementation for their core functionality. 
If a future dask supports `numpy.ma` I would be very keen to adopt that in short order for our library, this would be a great benefit, I feel.

I think that driving a conversation about how a future numpy and dask could handle the concept of missing data and processing in different ways is a really valuable conversation.  I think that this is a long term activity, which would be informed by dask adopting numpy's masked array as is.

I'd like to offer my encouragement for dask to adopt `numpy.ma`, with all its foibles and use this activity to explore what could be done better.  Perhaps this would encourage a conversation across the numpy community about having a fresh look at this topic; I'd certainly be interested in engaging with such a discussion and offering key use cases if I can generate them.

many thanks
mark 
jcrist: Ok, I've brought this PR back in sync with master.

I don't think this adds much complexity to `dask.array` internals, and think it'd be fine to merge. However, it might be nice to try it out on some simple problem to see if it solves real world needs. Specifically, the `fill_value` is not metadata on `dask.array`, and so can't be accessed statically (e.g. can't do `a.fill_value`). If this is needed for your work, then this PR doesn't satisfy that need. We could add a function to get it lazily (`ma.get_fill_value(x)`?), but having it as static metadata would require a redesign.

@marqh, @bjlittle do you have examples of what kind of operations you'd like to do with masked arrays? Do you have time/motivation to try this PR out?
pelson: @jcrist - great stuff! 👍 
jakirkham: Did you try with `nomask`?
jcrist: Thanks for the review @pelson. Responding to a few comments:

---

The errors you're seeing in `masked_where` and `masked_outside` have to do with lack of support for numpy arrays as arguments. Will fix.

One question is what do you expect to happen if you call a `da.ma.masked_*` function on *only* numpy arrays (as you did with `masked_outside`). Dask's convention so far is to call the equivalent numpy method in these cases. There's a question whether this should change (see https://github.com/dask/dask/issues/2103).

---

The `from_array` issue is a bug in our tokenizing function (will fix), but even then `from_array` won't work as it always calls `np.asarray` on all chunks. A few options here:

- Don't call `asarray` if the input is a masked array. I'm not as fond of this because it special cases an input type.

- Add a keyword to turn off `asarray`. NumPy seems to use `subok` for this in places, but I'm tempted to use `asarray=True` instead because I think it's clearer.

Thoughts?

---

> I was somewhat concerned about the numerical accuracy of things like standard deviation

We also work hard to implement numerically stable parallel algorithms. The one used here is from http://prod.sandia.gov/techlib/access-control.cgi/2008/086212.pdf and has been working well for us. It's used for computing all our moments (`var`, `skew`, `kurt`, etc...).

> Which, just like the biggus implementation has a (necessary) bottleneck at the sqrt

Just to clarify, depending on your chunking and reduction axis there may be multiple parallel calls to `sqrt`. In this case you're reducing down to a single chunk so there's only one call.
pelson: > Dask's convention so far is to call the equivalent numpy method in these cases. 

Personally, I'd expect to *always* have a lazy dask thing if I've called dask functions/methods. If I wanted immediate, I'd call numpy directly...

> bug in our tokenizing function (will fix)... Thoughts?

At this moment, I'm afraid I don't have enough background to possibly give an educated answer.

> We also work hard to implement numerically stable parallel algorithms.

Great. Biggus used a method described in ```Welford, BP (August 1962). "Note on a Method for Calculating Corrected Sums of Squares and Products``` and needed a particular implementation for masked arrays, so I think it is a good sign to see that your implementation of ``moment_chunk`` has held up well with this extension. 👍 
marqh: @jcrist many thanks for the updates

> @marqh, @bjlittle do you have examples of what kind of operations you'd like to do with masked arrays? 

key operations I know about with masks include:
 * statistical aggregations that are mask aware
 * statistical aggregations of float arrays using masks and NaNs with different semantics
 * masked integer arrays
 * conditional masks, calculated differently based on metadata
 * plotting of masked data

Are you interested in general problem statements like this, or are you more keen for actual examples of data processing which could form test cases?

Our current implementation converts masked arrays to NaN filled float arrays, which gives partial support, but there are edge cases which are causing concerns and lacks of capability for some cases. 
  
> Do you have time/motivation to try this PR out?

motivation, for sure, time is more of a challenge just now, i'll let you know if I make progress.
jcrist: I'm running into some hard-to-solve issues with masked arrays if the `fill_value` isn't a scalar. In these cases it becomes unclear what the intended behavior is when concatenating blocks, and is tricky enough that I'd like to just forbid it/not handle it. It seems to only come up in `masked_equals` and `masked_values` (there's a few numpy issues about how this is confusing and bad but needs to remain for backwards compat).

I'd like to not handle this case and only work with scalar `fill_value` - is that ok for your needs? I can make it work with arrays, but there'd probably be edge cases. If I can forbid it, what should the behavior be if an array fill_value shows up (say from some user function and `map_blocks`)?

- Error explicitly with "unsupported, can't be consistent here with numpy". This would be hard to catch in all cases, we'd really only see it when merging blocks.
- Silently reset `fill_value` to the default. This is what `np.ma.concatenate` does anyway, but seems like it could be unexpected.

If given no direction, I plan to error explicitly where easy/possible, and forbid operations that could result in this case (e.g. forbid `masked_equal` with an array `value`).
jcrist: > If given no direction, I plan to error explicitly where easy/possible, and forbid operations that could result in this case (e.g. forbid masked_equal with an array value).

This is what I ended up doing.

---

At this point this PR *may* be ready to go. Three things that may still need to be done:

- Improve the test suite a bit. I think we're good, but I want to do a once-over one more time with fresh eyes tomorrow.
- Add masked ufuncs? Currently we don't implement things like `np.ma.sin`. The top level ufuncs (e.g. `da.sin`) do work on masked arrays, *but* they are applied to the entire `.data` array, which the masked ufuncs avoid. I'm not sure how important this is in practice. The non-masked results are equivalent, but the data under the mask may be different.
- Add some docs? I'd prefer to wait on this until it sees some use, as we may want to change things after it's tried out. Happy to add docs now though if desired.
jcrist: Ok, I'm happy with this as is. I expect there will be future touch-ups as this feature gets use, but I think it'd be good to get it in so it's easier for others to try out. This could use a review if anyone has time (cc @mrocklin, @pelson, @bjlittle, or @marqh).
mrocklin: Just planning things out here, @pelson, @bjlittle, or @marqh do any of you have time to test this against Iris?
marqh: > I'd like to not handle this case and only work with scalar fill_value - is that ok for your needs? I can make it work with arrays, but there'd probably be edge cases. If I can forbid it, what should the behavior be if an array fill_value shows up (say from some user function and map_blocks)?
> 
>     Error explicitly with "unsupported, can't be consistent here with numpy". This would be hard to catch in all cases, we'd really only see it when merging blocks.
>     Silently reset fill_value to the default. This is what np.ma.concatenate does anyway, but seems like it could be unexpected.

@jcrist 

in my view, single fill values are the only case I can see any reason to support

I very much support failing with an exception if a case occurs where this is not the case, user of library intervention would be required at this stage

marqh: > Just planning things out here, @pelson, @bjlittle, or @marqh do any of you have time to test this against Iris?

@mrocklin I'd like to set up a testing framework to look at this

@bjlittle @pelson may I suggest that we create an iris branch somewhere and build dask from this branch in situ?  We can then apply changes to this iris branch to remove the work arounds for handling masked arrays and get the tests running in a place we can all see the results 

Would it make sense to have such a branch inside the scitools/iris project, where we all have eyes and permissions to change?
marqh: @pp-mo for your interest and information
lbdreyer: I have just created the iris branch [dask_mask_array](https://github.com/SciTools/iris/tree/dask_mask_array) that I will be targeting with my initial investigation into integrating this branch into Iris.


jakirkham: Did `compress` and/or `nonzero` get defined for Masked Arrays?
lbdreyer: Hi @jcrist 
We're sprinting this week with the goal of having a POC of using your branch in Iris. We'll feedback any questions we have.

We have come across one thing we wanted to check with you regarding fill value.
The following is using your most recent commit (d504e35)
```
>>> import dask.array.ma as dma
>>> 
>>> x = ma.masked_array([1,2,3,4], mask=[1,0,0,0])
>>> dx = da.from_array(x, chunks=(2,), asarray=False)
>>> dma.set_fill_value(dx, 8)
>>> print dx.compute().fill_value
8
```

However, if we call compute, then set the fill value, the next time we call compute the array fill value hasn't been changed, e.g.
```
>>> x = ma.masked_array([1,2,3,4], mask=[1,0,0,0])
>>> dx = da.from_array(x, chunks=(2,), asarray=False)
>>> print dx.compute().fill_value
999999
>>> 
>>> dma.set_fill_value(dx, 8)
>>> 
>>> print dx.compute().fill_value
999999
```

Is this as you would expect?
jcrist: Thanks for the issue report, this was a tricky bug that affected other parts of `dask.array` too. Fixed in #2572.

```python
In [1]: import numpy as np

In [2]: import dask.array as da

In [3]: x = np.ma.masked_array([1, 2, 3, 4], mask=[1, 0, 0, 0])

In [4]: dx = da.from_array(x, chunks=2, asarray=False)

In [5]: dx.compute().fill_value
Out[5]: 999999

In [6]: da.ma.set_fill_value(dx, 8)

In [7]: dx.compute().fill_value
Out[7]: 8
```
bjlittle: Hi @jcrist,

We just stumbled across the following behaviour ... apologies, it's your old friend `fill_value` again:
```python
>>> a = da.from_array(np.array([1, 2, 3]), chunks=(1,))
>>> from dask.array import ma as dma
>>> dma.set_fill_value(a, 123)
>>> a.compute()
array([1, 2, 3])
```
Silently ignores the `set_fill_value` on an `ndarray` ... seems reasonable, since we can't actually get a lazy `fill_value` anyways.

This is just an observation really, not quite sure what your expectation is on how this type of use pattern should behave ...

jcrist: I'm not sure what your expectation is here? In numpy, `set_fill_value` is a no-op on non-masked arrays:

```python
In [1]: import numpy as np

In [2]: x = np.array([1, 2, 3])

In [3]: np.ma.set_fill_value(x, 0.5)

In [4]: x
Out[4]: array([1, 2, 3])
```

We just do the same here. What were you expecting to happen?
bjlittle: Oh boy, numpy lets you set the value of a non-existent attribute and silently ignores the fact, of course it does ... I should've checked that first, apologies @jcrist. You're just honouring that behaviour, which if ideal.

I guess my beef/misunderstanding is with numpy. No worries.

btw @lbdreyer and I are actively re-engineering iris (in a feature branch) at this moment to use this PR #2301. If we discover any issues we'll be sure to let you know asap. We're hoping to be in a position by the end of this week to have a fully functional iris with all unit and integration tests passing. At that point, I guess we'll be keen to know "what's the blocker to #2301 being merged?" and "when is your next release?"
mrocklin: > btw @lbdreyer and I are actively re-engineering iris (in a feature branch) at this moment to use this PR #2301. If we discover any issues we'll be sure to let you know asap. We're hoping to be in a position by the end of this week to have a fully functional iris with all unit and integration tests passing. At that point, I guess we'll be keen to know "what's the blocker to #2301 being merged?" and "when is your next release scheduled?"

@bjlittle my understanding is that we're just waiting on you all for feedback.  I don't think that there are any blockers to merging this.  We can probably release within a week or two after that.
pp-mo: Hi @mrocklin @jcrist + people

Work is now progressing with our Iris "dask_mask_array" branch, using this dask to re-instate sensible masked operations.

All this is good, but there are a couple of points currently causing some problems ...

**One problem**
is that we don't have the handling of fill_value that we'd like, in that ...
  * (a) it is not an available thing on lazy arrays (i.e. without calling `compute()`), and 
  * (b) it is not a fixed value that applies to the whole of the array (e.g. `lazy[0].compute().fill_value` can be different from `lazy[0].compute().fill_value`)
We are in any case struggling right now to suggest what might work better, so I'm just putting that one out for your thoughts.

**A more obvious problem**
is that some basic statistical operations don't seem work yet on masked data.
For example:
```
>>> from numpy import ma
>>> import dask.array as da
>>> 
>>> masked_data = ma.masked_array([1, 2, 777, 3], mask=[0, 0, 1, 0])
>>> print 'masked mean = ', masked_data.mean()
masked mean =  2.0
>>> 
>>> for array_cast in (True, False):
...     print
...     print 'asarray={}:'.format(array_cast)
...     lazy_data = da.from_array(masked_data, chunks=masked_data.shape, asarray=array_cast)
...     print
...     print '   lazy data repr={!r}'.format(lazy_data)
...     print '   lazy computed-result={!r}'.format(type(lazy_data.compute()))
...     lazy_mean = lazy_data.mean()
...     print '   lazy mean repr={!r}'.format(lazy_mean)
...     print '   lazy mean computed-result={!r}'.format(lazy_mean.compute())
... 

asarray=True:

   lazy data repr=dask.array<array, shape=(4,), dtype=int64, chunksize=(4,)>
   lazy computed-result=<type 'numpy.ndarray'>
   lazy mean repr=dask.array<mean_agg-aggregate, shape=(), dtype=float64, chunksize=()>
   lazy mean computed-result=195.75

asarray=False:

   lazy data repr=dask.array<array, shape=(4,), dtype=int64, chunksize=(4,)>
   lazy computed-result=<class 'numpy.ma.core.MaskedArray'>
   lazy mean repr=dask.array<mean_agg-aggregate, shape=(), dtype=float64, chunksize=()>
Traceback (most recent call last):
[[ . . . ]]
  File "/data/local/itpp/git/dask/dask/array/reductions.py", line 261, in mean_agg
    return divide(pair['total'].sum(dtype=dtype, **kwargs),
  File "/opt/scitools/environments/experimental/2017_07_14/lib/python2.7/site-packages/numpy/core/_methods.py", line 32, in _sum
    return umr_sum(a, axis, dtype, out, keepdims)
ValueError: 'axis' entry is out of bounds
>>> 
```

Note that the 'asarray=True' case is not the correct result (!)

In our code, the more normal case is that the 'lazy_data' is wrapped around a "proxy" object that yields data on request (which happens to be masked) ...
The results are the same ...
```
>>> class DataProxy(object):
...     def __init__(self, shape, dtype, deferred_data_call):
...         self._fetch_call = deferred_data_call
...         self.shape = shape
...         self.dtype = dtype
...     def __getitem__(self, keys):
...         data = self._fetch_call()
...         assert data.shape == self.shape and data.dtype == self.dtype
...         return data
... 
>>> ma_proxy = DataProxy((4,), np.int64, lambda: masked_data)
>>> 
>>> for array_cast in (True, False):
...     print
...     print 'asarray={}:'.format(array_cast)
...     lazy_data = da.from_array(ma_proxy, chunks=(4,), asarray=array_cast)
...     print
...     print '   lazy data repr={!r}'.format(lazy_data)
...     print '   lazy computed-result={!r}'.format(type(lazy_data.compute()))
...     lazy_mean = lazy_data.mean()
...     print '   lazy mean repr={!r}'.format(lazy_mean)
...     print '   lazy mean computed-result={!r}'.format(lazy_mean.compute())
... 

asarray=True:

   lazy data repr=dask.array<array, shape=(4,), dtype=int64, chunksize=(4,)>
   lazy computed-result=<type 'numpy.ndarray'>
   lazy mean repr=dask.array<mean_agg-aggregate, shape=(), dtype=float64, chunksize=()>
   lazy mean computed-result=195.75

asarray=False:

   lazy data repr=dask.array<array, shape=(4,), dtype=int64, chunksize=(4,)>
   lazy computed-result=<class 'numpy.ma.core.MaskedArray'>
   lazy mean repr=dask.array<mean_agg-aggregate, shape=(), dtype=float64, chunksize=()>
Traceback (most recent call last):
[[ . . . ]]
ValueError: 'axis' entry is out of bounds
>>> 
```
jakirkham: The examples in `b` (e.g. `lazy[0].compute().fill_value`) don't appear to be different. Are they suppose to be?
jcrist: > A more obvious problem is that some basic statistical operations don't seem work yet on masked data...

A couple points here:

- When calling `from_array` on `np.ma.masked_array` objects, you *must* pass in `asarray=False`. Otherwise `np.asarray` is called on the masked array, which results in a numpy array. When you say that the "'asarray=True' case is not the correct result" this is because the mean is being taken on all the data, not just the masked part:

```python
In [1]: import numpy as np

In [2]: masked_data = np.ma.masked_array([1, 2, 777, 3], mask=[0, 0, 1, 0])

In [3]: x = np.asarray(masked_data)  # all the data, not just the masked bit

In [4]: x.mean()
Out[4]: 195.75
```

- I'm unable to produce the error you give above. What version of numpy are you using? I suspect that this is a bug in old versions of `numpy.ma`.

```python
In [5]: import dask.array as da

In [6]: x = da.from_array(masked_data, chunks=(4,), asarray=False)

In [7]: x.mean().compute()
Out[7]: 2.0

In [8]: x
Out[8]: dask.array<array, shape=(4,), dtype=int64, chunksize=(4,)>

In [9]: np.__version__
Out[9]: '1.13.1'
```

---

> One problem is that we don't have the handling of fill_value that we'd like...

This is something that I mentioned in a call with @pelson a few weeks ago. The current design doesn't keep track of the `fill_value` statically, so you have to compute something to query it. He mentioned that that might not matter for your code. If you *do* need it statically this would require a mild rewrite of the masked array support on our side (which is still doable if necessary).

If you don't need to know the `fill_value`, but do want to set it to something consistent across all chunks, you can call `da.ma.set_fill_value`. This will change the fill_value for every chunk. Note that this is fine to call on `np.ndarray`s as well, they're passed through unchanged (matching numpy).

```python
In [23]: m1 = np.ma.masked_array([1, 2, 3, 4], mask=[0, 1, 0, 1], fill_value=10)

In [24]: m2 = np.ma.masked_array([5, 6, 7, 8], mask=[0, 1, 0, 1], fill_value=20)  # different fill value

In [25]: dm1 = da.from_array(m1, chunks=2, asarray=False)

In [26]: dm2 = da.from_array(m2, chunks=2, asarray=False)

In [27]: dm3 = da.concatenate([dm1, dm2])  # each chunk has a different fill value

In [28]: dm3.compute()  # because of this, the fill_value is reset on compute, since they're not consistent
Out[28]:
masked_array(data = [1 -- 3 -- 5 -- 7 --],
             mask = [False  True False  True False  True False  True],
       fill_value = 999999)

In [29]: da.ma.set_fill_value(dm3, 25)  # if you set it though, things work fine

In [30]: dm3.compute()
Out[30]:
masked_array(data = [1 -- 3 -- 5 -- 7 --],
             mask = [False  True False  True False  True False  True],
       fill_value = 25)
```

pp-mo: > @jakirkham The examples in b (e.g. lazy[0].compute().fill_value) don't appear to be different (i.e. the code appears to be identical in both case). Are they suppose to be?

Apologies, meant to be just lazy[0].compute().fill_value and lazy[**1**].compute().fill_value
pp-mo: @jcrist thanks for your quick response.
> When calling from_array on np.ma.masked_array objects, you must pass in asarray=False. 

Thanks, I think I get that now !

> I'm unable to produce the error you give above. What version of numpy are you using?

I'm testing with 1.11.1
I'll look into testing against 1.13 instead.

> > One problem is that we don't have the handling of fill_value that we'd like...

> ... The current design doesn't keep track of the fill_value statically, so you have to compute something to query it ... If you do need it statically this would require a mild rewrite of the masked array support on our side (which is still doable if necessary).

Thanks for your views.  As I said, we still aren't sure what we *would* ideally like !
Latest thoughts are that we have some serious reservations about how this is managed in numpy.ma anyway, so we may need to go our own way to control it.  Our original [with-nonmasked-dask version of Iris](https://github.com/SciTools/iris/tree/7ec22fa3f14b1aea454c5a51ac350de190603d9c/lib/iris) was actually managing fill-values *separately* from the data arrays.




pp-mo: > I'm unable to produce the error you give above. What version of numpy are you using?

Aha, it does works properly with 1.13.1 !

pp-mo: > @jcrist it does work properly with 1.13.1 

Apologies for my panicking too soon :  Not being involved earlier, I'm still playing catch-up.

Digging further into what we need to re-factor the Iris usage, a thing has come up ...
There are as yet no dask operations allowing us to work with the data and mask components of arrays separately, for example ...
```
>>> import numpy.ma as ma
>>> masked = ma.masked_array([1, 2, 777, 3], mask=[0,0,1,0])
>>> da_mask = da.from_array(masked, chunks=masked.shape, asarray=False)
>>> 
>>> # WANTED: something like ma.getmaskarray()
>>> mask = masked.maskarray()
>>> mask.compute()
[False, False, True, False]
>>> 
>>> # WANTED: something like da.filled()
>>> filled = masked.filled(-99)
>>> filled.compute()
[1, 2, -99, 3]
>>> 
>>> # WANTED: build from separate data + mask, like da.masked_array()
>>> da_mask = da.from_array(np.array([0, 0, 0, 1]), chunks=4)
>>> remasked = da.masked_data(masked.filled(55), mask=da_mask)
>>> remasked.compute()
[1, 2, 55, --]
```

The one I have immediate use for is the "maskarray" operation :  For this we have previously been using `da.isnan()`.

**Or, maybe you can make these things for yourself ?**

Is it possible or wise to construct such operators yourself using something like `dask.array.core.elemwise` ?
I note that this is a "public" function, but it is not mentioned in the docs.

I have currently tried this,  ...
```
>>> def getmask(dask_array):
...   nd = dask_array.ndim
...   allinds = 'ijklmnopqr'[:nd]
...   result = da.atop(ma.getmaskarray, allinds, dask_array, allinds, dtype=dask_array.dtype)
...   return result
```
That seems to function, but I strongly suspect I haven't really understood what I am doing here.
Can anyone advise ?

pp-mo: > I have currently tried this, ... da.atop( ...

Update: using "elemwise" seems to work nicely.
See [here](https://github.com/SciTools/iris/pull/2743)

jcrist: Apologies for the late response, I've been away from my computer for a few days.

> There are as yet no dask operations allowing us to work with the data and mask components of arrays separately

`da.ma.getmaskarray`, `da.ma.getdata`, and `da.ma.filled` are already implemented. We don't have matching methods, because we don't have a `masked_array` class, but the function versions for those are already in the `da.ma` namespace. I'll add `da.ma.masked_array` as well. I suggest looking at what's in `da.ma` already.

> Is it possible or wise to construct such operators yourself using something like `dask.array.core.elemwise` ?

Yes. For functionality that matches those in `np.ma` though it'd be best to just add it to this PR. `elemwise`, `atop`, and `map_blocks` should all be considered public functions, feel free to use them.

pp-mo: > da.ma.getmaskarray, da.ma.getdata, and da.ma.filled 

Again, I hadn't spotted that !  Very useful.
Also, I see I should really be using `getdata`  in my code.

What about the "constructor" `ma.masked_array` method ?
I am still needing that to reconstruct my final result.
jcrist: > That would be great -- I am still needing that to reconstruct my final result.

Added in latest commit. Please let me know if you need anything else.
pp-mo: > Added in latest commit.

Making for [slightly neater code](https://github.com/SciTools/iris/pull/2743/commits/7c71140985b63d50bf3091250eaba67160198c85)

Great stuff, thanks @jcrist  !

pp-mo: Just to be _really_ clear, we like this + we want it !
We were already 100% committed to using dask anyway, but this is generally neater and much less of a wrench for existing code.
A lot of dust still to settle [here](https://github.com/SciTools/iris/projects/18), but basically 👍 👏 

mrocklin: @pp-mo @pelson we're likely to issue a small release soon.  We're now faced with the question of whether or not to include this or wait a bit for additional issues to crop up over the next week or two and release then. 

Do you have a sense for how much additional work you're likely to want on this?
jcrist: My personal preference is to do the small release and merge this after. I suspect we'll want to iterate on this a bit more, and releases are cheap to do. Up to y'all though.
pelson: Honestly, we have taken this to the point of where we are able to say with confidence that it has value and that we will be using it in iris 2.0. Additionally, we don't currently have any work planned for the next few weeks that will reveal any subsequent improvements/bugs.

Given the above, I'd love to see this in the release.
pp-mo: > I'd love to see this in the release.

👍 
mrocklin: There are likely to be a few other large array improvements in the next
release.  I'm inclined to follow @jcrist's intuition here and wait for all
of them to land before releasing a 0.X.0 release.  This would likely be in
a small number of weeks.

On Fri, Aug 25, 2017 at 1:36 PM, Patrick Peglar <notifications@github.com>
wrote:

> I'd love to see this in the release.
>
> 👍
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/dask/dask/pull/2301#issuecomment-324988108>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/AASszJB6SQClKEY-mAkKQwyiabMISEszks5sbwYIgaJpZM4NRcMo>
> .
>

mrocklin: Dask 0.15.2 has been released.. Should we merge this?  Is there additional work to do here?
pelson: > Is there additional work to do here?

We now have some confidence in saying that, even if it does transpire that there is more to do, it won't be a long way away from the functionality already provided by this PR.

Given this, I'm an advocate getting this in the hands of users ASAP.
jcrist: Fine by me. I've fixed the merge conflicts and updated some docs. I think this is good to go now.
mrocklin: +1 from me

Patch:
diff --git a/dask/array/__init__.py b/dask/array/__init__.py
index febd61182..679c66d99 100644
--- a/dask/array/__init__.py
+++ b/dask/array/__init__.py
@@ -32,6 +32,8 @@ from .reductions import (sum, prod, mean, std, var, any, all, min, max, vnorm,
 from .percentile import percentile
 with ignoring(ImportError):
     from .reductions import nanprod, nancumprod, nancumsum
+with ignoring(ImportError):
+    from . import ma
 from . import random, linalg, ghost, learn, fft
 from .wrap import ones, zeros, empty, full
 from .creation import ones_like, zeros_like, empty_like, full_like
diff --git a/dask/array/core.py b/dask/array/core.py
index 32d623d96..9cca27f99 100644
--- a/dask/array/core.py
+++ b/dask/array/core.py
@@ -33,7 +33,7 @@ from ..base import Base, tokenize, normalize_token
 from ..context import _globals
 from ..utils import (homogeneous_deepmap, ndeepmap, ignoring, concrete,
                      is_integer, IndexCallable, funcname, derived_from,
-                     SerializableLock, ensure_dict, package_of)
+                     SerializableLock, ensure_dict, Dispatch)
 from ..compatibility import unicode, long, getargspec, zip_longest, apply
 from ..delayed import to_task_dask
 from .. import threaded, core
@@ -41,6 +41,20 @@ from .. import sharedict
 from ..sharedict import ShareDict
 
 
+concatenate_lookup = Dispatch('concatenate')
+tensordot_lookup = Dispatch('tensordot')
+concatenate_lookup.register((object, np.ndarray), np.concatenate)
+tensordot_lookup.register((object, np.ndarray), np.tensordot)
+
+
+@tensordot_lookup.register_lazy('sparse')
+@concatenate_lookup.register_lazy('sparse')
+def register_sparse():
+    import sparse
+    concatenate_lookup.register(sparse.COO, sparse.concatenate)
+    tensordot_lookup.register(sparse.COO, sparse.tensordot)
+
+
 def getter(a, b, asarray=True, lock=None):
     if isinstance(b, tuple) and any(x is None for x in b):
         b2 = tuple(x for x in b if x is not None)
@@ -477,8 +491,8 @@ def _concatenate2(arrays, axes=[]):
         return arrays
     if len(axes) > 1:
         arrays = [_concatenate2(a, axes=axes[1:]) for a in arrays]
-    module = package_of(type(max(arrays, key=lambda x: x.__array_priority__))) or np
-    return module.concatenate(arrays, axis=axes[0])
+    concatenate = concatenate_lookup.dispatch(type(max(arrays, key=lambda x: x.__array_priority__)))
+    return concatenate(arrays, axis=axes[0])
 
 
 def apply_infer_dtype(func, args, kwargs, funcname, suggest_dtype=True):
@@ -2579,13 +2593,22 @@ def _enforce_dtype(*args, **kwargs):
     function = kwargs.pop('enforce_dtype_function')
 
     result = function(*args, **kwargs)
-    if dtype != object:
+    if dtype != result.dtype and dtype != object:
+        if not np.can_cast(result, dtype, casting='same_kind'):
+            raise ValueError("Inferred dtype from function %r was %r "
+                             "but got %r, which can't be cast using "
+                             "casting='same_kind'" %
+                             (funcname(function), str(dtype), str(result.dtype)))
         if np.isscalar(result):
             # scalar astype method doesn't take the keyword arguments, so
             # have to convert via 0-dimensional array and back.
-            result = result[...].astype(dtype, casting='same_kind', copy=False)[()]
+            result = result.astype(dtype)
         else:
-            result = result.astype(dtype, casting='same_kind', copy=False)
+            try:
+                result = result.astype(dtype, copy=False)
+            except TypeError:
+                # Missing copy kwarg
+                result = result.astype(dtype)
     return result
 
 
@@ -2832,8 +2855,7 @@ def concatenate3(arrays):
 
     advanced = max(core.flatten(arrays, container=(list, tuple)),
                    key=lambda x: getattr(x, '__array_priority__', 0))
-    module = package_of(type(advanced)) or np
-    if module is not np and hasattr(module, 'concatenate'):
+    if concatenate_lookup.dispatch(type(advanced)) is not np.concatenate:
         x = unpack_singleton(arrays)
         return _concatenate2(arrays, axes=list(range(x.ndim)))
 
diff --git a/dask/array/ma.py b/dask/array/ma.py
new file mode 100644
index 000000000..f656eaf97
--- /dev/null
+++ b/dask/array/ma.py
@@ -0,0 +1,251 @@
+from __future__ import absolute_import, division, print_function
+
+from functools import wraps
+from distutils.version import LooseVersion
+
+import numpy as np
+
+from ..base import normalize_token
+from .core import (concatenate_lookup, tensordot_lookup, map_blocks,
+                   asanyarray, atop)
+
+
+if LooseVersion(np.__version__) < '1.11.0':
+    raise ImportError("dask.array.ma requires numpy >= 1.11.0")
+
+
+@normalize_token.register(np.ma.masked_array)
+def normalize_masked_array(x):
+    data = normalize_token(x.data)
+    mask = normalize_token(x.mask)
+    fill_value = normalize_token(x.fill_value)
+    return (data, mask, fill_value)
+
+
+@concatenate_lookup.register(np.ma.masked_array)
+def _concatenate(arrays, axis=0):
+    out = np.ma.concatenate(arrays, axis=axis)
+    fill_values = [i.fill_value for i in arrays if hasattr(i, 'fill_value')]
+    if any(isinstance(f, np.ndarray) for f in fill_values):
+        raise ValueError("Dask doesn't support masked array's with "
+                         "non-scalar `fill_value`s")
+    if fill_values:
+        # If all the fill_values are the same copy over the fill value
+        fill_values = np.unique(fill_values)
+        if len(fill_values) == 1:
+            out.fill_value = fill_values[0]
+    return out
+
+
+@tensordot_lookup.register(np.ma.masked_array)
+def _tensordot(a, b, axes=2):
+    # Much of this is stolen from numpy/core/numeric.py::tensordot
+    # Please see license at https://github.com/numpy/numpy/blob/master/LICENSE.txt
+    try:
+        iter(axes)
+    except TypeError:
+        axes_a = list(range(-axes, 0))
+        axes_b = list(range(0, axes))
+    else:
+        axes_a, axes_b = axes
+    try:
+        na = len(axes_a)
+        axes_a = list(axes_a)
+    except TypeError:
+        axes_a = [axes_a]
+        na = 1
+    try:
+        nb = len(axes_b)
+        axes_b = list(axes_b)
+    except TypeError:
+        axes_b = [axes_b]
+        nb = 1
+
+    # a, b = asarray(a), asarray(b)  # <--- modified
+    as_ = a.shape
+    nda = a.ndim
+    bs = b.shape
+    ndb = b.ndim
+    equal = True
+    if na != nb:
+        equal = False
+    else:
+        for k in range(na):
+            if as_[axes_a[k]] != bs[axes_b[k]]:
+                equal = False
+                break
+            if axes_a[k] < 0:
+                axes_a[k] += nda
+            if axes_b[k] < 0:
+                axes_b[k] += ndb
+    if not equal:
+        raise ValueError("shape-mismatch for sum")
+
+    # Move the axes to sum over to the end of "a"
+    # and to the front of "b"
+    notin = [k for k in range(nda) if k not in axes_a]
+    newaxes_a = notin + axes_a
+    N2 = 1
+    for axis in axes_a:
+        N2 *= as_[axis]
+    newshape_a = (-1, N2)
+    olda = [as_[axis] for axis in notin]
+
+    notin = [k for k in range(ndb) if k not in axes_b]
+    newaxes_b = axes_b + notin
+    N2 = 1
+    for axis in axes_b:
+        N2 *= bs[axis]
+    newshape_b = (N2, -1)
+    oldb = [bs[axis] for axis in notin]
+
+    at = a.transpose(newaxes_a).reshape(newshape_a)
+    bt = b.transpose(newaxes_b).reshape(newshape_b)
+    res = np.ma.dot(at, bt)
+    return res.reshape(olda + oldb)
+
+
+@wraps(np.ma.filled)
+def filled(a, fill_value=None):
+    a = asanyarray(a)
+    return a.map_blocks(np.ma.filled, fill_value=fill_value)
+
+
+def _wrap_masked(f):
+
+    @wraps(f)
+    def _(a, value):
+        a = asanyarray(a)
+        value = asanyarray(value)
+        ainds = tuple(range(a.ndim))[::-1]
+        vinds = tuple(range(value.ndim))[::-1]
+        oinds = max(ainds, vinds, key=len)
+        return atop(f, oinds, a, ainds, value, vinds, dtype=a.dtype)
+
+    return _
+
+
+masked_greater = _wrap_masked(np.ma.masked_greater)
+masked_greater_equal = _wrap_masked(np.ma.masked_greater_equal)
+masked_less = _wrap_masked(np.ma.masked_less)
+masked_less_equal = _wrap_masked(np.ma.masked_less_equal)
+masked_not_equal = _wrap_masked(np.ma.masked_not_equal)
+
+
+@wraps(np.ma.masked_equal)
+def masked_equal(a, value):
+    a = asanyarray(a)
+    if getattr(value, 'shape', ()):
+        raise ValueError("da.ma.masked_equal doesn't support array `value`s")
+    inds = tuple(range(a.ndim))
+    return atop(np.ma.masked_equal, inds, a, inds, value, (), dtype=a.dtype)
+
+
+@wraps(np.ma.masked_invalid)
+def masked_invalid(a):
+    return asanyarray(a).map_blocks(np.ma.masked_invalid)
+
+
+@wraps(np.ma.masked_inside)
+def masked_inside(x, v1, v2):
+    x = asanyarray(x)
+    return x.map_blocks(np.ma.masked_inside, v1, v2)
+
+
+@wraps(np.ma.masked_outside)
+def masked_outside(x, v1, v2):
+    x = asanyarray(x)
+    return x.map_blocks(np.ma.masked_outside, v1, v2)
+
+
+@wraps(np.ma.masked_where)
+def masked_where(condition, a):
+    cshape = getattr(condition, 'shape', ())
+    if cshape and cshape != a.shape:
+        raise IndexError("Inconsistant shape between the condition and the "
+                         "input (got %s and %s)" % (cshape, a.shape))
+    condition = asanyarray(condition)
+    a = asanyarray(a)
+    ainds = tuple(range(a.ndim))
+    cinds = tuple(range(condition.ndim))
+    return atop(np.ma.masked_where, ainds, condition, cinds, a, ainds,
+                dtype=a.dtype)
+
+
+@wraps(np.ma.masked_values)
+def masked_values(x, value, rtol=1e-05, atol=1e-08, shrink=True):
+    x = asanyarray(x)
+    if getattr(value, 'shape', ()):
+        raise ValueError("da.ma.masked_values doesn't support array `value`s")
+    return map_blocks(np.ma.masked_values, x, value, rtol=rtol,
+                      atol=atol, shrink=shrink)
+
+
+@wraps(np.ma.fix_invalid)
+def fix_invalid(a, fill_value=None):
+    a = asanyarray(a)
+    return a.map_blocks(np.ma.fix_invalid, fill_value=fill_value)
+
+
+@wraps(np.ma.getdata)
+def getdata(a):
+    a = asanyarray(a)
+    return a.map_blocks(np.ma.getdata)
+
+
+@wraps(np.ma.getmaskarray)
+def getmaskarray(a):
+    a = asanyarray(a)
+    return a.map_blocks(np.ma.getmaskarray)
+
+
+def _masked_array(data, mask=np.ma.nomask, **kwargs):
+    dtype = kwargs.pop('masked_dtype', None)
+    return np.ma.masked_array(data, mask=mask, dtype=dtype, **kwargs)
+
+
+@wraps(np.ma.masked_array)
+def masked_array(data, mask=np.ma.nomask, fill_value=None,
+                 **kwargs):
+    data = asanyarray(data)
+    inds = tuple(range(data.ndim))
+    arginds = [inds, data, inds]
+
+    if getattr(fill_value, 'shape', ()):
+        raise ValueError("non-scalar fill_value not supported")
+    kwargs['fill_value'] = fill_value
+
+    if mask is not np.ma.nomask:
+        mask = asanyarray(mask)
+        if mask.size == 1:
+            mask = mask.reshape((1,) * data.ndim)
+        elif data.shape != mask.shape:
+            raise np.ma.MaskError("Mask and data not compatible: data shape "
+                                  "is %s, and mask shape is "
+                                  "%s." % (repr(data.shape), repr(mask.shape)))
+        arginds.extend([mask, inds])
+
+    if 'dtype' in kwargs:
+        kwargs['masked_dtype'] = kwargs['dtype']
+    else:
+        kwargs['dtype'] = data.dtype
+
+    return atop(_masked_array, *arginds, **kwargs)
+
+
+def _set_fill_value(x, fill_value):
+    if isinstance(x, np.ma.masked_array):
+        x = x.copy()
+        np.ma.set_fill_value(x, fill_value=fill_value)
+    return x
+
+
+@wraps(np.ma.set_fill_value)
+def set_fill_value(a, fill_value):
+    a = asanyarray(a)
+    if getattr(fill_value, 'shape', ()):
+        raise ValueError("da.ma.set_fill_value doesn't support array `value`s")
+    fill_value = np.ma.core._check_fill_value(fill_value, a.dtype)
+    res = a.map_blocks(_set_fill_value, fill_value)
+    a.dask = res.dask
+    a.name = res.name
diff --git a/dask/array/numpy_compat.py b/dask/array/numpy_compat.py
index 1ba08c302..7a289a88c 100644
--- a/dask/array/numpy_compat.py
+++ b/dask/array/numpy_compat.py
@@ -34,6 +34,7 @@ try:
             raise TypeError('Divide not working with dtype: '
                             'https://github.com/numpy/numpy/issues/3484')
         divide = np.divide
+        ma_divide = np.ma.divide
 
 except TypeError:
     # Divide with dtype doesn't work on Python 3
@@ -47,6 +48,10 @@ except TypeError:
             x = x.astype(dtype)
         return x
 
+    ma_divide = np.ma.core._DomainedBinaryOperation(divide,
+                                                    np.ma.core._DomainSafeDivide(),
+                                                    0, 1)
+
 # functions copied from numpy
 try:
     from numpy import broadcast_to, nanprod, nancumsum, nancumprod
diff --git a/dask/array/reductions.py b/dask/array/reductions.py
index b9e76a5b5..2a9def705 100644
--- a/dask/array/reductions.py
+++ b/dask/array/reductions.py
@@ -12,14 +12,29 @@ from . import chunk
 from .core import _concatenate2, Array, atop, lol_tuples, handle_out
 from .ufunc import sqrt
 from .wrap import zeros, ones
-from .numpy_compat import divide
+from .numpy_compat import ma_divide, divide as np_divide
 from ..compatibility import getargspec, builtins
 from ..base import tokenize
 from ..context import _globals
-from ..utils import ignoring, funcname
+from ..utils import ignoring, funcname, Dispatch
 from .. import sharedict
 
 
+# Generic functions to support chunks of different types
+empty_lookup = Dispatch('empty')
+empty_lookup.register((object, np.ndarray), np.empty)
+empty_lookup.register(np.ma.masked_array, np.ma.empty)
+divide_lookup = Dispatch('divide')
+divide_lookup.register((object, np.ndarray), np_divide)
+divide_lookup.register(np.ma.masked_array, ma_divide)
+
+
+def divide(a, b, dtype=None):
+    key = lambda x: getattr(x, '__array_priority__', float('-inf'))
+    f = divide_lookup.dispatch(type(builtins.max(a, b, key=key)))
+    return f(a, b, dtype=dtype)
+
+
 def reduction(x, chunk, aggregate, axis=None, keepdims=None, dtype=None,
               split_every=None, combine=None, name=None, out=None):
     """ General version of reductions
@@ -225,8 +240,8 @@ def nannumel(x, **kwargs):
 def mean_chunk(x, sum=chunk.sum, numel=numel, dtype='f8', **kwargs):
     n = numel(x, dtype=dtype, **kwargs)
     total = sum(x, dtype=dtype, **kwargs)
-    result = np.empty(shape=n.shape,
-                      dtype=[('total', total.dtype), ('n', n.dtype)])
+    empty = empty_lookup.dispatch(type(n))
+    result = empty(n.shape, dtype=[('total', total.dtype), ('n', n.dtype)])
     result['n'] = n
     result['total'] = total
     return result
@@ -235,7 +250,8 @@ def mean_chunk(x, sum=chunk.sum, numel=numel, dtype='f8', **kwargs):
 def mean_combine(pair, sum=chunk.sum, numel=numel, dtype='f8', **kwargs):
     n = sum(pair['n'], **kwargs)
     total = sum(pair['total'], **kwargs)
-    result = np.empty(shape=n.shape, dtype=pair.dtype)
+    empty = empty_lookup.dispatch(type(n))
+    result = empty(n.shape, dtype=pair.dtype)
     result['n'] = n
     result['total'] = total
     return result
@@ -275,14 +291,15 @@ with ignoring(AttributeError):
 
 def moment_chunk(A, order=2, sum=chunk.sum, numel=numel, dtype='f8', **kwargs):
     total = sum(A, dtype=dtype, **kwargs)
-    n = numel(A, **kwargs).astype(np.int64, copy=False)
+    n = numel(A, **kwargs).astype(np.int64)
     u = total / n
-    M = np.empty(shape=n.shape + (order - 1,), dtype=dtype)
+    empty = empty_lookup.dispatch(type(n))
+    M = empty(n.shape + (order - 1,), dtype=dtype)
     for i in range(2, order + 1):
         M[..., i - 2] = sum((A - u)**i, dtype=dtype, **kwargs)
-    result = np.empty(shape=n.shape, dtype=[('total', total.dtype),
-                                            ('n', n.dtype),
-                                            ('M', M.dtype, (order - 1,))])
+    result = empty(n.shape, dtype=[('total', total.dtype),
+                                   ('n', n.dtype),
+                                   ('M', M.dtype, (order - 1,))])
     result['total'] = total
     result['n'] = n
     result['M'] = M
@@ -308,14 +325,15 @@ def moment_combine(data, order=2, ddof=0, dtype='f8', sum=np.sum, **kwargs):
     n = sum(ns, **kwargs)
     mu = divide(total, n, dtype=dtype)
     inner_term = divide(totals, ns, dtype=dtype) - mu
-    M = np.empty(shape=n.shape + (order - 1,), dtype=dtype)
+    empty = empty_lookup.dispatch(type(n))
+    M = empty(n.shape + (order - 1,), dtype=dtype)
 
     for o in range(2, order + 1):
         M[..., o - 2] = _moment_helper(Ms, ns, inner_term, o, sum, kwargs)
 
-    result = np.zeros(shape=n.shape, dtype=[('total', total.dtype),
-                                            ('n', n.dtype),
-                                            ('M', Ms.dtype, (order - 1,))])
+    result = empty(n.shape, dtype=[('total', total.dtype),
+                                   ('n', n.dtype),
+                                   ('M', Ms.dtype, (order - 1,))])
     result['total'] = total
     result['n'] = n
     result['M'] = M
@@ -471,6 +489,13 @@ def arg_chunk(func, argfunc, x, axis, offset_info):
     else:
         arg += offset_info
 
+    if isinstance(vals, np.ma.masked_array):
+        if 'min' in argfunc.__name__:
+            fill_value = np.ma.minimum_fill_value(vals)
+        else:
+            fill_value = np.ma.maximum_fill_value(vals)
+        vals = np.ma.filled(vals, fill_value)
+
     result = np.empty(shape=vals.shape, dtype=[('vals', vals.dtype),
                                                ('arg', arg.dtype)])
     result['vals'] = vals
@@ -653,14 +678,28 @@ def cumreduction(func, binop, ident, x, axis, dtype=None, out=None):
     return handle_out(out, result)
 
 
+def _cumsum_merge(a, b):
+    if isinstance(a, np.ma.masked_array) or isinstance(b, np.ma.masked_array):
+        values = np.ma.getdata(a) + np.ma.getdata(b)
+        return np.ma.masked_array(values, mask=np.ma.getmaskarray(b))
+    return a + b
+
+
+def _cumprod_merge(a, b):
+    if isinstance(a, np.ma.masked_array) or isinstance(b, np.ma.masked_array):
+        values = np.ma.getdata(a) * np.ma.getdata(b)
+        return np.ma.masked_array(values, mask=np.ma.getmaskarray(b))
+    return a * b
+
+
 @wraps(np.cumsum)
 def cumsum(x, axis, dtype=None, out=None):
-    return cumreduction(np.cumsum, operator.add, 0, x, axis, dtype, out=out)
+    return cumreduction(np.cumsum, _cumsum_merge, 0, x, axis, dtype, out=out)
 
 
 @wraps(np.cumprod)
 def cumprod(x, axis, dtype=None, out=None):
-    return cumreduction(np.cumprod, operator.mul, 1, x, axis, dtype, out=out)
+    return cumreduction(np.cumprod, _cumprod_merge, 1, x, axis, dtype, out=out)
 
 
 def validate_axis(ndim, axis):
diff --git a/dask/array/routines.py b/dask/array/routines.py
index 53081d46a..e07082bb3 100644
--- a/dask/array/routines.py
+++ b/dask/array/routines.py
@@ -14,12 +14,11 @@ from toolz import concat, sliding_window, interleave
 from .. import sharedict
 from ..core import flatten
 from ..base import tokenize
-from ..utils import package_of
 from . import numpy_compat, chunk
 
 from .core import (Array, map_blocks, elemwise, from_array, asarray,
                    concatenate, stack, atop, broadcast_shapes,
-                   is_scalar_for_elemwise, broadcast_to)
+                   is_scalar_for_elemwise, broadcast_to, tensordot_lookup)
 
 
 @wraps(np.array)
@@ -111,8 +110,8 @@ ALPHABET = alphabet.upper()
 
 def _tensordot(a, b, axes):
     x = max([a, b], key=lambda x: x.__array_priority__)
-    module = package_of(type(x)) or np
-    x = module.tensordot(a, b, axes=axes)
+    tensordot = tensordot_lookup.dispatch(type(x))
+    x = tensordot(a, b, axes=axes)
     ind = [slice(None, None)] * x.ndim
     for a in sorted(axes[0]):
         ind.insert(a, None)
diff --git a/dask/utils.py b/dask/utils.py
index 51f703f2b..f16a3c283 100644
--- a/dask/utils.py
+++ b/dask/utils.py
@@ -824,28 +824,6 @@ def ensure_dict(d):
     return dict(d)
 
 
-_packages = {}
-
-
-def package_of(typ):
-    """ Return package containing type's definition
-
-    Or return None if not found
-    """
-    try:
-        return _packages[typ]
-    except KeyError:
-        # http://stackoverflow.com/questions/43462701/get-package-of-python-object/43462865#43462865
-        mod = inspect.getmodule(typ)
-        if not mod:
-            result = None
-        else:
-            base, _sep, _stem = mod.__name__.partition('.')
-            result = sys.modules[base]
-        _packages[typ] = result
-        return result
-
-
 # XXX: Kept to keep old versions of distributed/dask in sync. After
 # distributed's dask requirement is updated to > this commit, this function can
 # be moved to dask.bytes.utils.
diff --git a/docs/source/array-api.rst b/docs/source/array-api.rst
index 2091bd746..98ff6fa02 100644
--- a/docs/source/array-api.rst
+++ b/docs/source/array-api.rst
@@ -192,6 +192,28 @@ Linear Algebra
    linalg.svd_compressed
    linalg.tsqr
 
+Masked Arrays
+~~~~~~~~~~~~~
+
+.. autosummary::
+   ma.filled
+   ma.fix_invalid
+   ma.getdata
+   ma.getmaskarray
+   ma.masked_array
+   ma.masked_equal
+   ma.masked_greater
+   ma.masked_greater_equal
+   ma.masked_inside
+   ma.masked_invalid
+   ma.masked_less
+   ma.masked_less_equal
+   ma.masked_not_equal
+   ma.masked_outside
+   ma.masked_values
+   ma.masked_where
+   ma.set_fill_value
+
 Random
 ~~~~~~
 
@@ -458,6 +480,25 @@ Other functions
 .. autofunction:: svd_compressed
 .. autofunction:: tsqr
 
+.. currentmodule:: dask.array.ma
+.. autofunction:: filled
+.. autofunction:: fix_invalid
+.. autofunction:: getdata
+.. autofunction:: getmaskarray
+.. autofunction:: masked_array
+.. autofunction:: masked_equal
+.. autofunction:: masked_greater
+.. autofunction:: masked_greater_equal
+.. autofunction:: masked_inside
+.. autofunction:: masked_invalid
+.. autofunction:: masked_less
+.. autofunction:: masked_less_equal
+.. autofunction:: masked_not_equal
+.. autofunction:: masked_outside
+.. autofunction:: masked_values
+.. autofunction:: masked_where
+.. autofunction:: set_fill_value
+
 .. currentmodule:: dask.array.ghost
 
 .. autofunction:: ghost
diff --git a/docs/source/array-sparse.rst b/docs/source/array-sparse.rst
index df2a2c8c7..952c9d591 100644
--- a/docs/source/array-sparse.rst
+++ b/docs/source/array-sparse.rst
@@ -55,18 +55,21 @@ following operations:
 
 1.  Simple slicing with slices, lists, and elements (for slicing, rechunking,
     reshaping, etc).
-2.  A ``concatenate`` function at the top level of the package (for assembling
-    results)
+2.  A ``concatenate`` function matching the interface of ``np.concatenate``.
+    This must be registered in ``dask.array.core.concatenate_lookup``.
 3.  All ufuncs must support the full ufunc interface, including ``dtype=`` and
     ``out=`` parameters (even if they don't function properly)
 4.  All reductions must support the full ``axis=`` and ``keepdims=`` keywords
     and behave like numpy in this respect
 5.  The array class should follow the ``__array_priority__`` protocol and be
     prepared to respond to other arrays of lower priority.
+6.  If ``dot`` support is desired, a ``tensordot`` function matching the
+    interface of ``np.tensordot`` should be registered in
+    ``dask.array.core.tensordot_lookup``.
 
-The implementation of other operations like tensordot, reshape, transpose, etc.
+The implementation of other operations like reshape, transpose, etc.
 should follow standard NumPy conventions regarding shape and dtype.  Not
-implementing these is fine; the parallel dask.array will err at runtime if
+implementing these is fine; the parallel ``dask.array`` will err at runtime if
 these operations are attempted.
 
 


Test patch:
diff --git a/dask/array/tests/test_masked.py b/dask/array/tests/test_masked.py
new file mode 100644
index 000000000..a594bbfb9
--- /dev/null
+++ b/dask/array/tests/test_masked.py
@@ -0,0 +1,315 @@
+import random
+from itertools import product
+
+import numpy as np
+import pytest
+
+import dask.array as da
+from dask.base import tokenize
+from dask.array.utils import assert_eq
+
+pytest.importorskip("dask.array.ma")
+
+
+def test_tokenize_masked_array():
+    m = np.ma.masked_array([1, 2, 3], mask=[True, True, False], fill_value=10)
+    m2 = np.ma.masked_array([1, 2, 3], mask=[True, True, False], fill_value=0)
+    m3 = np.ma.masked_array([1, 2, 3], mask=False, fill_value=10)
+    assert tokenize(m) == tokenize(m)
+    assert tokenize(m2) == tokenize(m2)
+    assert tokenize(m3) == tokenize(m3)
+    assert tokenize(m) != tokenize(m2)
+    assert tokenize(m) != tokenize(m3)
+
+
+def test_from_array_masked_array():
+    m = np.ma.masked_array([1, 2, 3], mask=[True, True, False], fill_value=10)
+    dm = da.from_array(m, chunks=(2,), asarray=False)
+    assert_eq(dm, m)
+
+
+functions = [
+    lambda x: x,
+    lambda x: da.expm1(x),
+    lambda x: 2 * x,
+    lambda x: x / 2,
+    lambda x: x**2,
+    lambda x: x + x,
+    lambda x: x * x,
+    lambda x: x[0],
+    lambda x: x[:, 1],
+    lambda x: x[:1, None, 1:3],
+    lambda x: x.T,
+    lambda x: da.transpose(x, (1, 2, 0)),
+    lambda x: x.sum(),
+    lambda x: x.dot(np.arange(x.shape[-1])),
+    lambda x: x.dot(np.eye(x.shape[-1])),
+    lambda x: da.tensordot(x, np.ones(x.shape[:2]), axes=[(0, 1), (0, 1)]),
+    lambda x: x.sum(axis=0),
+    lambda x: x.max(axis=0),
+    lambda x: x.sum(axis=(1, 2)),
+    lambda x: x.astype(np.complex128),
+    lambda x: x.map_blocks(lambda x: x * 2),
+    lambda x: x.round(1),
+    lambda x: x.reshape((x.shape[0] * x.shape[1], x.shape[2])),
+    lambda x: abs(x),
+    lambda x: x > 0.5,
+    lambda x: x.rechunk((4, 4, 4)),
+    lambda x: x.rechunk((2, 2, 1)),
+]
+
+
+@pytest.mark.parametrize('func', functions)
+def test_basic(func):
+    x = da.random.random((2, 3, 4), chunks=(1, 2, 2))
+    x[x < 0.4] = 0
+
+    y = da.ma.masked_equal(x, 0)
+
+    xx = func(x)
+    yy = func(y)
+
+    assert_eq(xx, da.ma.filled(yy, 0))
+
+    if yy.shape:
+        zz = yy.compute()
+        assert isinstance(zz, np.ma.masked_array)
+
+
+def test_tensordot():
+    x = da.random.random((2, 3, 4), chunks=(1, 2, 2))
+    x[x < 0.4] = 0
+    y = da.random.random((4, 3, 2), chunks=(2, 2, 1))
+    y[y < 0.4] = 0
+
+    xx = da.ma.masked_equal(x, 0)
+    yy = da.ma.masked_equal(y, 0)
+
+    assert_eq(da.tensordot(x, y, axes=(2, 0)),
+              da.ma.filled(da.tensordot(xx, yy, axes=(2, 0)), 0))
+    assert_eq(da.tensordot(x, y, axes=(1, 1)),
+              da.ma.filled(da.tensordot(xx, yy, axes=(1, 1)), 0))
+    assert_eq(da.tensordot(x, y, axes=((1, 2), (1, 0))),
+              da.ma.filled(da.tensordot(xx, yy, axes=((1, 2), (1, 0))), 0))
+
+
+@pytest.mark.parametrize('func', functions)
+def test_mixed_concatenate(func):
+    x = da.random.random((2, 3, 4), chunks=(1, 2, 2))
+    y = da.random.random((2, 3, 4), chunks=(1, 2, 2))
+
+    y[y < 0.4] = 0
+    yy = da.ma.masked_equal(y, 0)
+
+    d = da.concatenate([x, y], axis=0)
+    s = da.concatenate([x, yy], axis=0)
+
+    dd = func(d)
+    ss = func(s)
+    assert_eq(dd, ss)
+
+
+@pytest.mark.parametrize('func', functions)
+def test_mixed_random(func):
+    d = da.random.random((4, 3, 4), chunks=(1, 2, 2))
+    d[d < 0.4] = 0
+
+    fn = lambda x: np.ma.masked_equal(x, 0) if random.random() < 0.5 else x
+    s = d.map_blocks(fn)
+
+    dd = func(d)
+    ss = func(s)
+
+    assert_eq(dd, ss)
+
+
+def test_mixed_output_type():
+    y = da.random.random((10, 10), chunks=(5, 5))
+    y[y < 0.4] = 0
+
+    y = da.ma.masked_equal(y, 0)
+    x = da.zeros((10, 1), chunks=(5, 1))
+
+    z = da.concatenate([x, y], axis=1)
+    assert z.shape == (10, 11)
+    zz = z.compute()
+    assert isinstance(zz, np.ma.masked_array)
+
+
+def test_creation_functions():
+    x = np.array([-2, -1, 0, 1, 2] * 20).reshape((10, 10))
+    y = np.array([-2, 0, 1, 1, 0] * 2)
+    dx = da.from_array(x, chunks=5)
+    dy = da.from_array(y, chunks=4)
+
+    sol = np.ma.masked_greater(x, y)
+    for (a, b) in product([dx, x], [dy, y]):
+        assert_eq(da.ma.masked_greater(a, b), sol)
+
+    # These are all the same as masked_greater, just check for correct op
+    assert_eq(da.ma.masked_greater(dx, 0), np.ma.masked_greater(x, 0))
+    assert_eq(da.ma.masked_greater_equal(dx, 0), np.ma.masked_greater_equal(x, 0))
+    assert_eq(da.ma.masked_less(dx, 0), np.ma.masked_less(x, 0))
+    assert_eq(da.ma.masked_less_equal(dx, 0), np.ma.masked_less_equal(x, 0))
+    assert_eq(da.ma.masked_equal(dx, 0), np.ma.masked_equal(x, 0))
+    assert_eq(da.ma.masked_not_equal(dx, 0), np.ma.masked_not_equal(x, 0))
+
+    # masked_where
+    assert_eq(da.ma.masked_where(False, dx), np.ma.masked_where(False, x))
+    assert_eq(da.ma.masked_where(dx > 2, dx), np.ma.masked_where(x > 2, x))
+
+    with pytest.raises(IndexError):
+        da.ma.masked_where((dx > 2)[:, 0], dx)
+
+    assert_eq(da.ma.masked_inside(dx, -1, 1), np.ma.masked_inside(x, -1, 1))
+    assert_eq(da.ma.masked_outside(dx, -1, 1), np.ma.masked_outside(x, -1, 1))
+    assert_eq(da.ma.masked_values(dx, -1), np.ma.masked_values(x, -1))
+
+    # masked_equal and masked_values in numpy sets the fill_value to `value`,
+    # which can sometimes be an array. This is hard to support in dask, so we
+    # forbid it. Check that this isn't supported:
+    with pytest.raises(ValueError):
+        da.ma.masked_equal(dx, dy)
+
+    with pytest.raises(ValueError):
+        da.ma.masked_values(dx, dy)
+
+    y = x.astype('f8')
+    y[0, 0] = y[7, 5] = np.nan
+    dy = da.from_array(y, chunks=5)
+
+    assert_eq(da.ma.masked_invalid(dy), np.ma.masked_invalid(y))
+
+    my = np.ma.masked_greater(y, 0)
+    dmy = da.ma.masked_greater(dy, 0)
+
+    assert_eq(da.ma.fix_invalid(dmy, fill_value=0),
+              np.ma.fix_invalid(my, fill_value=0))
+
+
+def test_filled():
+    x = np.array([-2, -1, 0, 1, 2] * 20).reshape((10, 10))
+    dx = da.from_array(x, chunks=5)
+
+    mx = np.ma.masked_equal(x, 0)
+    mdx = da.ma.masked_equal(dx, 0)
+
+    assert_eq(da.ma.filled(mdx), np.ma.filled(mx))
+    assert_eq(da.ma.filled(mdx, -5), np.ma.filled(mx, -5))
+
+
+def assert_eq_ma(a, b):
+    res = a.compute()
+    assert type(res) == type(b)
+    if hasattr(res, 'mask'):
+        np.testing.assert_equal(res.mask, b.mask)
+        a = da.ma.filled(a)
+        b = np.ma.filled(b)
+    assert_eq(a, b, equal_nan=True)
+
+
+@pytest.mark.parametrize('dtype', ('i8', 'f8'))
+@pytest.mark.parametrize('reduction', ['sum', 'prod', 'mean', 'var', 'std',
+                                       'min', 'max', 'any', 'all'])
+def test_reductions(dtype, reduction):
+    x = (np.random.RandomState(42).rand(11, 11) * 10).astype(dtype)
+    dx = da.from_array(x, chunks=(4, 4))
+    mx = np.ma.masked_greater(x, 5)
+    mdx = da.ma.masked_greater(dx, 5)
+
+    dfunc = getattr(da, reduction)
+    func = getattr(np, reduction)
+
+    assert_eq_ma(dfunc(mdx), func(mx))
+    assert_eq_ma(dfunc(mdx, axis=0), func(mx, axis=0))
+    assert_eq_ma(dfunc(mdx, keepdims=True, split_every=4),
+                 func(mx, keepdims=True))
+    assert_eq_ma(dfunc(mdx, axis=0, split_every=2), func(mx, axis=0))
+    assert_eq_ma(dfunc(mdx, axis=0, keepdims=True, split_every=2),
+                 func(mx, axis=0, keepdims=True))
+    assert_eq_ma(dfunc(mdx, axis=1, split_every=2), func(mx, axis=1))
+    assert_eq_ma(dfunc(mdx, axis=1, keepdims=True, split_every=2),
+                 func(mx, axis=1, keepdims=True))
+
+
+@pytest.mark.parametrize('reduction', ['argmin', 'argmax'])
+def test_arg_reductions(reduction):
+    x = np.random.random((10, 10, 10))
+    dx = da.from_array(x, chunks=(3, 4, 5))
+    mx = np.ma.masked_greater(x, 0.4)
+    dmx = da.ma.masked_greater(dx, 0.4)
+
+    dfunc = getattr(da, reduction)
+    func = getattr(np, reduction)
+
+    assert_eq_ma(dfunc(dmx), func(mx))
+    assert_eq_ma(dfunc(dmx, 0), func(mx, 0))
+    assert_eq_ma(dfunc(dmx, 1), func(mx, 1))
+    assert_eq_ma(dfunc(dmx, 2), func(mx, 2))
+
+
+def test_cumulative():
+    x = np.random.RandomState(0).rand(20, 24, 13)
+    dx = da.from_array(x, chunks=(6, 5, 4))
+    mx = np.ma.masked_greater(x, 0.4)
+    dmx = da.ma.masked_greater(dx, 0.4)
+
+    for axis in [0, 1, 2]:
+        assert_eq_ma(dmx.cumsum(axis=axis), mx.cumsum(axis=axis))
+        assert_eq_ma(dmx.cumprod(axis=axis), mx.cumprod(axis=axis))
+
+
+def test_accessors():
+    x = np.random.random((10, 10))
+    dx = da.from_array(x, chunks=(3, 4))
+    mx = np.ma.masked_greater(x, 0.4)
+    dmx = da.ma.masked_greater(dx, 0.4)
+
+    assert_eq(da.ma.getmaskarray(dmx), np.ma.getmaskarray(mx))
+    assert_eq(da.ma.getmaskarray(dx), np.ma.getmaskarray(x))
+    assert_eq(da.ma.getdata(dmx), np.ma.getdata(mx))
+    assert_eq(da.ma.getdata(dx), np.ma.getdata(x))
+
+
+def test_masked_array():
+    x = np.random.random((10, 10)).astype('f4')
+    dx = da.from_array(x, chunks=(3, 4))
+    f1 = da.from_array(np.array(1), chunks=())
+
+    fill_values = [(None, None), (0.5, 0.5), (1, f1)]
+    for data, (df, f) in product([x, dx], fill_values):
+        assert_eq(da.ma.masked_array(data, fill_value=df),
+                  np.ma.masked_array(x, fill_value=f))
+        assert_eq(da.ma.masked_array(data, mask=data > 0.4, fill_value=df),
+                  np.ma.masked_array(x, mask=x > 0.4, fill_value=f))
+        assert_eq(da.ma.masked_array(data, mask=data > 0.4, fill_value=df),
+                  np.ma.masked_array(x, mask=x > 0.4, fill_value=f))
+        assert_eq(da.ma.masked_array(data, fill_value=df, dtype='f8'),
+                  np.ma.masked_array(x, fill_value=f, dtype='f8'))
+
+    with pytest.raises(ValueError):
+        da.ma.masked_array(dx, fill_value=dx)
+
+    with pytest.raises(np.ma.MaskError):
+        da.ma.masked_array(dx, mask=dx[:3, :3])
+
+
+def test_set_fill_value():
+    x = np.random.randint(0, 10, (10, 10))
+    dx = da.from_array(x, chunks=(3, 4))
+    mx = np.ma.masked_greater(x, 3)
+    dmx = da.ma.masked_greater(dx, 3)
+
+    da.ma.set_fill_value(dmx, -10)
+    np.ma.set_fill_value(mx, -10)
+    assert_eq_ma(dmx, mx)
+
+    da.ma.set_fill_value(dx, -10)
+    np.ma.set_fill_value(x, -10)
+    assert_eq_ma(dx, x)
+
+    with pytest.raises(TypeError):
+        da.ma.set_fill_value(dmx, 1e20)
+
+    with pytest.raises(ValueError):
+        da.ma.set_fill_value(dmx, dx)
diff --git a/dask/tests/test_utils.py b/dask/tests/test_utils.py
index a7780279a..5258cccd4 100644
--- a/dask/tests/test_utils.py
+++ b/dask/tests/test_utils.py
@@ -8,7 +8,7 @@ from dask.sharedict import ShareDict
 from dask.utils import (takes_multiple_arguments, Dispatch, random_state_data,
                         memory_repr, methodcaller, M, skip_doctest,
                         SerializableLock, funcname, ndeepmap, ensure_dict,
-                        package_of, extra_titles, asciitable, itemgetter)
+                        extra_titles, asciitable, itemgetter)
 from dask.utils_test import inc
 
 
@@ -304,17 +304,6 @@ def test_ensure_dict():
     assert ensure_dict(md) == d
 
 
-def test_package_of():
-    import math
-    assert package_of(math.sin) is math
-    try:
-        import numpy
-    except ImportError:
-        pass
-    else:
-        assert package_of(numpy.memmap) is numpy
-
-
 def test_itemgetter():
     data = [1, 2, 3]
     g = itemgetter(1)