ENH: add function to get broadcast shape from a given set of shapes.

[madhulikajc]

Add new function numpy.broadcast_shape which takes tuples
for the shapes to be broadcast against each other.
Return the broadcasted shape as a tuple.
See #17217

[tillahoffmann]

I was incidentally having a look at the same thing and ended up with the following implementation. It doesn't use the builtin functionality of numpy, but it also doesn't create any temporary arrays. Might be another option?

def broadcast_shapes(*shapes):
    """
    Evaluate the shape due to broadcasting any number of shapes against each other.

    Parameters
    ----------
    `*shapes` : tuples
        The shapes to broadcast.

    Returns
    -------
    broadcasted_shape : tuple
        Shape of the array that would be obtained when broadcasting the shapes against each other.

    Examples
    --------
    >>> broadcast_shapes((2, 1), (1, 3))
    (2, 3)
    """
    shapes = [shape if isinstance(shape, tuple) else (shape,) for shape in shapes
              if shape is not None]
    maxdim = max(map(len, shapes))
    broadcast_shape = []
    for i, sizes in enumerate(it.zip_longest(*map(reversed, shapes), fillvalue=1)):
        sizes = set(sizes) | {1}
        if len(sizes) > 2:
            raise ValueError(f'unmatched dimensions at dimension {maxdim - i}')
        broadcast_shape.append(max(sizes))
    return tuple(reversed(broadcast_shape))

[eric-wieser]

@tillahoffmann, with #17535 (comment) the temporary arrays all have 0 bytes of data, so don't really cost anything

[tillahoffmann]

Excellent, let's avoid the duplication of code in that case. Looking forward to having this functionality available natively in numpy. @madhulikajc, thank you.

[madhulikajc]

Thank you @tillahoffmann and @eric-wieser. One thing I notice from @tillahoffmann's code is that an integer is accepted as an array shape. My code will also do this as np.empty() accepts integers and converts them into rank one arrays. But I do not explicitly test for this in my tests. I will add tests to check for this functionality.

[eric-wieser]

My only concern is that it might make sense to expose this more directly via a C API rather than with the strange workaround we've settled on here. That's always something that could be explored more in a future PR though.

[madhulikajc]

@eric-wieser, happy to explore exposing this functionality via a C API instead (i.e. the logic in PyArray_Broadcast). I had originally thought I would do this via a C API as well, but this ended up seeming simpler and cleaner (even though I agree it's a bit of an odd workaround). Let me know if you think we should proceed with this PR as-implemented and explore other options in future. If so, I can put together release notes and any other final housekeeping items. Thanks again for all your guidance here.

[seberg]

I like the plan to explore C (although I agree it is not absolutely necessary). Especially if we can consolidate a bit of code in PyArray_Broadcast to use the same in the end.

It is a bit unfortunate that I guess np.nditer is too complicated to consolidate here (broadcasting arrays can also be written as np.nditer(arrays, flags=["multi_index"]).shape). But that one supports too many other things we probably do not need (e.g. not broadcasting certain arrays, and nowadays, even broadcasting only specified axes).

[charris]

Needs a release note. Look in doc/release/upcoming_changes to see how it is done.

[mhvk]

On C vs python, I think it makes sense to merge this as is, since it is a helpful addition and any C implementation can just be slotted in anyway, and would then have the tests all ready.

[seberg]

We should probably just put this in before it gets stalled for no good reason. Since it is an API addition, might want to ping the mailing list, but honestly, it is small enough to do that after the fact.

Unfortunately, I can't help but bikeshed the name: Do we have a preference for broadcast_shapes vs. broadcast_shape? For some reason I feel the first may be better, so just in case everyone leans that way. If not, either name is completely fine to be honest.

[madhulikajc]

I prefer broadcast_shapes too actually. Curious how others feel. I guess it's a question of whether it's an implicit "get" as in "get_broadcast_shape" or "broadcast these input shapes". Using the plural feels consistent with other numpy functions like expand_dims etc.

[mhvk]

👍 to broadcast_shapes. After all, it is broadcast_arrays too.

[seberg]

I suppose the question of plural vs. singular is whether we are referring to the return value or the input value(s). Btw. noticed one possible addition to the test: Making sure that no input also returns () as a shape.

[mattip]

LGTM. One tiny nit, but can be ignored.

As for moving to C: definitely should go to another PR. The chain goes

• add to numpy/core/multiarray/all
• add a C function to array_module_methods. There are lots of examples there on different arg parsing and signatures. If you are new to this, you might want to look at the Python docs together with the code examples.
• rummage around the code base looking at things like PyArray_Broadcase to see if it can be refactored

Maybe the correct way to reason about this is in the opposite order: there is no reason to move things into C just for fun if no real advantage can be realized.

[eric-wieser]

Maybe the correct way to reason about this is in the opposite order: there is no reason to move things into C just for fun if no real advantage can be realized.

My thinking is we probably already have the code in C somewhere, and if we don't we might be able to refactor it out and use it elsewhere. It's certainly not worth a standalone C API if all it's used for is this function.

[asmeurer]

Thanks for adding this. This looks good to me. The dtype=[] trick looks good (taking Eric's word that it uses no memory), as it will keep things like error messages the same. Though perhaps if the logic in C is ever factored out, it would make more sense for this to call the C function directly.

[mattip]

One last thing. This should be indexed in the reference documents somewhere so it shows up when searching for broadcast_shapes in the rendered documentation (available by looking at the last CI run -> ci/circleci: build artifact -> details). The documentation team might have an idea where the best place for this would be, maybe under utility functions, the source of which is routines.other.rst?

[eric-wieser]

taking Eric's word that it uses no memory

In case you want more than my word:

>>> np.zeros(1000000, dtype=[]).nbytes
0

[rgommers]

The documentation team might have an idea where the best place for this would be, maybe under utility functions, the source of which is routines.other.rst?

Good point. That seems like a decent place. To make it more discoverable, adding an entry to the See Also sections in the docstrings of broadcast_to, broadcast_arrays and perhaps one or two more of the most related functions would be good. And vice versa, the docstring of this new function could also use a See Also.

[rgommers]

Docs also look good now. This is very nice to have, thanks @madhulikajc. Merging.