Monkey-patch non-ufunc Numpy functions

@iguananaut

Once #929 is merged, we should try and monkey-patch a few of the common Numpy functions that are not ufuncs, such as np.dot.

@iguananaut - maybe you'll be able to take a stab at this since you had some ideas?

Just for the records, there's already a file to put tests for these - astropy/units/tests/test_quantity_non_ufuncs.py

For example,

In [5]: w1 = u.Quantity(1000, u.AA)

In [6]: w2 = u.Quantity(9000, u.AA)

In [7]: np.hstack([w1, w2])
Out[7]: <Quantity [1000 9000] (Unit not initialised)>

does not work properly.

Just so searches for this issue also lead to the current work-around (possibly even the preferred way to do it, though I do hope we can get np.hstack to just work too) [1]:

u.Quantity([w1, w2])
<Quantity [1000,9000] Angstrom>

EDIT: in terms of other work-arounds: while np.dot(q1,q2) does not work q1.dot(q2) is fine; see #1930.

EDIT: another case to solve is np.linalg.norm (see #1930)
And another is scipy.integrate.odeint (maybe rather hard to get right)

[1] https://groups.google.com/forum/#!topic/astropy-dev/zH98DFp0YKM

@astrofrog

Just so we do not forget: Tests of numpy 1.8.0rc1 showed a similar problem with np.median [1]; following @astrofrog's posting to the numpy issue list [2], pointers were given to the planned appearance of a __numpy_ufunc__ method [3], and the suggestion was made that the best solution for quantities was to use a separate dtype rather than to subclass ndarray -- likely similar to the hotel-california dtype discussed in #929.

[1] https://groups.google.com/forum/#!topic/astropy-dev/-heqFIhQHP0
[2] numpy/numpy#3846
[3] e.g., http://docs.scipy.org/doc/numpy-dev/reference/arrays.classes.html

Hopefully, we'll solve this at some point; in the meantime, I added this to the "known issues" (see #1937). One PR I hope can be undone soon!

I've made some initial trials with __numpy_ufunc__ for the development version of numpy 1.9 [1]. Most tests for units work (but not yet the tests involving Quantity subclasses), and it solves the problem for np.dot (but not for np.hstack, np.inner, np.outer). Unfortunately, there is a numpy bug preventing the inplace use of ufuncs with two arguments [2].

[1] https://github.com/mhvk/astropy/compare/quantity-numpy-ufunc?expand=1
[2] numpy/numpy#4160

I take it you'd vote for numpy/numpy#4164 :)

@embray

@embray @mhvk - at this point I'm thinking that maybe we should at least provide equivalents to Numpy functions in astropy.units, e.g.

d = u.linspace(3 * u.kpc, 10 * u.kpc, 10)

and maybe not monkey-patch. Monkey patching is always tricky and we don't really want to mess with people's numpy, right?

Maybe the functions should be in a separate namespace though, e.g.

from astropy.units import numpy_wrappers as unp
unp.linspace(...)

or similar?

In any case, if we did decide to monkey patch, we could then use these functions as a starting point.

@astrofrog

@astrofrog - I like the numpy_wrappers idea. Though perhaps we should try to ensure that with each wrapper we also submit a bug report or make a PR to numpy to fix the particular issue (I'm still wondering whether I could tackle numpy/numpy#4164, but as so often time is in short supply).

@mhvk

@mhvk - the problem is that this is not really a bug in numpy, it is just that by design non-ufuncs don't have to call any methods on the object like ufuncs do. Of course numpy/numpy#4164 would fix it but until then, I don't think it can be treated as a bug?

I think this can just be closed now that the consensus has been against monkey-patching.

@mhvk

@mhvk In case it's relevant to your efforts to improve Quantity / numpy compatibility, I wanted to mention that with the update to Numpy 1.9 in gammapy/gammapy#200 for me code that was using Quantity broke.

In Numpy 1.8 np.where returns an array and we had code like this:

>>> import numpy as np
>>> from astropy.units import Quantity
>>> a = np.where([True, False], Quantity([1,2], 'cm'), Quantity([10, 20], 'cm'))
>>> b = Quantity(a, 'cm')
>>> a
array([  1.,  20.])
>>> b
<Quantity [  1., 20.] cm>

With numpy 1.9 np.where retuns an uninitialised Quantity and our code broke like this:

>>> import numpy as np
>>> from astropy.units import Quantity
>>> a = np.where([True, False], Quantity([1,2], 'cm'), Quantity([10, 20], 'cm'))
>>> b = Quantity(a, 'cm')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.0.dev9928-py3.4-macosx-10.9-x86_64.egg/astropy/units/quantity.py", line 178, in __new__
    value = value.to(unit)
  File "/Users/deil/Library/Python/3.4/lib/python/site-packages/astropy-1.0.dev9928-py3.4-macosx-10.9-x86_64.egg/astropy/units/quantity.py", line 579, in to
    self.unit.to(unit, self.value, equivalencies=equivalencies))
AttributeError: 'NoneType' object has no attribute 'to'
>>> a
<Quantity [  1., 20.] (Unit not initialised)>

@cdeil

@cdeil - yes, this error also is biting us in #2958. In that issue, it was suggested to also report this with numpy, since it is indeed annoying. In thinking about it, though, I wasn't quite sure what the preferred behaviour should be. What do you think? One option would be that the output takes its class from the first item in np.where(condition, first, second), but that somehow it is tested whether second is compatible with first. E.g., for Quantity, one would need to check units and possibly convert... Possibly, most of this could be handled if np.where were treated as a ufunc.

@mhvk

@mhvk I'm not sure what the behaviour for np.where should / can be.

But we'll have to support Numpy 1.8 and 1.9 for many years to come, so a more concrete thing that could be done now is improve the documentation here.
How many of the numpy or scipy functions work with quantities. Is it the case that most numpy function do work and most scipy functions don't? It should be possible to check this mostly with an automatic script?
If so, I would find it helpful to state this in the docs and give a full list of numpy functions somewhere that don't ... these issues will still bite many Astropy devs and users in any case, but having the issue well-documented could help reduce the number of bites.

@cdeil

@cdeil - agreed, we should have a more complete list of things that work and do not work. @astrofrog added some specific tests in units that currently xfail, but I do not think this is complete. Unfortunately, making this list is not an entirely trivial endeavour (at least, I do not know how to do it in the bits of time that I have... Since every initiative is punishable, could you do a PR on adding the np.where issue to our documentation?). In the meantime, I think it is very useful to report issues we find upstream (ideally, of course, with solutions). For this particular one, see numpy/numpy#5095.

Navigation Menu

Monkey-patch non-ufunc Numpy functions #1274

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Monkey-patch non-ufunc Numpy functions #1274

Description

Activity

astrofrog commented on Sep 28, 2013

mhvk commented on Sep 28, 2013

mhvk commented on Oct 1, 2013

mhvk commented on Dec 31, 2013

mhvk commented on Dec 31, 2013

pv commented on Jan 2, 2014

astrofrog commented on Mar 11, 2014

astrofrog commented on Mar 11, 2014

astrofrog commented on Mar 11, 2014

mhvk commented on Mar 11, 2014

astrofrog commented on Mar 11, 2014

embray commented on May 2, 2014

cdeil commented on Sep 20, 2014

mhvk commented on Sep 21, 2014

cdeil commented on Sep 21, 2014

mhvk commented on Sep 21, 2014

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions