nd-rolling

[fujiisoup]

• Closes Convolution operation #4196
• Tests added
• Passes isort -rc . && black . && mypy . && flake8
• User visible changes (including notable bug fixes) are documented in whats-new.rst
• New functions/methods are listed in api.rst

I noticed that the implementation of nd-rolling is straightforward.
The core part is implemented but I am wondering what the best API is, with keeping it backward-compatible.

Obviously, it is basically should look like

da.rolling(x=3, y=3).mean()

A problem is other parameters, centers and min_periods.
In principle, they can depend on dimension.
For example, we can have center=True only for x but not for y.

So, maybe we allow dictionary for them?

da.rolling(x=3, y=3, center={'x': True, 'y': False}, min_periods={'x': 1, 'y': None}).mean()

The same thing happens for .construct method.

da.rolling(x=3, y=3).construct(x='x_window', y='y_window', stride={'x': 2, 'y': 1})

I'm afraid if this dictionary argument was a bit too redundant.

Does anyone have another idea?

[fujiisoup]

This is a basic test for the new nd-rolling

[max-sixty]

This looks very promising; I'm surprised it was possible without more code. I'm being slow, but where is the nd-rolling algo? I had thought bottleneck didn't support more than one dimension? https://bottleneck.readthedocs.io/en/latest/bottleneck.move.html, and that we'd have to implement our own in numbagg (which would be very possible)

[max-sixty]

Re the API, I think the dict is probably the best option, although it does complicate as the arguments become differently typed depending on one vs multiple dimensions.

One alternative is to allow fluent args, like:

(
    da
    .rolling(x=3, center=True)
    .rolling(y=5, min_periods=2)
    .mean()
)

...but does that then seem like the second rolling is operating on the result of the first?

[fujiisoup]

Hi @max-sixty

One alternative is to allow fluent args, like:
...but does that then seem like the second rolling is operating on the result of the first?

I couldn't think of it until just now. But yes, it sounds to me like a repeated rolling operation.

I'm being slow, but where is the nd-rolling algo? I had thought bottleneck didn't support more than one dimension?

No. With nd-rolling, we need to use numpy reductions.
Its skipna=True operation is currently slow, but it can be improved replacing nan before the stride-trick.

[fujiisoup]

Another API concern. We now use min_periods, in which we implicitly assume one-dimension cases.

With nd-dimension, I think min_counts argument is more appropriate like bottleneck, which will limit the lower bound of the number of missing entries in the n-dimensional window.

Even if we leave it, we may disallow nd-argument of min_periods but keep it a scalar.

[pep8speaks]

Hello @fujiisoup! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2020-08-07 22:59:00 UTC

[fujiisoup]

I hope it is true, though I didn't check it.
A suspicious part is da.concatenate, which I expect does not copy the original (strided-)array.

[fujiisoup]

I noticed that ghosting breaks my expectation. I need to update the algo here.

[dcherian]

I think the dictionary is OK. We can allow scalars for center, stride etc. in which case the scalar applies to all dimensions. This would preserve backward compatibility.

scalar for min_counts also seems OK to me. We can update it if someone requests it.

[fujiisoup]

I think now it is ready for review, though I'm sure tests miss a lot of edge cases.
Maybe we can fix them if pointed out.

[fujiisoup]

A possible improvement will be nan-reduction methods for nd-rolling.
Currently, we just use numpy nan-reductions, which is memory consuming for strided arrays.

This issue can be solved by replacing nan by appropriate values and applying nonnan-reduction methods,
e.g.,

da.rolling(x=2, y=3).construct(x='xw', y='yw').sum(['xw', 'yw'])

should be the same with

da.rolling(x=2, y=3).construct(x='xw', y='yw', fill_value=0).sum(['xw', 'yw'], skipna=False)

and the latter is much more memory efficient.

I'd like to leave this improvement to future PR.

[fujiisoup]

I got an error for typechecking, only in CI but not in local, from the code that I didn't change.

[fujiisoup]

Could anyone kindly review this?

[max-sixty]

Looks really good, great implementation imo

[max-sixty]

Types on the signature should change too (indeed, I'm surprised mypy doesn't fail)

[max-sixty]

Isn't this a list of ints? Is it going to do getattr(da, 3)?

[max-sixty]

This is very clever! I spent a while trying to figure out how it works...

[max-sixty]

Will this test will fail for something like std?
IIUC iteratively applying the function to the array for each dim only works for some functions.

Which is probably fine, no need to write more tests imo, though maybe worth making a note

[fujiisoup]

In this case, many other functions should have the same result.
If there are some nans in random positions, the result should be different.
I added a few more tests and include function std, but I don't think it is necessary to test all nanops.

[fujiisoup]

Thanks @max-sixty for the review ;)
I'll work for the update in a few days.

[max-sixty]

Suggested change

	# Do it in the oposite order
	# Do it in the opposite order

[max-sixty]

I'm still a bit confused. How does the "roll over each dimension in turn" approach equal the "roll over both dimension together" approach with a function like std? Here's a proposed counter example:

import xarray as xr
import numpy as np

da = xr.DataArray(np.asarray([[0,10,0],[0,10,0], [0,10,0]]), dims=list('xy'))
print(da)

<xarray.DataArray (x: 3, y: 3)>
array([[ 0, 10,  0],
       [ 0, 10,  0],
       [ 0, 10,  0]])
Dimensions without coordinates: x, y

x_std = da.rolling(dict(x=2)).std()
print(x_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
Dimensions without coordinates: x, y

x_then_y_std = x_std.rolling(dict(y=2)).std()
print(x_then_y_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [nan,  0.,  0.],
       [nan,  0.,  0.]])
Dimensions without coordinates: x, y

combined_std = da.rolling(dict(x=2, y=2)).std()
print(combined_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [nan,  5.,  5.],
       [nan,  5.,  5.]])
Dimensions without coordinates: x, y

[dcherian]

I think it's more like

da.rolling(x=2).construct("xr").rolling(y=2).construct("yr").std(["xr", "yr"])

but that gives

array([[0., 5., 5.],
       [0., 5., 5.],
       [0., 5., 5.]])

[fujiisoup]

Thanks @max-sixty .

You are completely correct.
As the test pass, I was fooling myself.

The reason was that the dataset I was using for the test does not have time and x simultaneously.
So I was not testing the 2d-rolling but just 1d-rolling.

Fixed. Now it correctly fails for mean and std, but passes with max and sum

[max-sixty]

Great! Thanks @fujiisoup !

Ready to go from my POV

[fujiisoup]

@max-sixty
thanks for the review.
merged