New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some NumPy and SciPy operations drop units from arrays #59

Open
namurphy opened this Issue Jul 20, 2017 · 11 comments

Comments

Projects
None yet
4 participants
@namurphy
Copy link
Member

namurphy commented Jul 20, 2017

A known issue with astropy.units is that some NumPy operations take in arguments that have units, but return arguments that do not have units even though they should. This issue requires an upstream fix in NumPy (e.g., numpy/numpy#4164). The SciPy 2017 talk on MetPy's Choice of Unit Support delves into this in great detail. An example from the talk:

>>> import numpy as np
>>> from astropy import units
>>> a = np.array([3.]) * units.m
>>> b = np.array([4.]) * units.m
>>> np.concatenate((a, b))
<Quantity [ 3., 4.] (Unit not initialised)>

There are some workarounds to some of these problems.

>>> np.isclose(500* u.km/u.s, 300 * u.km / u.s
UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)
>>> np.isclose(500* u.km/u.s, 300 * u.km / u.s, atol=1e-8 * u.mm / u.s)
array([False], dtype=bool)

We cannot do anything from within PlasmaPy, but will need to wait for fixes within NumPy (e.g., new support for __array_ufunc__ will improve overriding behavior for ufuncs).

@namurphy

This comment has been minimized.

Copy link
Member

namurphy commented Jul 20, 2017

A takeaway point: when writing functions that take Astropy Quantities as arguments, always include a unit test to test units.

>>> density_of_wombats = 18.7 * u.parsec**-3
>>> velocity_of_wombats = 0.04 * c
>>> flux_of_wombats = density_of_wombats * velocity_of_wombats  # second golden rule of astrophysics
>>> assert flux_of_wombats.si.unit == units.m**-2 * units.s**-1
@Cadair

This comment has been minimized.

Copy link
Contributor

Cadair commented Jul 21, 2017

iirc the __array_ufunc__ stuff is in astropy 2.0 / numpy 1.13 so you could require these versions to avoid some issues.

There is a quantity aware version of isclose and assert_all_close in astropy.test.helpers I think.

Also a good way to make sure your functions return the correct units is to use the return annotations and u.quantity_input in astropy 2.0 which now converts the output of a function to match the return annotation.

@SolarDrew

This comment has been minimized.

Copy link
Contributor

SolarDrew commented Jul 21, 2017

Given that extra functionality, I propose again that we only support astropy 2.x.

@namurphy

This comment has been minimized.

Copy link
Member

namurphy commented Jul 24, 2017

It was decided during our last telecon that we shall support only Astropy 2.0 and later, and NumPy 1.13 and later. It was also decided that Nick's puns are hilarious and that there should be more of them in the code base.

@SolarDrew

This comment has been minimized.

Copy link
Contributor

SolarDrew commented Jul 24, 2017

For the record, I whole-heartedly support both of these decisions.

@namurphy

This comment has been minimized.

Copy link
Member

namurphy commented Jul 24, 2017

Once Astropy Quantities can work with scipy.special.erf, then we should:

  • Refactor plasma_dispersion_function so that it works with Quantities directly and we don't need to use .value to get the numerical value(s) as follows:
    if isinstance(zeta, u.Quantity):
        if zeta.unit == u.dimensionless_unscaled:
            zeta = zeta.value

See: astropy/astropy#6390

@namurphy namurphy changed the title Some NumPy operations drop units from arrays Some NumPy and SciPy operations drop units from arrays Jul 24, 2017

@blalterman

This comment has been minimized.

Copy link

blalterman commented Aug 8, 2017

Would you mind if I added my 2 cents to the discussion?

I happen to be visiting @namurphy and we're discussing several piles of code I'd like to contribute to PlasmaPy. There are two concerns that have come up in our discussions. (1) How will data analysts and experimentalists use PlasmaPy? (2) How can we leverage already existing standards? I bring up (1) in #86. I'd like to focus on (2) here.

There is a lot of development in projects like pandas and xarray that handle a large number of proposed developments I've seen in the GitHub PlasmaPy e-mails I've received. Interpolation is just one example. Why re-implement an interpolation method in a new class that has to handle a whole bunch of things that others have already written into their subclasses. Pandas' DataFrame.interpolate is a great example.

I bring up interpolation as it relates to units. There's a big discussion about how to implement units in PlasmaPy. Right now, we're busy creating a brand new Plasma object and developing it. At some point, we're going to redo a huge pile of development that already exists in Pandas or Xarray. Why not spend the time we'd be using to reinvent the wheel and just use it to implement units in a functional way for Pandas, Xarray, or even Numpy? That way, we can leverage the development they have done and not worry about writing our own test cases for things they have. After all, isn't PlasmaPy part of the growing set of open-source python projects?

@SolarDrew

This comment has been minimized.

Copy link
Contributor

SolarDrew commented Aug 9, 2017

Would you mind if I added my 2 cents to the discussion?

Please always feel free, it's an open forum :)

Why not spend the time we'd be using to reinvent the wheel and just use it to implement units in a functional way for Pandas, Xarray, or even Numpy?

My understanding is that none of those things is trivial. They may be slightly more efficient in the long term, but we also need to remember the short-term goal of providing useful funcitonality to users now.

After all, isn't PlasmaPy part of the growing set of open-source python projects?

I don't see that it's any less so if we decide existing tools don't serve our purposes and implement new ones that do.

To get back to the main topic of the issue, did we decide that using only Astropy >2.0 and Numpy >1.13 solves the problem, or are there still operations that don't output the untis correctly?

@Cadair

This comment has been minimized.

Copy link
Contributor

Cadair commented Aug 9, 2017

There are still issues with numpy 1.13 and astropy 2.0 but at least there are less.

@SolarDrew

This comment has been minimized.

Copy link
Contributor

SolarDrew commented Aug 9, 2017

Fair enough. I think we should log these as individual issues as and when we find them (so we can tick them off) and refer back to this one.

@namurphy

This comment has been minimized.

Copy link
Member

namurphy commented Sep 17, 2017

If anyone happens to run into this problem, the first thing I would suggest is making sure you're using the most recent versions of NumPy and Astropy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment