Unit Test Results

See test report for an extended history of previous test failures. This is useful for diagnosing flaky tests.

      9 files  ±0        9 suites  ±0   3h 14m 2s ⏱️ -31s
 18 059 tests ±0   16 833 ✅  -  9   1 217 💤 ±0   9 ❌ + 9 
161 605 runs  ±0  149 403 ✅  - 81  12 121 💤 ±0  81 ❌ +81 

For more details on these failures, see this check.

Results for commit 196b85e. ± Comparison against base commit 37835c4.

Hi, thanks for digging so deeply in to this.

I had trouble understanding the description so here's what I have understood so far. In your example problem, indices is a list of 3 ndarrays with shapes shape=(1000, 61, 1, 1), shape=(1000, 1, 61, 1), shape=(1000, 1, 1, 61) so it's the case highlighted in this comment: https://github.com/dask/dask/pull/11625/files#r1901046262 . THe current algorithm operates on the full broadcasted versions of these 3 arrays, so it scales badly.

It would certainly be good to come up with a better algorithm for this case; I spent a bunch of time trying to do it but didn't succeed (at least not cleanly). However I don't understand what you've changed. Can you add some comments to the PR please explaining the logic.

FYI I developed #11625 looking at the benchmark in #10237. It would be good to not regress here.

import dask.array as da
import numpy as np
arr_1 = da.ones([6000, 6000], chunks=[3000, 3000])  # 36M elements in 4 chunks
idx = np.repeat(np.arange(0, 6000, 6), 1000)  # 1M index points
arr_2 = arr_1.vindex[idx, idx[::-1]]

In particular, writing code like this: mask = block_coords[dim] == coord inside a hot loop is not great for performance, which is why the current code, does an argsort outside the loop.

Not sure if he didn't leave for the week-end already, but @GwenaelCaer has created a documentation website with a lot of details on the optimization: https://odatis-public.gitlab-pages.ifremer.fr/fast_vindex/index.html

Does that answer your questions?

Edit: in particular, this might be interesting: https://odatis-public.gitlab-pages.ifremer.fr/fast_vindex/development/algorithm_step_by_step.html

ah thanks, that does help!

Conceptually, the trick is to completely avoid broadcasting, which I agree with. However, I'm not sure this implementation works well for the pointwise case (which is the case I heavily optimized for). Note that there is no advantage to be gained by not broadcasting in that case.

My conclusion earlier was that we'd need two algorithms to handle the two cases; in fact i remember thinking it's probably better to ask the user to use oindex directly. IIRC that algo is already decent (but I may be wrong)

My conclusion earlier was that we'd need two algorithms to handle the two cases

If I remember correctly, the proposed algorithm doesn't try to work with point-wise indexing, it is specialized on extracting regions (conceptually similar to what "vectorized multi-slices" would do). Which is why I suggested to add the proposed algorithm as a optimization for that case, and keep the existing algorithm for the other cases (not sure if the code already does that).

in fact i remember thinking it's probably better to ask the user to use oindex directly.

Not sure oindex would work for this particular use-case (cubes around potentially randomly distributed points), as the indices are only locally orthogonal.

it is specialized on extracting region

as the indices are only locally orthogonal.

Hmmm I don't understand --- the indices going in are broadcastable against each other. I'll need to spend some time thinking about it. This feels like a different indecing type almost

Which is why I suggested to add the proposed algorithm as a optimization for that case,

In that case the first step may be to clearly write out a condition that triggers the optimization after checking a condition. and then we could explore how to share code.

Hmmm I don't understand

It appears I might have confused myself. After clearing my head a bit I think you're right in that the arrays could be broadcasted against each other (to get three (n_points, nx, ny, nz) arrays, which we'd like to avoid to reduce the size in memory).

Fundamentally, the problem can be decomposed into n_points orthogonal indexing operations that we're trying to do at the same time to reduce the overhead.

In that case the first step may be to clearly write out a condition that triggers the optimization

Agreed, that makes sense to me.

Fundamentally, the problem can be decomposed into n_points orthogonal indexing operations that we're trying to do at the same time to reduce the overhead.

Yes i agree here. The challenge is figuring out when to do it, which is why maybe .oindex is a better place to put this kind of work (though i'm not sure if this counts as oindexing, i think it does)

I believe it's a hybrid between vectorized and orthogonal indexing: if we have the three arrays of shape (n_points, index_size), we want the last dimension to be orthogonal and the first dimension vectorized.

Thank you @keewis and @dcherian for the very interesting discussion. My algorithm is designed for a specific case with a hybrid of vectorized and orthogonal indexing, and it is therefore not suitable for point-wise indexing. It would therefore be useful to identify a condition that indicates when it could actually be applied.