Investigate enabling array creation routines for quantities

@maxnoe

What is the problem this feature will solve?

Currently, Numpy array creation functions cannot be used to initialize Quantity. However, @maxnoe speculates:

I think it could be implemented now using the like option. Passing an empty quantity to like now creates an error inside astropy:

In [12]: np.arange(0 * u.m, 0.1 * u.m, 0.01 * u.m, like=u.Quantity([], u.m))
WARNING: function 'arange' is not known to astropy's Quantity. Will run it anyway, hoping it will treat ndarray subclasses correctly. Please raise an issue at https://github.com/astropy/astropy/issues. [astropy.units.quantity]
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[12], line 1
----> 1 np.arange(0 * u.m, 0.1 * u.m, 0.01 * u.m, like=u.Quantity([], u.m))

File ~/.local/conda/envs/cta-dev/lib/python3.11/site-packages/astropy/units/quantity.py:1888, in Quantity.__array_function__(self, function, types, args, kwargs)
   1880 else:
   1881     warnings.warn(
   1882         f"function '{function.__name__}' is not known to astropy's Quantity."
   1883         " Will run it anyway, hoping it will treat ndarray subclasses"
   (...)
   1886         AstropyWarning,
   1887     )
-> 1888     return super().__array_function__(function, types, args, kwargs)
   1890 # If unit is None, a plain array is expected (e.g., boolean), which
   1891 # means we're done.
   1892 # We're also done if the result was NotImplemented, which can happen
   1893 # if other inputs/outputs override __array_function__;
   1894 # hopefully, they can then deal with us.
   1895 if unit is None or result is NotImplemented:

AttributeError: 'builtin_function_or_method' object has no attribute '_implementation'

Describe the desired outcome

Enable the Numpy array creation routines for use with quantities: full, zeros, ones, empty, and arange.

Additional context

CC @neutrinoceros.

@mhvk

Thanks for opening this issue. I think this is a reasonable expectation and the way forward is simple but requires a bit of work. Specifically, we need to add implementations like

@dispatched_function
def arange(*args, **kwargs):
    ...

to astropy/units/quantity_helper.py, together with a bunch of tests. I'm game to tackle this, just want to ask @mhvk for his blessing before I get started :)

Thanks! Yes, this has been a hole in the implementations, unfilled mostly because none of the routines are particularly useful for Quantity -- zeros_like, etc., are better and arange doesn't make much sense - and because it is a bit of a pain to use with the like= option. That said, we might as well fill the hole, both for Quantity and Masked!

arange doesn't make much sense

I guess would be nice to be able to do something like: np.arange(0 * u.day, 1 * u.day, 1 * u.s) or np.arange(0 * u.m, 10 * u.m, 10 * u.cm)

arange doesn't make much sense

I guess would be nice to be able to do something like: np.arange(0 * u.day, 1 * u.day, 1 * u.s) or np.arange(0 * u.m, 10 * u.m, 10 * u.cm)

That is just the kind of usage that I have worked around for now: https://m4opt.readthedocs.io/en/latest/api/m4opt.utils.numpy.arange_with_units.html

np.arange(0 * u.day, 1 * u.day, 1 * u.s)

But what is the expected output here?

Also how do you pass in equivalencies?

Note that it would have to be np.arange(0*u.day, 1*u.day, 1*u.s, like=<some-quantity>), since the inputs cannot be numpy arrays and hence numpy does not override on those being Quantity.

I guess the output is in principle fairly obvious (here, effectively the same as np.arange(0., 1., 1./86400) << u.day), but np.arange is really not ideal for floats, because rounding errors can make it unpredictable whether or not the final point is included. So, generally something like np.linspace is better (and does recognize quantity input - indeed, can even do it on array-shaped start and stop!).

Anyway, as said, no reason not to support the __array_function__ interface fully.

@lpsinger

One question: @lpsinger - in your implementation, you use the unit of the start point, which is certainly what I thought one should do initially too. But it may be good to think through what one should really expect, which perhaps depends on what one is trying to solve with np.arange with units. For that purpose, let's take a different example: np.arange(0.*u.yr, 1*u.yr, 1*u.day). Here, there is not an integer number of points. Might it be more logical to take the unit of step? Also, what should np.arange(10*u.day) do? Assume a start of 0 and step=1*stop.unit? Though maybe best to just raise a UnitsError since the default step=1 implies dimensionless.

I guess this brings up more generally what one wants: samples spread over a year separated by exactly one day but not necessarily covering exactly one year, or samples covering exactly one year with a spacing that may not be exactly one day. For arange, I guess it is the first (which is why perhaps the unit of step should be used), but perhaps it would be nice to support the second usage. Logically, that is more like np.linspace, but with a step option. Though then one has to think what step means... Checking numpy, I see there was once a PR for that (numpy/numpy#2945), though it ended because it was not clear how strict one should be in (stop-start)/step being close to an integer.

Anyway, perhaps overthinking things here...

Anyway, perhaps overthinking things here...

Maybe it just shows it's complex and the real solution would be just to raise a nicer error message in case someone tries arange(..., like=u.Quantity(...)

(in passing: I'm listening, and I'll wait until a consensus is reached before I start implementing the agreed solution)

@lpsinger

One question: @lpsinger - in your implementation, you use the unit of the start point, which is certainly what I thought one should do initially too. But it may be good to think through what one should really expect, which perhaps depends on what one is trying to solve with np.arange with units.

I would just take the units of the first argument.

Also, what should np.arange(10*u.day) do? Assume a start of 0 and step=1*stop.unit?

Yes.

OK, I'm fine to go with that.

@mhvk

Sorry for the delay, I'm going to start working on this now.
Early notes:

full, zeros, ones, empty all have x_like pendants that we already support, so there should be some overlap in the implementation (avoid duplicating effort)
arange is the only truly unique case, and since it's also the most discussed one here, I think I'll start with it.
As previously mentionned by @mhvk, we have the same issue between astropy.units and astropy.utils.masked, but for now I'll focus on units.
maybe a preliminary problem to figure out is adjusting our existing tests checking for implementation completion: these functions are not currently explicitly registered anywhere in astropy.units.quantity_helper.function_helpers (not even in UNSUPPORTED_FUNCTIONS), but our tests still pass

I just opened an issue upstream about that last point numpy/numpy#27451

Another note that a subset of NEP 35 supporting functions can be discovered as

diff --git a/astropy/units/tests/test_quantity_non_ufuncs.py b/astropy/units/tests/test_quantity_non_ufuncs.py
index c0e9c7c31e..612ce7f7ee 100644
--- a/astropy/units/tests/test_quantity_non_ufuncs.py
+++ b/astropy/units/tests/test_quantity_non_ufuncs.py
@@ -59,11 +59,25 @@ def get_wrapped_functions(*modules: ModuleType) -> set[FunctionType]:
     else:
         from numpy.testing.overrides import allows_array_function_override
 
+    def has_like_keyword(f) -> bool:
+        # temporary solution, waiting for a stable upstream API
+        # https://github.com/numpy/numpy/issues/27451
+        try:
+            sig = inspect.signature(f)
+        except ValueError:
+            # uninspectable builtin function
+            return False
+
+        if "like" not in sig.parameters:
+            return False
+
+        return sig.parameters["like"].kind is KEYWORD_ONLY
+
     return {
         f
         for module in modules
         for f in module.__dict__.values()
-        if callable(f) and allows_array_function_override(f)
+        if callable(f) and (allows_array_function_override(f) or has_like_keyword(f))
     }

however this only works for Python function, not C functions, so it's not a satisfactory replacement for the missing upstream API. Anyway, I'll leave this part of the problem aside, and focus on implementing the functions we already know we care about for now !

progress report on np.arange:
this function is a bit complicated to implement correctly because of how flexible its signature is, but Python' structural pattern matching (i.e. match/case) comes to the rescue and makes this actually doable. One thing that surprised me in my initial experimentation is that numpy doesn't seem to pass down the exact value of the like argument: we can't access it within our dispatched implementation of np.arange. I don't think it's blocking in any way, but it's something to keep in mind.

I think we do get the class, though, right? So, we can create an Angle or a Quantity?

If I'm not mistaken, we don't get to see anything about the like argument from where we stand. If Angle merely gets its __array_function__ implementation from Quantity by inheritance, then we cannot tell from within our dispatch whether the expected return type is Quantity or Angle.

my implementation is at: #17059
The PR isn't finished but early feedback is welcome.

That cannot be quite true: in then, the __array_function__ implementation is on the class, but it is the instance's version that is called, and its type is passed on as well:

In [2]: class A:
   ...:     def __array_function__(self, *args, **kwargs):
   ...:         print(repr(self), *args, **kwargs)
   ...: 

In [3]: np.arange(0., 10., 2., like=A())
<__main__.A object at 0x7f2a221768d0> <built-in function arange> (<class '__main__.A'>,) (0.0, 10.0, 2.0) {}

In [4]: class B(A):
   ...:     pass
   ...: 

In [5]: np.arange(0., 10., 2., like=B())
<__main__.B object at 0x7f2a21faf050> <built-in function arange> (<class '__main__.B'>,) (0.0, 10.0, 2.0) {}

I.e., basically, we have self and its type. But I don't think so far we have used either in any of the implementations, since it wasn't needed.

p.s. Did not have a change to look at the implementation yet...

I.e., basically, we have self and its type.
oh, right.

But I don't think so far we have used either in any of the implementations, since it wasn't needed.

yes, NEP 35 != NEP 18 after all, but let's see how far we can stretch the existing facilities before we need to build more !

#17059 is now ready for review !

With #17059, I'll make it my top priority to implement the rest of NEP 35 functions for Quantity, and I'm hoping I'll also be able to tackle astropy.utils.masked, though that is a secondary goal at this stage. Right now I'll make an attempt at addressing the upstream issue (numpy/numpy#27451) first so the goal can be more clearly defined for us all (including other libraries).

@lpsinger

Here are the rest of the functions requested by @lpsinger : #17120

Navigation Menu

Repository navigation

Investigate enabling array creation routines for quantities #17001

What is the problem this feature will solve?

Describe the desired outcome

Additional context

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Investigate enabling array creation routines for quantities #17001

Description

What is the problem this feature will solve?

Describe the desired outcome

Additional context

Activity

neutrinoceros commented on Sep 12, 2024

mhvk commented on Sep 12, 2024

maxnoe commented on Sep 12, 2024

lpsinger commented on Sep 12, 2024

pllim commented on Sep 12, 2024

pllim commented on Sep 12, 2024

mhvk commented on Sep 12, 2024

mhvk commented on Sep 12, 2024

maxnoe commented on Sep 12, 2024

neutrinoceros commented on Sep 13, 2024

lpsinger commented on Sep 13, 2024

mhvk commented on Sep 13, 2024

neutrinoceros commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

mhvk commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

mhvk commented on Sep 23, 2024

neutrinoceros commented on Sep 23, 2024

neutrinoceros commented on Sep 25, 2024

neutrinoceros commented on Oct 2, 2024

neutrinoceros commented on Oct 3, 2024

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions