Indexing with arrays only seems to work with 1D arrays

[jamesavery]

Small example:

a  = np.arange(10);
ix = a[::-1];
print(a[ix])    # works fine
b  = a[:,None]*a[None,:];

print(b[:,ix])  # Explodes in array_create.py, line 147
print(b[ix,:])  # Explodes in array_create.py, line 147
print(b[ix])    # Handed off to NumPy
print(b[:][ix]) # Handed off to NumPy
print(b.flatten()[b]) # 2D-index on flat array works fine

[jamesavery]

@madsbk I've written a user kernel that implements numpy.take for any axis (and any ndim for the index matrix too). How do I link it up with the actual indexing, so I can get A[:,ix,:], etc., to work in Bohrium?

[dionhaefner]

print(b.flatten()[b]) # 2D-index on flat array works fine

If that's the case, the solution could be in pure Python; just use the strides to convert the indices into absolute offsets, then index into the flat array. But I'm afraid any solution will have to be implemented inside Bohrium, as it's not possible to monkey-patch __getitem__ of bh.ndarray since it's a C extension type.

[jamesavery]

Same solution as for #586. New GATHER instruction that allows indexing over any axis.