Several modifications to fix slicing bug - #1776

I don't think h5py actually supports outer indexing -- it raises an error if you attempt to index along multiple dimensions with lists. It might be more accurate to say that it supports the "intuitive" part of NumPy's mixed slice/array fancy indexing, where axes end up in the "obvious" order. (I use scare quotes here because this isn't really obvious when you consider fancy indexing in full generality, e.g., as .vindex for NumPy)

Likewise, netCDF4 doesn't support "outer" slicing in an efficient manner -- under the hood it basically generates a loop (in Python) of calls to the underlying HDF5 library.

It's not immediately obvious what the right solution here is, but maybe this provide some valuable context?

You're right, h5py does not support outer indexing! I guess I was fooled by the behaviour of h5netcdf.legacyapi. Nevertheless, in h5py, unlike in numpy,

a[0, :, [0, 2]] == a[0, :, :][:, [0, 2]]

following the example described in #1770.

From my testing, it seems clear that such slices should not generally be fused, otherwise issue #1770 could pop-up.
One option, is to prevent these slices in those cases. However, in the case of a netCDF4/h5py file, no matter the efficiency of the indexing implementation, computing a[0, :, [0, 2]] will be more memory-efficient that computing a[0, :, :][:, [0, 2]], especially if the last axis is large. With the current fix to #1770 but without the distinction between getarray_outer and getarray_vector, the former would be computed.
Another option, as proposed in this PR, is to allow users to specify whether this getitem method works on the underlying array.

I agree with you that it is not completely clear whether the right course of action is either one of these options.

Nevertheless, in h5py, unlike in numpy,

a[0, :, [0, 2]] == a[0, :, :][:, [0, 2]]

following the example described in #1770.

Maybe raise an issue about this in the h5py issue tracker? It seems like something that they should at least know about.

In fact it's in issue h5py/h5py/issues/652. It has been raised about a year ago.

Another option, would be to revert to fancy being a boolean and eagerly check the behaviour of the array on a minimal subset of the array in da.from_array. Then one could decide whether fancy=True should use getarray_outer or getarray_vector based on this check.
In the same way, we could always check fancy indexing rules on the underlying and essentially forgo the fancy keyword...

Another option, would be to revert to fancy being a boolean and eagerly check the behaviour of the array on a minimal subset of the array in da.from_array. Then one could decide whether fancy=True should use getarray_outer or getarray_vector based on this check.

This seems like a bad idea. Dask should not be evaluating pieces of arrays values. I like explicit options... just need to think through this a little more.

I agree with you that explicit options are usually better!

One approach is to default optimizations to vector indexing when fancy=True unless method .oindex exists and method .vindex does not. Then we could create a simple wrapper that would allow users to expose an .oindex method on their objects that support outer indexing. In summary:

1. da.from_array(x, fancy=False) -> current optimization
2. da.from_array(x) -> current optimization (but maybe done through getarray_vector to allow fixing of slicing bug) if (neither .vindex nor .oindex are methods of x) OR if (.vindex is a method of x).
3. da.from_array(x) -> new optimization possibilities (done through getarray_outer or equivalent) if .oindex is a method of x but .vindex is not.
4. Create a wrapper e.g. da.outer_array(x) that exposes a .oindex method on x and hides the .vindex method (if it exists) on x.

This would seem to impose minimal changes to the user interface while allowing the optimization to be aware of the allowed optimal (in the sense of minimal amount of data brought to memory) underlying slicing methods.

One advantage of this approach is that with the current PR, I could remove the possibility of specifying fancy='outer' in da.from_array while keeping the distinction between getarray_vector and getarray_outer for optimizations purposes. This would allow to move pass #1770 while making sure that if a distinction is made between vector and outer indexing in the future this would not result in suboptimal slicing (again in the sense of minimal data brought to memory) for outer indexing.

What is the status of this?

This is how I would sum it up:

1. vector indexing (fancy=True) can break dask in special cases. The easiest fix is to turn off optimization in those cases.
2. outer indexing would however work well in those special cases and benefit from much needed optimization for some problems. The easy fix in (1) would prevent this desired optimization in arrays supporting outer indexing, so it would be kind of important to allow outer indexing-specific optimizations.

The current PR proposes a simple solution that allows to specify fancy={'vector', 'outer', True} where fancy=True is the same as fancy='vector'. That was the simplest solution I could come with...

I encountered and fixed this issue in #2313 without realizing that this was already raised (and a fix proposed). Sorry about that.

In #2313 I prevent fusing these indices in all cases. I like this because it's simple and correct for all backends, but it does prevent the optimizations coming into effect for hdf5 where they're valid(?) (although maybe they shouldn't be). I'm not sure what's best here.

@jcrist
Sorry I did not reply to your comment in #2313 more promptly. I think your approach is the right one, at least for the time being. I'm going to think of possible optimizations that build upon your approach and see if I can come up with something that makes sense.

No worries. Do you want to still keep this PR open then?

No, not really. I'm closing it!