Pinging @dask/array folks since I suspect someone in there will be able to review.

Is there a helper function to validate that a valid array graph is being constructed?

For e.g., with this branch, block sizes change after compute. This is resulting in some hard to debug errors.

before compute: (20, 3) | after compute: (17, 3)
before compute: (8, 3) | after compute: (11, 3)

Earlier, I had an issue where the shape was changing before/after compute.

It'd be nice if there was a function to check these things.

Is there a helper function to validate that a valid array graph is being constructed?

assert_eq should check the shape before and after compute.

assert_eq should check the shape before and after compute.

I don't think that will check that the individual chunks have the sizes that they're supposed to though (and which I also have my doubts about after reading the code), just that the pasted-together array does. I'm an infrequent dask developer so I don't know if there are any tools to do that, but if you wanted to write your own, you could probably do something with map_blocks and block_info to assert that the expected shape from block_info makes the actual shape for each chunk.

I've pushed an experiment to https://github.com/bmerry/sliding-window which greatly simplifies the implementation, and I think fixes a few bugs too (in particular, it should handle repeated axes). It passes the tests you wrote, but I can't say for sure whether it gets the newchunks correct. You're welcome to cherry-pick it.

Thanks @bmerry.

I have a simpler implementation too now by using symmetric overlaps and boundary padding.

I'll take a closer look at yours soon, it looks cleaner!

I don't think that will check that the individual chunks have the sizes that they're supposed to though (and which I also have my doubts about after reading the code), just that the pasted-together array does. I'm an infrequent dask developer so I don't know if there are any tools to do that, but if you wanted to write your own, you could probably do something with map_blocks and block_info to assert that the expected shape from block_info makes the actual shape for each chunk.

#7277 should help with this.

Very nice implementation @bmerry !

this version passes all of xarray's tests too so it should be ready to go.

It's looking good. The only thing I think still left is a test for a window size of 0, which is legal in numpy and has some slightly odd behaviour that I suspect might not work as is. It's also pretty useless (you get an array with a length-zero dimension) so if it's difficult to make work I don't think it would be unreasonable to just raise an exception stating that it's not supported by Dask.

The only thing I think still left is a test for a window size of 0, which is legal in numpy and has some slightly odd behaviour that I suspect might not work as is

I think we can leave that for someone who really needs it :) . I am raising ValueError now. I could change it to NotImplementedError if some one felt strongly about it.

FWIW I have been working on something similar, but only came across this PR today as I was cleaning up my code. I have uploaded it here https://github.com/dtpc/dask-windows

This PR is obviously must cleaner. I am only posting my version as it is based on skimage.util.view_as_windows and so supports a step argument as well, in case anybody is looking for that. :)

It seems like this is very close to being ready. I am not quite sure if the failures are known ones or not (ping @jrbourbeau who will know), but it's probably worth merging main in to see if this passes.

yay. tests pass!

Thanks for persevering! :)

Thanks @jsignell and @bmerry

Just testing out this new function. Are you aware that it can easily result in humongous chunks?
I got three years worth of wind-velocity measurements as lat/long images:

dataset:

Transforming it to a one-year sliding window results in three 500gb chunks.

sliding_window_view(dataset, 2920, axis=0):

safe_chunks are ((2920, 2928, 2920), (64,), (128,))

This is misleading unfortunately. it does rechunk to your window size, but after that it uses numpy's sliding_window_view to create a view that looks like a giant array 2928x64x128x2920 but is really using memory for 2928x64x128 = 200MB

@dcherian True, the last dimension is fake. Yet, if I try to compute anything with these "mega" chunks, serialization fails (probably because dask realizes early it can't hold that much memory). This effectively prevents me from using the sliding window function.

Anyways, I believe it should be possible to circumvent the rechunk() entirely (which also incurs additional overhead). I thought I had a working solution yesterday (by reusing the tasks from map_overlay and adjusting chunks afterwards by "creating" a new da.Array with identical name), however on closer inspection it fails on some edge cases.

My current approach is to call da.overlap.overlap and da.map_blocks directly. I will report back once I got a working proof of concept.

While it is possible to circumvent the additional rechunk(), the final chunks will of course be as big as the current solution. Yet, if you are interested in improving the current solution by removing the need for the rechunk, then here is a rough concept code for a single axis only:

org =  np.broadcast_to(np.arange(8)[:,None],(8,3))
X = da.from_array(org, chunks=(3, 3))

>>> org
array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2],
       [3, 3, 3],
       [4, 4, 4],
       [5, 5, 5],
       [6, 6, 6],
       [7, 7, 7]])

window_size = 3
X_overlap = da.overlap.overlap(X, depth={0: (0, window_size - 1)}, boundary='none')

>>> X_overlap.chunks
((5, 5, 2), (3,))

# this has to be calculated algorithmically in a full impl:
# subtract (window_size-1) from each of X_overlap's chunks along window axis, and set negative ones to zero afterwards
new_chunks = ((3,3,0), (3,), (window_size,))

mapped = da.map_blocks(
    sliding_window_view, 
    X_overlap,
    chunks=new_chunks,
    new_axis=2,
    meta=X_overlap,

    window_shape=window_size,
    axis=0
)
>>> mapped
dask.array<sliding_window_view, shape=(6, 3, 3), dtype=int64, chunksize=(3, 3, 3), chunktype=numpy.ndarray>

>>> mapped.chunks
((3, 3, 0), (3,), (3,))

# mapped has correct shape, but zero chunks will get computed anyways and will produce errors when computed
# in the simple / usual case where all zero chunks are at the end, we can just opt to not include them

fixed = da.Array(mapped.dask, mapped.name, ((3,3), (3,), (3,)), meta=mapped)
>>> fixed
dask.array<sliding_window_view, shape=(6, 3, 3), dtype=int64, chunksize=(3, 3, 3), chunktype=numpy.ndarray>

computed = fixed.compute()
>>> computed.shape
(6, 3, 3)
>>> computed
array([[[0, 1, 2],
        [0, 1, 2],
        [0, 1, 2]],

       [[1, 2, 3],
        [1, 2, 3],
        [1, 2, 3]],

       [[2, 3, 4],
        [2, 3, 4],
        [2, 3, 4]],

       [[3, 4, 5],
        [3, 4, 5],
        [3, 4, 5]],

       [[4, 5, 6],
        [4, 5, 6],
        [4, 5, 6]],

       [[5, 6, 7],
        [5, 6, 7],
        [5, 6, 7]]])

# for a robust solution also working on irregular chunks,  we have to "reindex" the chunks by creating a new task graph
mapped_irregular_chunks = ((3, 0, 3, 0), (3,), (3,))
fixed2_chunks = ((3,3), (3,), (3,))

new_tasks = {
('fixed2', 0,0,0): (mapped.name, 0,0,0),
('fixed2', 1,0,0): (mapped.name, 2,0,0),
}

fixed2 = da.Array(new_tasks + dict(mapped.dask), 'fixed2', fixed2_chunks, meta=mapped)