ENH: (NEP 18) implement and test cumulative functions

[neutrinoceros]

cumprod is a bit tricky, I don't think there's a right way to handle it, so I'm just adding a warning to raise awareness

[ngoldbaum]

I think it would be more useful to raise an error. It looks like astropy does this generically for the accumulate method on the multiply ufunc (which cumprod calls internally):

In [1]: import numpy as np

In [2]: from astropy.units import gram

In [3]: arr = [1, 2, 3]*gram

In [4]: arr
Out[4]: <Quantity [1., 2., 3.] g>

In [5]: np.cumprod(arr)
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
Cell In [5], line 1
----> 1 np.cumprod(arr)

File <__array_function__ internals>:185, in cumprod(*args, **kwargs)

File ~/.pyenv/versions/3.10.8-debug/lib/python3.10/site-packages/astropy/units/quantity.py:1715, in Quantity.__array_function__(self, function, types, args, kwargs)
   1702 # A function should be in one of the following sets or dicts:
   1703 # 1. SUBCLASS_SAFE_FUNCTIONS (set), if the numpy implementation
   1704 #    supports Quantity; we pass on to ndarray.__array_function__.
   (...)
   1712 # function is in none of the above, we simply call the numpy
   1713 # implementation.
   1714 if function in SUBCLASS_SAFE_FUNCTIONS:
-> 1715     return super().__array_function__(function, types, args, kwargs)
   1717 elif function in FUNCTION_HELPERS:
   1718     function_helper = FUNCTION_HELPERS[function]

File ~/.pyenv/versions/3.10.8-debug/lib/python3.10/site-packages/numpy/core/fromnumeric.py:3145, in cumprod(a, axis, dtype, out)
   3084 @array_function_dispatch(_cumprod_dispatcher)
   3085 def cumprod(a, axis=None, dtype=None, out=None):
   3086     """
   3087     Return the cumulative product of elements along a given axis.
   3088 
   (...)
   3143 
   3144     """
-> 3145     return _wrapfunc(a, 'cumprod', axis=axis, dtype=dtype, out=out)

File ~/.pyenv/versions/3.10.8-debug/lib/python3.10/site-packages/numpy/core/fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds)
     54     return _wrapit(obj, method, *args, **kwds)
     56 try:
---> 57     return bound(*args, **kwds)
     58 except TypeError:
     59     # A TypeError occurs if the object does have such a method in its
     60     # class, but its signature is not identical to that of NumPy's. This
   (...)
     64     # Call _wrapit from within the except clause to ensure a potential
     65     # exception has a traceback chain.
     66     return _wrapit(obj, method, *args, **kwds)

File ~/.pyenv/versions/3.10.8-debug/lib/python3.10/site-packages/astropy/units/quantity.py:591, in Quantity.__array_ufunc__(self, function, method, *inputs, **kwargs)
    569 """Wrap numpy ufuncs, taking care of units.
    570 
    571 Parameters
   (...)
    585     Results of the ufunc, with the unit set properly.
    586 """
    587 # Determine required conversion functions -- to bring the unit of the
    588 # input to that expected (e.g., radian for np.sin), or to get
    589 # consistent units between two inputs (e.g., in np.add) --
    590 # and the unit of the result (or tuple of units for nout > 1).
--> 591 converters, unit = converters_and_unit(function, method, *inputs)
    593 out = kwargs.get('out', None)
    594 # Avoid loop back by turning any Quantity output into array views.

File ~/.pyenv/versions/3.10.8-debug/lib/python3.10/site-packages/astropy/units/quantity_helper/converters.py:276, in converters_and_unit(function, method, *args)
    265         raise TypeError("Cannot use '{1}' method on ufunc {0} with a "
    266                         "Quantity instance as the result is not a "
    267                         "Quantity.".format(function.__name__, method))
    269     if (converters[0] is not None or
    270         (unit is not None and unit is not result_unit and
    271          (not result_unit.is_equivalent(unit) or
   (...)
    274         # then things like np.cumprod will not longer fail (they check
    275         # for TypeError).
--> 276         raise UnitsError("Cannot use '{1}' method on ufunc {0} with a "
    277                          "Quantity instance as it would change the unit."
    278                          .format(function.__name__, method))
    280 return converters, result_unit

UnitsError: Cannot use 'accumulate' method on ufunc multiply with a Quantity instance as it would change the unit.

[ngoldbaum]

although looking at the traceback it seems they do have a cumprod wrapper, it just always (?) raises an error.

[neutrinoceros]

It seems that it does error 100% of the time yes, including for scalar arrays (to which cumprod, while meaningless, is also arguably harmless).
I guess it's best to align with astropy here and raise an error too, thank you for pointing this out.

[neutrinoceros]

I'm not sure which existing exception I should raise however, which makes me wonder wether it'd be a good idea to have a generic UnytError exception type (that all other unyt exceptions would possibly inherit from ?).

I'll do that here while I wait for your feedback, but I'll leave the class hierarchy alone for now.

[neutrinoceros]

@ngoldbaum any idea when you might be able to re-iterate on this one ? I usually wait for the current NEP 18 PR to go through before I start working on another one but I could adapt my strategy if appropriate. So I absolutely don't mean to pressure you, just let me know if you can.

[ngoldbaum]

Hi sorry for taking so long to look at this. I took last week off from coding and have been focusing on other stuff at work.

I'm not really a fan of UnytError but I also don't think that should block getting the __array_function__ stuff working. I wish this was raising UnitOperationError, or we somehow made UnitOperationError more general since I would guess that's the most common type of exception people would be catching for this sort of thing and it irks me a bit that they'd need to catch more than one kind of exception for different corner cases.

We probably need to more carefully look at how exceptions work in unyt in general since right now the situation is kind of a hodgepodge, although that might need a deprecation cycle since we'd be doing an API break.

For now I'm just going to merge this, but I'd like to have a discussion about how to handle exceptions, whether we need to do some sort of deprecation cycle, and how we can make it simpler to deal with exceptions raised by unyt before we do the final release.

[neutrinoceros]

Thanks. I agree it'd be good to discuss exceptions in more details. I opened an issue from your comment so it doesn't get lost. I'll give it some thoughts tomorrow on a rested mind.