Modified astropy.units.Quantity.__array_ufunc__() to return NotImplemented if the inputs are incompatible

[byrdie]

Description

This pull request is to address the current behavior of astropy.units.Quantity.__array_ufunc__(), which raises a ValueError instead of returning NotImplemented if the inputs are incompatible.

Fixes #13976

Checklist for package maintainer(s)

This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see "When to rebase and squash commits".
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the Extra CI label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the no-changelog-entry-needed label. If this is a manual backport, use the skip-changelog-checks label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• Is a milestone set? Milestone must be set but astropy-bot check might be missing; do not let the green checkmark fool you.
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate backport-X.Y.x label(s) before merge.

[pllim]

This is a behavior change, so I don't think we should backport.

[mhvk]

Really nice that this works! I have to think about it a little more, so just some initial comments for now.

[pllim]

Just keep in mind that there might be a bunch of conflicts you need to resolve when we apply black to the main branch, which is going to happen over the next 2 weeks. FYI.

[byrdie]

@mhvk, I've been thinking about this more and I've been wondering if it would really be that bad of a performance hit if we just checked the types at the top of the function? It just seems that the way we have it now, duck types will have mildly poor performance since it has to evaluate most of the try block before deciding to return NotImplemented.

[mhvk]

@byrdie - I think we should optimize for Quantity, since that is mostly what the function will get; it should be OK to be a but slow in returning NotImplemented.

[byrdie]

@mhvk, ok I just thought I'd ask since it would allow us to avoid this problem of having such a large try block, and the speed increase for duck types would be considerable compared to the penalty for pure Quantities.

FWIW I just timed it on my machine to check, and it seems pretty negligible to me.

inputs = (np.empty((3, 4, 5)) * u.mm, np.empty((3, 4, 5)) * u.mm)
outputs = (np.empty((3, 4, 5)) * u.mm, )
%time not all(isinstance(io, (int, float, tuple, list, np.ndarray, u.Quantity, astropy.table.Column)) for io in inputs + outputs)
CPU times: total: 0 ns
Wall time: 0 ns

Is this not a very accurate test?

[mhvk]

Use %timeit instead -- I get 1.69 micro-sec, which is not quite negligible.

I think we can put the try/except around the loop that converts inputs instead (note that these days try/except has zero costs if no exception is raised).

[mhvk]

p.s. Sorry to be so picky, but I try to be super careful to ensure use of Quantity is as fast as possible (e.g., we made quite a few micro-optimizations in the initializer). This since many have complained it is too slow and one should just use regular arrays -- thus loosing the benefit of having correct units and potentially giving a few more incorrect scientific results.

[byrdie]

I don't think you're being picky, I'm totally on-board with the need for speed and I agree that 1.69 microseconds is too slow.

What do you mean by putting the try/except around the loop that converts inputs? Is it the converters_and_unit() function? If so, see the issues I had with that in our conversation above.

[mhvk]

(Without further checking) I meant

astropy/astropy/units/quantity.py

Lines 666 to 670 in a93dd3d

	# Same for inputs, but here also convert if necessary.
	arrays = []
	for input_, converter in zip(inputs, converters):
	input_ = getattr(input_, "value", input_)
	arrays.append(converter(input_) if converter else input_)

[byrdie]

@mhvk, I think that is what I originally tried. If I do that, I get the error:

function = <ufunc 'add'>, method = '__call__'
args = (<Quantity 1. m>, DuckQuantity1(data=<Quantity 1. mm>))
ufunc_helper = <function get_converters_and_unit at 0x00000231DAD5DC60>
units = [Unit("m"), None], converters = [None, False]

    def converters_and_unit(function, method, *args):
        """Determine the required converters and the unit of the ufunc result.
    
        Converters are functions required to convert to a ufunc's expected unit,
        e.g., radian for np.sin; or to ensure units of two inputs are consistent,
        e.g., for np.add.  In these examples, the unit of the result would be
        dimensionless_unscaled for np.sin, and the same consistent unit for np.add.
    
        Parameters
        ----------
        function : `~numpy.ufunc`
            Numpy universal function
        method : str
            Method with which the function is evaluated, e.g.,
            '__call__', 'reduce', etc.
        *args :  `~astropy.units.Quantity` or ndarray subclass
            Input arguments to the function
    
        Raises
        ------
        TypeError : when the specified function cannot be used with Quantities
            (e.g., np.logical_or), or when the routine does not know how to handle
            the specified function (in which case an issue should be raised on
            https://github.com/astropy/astropy).
        UnitTypeError : when the conversion to the required (or consistent) units
            is not possible.
        """
    
        # Check whether we support this ufunc, by getting the helper function
        # (defined in helpers) which returns a list of function(s) that convert the
        # input(s) to the unit required for the ufunc, as well as the unit the
        # result will have (a tuple of units if there are multiple outputs).
        ufunc_helper = UFUNC_HELPERS[function]
    
        if method == "__call__" or (method == "outer" and function.nin == 2):
            # Find out the units of the arguments passed to the ufunc; usually,
            # at least one is a quantity, but for two-argument ufuncs, the second
            # could also be a Numpy array, etc.  These are given unit=None.
            units = [getattr(arg, "unit", None) for arg in args]
    
            # Determine possible conversion functions, and the result unit.
            converters, result_unit = ufunc_helper(function, *units)
    
            if any(converter is False for converter in converters):
                # for multi-argument ufuncs with a quantity and a non-quantity,
                # the quantity normally needs to be dimensionless, *except*
                # if the non-quantity can have arbitrary unit, i.e., when it
                # is all zero, infinity or NaN.  In that case, the non-quantity
                # can just have the unit of the quantity
                # (this allows, e.g., `q > 0.` independent of unit)
                try:
                    # Don't fold this loop in the test above: this rare case
                    # should not make the common case slower.
                    for i, converter in enumerate(converters):
                        if converter is not False:
                            continue
                        if can_have_arbitrary_unit(args[i]):
                            converters[i] = None
                        else:
                            raise UnitConversionError(
                                f"Can only apply '{function.__name__}' function to "
                                "dimensionless quantities when other argument is not "
                                "a quantity (unless the latter is all zero/infinity/nan)."
                            )
                except TypeError:
                    # _can_have_arbitrary_unit failed: arg could not be compared
                    # with zero or checked to be finite. Then, ufunc will fail too.
>                   raise TypeError(
                        "Unsupported operand type(s) for ufunc {}: '{}'".format(
                            function.__name__,
                            ",".join([arg.__class__.__name__ for arg in args]),
                        )
                    )
E                   TypeError: Unsupported operand type(s) for ufunc add: 'Quantity,DuckQuantity1'

..\quantity_helper\converters.py:208: TypeError

inside converters_and_unit. Do you think I should have two try/except blocks or should we just leave it as one large block?

[mhvk]

OK, I'll try to think a bit about this -- you are right that a super-big try/except is really not that nice (even if it does the job!).

[byrdie]

At this point, I'm tempted to just call it good the way it is, since it works and I've been testing it with my code for a while. The only other thing I can think of is having two try/execpt blocks with a separate function for the body of the except block to avoid repetition.

[mhvk]

Two small little clean-ups found upon another read. With that, ready to go in!

[mhvk]

All OK now, thanks!

[mhvk]

Hmm, the docs build failed, which I really do not understand since the PR that goes before this one (#14140) passed. Yet it seems related to that PR, not to what you have here. I'll need to investigate...

[byrdie]

My other PR is failing the docs build as well which has been confusing me.

[mhvk]

@byrdie - those failures are definitely a leftover of #14140 -- no idea how that passed with those omissions! Anyway, hopefully that's fixed in #14171.

For here, thinking about it more, really we should return NotImplemented only if it is not one of the list you have, or a subclass of ndarray that overrides __array_ufunc__. More generally, for anything to happen with NotImplemented, some of the other arguments have to actually define __array_ufunc__, so I guess one could do something like

if not all(getattr(obj.__class__, '__array_ufunc__', None) in (None, ndarray.__array_ufunc__, self.__class__.__array_ufunc__)
    for obj in ...):

What do you think?

[byrdie]

@mhvk, that looks good to me, I'll push an update shortly.

[mhvk]

@byrdie - sorry, but on yet another look, I found yet another really minor issue. With that out of the way, this really can go in...

[mhvk]

All OK now, let's get it in. Thanks, @byrdie!

[byrdie]

Thanks for all your help @mhvk!