Hello,

I have been looking into the feasibility of replacing the internals of the data object in cf-python with dask. I am very impressed by dask's functionality, and it is promising that we can use it.

However, the lack of a fully featured __setitem__ is a problem for us - particulary becuase we need to assign to larger-than-memory arrays (possibly with larger-than-memory indices!).

Having read some of the previous discussions on this topic in this issue tracker, it has occured to me that the cf-python implementation of __setitem__ could actually be applied to dask: The cf-python implementation of __setitem__ is domain-distributed with nearly all of the numpy API, but does not have a lazy functionality.

I have found that by combining the cf-python approach with map_blocks I have (it seems!) implemented a much more richly featured lazy __setitem__ in dask.

I would be happy to submit a pull request with the new code (which only modifies Array.__setitem__) that implements this.

With my new code I can do things like:

>>> import dask.array as da  # from a branch in my fork
>>> import numpy
>>> x = da.from_array(numpy.arange(40.).reshape(5, 8), chunks=(2, 3))
>>> x.compute() 
array([[ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.],
       [ 8.,  9., 10., 11., 12., 13., 14., 15.],
       [16., 17., 18., 19., 20., 21., 22., 23.],
       [24., 25., 26., 27., 28., 29., 30., 31.],
       [32., 33., 34., 35., 36., 37., 38., 39.]])

>>> x[..., 3:1:-1] = -x[2, 4:6]
>>> x.compute()
array([[  0.,   1., -21., -20.,   4.,   5.,   6.,   7.],
       [  8.,   9., -21., -20.,  12.,  13.,  14.,  15.],
       [ 16.,  17., -21., -20.,  20.,  21.,  22.,  23.],
       [ 24.,  25., -21., -20.,  28.,  29.,  30.,  31.],
       [ 32.,  33., -21., -20.,  36.,  37.,  38.,  39.]])

which I hope demonstrates that it is "safe".

Perhaps a tougher test is the situation described in #2000 (comment). This works for me in the "documented gotcha" sense.

>>> import dask.array as da
>>> import numpy
>>> x = da.from_array(numpy.arange(40.).reshape(5, 8), chunks=(2, 3))
>>> y = x[0, :]
>>> y[:] = 0
>>> x[0, :].compute()
array([0. 1. 2. 3. 4. 5. 6. 7.])
>>> y.compute()
array([0., 0., 0., 0., 0., 0., 0., 0.])
>>> x[0, :].compute()
array([0., 0., 0., 0., 0., 0., 0., 0.])

If we did the same thing in numpy:

>>> import numpy
>>> x = numpy.arange(40.).reshape(5, 8)
>>> y = x[0, :]
>>> y[:] = 0
>>> x[0, :]
array([0., 0., 0., 0., 0., 0., 0., 0.])

the row in x is zero-ed straight away, because numpy is non-lazy. In both cases y is a view of part of x, but in the dask case we have to factor in lazy evaluation in out understanding of views. How does that sound (I'm aware that I surely don't have enough experience and knowledge on how dask behaves in these areas)?

The full numpy broadcasting rules are implemented.

The RHS assignement value can be any array-like as accepted by dask.array.asanyarray.

The existing 1-d boolean array functionality is retained (but not implemented with where anymore).

I think that the parts of the numpy API that are missing are a) the ability provide two or more indices that are sequences of integers or booleans, (such as x[[1, 2], [2, 3]]) and b) the ability to provide a non-strictly monotonically increasing or decreasing sequence (such as x[[3, 6, 1]] or x[[3, 3, 1]]). I suspect that it the code could be extended to include these cases - the main reason they are not already there is because these use cases are not allowed by cf-python.

Everything else is allowed, such as:

x[3, ..., -2:-1] = [[-99]]
x[[True False, True, False, False], -2:1:-3] = y

etc.

I well appreciate that reviewing pull requests can be time consuming, but it seems that there is some general interest in this feature if people like the functionality, and what I have done could well be improved by a dask expert.

Please let me know if a documented and commented PR would be welcome,
Many thanks,

David Hassell