issue_comments

You are right about the behavior of the code. I don't see any way to enhance that in the general case.

Maybe, in your case, rechunking before interpolating might be a good idea

Yep this look good to me too. Thank you

I thing this has been solve with v0.16.0 and https://github.com/pydata/xarray/pull/4155

And I forgot to take into account that your interpolation only need 48² points of the input array, so the input array will be reduced at the start of the process (you can replace every 100 by 48 in my previous answers)

I forgot to take into account that the interpolations are orthogonal So in sequential we are doing 2 interpolation first x then y In parallel we do the same: The fist interpolation will have 20 000 tasks, each task will have the totality of the input array, and compute an interpolation of 5 point of the output (x) producing an array of 5x100 per task or 100 000x100 full result as an intermediate array. The second interpolation will have 20 000² tasks each task will have a block of 5x100 point of the intermediate array and compute an interpolation on 5 point of the output (y) resulting in a 5² array per task and the 100 000² full result.

So plenty of room for overhead...

In your case, each task (20 000²) will have the entire input (100²), and interpolate a few points (5²).

Maybe the overhead comes with duplicating the input array 20 000² times, maybe it comes with the fact that you are doing 20 000² small interpolation instead of 1 big interpolation

Do this answer your question?

If the input array is chunked in the interpolated dimension, the chunks will be merged during the interpolation.

This may induce a large memory cost at some point, but I do not know how to avoid it...

@cyhsu I can answer this question.

For best performance you should chunk the input array on the non interpolated dimensions and chunk the destination. Aka :

``` datax = xr.DataArray(data=np.arange(0, 4), coords={"x": np.linspace(0, 1, 4)}, dims="x") datay = xr.DataArray(data=da.from_array(np.arange(0, 4), chunks=2), coords={"y": np.linspace(0, 1, 4)}, dims="y") data = datax * datay

x = xr.DataArray(data = da.from_array(np.linspace(0,1), chunks=2), dims='x')

res = data.interp(x=x) ```

You're welcome :)

Hi.

I agree, part of this work might belong in dask. But I don't know dask internals enough to go there. In this case, everything was already in place.

Moreover I do think that there is room for optimization. In particular, in this implementation, the work is distributed along chunks corresponding to destination. This means that one may have big intermediate array. For example interpolating one value in a chunked vector will load the full vector in memory (first localization aside). In my previous implementation (and uglier), the interpolation was done with the chunks of the starting array. This might be a better choice sometimes.

OK, I'm happy with the results now (better than my first submission of course).

I did not add so much tests since the result replace what was done before, thus the previous tests applies.

I'm going for some holidays so I won't work that much for the time being. But I'll be able to answer any questions.

Thanks for the reviewing and pushing me into doing a much better job.

FYI, don't merge yet. I fixed a bug today, but did not push it. And there is some work to do on the testing side.

Since I was on it, I extended the decomposition of orthogonal interpolation. If you want I can break this into two PR.

Guys, I got it. I managed to use da.blockwise which allows me to overcome all the previous limitations.

The result is much more simple, much more reliable.

While at it, I added the missing bit to make it work with cubic or quadratic method. I'm not touching the code anymore, waiting for review.

@fujiisoup I managed to implement the support of unsorted interpolation.

Also, I reworked the tests, I now test for much more situations.

Thanks, @pums974 for this update. I left some comments.

Can you add some tests for more edge cases? Something we may want to check would be

* scalar interpolation
* interpolation into an unsorted dimension (e.g., `da.interp(x=[0, 3, 2])`)

Your welcome, thanks for the feedback

• scalar interpolation: you mean a test like test_interpolate_nd_scalar but between chunks ?
• Unsorted interpolation: As I said, I did not looked into it, This need some work, presumably an argsort at the begining in order to interpolate in a sorted dimension and reorder the result into the requested order. I'm trying to implement it, but this seems a bit more challenging than I thought.

@fujiisoup I managed to solve the issues you raised about AttributeError: 'memoryview' object has no attribute 'dtype' This was due to datetime index. Using an IndexVariable seems to solve it

Also I realize that for 1d interpolation cubic and quadratic method are allowed which may not give the same result with chunked data (or even crash if there is not enough data in the chunked direction). Now they are forbidden

I mean, in this case you have to interpolate in another direction. You cannot consider having a 1d function.

Thanks, That's weird, I have no problem in mine... What are your versions of dask and numpy ?

As for implementing this in dask, you may be right, it probably belong there, But I am even less use to their code base, and have no clue where to put it.

And for unsorted destination, that's something I didn't think about. maybe we can add an argsort at the beggining.

ok, but what about python res = data.interp(y=0.5)

Hum, ok, but I don't see how it would work if all points are between chunks (see my second example)

No problem, we are all very busy. But thanks for your message.

On my computer it passes pytest: ``` $> pytest . ======================= test session starts ================================= platform linux -- Python 3.8.3, pytest-5.4.3, py-1.8.2, pluggy-0.13.1 [...] ===== 3822 passed, 2710 skipped, 77 xfailed, 24 xpassed, 32 warnings in 48.25s ========

$> pip freeze appdirs==1.4.4 attrs==19.3.0 black==19.10b0 click==7.1.2 dask==2.18.1 flake8==3.8.3 isort==4.3.21 mccabe==0.6.1 more-itertools==8.4.0 numpy==1.18.5 packaging==20.4 pandas==1.0.4 pathspec==0.8.0 pluggy==0.13.1 py==1.8.2 pycodestyle==2.6.0 pyflakes==2.2.0 pyparsing==2.4.7 pytest==5.4.3 python-dateutil==2.8.1 pytz==2020.1 PyYAML==5.3.1 regex==2020.6.8 scipy==1.4.1 six==1.15.0 toml==0.10.1 toolz==0.10.0 typed-ast==1.4.1 wcwidth==0.2.4 -e git+git@github.com:pums974/xarray.git@c47a1d5d8fd7ca401a0dddea67574af00c4d8e3b#egg=xarray ```

I also want to alert that my version does not work with "advanced interpolation" (as shown in the xarray documentation) Also, my version cannot be used to make interpolate_na work with chunked data

When using the current official interp function (xarray v0.15.1), the code: ```python datax = xr.DataArray(data=da.from_array(np.arange(0, 4), chunks=2), coords={"x": np.linspace(0, 1, 4)}, dims="x") datay = xr.DataArray(data=da.from_array(np.arange(0, 4), chunks=2), coords={"y": np.linspace(0, 1, 4)}, dims="y") data = datax * datay

# both of these interp call fails
res = datax.interp(x=np.linspace(0, 1))
print(res.load())

res = data.interp(x=np.linspace(0, 1), y=0.5)
print(res.load())

`` fails withNotImplementedError: Chunking along the dimension to be interpolated (0) is not yet supported.`

but succeed with the monkey patched version

EDIT : added the second interp in order to show more general use

Thanks, I'll try to make you tests and example quickly.

As for map_overlap, I tried, but something went wrong (I don't remember what though) . Maybe I didn't try enough.

EDIT: I cannot use map_overlap because it doesn't pass *args through to map_blocks

Any feedback ? Is my issue usefull ? Should I write a merge request ?