Hrm, what if slicing changes the type of the array?

For example, I could imagine that slicing some types of sparse matrices might result in other types of sparse matrices, or possibly even numpy arrays. Ideally we transform meta as we perform the operation, rather than just changing around the dtype. In this way our dask array's meta type will track what happens in the lazy chunks.

So for an operation like __getitem__ we might make a new index that was just like the input, but had a bunch of elements that were safe to apply against an empty array.

In [16]: index = (slice(0, 5), [1, 2, 3], None, slice(0, 2))

In [17]: x = np.empty((0, 0, 0))

In [18]: index2 = (slice(0, 0), [], None, slice(0, 0))  # like index, but empty

In [19]: x[index2]
Out[19]: array([], shape=(0, 0, 1, 0), dtype=float64)

However I'm not sure how best to handle integer elements in the index. I'm not sure that there is an empty-safe way to handle this.

Yes, I think that will need some more careful thinking. We were aware of this kind of situation, unless I misunderstood, this case is analogous to the case of a Dask array containing multiple arrays of different types, which was mentioned by you in #2977:

The plan I've proposed above doesn't support the mixed-array case (where we have a mix of dense and sparse arrays in the same dask array for example), but I suggest that we start there. I think that it will be a good and productive first step.

I think that this concern is independent of mixed chunks.

Lets say that some CSR matrix type exists, and I have a dask array full of them.

When I do something like x[0, :] I pull out a row. It could be that the return type of that object is no longer a CSR matrix, but some other form (maybe COO). The output still has a uniform chunk type, it's just a new type.

I think I get it now. I'll try to think how to handle it.

The best I could do so far was:

meta = self._meta[index].astype(self.dtype)
meta.shape = (0,) * len(chunks)

That second line probably suffices to keep the shape as we need it, unless the result should be a 0d array (for such case, I don't have a solution yet). But if the index could change the type of the matrix, isn't there a case that the actual array size will matter and not change it if the matrix is empty, i.e., after slicing, data and _meta be different anyway?

To be honest, at this point I'm starting to ask myself if this is really a good design. There are so many ifs and buts that I think it may immediately render the code unmaintainable.

Oh, I forgot to mention that the code above works only for slicing, indexing would need to be a separate case.

To be honest, at this point I'm starting to ask myself if this is really a good design.

It's a good question to ask. And thank you for going through the experiment that helps us to understand whether or not this is the right approach.

There are so many ifs and buts that I think it may immediately render the code unmaintainable.

We've gone through this process before both with dtype and with shape. It's definitely a challenge.

My guess is that we'll be able to find a solution for most cases. For the more difficult cases there are people watching this PR that have deeper NumPy experience than you or I do (people like @shoyer @jakirkham @hameerabbasi and @TomAugspurger ) and they can help once we identify the difficult issues.

At this point we've identified that slicing is going to be difficult. Lets move on to other operations before diving too deeply here, and maybe one of the people above can offer a suggestion based on the description above.

We've gone through this process before both with dtype and with shape. It's definitely a challenge.

I understand this, and I'm not even worried about NumPy, but answering the question of what could 3rd-party libraries do to their own types in every possible case may be untenable. I'm thinking that waiting for empty_like() to be available may be the right choice, since we could then just recreate _meta with the right size every time something may have changed, rather than spending so much time in handling all possible cases.

At this point we've identified that slicing is going to be difficult. Lets move on to other operations before diving too deeply here, and maybe one of the people above can offer a suggestion based on the description above.

I ended up on slicing, because fixing the issue below (in partial_reduce()) depends on slicing. At some point things fail because the _meta array has the wrong shape. I suspect this will be a requirement for other pieces of the code as well, since it's one of the most basic operations.

@mrocklin I pushed some new commits, one of them attempts to handle slicing of _meta. This has been capable of handling all the cases that previously failed in existing tests, but there's always a chance that I've missed cases that I couldn't test for.

However I'm not sure how best to handle integer elements in the index. I'm not sure that there is an empty-safe way to handle this.

There is, just prepend Ellipsis to the list of indexes.

>>> import numpy as np
>>> np.float64(0)[..., None, None]
array([[0.]])

And if you want to guarantee that the output is an array, even if 0-D, postpend with Ellipsis:

>>> np.array([[0]])[0, 0, ...]
array(0)

I think that ideally we would call

func(x._meta, keepdims=keepdims)

That way if the shape or type of the array changes we'll learn about it. Changing the dtype alone probably isn't enough.

Good catch, I'll try to do that.

It ended up being much more complicated to handle reductions, but just like slicing, the commit that I pushed to handle this case, seems to handle all cases I've seen in existing tests.

This function probably needs all of the inputs, not just the arrays. Maybe something like the following:

args_meta = [arg._meta if hasattr(arg, '_meta') else arg for arg in args]
kwargs_meta = {k: v._meta if hasattr(v, '_meta') else v for k, v in kwargs.items()}
func = func(*args_meta, **kwargs_meta)

This happens because sometimes people put non-arrays in map_blocks or other blockwise operations like the following:

x.map_blocks(func, y, 123, option='foo')

In these cases it's important to keep the extra args/keyword args

I was aware my handling wasn't complete, but I wasn't sure what else was missing there, thanks for pointing that out. And to be honest, I wouldn't be surprised if your changes don't make it complete either, there are really many things that can happen in this blockwise() function.

Yes, it's at the core of many calls and is the frequent target of many many feature requests. Keeping it as simple as it is simple has been a challenge.

(not that I think it's simple)

Regardless, thanks for being willing to dive into it. I acknowledge that it's not trivial.

I did the changes you suggested, it mostly works but breaks test_argwhere_str and test_count_nonzero_str. The issue is with passing **kwargs, after a lot of unsuccessful digging, I think what happens is the 'concatenate' argument is passed to numpy.vectorize() that doesn't expect it. I will come back to check this later.

I'll try to take a look at this later today as well

It's still difficult to handle all cases when we pass kwargs, I'm still working on fixing this.

I would expect this to be else arg instead of else np.empty((0,)). What if the user provides an argument that is not an array, like a string?

def f(a, b):
    if b == 'add':
        return a + 1
    elif b == 'sub':
        return a - 1
    else:
        raise ValueError(b)
    
blockwise(f, 'i', x, 'i', 'add', None)

At first, I'm trying to get current tests to work, but yes, this is a valid concern. In fact, I've already seen some issues with that part of the code, I'll try to fix what I can before moving on to more complicated cases.

In my latest commit I fixed the condition to else arg, as @mrocklin suggested.

I've added a second try/except to attempt passing **kwargs and if there's an exception, don't. I noticed the only case this happens is in the test that @mrocklin added here: https://github.com/dask/dask/pull/4543/files#diff-8f3045290928cb33f56d0b894e765b88R141

This happens because blockwise() is always called passing a concatenate=True in map_blocks(), as seen here https://github.com/dask/dask/blob/master/dask/array/core.py#L484, and sparse.COO.from_numpy() doesn't have this argument. What would be the correct solution here?

I'm confused about the objective of reduced_meta here. Can you help me understand this?

Also, do we need to call blockwise? Or are you just calling blockwise in order to get the meta object?

I'm trying to ensure that _meta goes through all the same steps as the main array, since I assume it's not possible to predict how every single step may change the arrays and their data type.

Can we use tmp._meta instead here perhaps?

In general calling compute can be quite complex to do. We engage the full optimization system. If we're using the distributed scheudler then we trigger network communications. If we're using Dask in an asynchronous mode then things will probably block in bad ways. If we need to call compute then we can turn all of this stuff off during that compute call, but I'd prefer to avoid it if we can.

I totally get your point now @mrocklin. However, I don't immediately see a different way here. What I could do is not "reduce" _meta, but I presume that could lead to different results for complex types. Do you have any suggestions on how to better handle it here?

Calling .compute() anywhere is somewhat concerning within dask array code. Ideally things would stay lazy.

reduced_meta is _meta wrapped in a Dask array, the one returned by blockwise(). So far, this is how I managed to ensure we go through all steps.

Can I ask for information about this? Why did we need this?

Or, maybe more directly, the fact that we need this is somewhat concerning. I'm curious about what the problem is so that maybe someone can find another approach.

I don't like this either, but what happens here is that in partial_reduce() I'm calling func() (or meta_func() which doesn't _concatenate2()), and depending on the function in question, it may accept different arguments. In particular the meta argument had to be added because functions chunk/combine/aggregate functions may return dicts, which aren't valid arrays to reduce.

Yeah I agree with both of you, we should find another approach here if possible. Do you have a couple of examples where this is needed so we can formulate and try out some alternatives?

Alright, I pushed commit ce7e5dd to get rid of (part of) that. We still need some exception handling for some of the reduction functions (e.g., mean_chunk, mean_combine, etc.) for which I added a meta= keyword argument to prevent them from returning dicts that we don't want for _meta construction, but most of them won't have such keyword. The second case is for ValueError exception, like ufuncs that can't be executed on empty arrays.

I now honestly don't know if I had a legitimate case where I've failure due to dtype= and keepdims=. The two possibilities I see now for that are:

1. I had some bug in another part of the code during my development process; or
2. I'm really going insane.

Please let me know if the solution on the aforementioned commit seems acceptable. But IMHO, it seems acceptable.

I'm going to try to defer to @jakirkham on this :)

Sorry I'm not totally following what is happening in this block. Could you please provide a summary, @pentschev?

Sure, what happens here is very simple, basically creating a new_index in this order:

1. Check all elements of the index2 tuple, if it's not a tuple, make it a slice(0, 0);
2. For elements that are not tuples (e.g., a single index), make it an Ellipsis if the index is valid, otherwise make it None.

This way, I convert all cases to 0-sized slices. For example:

• [1, slice(1:2), :] becomes shape (0, 0, ...)
• [None, 1, None] becomes shape (1, 0, 1), the 1s that arise from None should probably still be fixed, even though it may not be a problem generally

Thanks for providing more info. I'm still confused unfortunately. Have a couple of follow-up questions.

What is the difference between using slice(0, 0) and Ellipsis on _meta?

Also, I thought None would add singleton dimensions. Why do we want that for this use case?

What is the difference between using slice(0, 0) and Ellipsis on _meta?

I think they're equivalent, I probably started using Ellipsis and forgot to fix that. Plus, having two consecutive Ellipsis isn't allowed, so that's wrong for sure, I'll fix it.

Also, I thought None would add singleton dimensions. Why do we want that for this use case?

It does add a singleton dimension. We would want one if the user requests so, as per @mrocklin's previous example here.

Just to be clear here, I was attempting to cover all cases I could think of. Maybe not all of them would make sense in a realistic context, but are technically allowed, so we should cover them.

I have a slight preference to inline not NEP_18_ACTIVE here (and maybe the reason too). It's shorter and doesn't requrie the reader to go find the definition of missing_arrfunc_cond.

Just a style preference though. Not a big deal.

It definitely isn't shorter, due to the long reason string. But I have no strong preference, so I'm just gonna go ahead and inline things.

Assuming that assert_eq checks that type(x.compute()) == type(x._meta) then would it work to use assert_eq(a, a) here instead?

Yes, I feared that it would look ridiculous and someone would end up removing that at one point. I'm going to do that and add a comment that it's done on purpose.

We use assert_eq(a, a) in other places as well. I recommend just doing that. I think that most of the maintainers know that this is common, just because of all of the other checks that assert_eq runs.