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

Add units to gammapy.maps #1374

Merged
merged 9 commits into from Apr 16, 2018

Conversation

2 participants
@registerrier
Contributor

registerrier commented Apr 13, 2018

This is a first commit to add units to map. This contains:

  • addition of Map.unit as a string in the constructor and in the relevant Map.createmethods
  • addition of Map.quantityproperty and setter
  • addition of Map.unit in the header exported and imported respectively in Map.to_hdu() and in Map.from_hdu()
  • Tests for creation and hdu read/write

Still need to add unit in Map.fill_map_idx.

Note that the Map.quantity setter forces the unit so one can change the unit with Map.quantity = some_quantity. If this is not wanted one can add a check on the unit consistency. Or even remove the setter.

EDIT : this PR implements #1206

@registerrier registerrier added this to the 0.8 milestone Apr 13, 2018

@registerrier registerrier self-assigned this Apr 13, 2018

@registerrier registerrier requested review from cdeil and woodmd Apr 13, 2018

@cdeil

@registerrier - Thanks!

I left some inline comments.

You're trying to support both unit='' and unit=None, right?
I think that's a bad idea, you can already see in this PR how it increases the code needed for input handling and tests.
Suggest to just pick one and only support that.

map_out = cls(geom, meta=meta)
# Read unit
if 'UNIT' in hdu.header:

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

unit = hdu.header.get('UNIT', None)

This comment has been minimized.

@registerrier
@@ -102,7 +104,14 @@ def from_hdu(cls, hdu, hdu_bands=None):
shape_wcs = tuple([np.max(geom.npix[0]),
np.max(geom.npix[1])])
meta = cls._get_meta_from_header(hdu.header)
map_out = cls(geom, meta=meta)
# Read unit

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

Remove comment "read unit".
Comments that just say the same thing as the code are useless, even harmful, because 2x as many lines and very hard to keep comments updated.
So IMO we should strive for clean code and only add comments where they add extra value over the code.

This comment has been minimized.

@registerrier

registerrier Apr 13, 2018

Contributor

OK, indeed this one is pretty obvious

@@ -78,7 +80,12 @@ def from_hdu(cls, hdu, hdu_bands=None):
# TODO: Should we support extracting slices?
meta = cls._get_meta_from_header(hdu.header)
map_out = cls(hpx, None, meta=meta)
if 'UNIT' in hdu.header:

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

unit = hdu.header.get('UNIT', None)

This comment has been minimized.

@registerrier

registerrier Apr 13, 2018

Contributor

Done

self._geom = geom
self._data = data
if isinstance(data,Quantity):
self.unit = Quantity.unit.to_string()

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

I would suggest to call the data member self._unit to be consistent with the rest of the gammapy.maps codebase

This comment has been minimized.

@registerrier

registerrier Apr 13, 2018

Contributor

OK. Done.
I added a Map.unit property that return a ~astropy.units.Unit

@registerrier registerrier changed the title from First commit for maps.unit to Add units to gammapy.maps Apr 13, 2018

@cdeil

@registerrier - What about the None / empty string?
Otherwise looks good to me.

@@ -33,11 +34,24 @@ class Map(object):
Data array
meta : `~collections.OrderedDict`
Dictionary to store meta data.
unit : str or `~astropy.unit.Unit`

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

astropy.unit -> astropy.units

This comment has been minimized.

@registerrier
@property
def unit(self):
"""Return map unit as `~astropy.unit.Unit`"""

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

astropy.unit -> astropy.units

(I always check HTML docs locally for things I add)

self._data = data.value
else:
self._data = data
if unit is None:

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

You still have a mix of unit=None and unit=''.

Suggest to change to unit='' every to keep it simple.
Remember the Python zen "there should be one and only one way to do it".
Allowing multiple ways to do the same thing always leads to problems (only in coding, in life of course diversity is good).

This comment has been minimized.

@registerrier

registerrier Apr 13, 2018

Contributor

OK, right. Will do now. I have to adapt the tests then.

@registerrier

This comment has been minimized.

Contributor

registerrier commented Apr 13, 2018

This new version has a treatment of units for Map.fill_by_idx(idx, weights) with some test.
Now we check if weights is consistent with the map unit.
Similarly Map.coaddchecks that units are equivalent and convert to the unit of the target map.
Tests have been added for WcsNDMap and HpxNDMap

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2018

@registerrier - I've done some more small edits in 3f7730c

@woodmd - Do you have time to review this soon?
If no, I'd like to merge it in now, and then comments will still be welcome / addressed any time later.
(we have a presentation of Gammapy at the CTA meeting next week, and want to have a working 3D analysis in Gammapy master early next week, that's why we're rushing a bit)

@@ -410,6 +415,16 @@ def fill_by_idx(self, idx, weights=None):
msk = idx_local[0] >= 0
idx_local = [t[msk] for t in idx_local]
if weights is not None:
if isinstance(weights, Quantity):
if self.unit.is_equivalent(weights.unit):
weights = np.asarray(weights.to(self.unit).value, dtype=self.data.dtype)

This comment has been minimized.

@cdeil

cdeil Apr 13, 2018

Member

I think we should try to simplify these checks and handling for weights.

@registerrier - Why are you calling np.asarray and forcing a dtype?
That makes a copy, i.e. is inefficient.
There is no reason to introduce a copy and dtype change for weights, no?

Another thing I note is that these 10 complex lines of type checking are duplicated for HPX and WCS.

@registerrier - Can you try to remove the duplication (few line utility function?) and also try to simplify the logic to have less code here?

I haven't thought about it much, but shouldn't Astropy do the quantity handling for you and through an error if they don't match? Why do you have to do quantity -> array -> quantity and raise an error yourself here?

This comment has been minimized.

@registerrier

registerrier Apr 13, 2018

Contributor

I simply reused here the code as initially written in gammapy.maps, to ensure stable behavior at least for empty units. I think for the sake of stability I would leave as is for the moment.
Note BTW that from np.asarray documentation, no copy is performed if the input is already an ndarray with matching dtype and order. So the inefficiency is likely minimal.

Regarding the quantity-array conversion. The code is making some check for multiple entries to fill in a given bin. Hence the np.bincount, later in the code. I don"t know if it works well with weights of type Quantity.

Again if we are to make rapid changes to have 3D working, we should postpone this type of improvement for later, no?

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2018

My main remaining question here is this: we always store unit as string, but we only expose it as a Unit object? Is that a good idea? Or should we change that somehow?

Note that units can store a scale:

In [1]: from astropy import units as u

In [2]: u.Unit('42 second')
Out[2]: Unit("42 s")

In [3]: u.Unit('42 second').scale
Out[3]: 42.0

In [4]: str(u.Unit('42 second'))
Out[4]: '42 s'

In [5]: u.Unit('42 second').to_string('fits')
---------------------------------------------------------------------------
UnitScaleError                            Traceback (most recent call last)
<ipython-input-5-0a151fedd8f0> in <module>()
----> 1 u.Unit('42 second').to_string('fits')

~/Library/Python/3.6/lib/python/site-packages/astropy/units/core.py in to_string(self, format)
    601 
    602         f = unit_format.get_format(format)
--> 603         return f.to_string(self)
    604 
    605     def __format__(self, format_spec):

~/Library/Python/3.6/lib/python/site-packages/astropy/units/format/fits.py in to_string(cls, unit)
    122                     "The FITS unit format is not able to represent scales "
    123                     "that are not powers of 10.  Multiply your data by "
--> 124                     "{0:e}.".format(unit.scale))
    125             elif unit.scale != 1.0:
    126                 parts.append('10**{0}'.format(int(base)))

UnitScaleError: The FITS unit format is not able to represent scales that are not powers of 10.  Multiply your data by 4.200000e+01.

I think we had problems and confusion with that before.
So I'd prefer we invent some fool-safe way to only have units with the default scale=1, everywhere in Gammapy.

The FITS I/O should probably be changed to call to_string('fits')?
Also in the initialisers?

I'll leave this PR open for a day or two and would like to resolve those points; so @woodmd @joleroi @registerrier - please comment or make changes by attaching commits here.
But Sunday evening latest we need to merge this so that we can continue using this.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 16, 2018

There's a fail now in test_make_map_exposure_true_energy and test_make_map_fov_background :
https://travis-ci.org/gammapy/gammapy/jobs/366198442#L3828
I'll have a look now, and then merge this in.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 16, 2018

I rebased to get the gammapy/cube/new.py into this branch, and fixed the fails there with some edits in 9517cb8 .

@cdeil

cdeil approved these changes Apr 16, 2018

I've added one missing requires scipy in a test in dc48302 .

And I've simplified the weights handling in map fill in 609ebfb .
I think it's pretty much equivalent with the old version (at least all tests still pass), but simpler.
@registerrier - It's very well possible that there's a problem I don't understand here yet, and with the new implementation. If so, please open a new PR adding a test case that shows the problem, and then update the implementation again.

Merging now.

@cdeil cdeil merged commit 318ffd7 into gammapy:master Apr 16, 2018

0 of 2 checks passed

continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details

@cdeil cdeil added this to Done in Map analysis via automation Apr 16, 2018

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