Alternative to #929 to solve issues with array Quantity objects - #1139

Amazingly, this works without modifications:

In [13]: np.mean(q)
Out[13]: <Quantity 6.0 m>

In [14]: np.median(q)
Out[14]: <Quantity 6.0 m>

What does just np.array(q) do with this? My suspicion is it should just work. What about things like __radd__ and __rsub__?

@iguananaut - np.array(q) shouldn't convert to a unitless array in my opinion because the units are uncertain. However, note that it doesn't crash, and even returns sensible results for complicated methods like cumprod because it makes it into an array of quantities:

In [15]: np.array(np.array([1,2,3]) * u.m / u.cm)
Out[15]: array([<Quantity 1 m / cm>, <Quantity 2 m / cm>, <Quantity 3 m / cm>], dtype=object)

In [16]: np.array(np.array([1,2,3]) * u.m / u.cm).cumprod()
Out[16]: array([<Quantity 1 m / cm>, <Quantity 2 m2 / cm2>, <Quantity 6 m3 / cm3>], dtype=object)

Anyway, I think this should be discussed in the mailing list at this stage.

@iguananaut - subtraction works as expected:

In [17]: np.array([1,2,3]) * u.m + np.array([4,5,6]) * u.cm
Out[17]: <Quantity [ 1.04  2.05  3.06] m>

In [18]: np.array([1,2,3]) * u.m + np.array([4,5,6])
ERROR: TypeError: Object of type '<type 'numpy.ndarray'>' cannot be added with a Quantity object. Addition is only supported between Quantity objects with compatible units. [astropy.units.quantity]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-18-5d161d68162d> in <module>()
----> 1 np.array([1,2,3]) * u.m + np.array([4,5,6])

/Users/tom/Library/Python/2.7/lib/python/site-packages/astropy-0.3.dev4423-py2.7-macosx-10.6-x86_64.egg/astropy/units/quantity.pyc in __add__(self, other)
    268                 "Object of type '{0}' cannot be added with a Quantity "
    269                 "object. Addition is only supported between Quantity "
--> 270                 "objects with compatible units.".format(other.__class__))
    271 
    272     def __sub__(self, other):

TypeError: Object of type '<type 'numpy.ndarray'>' cannot be added with a Quantity object. Addition is only supported between Quantity objects with compatible units.

That's not radd and rsub. I mean like:

In [6]: (np.array([4, 5, 6]) * u.m) + np.array([1, 2, 3])
ERROR: TypeError: Object of type '<type 'numpy.ndarray'>' cannot be added with a Quantity object. Addition is only supported between Quantity objects with compatible units. [astropy.units.quantity]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-f1bae0925764> in <module>()
----> 1 (np.array([4, 5, 6]) * u.m) + np.array([1, 2, 3])

/internal/1/root/src/astropy/astropy/astropy/units/quantity.py in __add__(self, other)
    268                 "Object of type '{0}' cannot be added with a Quantity "
    269                 "object. Addition is only supported between Quantity "
--> 270                 "objects with compatible units.".format(other.__class__))
    271 
    272     def __sub__(self, other):

TypeError: Object of type '<type 'numpy.ndarray'>' cannot be added with a Quantity object. Addition is only supported between Quantity objects with compatible units.

In [7]: np.array([1, 2, 3]) + (np.array([4, 5, 6]) * u.m)
ERROR: TypeError: unsupported operand type(s) for +: 'int' and 'Quantity' [IPython.core.interactiveshell]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-7-288d9af21ded> in <module>()
----> 1 np.array([1, 2, 3]) + (np.array([4, 5, 6]) * u.m)

TypeError: unsupported operand type(s) for +: 'int' and 'Quantity'

At least in the latter case Numpy rejects it, albeit with a different exception. But as I wrote in the mailing list I think this can be handled in __array_prepare__ if need be.

Also I don't really like the np.array(np.array([1, 2, 3]) * u.m) returns an object array. This means you can replace values in the array with quantities of different units, for example, or even any other arbitrary object.

So I don't really see the advantage of removing __array__ so long as we add __array_prepare__ and __array_wrap__ methods that do the appropriate things.

@iguananaut - here's what bothers me about keeping __array__ - if we go with your suggestion of keeping it, I don't know how we can get rid of the following warning:

In [3]: np.array([1,2,3]) * (0.01 * u.Unit(0.1 * u.m))
WARNING: Converting Quantity object in units '0.1 m' to a Numpy array [astropy.units.quantity]
Out[3]: <Quantity [ 0.01  0.02  0.03] 0.1 m>

Can you see how this could be done?

I also am now against having implicit conversions that drop units, because they simply don't have predictable results to the users - this is another reason why I would prefer to drop array, though I'd also like to make sure that np.array(q) doesn't work either since as you pointed out it could be confusing.

I'll have to think more about what to do with the warning, but for the other stuff I refer you back to every single point I made last time we talked about this.

So far all I see is a desire to remove incredibly useful functionality for the sake a hypothetical problem that I've yet to see.

@iguananaut - so just so I understand your example with:

np.array([1, 2, 3]) + (np.array([4, 5, 6]) * u.m)

are you saying the exception is wrong and that we should change the exception?

If you can find a way to remove that warning, then I agree that there's no harm in keeping the implicit conversion that drops units provided that it does emit a warning in that case. I do agree it's useful functionality as long as the warning is there. What bothers me at the moment is just the presence of the warning in cases where units get added back, and if the only solution is to completely disable the warning, then I would want to disable implicit conversion to unitless objects. Do you see what I mean?

Ultimately, the real problem is that Numpy is buggy and that it shouldn't be calling __array__ in the first place when doing:

np.array([1,2,3]) * (0.01 * u.Unit(0.1 * u.m))

After all, it doesn't do it for:

(0.01 * u.Unit(0.1 * u.m)) * np.array([1,2,3])

So maybe the solution is to find a workaround to the Numpy bug rather than trying to selectively disable the warning. What do you think?

Okay, I've spent a little time debugging this to understand exactly what's going on in Numpy's head here.

Upon entering a ufunc one of the first things Numpy does is to attempt to create an ndarray from each argument to the function using every way it knows how, including calling that argument's __array__ if it has one.

When we remove __array__ what it ends up doing with a Quantity is treating it as an iterable, and building an object array from each scalar Quantity.

A little later down it does the check to see if the right-hand argument has an __r*__ method such as __rmul__, but for some reason it only performs this check if the dtype of array created from the right-hand side is NPY_OBJECT. That's what leads to the seemingly buggy behavior, and yet it seems too deliberate to be called a "bug" (though no justification is given in the code for this behavior).

So when we remove __array__ a Quantity converts to an array in the form like array([<Quantity 1 m / cm>, <Quantity 2 m2 / cm2>, <Quantity 6 m3 / cm3>], dtype=object) which is confusing. But because it's an object array Numpy remembers to check the original Quantity argument for an __rmul__ method, and calls that instead of proceeding with the ufunc as normal.

So really you're damned either way--if we drop __array__ then we can't issue any sort of warning about dropping the units. If we keep it we get the desired behavior with respect to operations, but get superfluous warnings on certain operations.

I still feel like there must be a workaround, but I don't have anything specific to give you yet. I still wonder if making Quantity an ndarray subclass isn't still the best option. Now that I have a better understanding how __array_prepare__ and __array_wrap__ work we can work around most of the issues with that solution.

@iguananaut - thanks for digging deeper into this - this is very helpful!

This still needs a lot of tests to cover the new functionality, but this implements all the functionality in http://mpia.de/~robitaille/astropy/quantity.html

@iguananaut - by the way, the behavior for rsub and radd should make more sense now. In all cases, the error is:

UnitsException: '' (dimensionless) and 'm' (length) are not convertible

@astrofrog - this is a wonderful start. Here, a specific 1+v/c test, setting up with:

import astropy.units as u
import numpy as np
from astropy.constants import c
v = np.arange(3)* 1e5 * u.km/ u.s

First test:

In [8]: 1+v/c, np.mean(1+v/c), np.mean(1+v/c).unit, np.mean(1+v/c).value
Out[8]: 
(<Quantity [ 0.001       0.00133356  0.00166713] km / m>,
 <Quantity 0.0013335640952 km / m km / m>,
 Unit("km / m"),
 array(<Quantity 0.0013335640952 km / m>, dtype=object))

So, something odd happens because the array of Quantities that is formed. This would seem one argument to disallow that from happening. Or perhaps this just requires one to define a "mean" attribute. I'll further comment on astropy-dev.

(I deleted a previous comment I made because I was mistaken)

@mhvk - I'm not sure I understand what is wrong in the example you sent. If I use this branch, I get:

In [9]: 1 + v/c
Out[9]: <Quantity [ 0.001       0.00133356  0.00166713] km / m>

In [10]: np.mean(1 + v/c)
Out[10]: <Quantity 0.0013335640952 km / m>

In [11]: np.mean(1 + v/c).unit
Out[11]: Unit("km / m")

In [12]: np.mean(1 + v/c).value
Out[12]: 0.0013335640951981521

Similar to what you get. This is correct, right? I mean ideally it would figure out it has to decompose, which you can do manually:

In [13]: (1 + v/c).decompose()
Out[13]: <Quantity [ 1.          1.3335641   1.66712819] >

but apart from that there's nothing fundamentally wrong with the results, right?

@astrofrog - very weird, I just redid what I had before, and I still get the same I had before, i.e., np.mean(1+v/c).value gives an array that contains a Quantity. Could it be a numpy issue? I'm only at 1.6 (which is what Debian has in testing and unstable; it has 1.7 in experimental). Can you easily test this? (I can too, but not easily). I guess one option would be to add a test and run it through travis.

@mhvk - I added a test, let's see what Travis says.

@mhvk - bingo, it's Numpy 1.6 that is causing issues. Will look into it locally.

@mhvk - this should now be fixed - could you try it out with your installation?

@astrofrog - yes, fixed indeed. Great!

@astrofrog - since it fairly directly affects the pull request, I thought I should repeat a comment also made on astropy-dev: Shouldn't all functions return quantities instead of arrays? This will need to be the case anyway for arcsin, etc., but also seems more logical in that one should not have unexpected type changes (and will help eventual inclusion of uncertainties 😈)

@mhvk - that is a good point, and would be easy enough to do. This reminds me I should implement the inverse trigonometric functions!

Numpy is driving me nuts. It turns out that things are made a lot harder by the fact that Numpy can convert Quantity arrays to arrays of Quantities:

In [5]: np.array(np.arange(3) * u.m)
Out[5]: array([<Quantity 0 m>, <Quantity 1 m>, <Quantity 2 m>], dtype=object)

When one does np.sqrt(q), it first looks for the sin method and runs that (which returns a Quantity) then it decides for some reason to convert it to an array, giving an array of Quantity objects as above, then it passes it to __array_wrap__. This led me to do the kind of hack you can see in __array_wrap__ for sqrt. Now I have to do that for all the trigonometric functions too if we want conversions from/to dimensionless Quantities, and it seems a bit silly. But it turns out that if one removes __len__ then Numpy suddenly doesn't know how to convert to an array, and this then actually allows everything to work properly. This seems like a whack-a-mole - if we want to be able to cleanly implement ufuncs, we have to get rid of __len__... A bit silly really...

@mhvk @iguananaut - the ufuncs now return dimensionless quantities (where appropriate) and the inverse trigonometric functions now work, so things like np.arccos(np.cos(q)) will round-trip.

@iguananaut - as you'll see, __array_wrap__ actually re-does the calculations as I mentioned above, which is not idea, but at least works for now. Even if we move to a better system, I'm hoping that at least the tests represent the desired behavior and can be re-used. Let me know if you can think of a better way to do it though!

@astrofrog - it is really nice to see this all working! What I personally still would like is to drop the test for whether scale is unity (this was something you raised on astropy-dev, but didn't get many replies). Since just converting a quantity with a dimension to numbers is no longer possible without .value, I think there no longer remains much ambiguity about what the user wants. Hence, I think for all those functions, if a quantity can be made dimensionless, it should be done, i.e., for

import astropy.units as u
import numpy as np
from astropy.constants import c
v = np.arange(3)*1.e5*u.km/u.s
np.log10(1+v.si/c)
np.log10(1+v/c)

both np.log10 should work. I think this would also be expected behaviour given that

a = np.arange(3)*u.deg
np.sin(a)
np.sin(a.si)

do both work (as does the somewhat nonsensical np.cos((1+v/c)*360.*u.deg) I mentioned on astropy-dev).
(Hey, with this implemented, my oldest still open issue - #862 - could be closed ;-)

p.s. of course this pull request would not just close #862, but also @mdboom's #955 and your own #669... (mentioning these in part to help remind what can be cleaned up).

If I am understanding correctly, @mhvk's last comment is going back to @astrofrog's question on astropy-dev, right? That is:

should implicit float conversion for dimensionless but *not*
unscaled quantities (such as m/km) raise an error (as it currently
does in my branch) or should it automatically decompose the values? If
we go with the latter, it could always raise an INFO message?

I'm fairly neutral on this, but I think @mhvk's example demonstrates that we would want this. I can't really think of any cases where this would be ambiguous, and it allows a nice way of doing the log of some unitful quantity:

v = somenum * u.km / u.s
r = somenum * u.kpc
m = v**2 * r / cnst.G # G is the gravitational G
logm = np.log10(m/u.solMass)

The alternative is to do np.log10(m.to(u.solMass).value), and I think the above form is more "right" the way I think about units.

Aside from that, this looks good to me. @astrofrog, is there pending work here, or is this all currently awaiting comment/endorsement?

@eteq - while the current PR works ok at first glance, it's horribly inefficient, because if you compute e.g. np.sin(q) it will first return a Quantity with an array inside, which Numpy then converts to an array of Quantities (very bad if there are e.g. millions of values!) and the __array_wrap__ recomputes the Quantity from scratch. This is all due to the way Numpy works. So I'll be organizing a small hangout to talk about these issues.

We decided to go with #929 instead of this, so closing!