Support Everything that XArray Expects

single argument ufuncs (sin, exp, etc.) and ufunc like functions (pd.isnull, notnull, astype, around, isclose)
broadcasting binary operations (e.g., ndarray arithmetic)
three argument version of where (preferably with broadcasting)
aggregation methods
- max/min/sum/prod
- argmax/argmin/std/var
- nan-skipping aggregations (e.g., nanmin, nanmax, nanprod, nansum)
indexing with integer arrays, booleans, slices, None
transpose
indexing:
- basic indexing (int/slice)
- outer indexing for a single array
- outer indexing (int, slice and 1d integer arrays separately applied to each axis)
- vectorized indexing (integer arrays with broadcasting, like NumPy)
broadcast_to (NumPy 1.10)
concatenate and stack (NumPy 1.10)

I updated the check-list (apparently I'm an "Owner" for pydata) to drop "insert/fancy indexing" and add a three-item checklist for indexing instead.

Do we have partial XArray compatibility or does XArray fail with sparse arrays outright?

@hameerabbasi

@hameerabbasi Not yet, see pydata/xarray#1375. But sparse is close enough that it might be worth trying soon.

From a usability perspective, the most helpful additional feature would be nan-skipping aggregations. I know sparse will probably rarely be used with NaNs, but xarray's aggregation methods default to skipna=True so this would be valuable for consistency.

In addition to the direct uses, xarray uses where and outer/vectorized indexing internally primarily for alignment with join='outer' (e.g., in concat or merge). This sort of alignment would require making sparse arrays dense, unless they have a fill value of NaN rather than 0. So these operations are less essential than they might otherwise seem.

I just tested with your example code, NaN-skipping aggregations use ufuncs with reduce, so they should already be supported in theory.

>>> class A(np.ndarray):
...     def __array_ufunc__(self, *args, **kwargs):
...         print('ufunc', args, kwargs)
...         
>>> np.nansum(np.arange(5).view(A))
ufunc (<ufunc 'add'>, 'reduce', A([0, 1, 2, 3, 4])) {'axis': None, 'dtype': None, 'keepdims': False}

Although I have absolutely now idea how to discern that this was NaN skipping from the args/kwargs.

Ah, it was a dtype issue.

>>> np.nansum(np.arange(5, dtype=np.float).view(A))
ufunc (<ufunc 'isnan'>, '__call__', A([0., 1., 2., 3., 4.])) {}
ufunc (<ufunc 'add'>, 'reduce', A([0., 1., 2., 3., 4.])) {'axis': None, 'dtype': None, 'keepdims': False}

So we're good to go on those already. :-)

Edit: Concrete example:

>>> ar = sparse.DOK((2, 2))
>>> ar[1, 1] = np.nan
>>> s = sparse.COO(ar)
>>> np.nansum(s)
0.0

@njsmith

FWIW this is a great demonstration of the value of NumPy's use of protocols. It would be great if we could get an alternative ndarray implementation, sparse, to work with a complex downstream user of NumPy code, XArray, without either library having to explicitly know about the other. cc @njsmith

I was also considering sparse arrays with different fill values. However, some operations (such as dot) blow up in complexity, and on others, it's fairly easy to support.

NumPy's nan-aggregations work, but only by converting sparse arrays into numpy arrays, inside np.lib.nanfunctions._replace_nan:

import numpy as np
import sparse

class NoncoercibleCOO(sparse.COO):
  def __array__(self, *args, **kwargs):
    raise Exception('cannot convert to numpy')
  
ar = sparse.DOK((2, 2))
ar[1, 1] = np.nan
s = NoncoercibleCOO(ar)
np.nansum(s)

This results in:

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-13-481c3bc38f6d> in <module>()
      9 ar[1, 1] = np.nan
     10 s = NoncoercibleCOO(ar)
---> 11 np.nansum(s)

/usr/local/lib/python3.6/dist-packages/numpy/lib/nanfunctions.py in nansum(a, axis, dtype, out, keepdims)
    579 
    580     """
--> 581     a, mask = _replace_nan(a, 0)
    582     return np.sum(a, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
    583 

/usr/local/lib/python3.6/dist-packages/numpy/lib/nanfunctions.py in _replace_nan(a, val)
     62 
     63     """
---> 64     a = np.array(a, subok=True, copy=True)
     65 
     66     if a.dtype == np.object_:

<ipython-input-13-481c3bc38f6d> in __array__(self, *args, **kwargs)
      4 class NoncoercibleCOO(sparse.COO):
      5   def __array__(self, *args, **kwargs):
----> 6     raise Exception('cannot convert to numpy')
      7 
      8 ar = sparse.DOK((2, 2))

Exception: cannot convert to numpy

So this really isn't a good solution yet :).

You might actually view this as the flip side of implementing protocols (like __array__): you get default implementations, but they aren't always good ones and it's not always obvious when that's the case.

Is this already reported upstream in NumPy?

Hmm. So for this what we really want is something like: implement nansum using which=, and make reductions support which=?

I was also considering sparse arrays with different fill values. However, some operations (such as dot) blow up in complexity, and on others, it's fairly easy to support.

I think this could be quite interesting indeed. NaN (which pandas's SparseArray uses) is the most obvious other default value to support. Even dot could be make to work if you make a nan-skipping version (nandot, nanmatmul?).

Is this already reported upstream in NumPy?

I don't know if this really qualifies as a bug, so much as an unintended consequence. I think the easiest fix would be to replace np.copyto() in _replace_nan with a protocol dispatching version of the three argument variant of np.where().

A simpler option would be to add skipnan= in ufunc reduce and accumulate methods.

Edit: And reduceat.

Edit: Or maybe in ufuncs themselves. ufunc(NaN, NaN) = NaN, otherwise, return the non-NaN item.

I think this could be quite interesting indeed. NaN (which pandas's SparseArray uses) is the most obvious other default value to support. Even dot could be make to work if you make a nan-skipping version (nandot, nanmatmul?).

I meant supporting arbitrary fill values. This would make many dense operations automatically possible, e.g. y = x + 5, y would have a fill value of 5.

You might, for example, do x[x == 0] = np.nan (please note, multidimensional boolean indexing not currently supported) and it would work.

I meant supporting arbitrary fill values. This would make many dense operations automatically possible, e.g. y = x + 5, y would have a fill value of 5.

Yes, I see the elegance in that.

NaN has a nice property (like zero) that NaN * anything -> NaN. This could make some operations easier to implement if the general fill value is not viable.

I do think that some static information like this is likely to be desirable when we start profiling. Also, I suspect that most use of elemwise will be by expert users (like us when we're using elemwise internally), and so providing this information might be easy.

I agree 100% with this. But (at the risk of repeating myself), this ignores the NaN/inf problem with mul. While doing this, we would have to filter/check for those, or risk inconsistency in these cases as compared to Numpy. Also, it doesn't work for arbitrary fill values. Also that refactor I did in #84 will need to be partially reversed.

Alternativley, what is our current support for operations that have non-zero return values when one of the inputs is zero? My guess was that we erred. If this is the case then can we just ignore the situation where we produce values on zero-valued inputs? Am I understanding things correctly?

Currently, we err iff func(all_zeros) will produce a nonzero. If zeros/nonzeros are mixed (like in add), we don't err but compute func(zero, other.data) etc. Here, I've just generalized that behavior to multiple operands. The plan is to make it so that func(all_fillvalues) becomes a fill-value for the return array.

@shoyer

@shoyer I don't suppose you happen to have a TestCase for what an XArray container class should support do you? If so it would be good to import that into our test suite to track progress in the future between the projects.

No, not yet, unfortunately. Historically, the API for xarray container classes has evolved in a dynamic way based on what we needed for numpy + dask. We'll need to do some internal clean-up before we can support sparse arrays well.

…

On Wed, Feb 21, 2018 at 6:24 AM Matthew Rocklin ***@***.***> wrote: @shoyer <https://github.com/shoyer> I don't suppose you happen to have a TestCase for what an XArray container class should support do you? If so it would be good to import that into our test suite to track progress in the future between the projects. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABKS1nYsOS6VNgx8ninOXeEf-QPNHwL-ks5tXCc5gaJpZM4M-1J4> .

Feel free to ping me or involve me in the project. :-)

Three-argument version of where and NaN-skipping aggregations are in!

@shoyer

@shoyer what is the right thing for folks to push on here if they have time? The multipledispatch VarArgs conversation? Internal fixes to XArray? A test suite?

My intuition would be the multipledispatch VarArgs conversation. We should do things properly. 😄 I don't mind it taking time, it'll save time and headache down the road. Of course, if the consensus is on another solution, I'd happily work on that, too.

VarArgs will be needed for finishing this in xarray, but really only for a few functions (concatenate and stack). So I think we could simultaneously work on:

Getting VarArgs of some form into multipledispatch.
Porting as much of the existing xarray dispatch machinery for NumPy and Dask to use multipledispatch. This has at least two components: the wrappers in xarray/core/duck_array_ops.py, and implementing __array_ufunc__ (__array_ufunc__ for xarray xarray#1617).

@mrocklin

Since the final review deadline for the papers in SciPy 2018 is June 18, I thought I'd revive this thread. I'd love to show off Dask/Sparse and XArray/Sparse integration in the paper and at SciPy 2018.

cc @mrocklin @shoyer If you need anything done on my end, or anything at all you think I'm capable of that'd speed this along (Including contributing to XArray or Dask, I'd love to dip my feet in), don't hesitate to say so.

I looked into multipledispatch but VarArgs seems like something better done by someone else. arrayish is already up, if any changes are required there, let me know.

@shoyer

Experiments with dask.array and sparse would be welcome. I think that this should work today (we do include some tests here in dask.array's CI). It might be interesting to do profile a few workflows. I'll think a bit about what these might be.

For XArray my guess is that the right thing to do is to follow @shoyer 's two paths listed just above. Looking at duck_array_ops.py and how that functionality gets used within XArray seems like a sensible place to start, though @shoyer may have other suggestions (I don't know this codebase as well). My guess is that there some work to get through here, but probably nothing that is too controversial.

It also looks like the conversation in mrocklin/multipledispatch#72 ended in a nice place. Again, it seems like consensus was reached and now this requires some hopefully straightforward work.

I’ve implemented outer indexing for a single array, but not multiple. When NumPy implements oindex and vindex, we’ll follow suit. I somehow don’t feel like implementing every single fancy index case with its weird behaviour.

Not sure if this is the right place to mention, but it doesn't appear elsewhere in issues/PRs on the repo.

I ran into the following trying to call .shift() on a sparse.COO-backed DataArray:

../../.local/lib/python3.8/site-packages/xarray/core/dataarray.py:3077: in shift
    variable = self.variable.shift(
../../.local/lib/python3.8/site-packages/xarray/core/variable.py:1205: in shift
    result = result._shift_one_dim(dim, count, fill_value=fill_value)
../../.local/lib/python3.8/site-packages/xarray/core/variable.py:1166: in _shift_one_dim
    data = duck_array_ops.pad(
../../.local/lib/python3.8/site-packages/xarray/core/duck_array_ops.py:56: in f
    return wrapped(*args, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

args = (<COO: shape=(2, 5), dtype=float64, nnz=4, fill_value=nan>, [(1, 0), (0, 0)]), kwargs = {'constant_values': nan, 'mode': 'constant'}
relevant_args = (<COO: shape=(2, 5), dtype=float64, nnz=4, fill_value=nan>,)

>   ???
E   TypeError: no implementation found for 'numpy.pad' on types that implement __array_function__: [<class 'sparse._coo.core.COO'>]

<__array_function__ internals>:5: TypeError

There were other similar cases. So far I've noticed:

numpy.flip for DataArray.bfill
numpy.nancumprod for DataArray.cumprod
numpy.pad for DataArray.shift

@khaeru

@khaeru Do you have use-cases for these for sparse arrays? If so, please raise separate issues for each missing function.

In a lot of cases, one can’t achieve better performance on these than dense arrays.

@hameerabbasi

@hameerabbasi the concern is less performance than being able to use the DataArray API as it is defined without a bunch of extra wrapper code—what I understood as the topic of this issue, given the title.

My current mitigation is to have a subclass that allows some of these calls by converting to dense/numpy.ndarray-backed DataArray, performing the operation, then re-converting to sparse.COO-backed. I know that subclassing the xarray classes is not encouraged, I would be very happy to not have to do that :)

I'll do as you suggest and open an issue for one example.

Navigation Menu

Repository navigation

Support Everything that XArray Expects #1

20 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Support Everything that XArray Expects #1

Description

Activity

shoyer commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

shoyer commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

mrocklin commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

shoyer commented on Jan 30, 2018

mrocklin commented on Jan 30, 2018

njsmith commented on Jan 30, 2018

shoyer commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

hameerabbasi commented on Jan 30, 2018

shoyer commented on Jan 30, 2018

20 remaining items

hameerabbasi commented on Feb 10, 2018

mrocklin commented on Feb 21, 2018

shoyer commented on Feb 21, 2018

hameerabbasi commented on Feb 21, 2018

hameerabbasi commented on Feb 22, 2018

mrocklin commented on Mar 4, 2018

hameerabbasi commented on Mar 4, 2018

shoyer commented on Mar 4, 2018

hameerabbasi commented on May 18, 2018

mrocklin commented on May 18, 2018

mrocklin commented on May 18, 2018

hameerabbasi commented on Aug 4, 2018

khaeru commented on Mar 7, 2021

hameerabbasi commented on Mar 7, 2021

khaeru commented on Mar 7, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions