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

Fix MapAxis interpolation FITS serialisation #1926

Merged
merged 5 commits into from Nov 15, 2018

Conversation

3 participants
@AtreyeeS
Contributor

AtreyeeS commented Nov 14, 2018

This PR addresses #1887

  1. Removed the special column names for the time.
  2. Added an INTERP%i (where i is the axis number) keyword to the header. If this exists, it is read appropriately. Otherwise, it defaults to log for energy and lin for other axes.

@AtreyeeS AtreyeeS requested review from registerrier and cdeil Nov 14, 2018

@cdeil

@AtreyeeS - thanks!

There's several nitpicky comments inline, but all really minor, overall this looks good to me.

@registerrier - Could you please also have a look?

energy_axis = MapAxis(nodes=np.logspace(-2., 2, 5), unit='TeV', name='energy', node_type=node_type, interp=interp)
time_axis = MapAxis(nodes=np.linspace(0.0, 1000.0, 4), unit='s', name='time', node_type=node_type, interp=interp)
geom = WcsGeom.create(skydir=(0, 0), binsz=1, npix=10, axes=[energy_axis, time_axis])
m = Map.from_geom(geom, unit='m2')

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

Can't you call Map.create directly?

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Changed.
Since I was checking just the geoms and not the full maps, Map.from_geom() was just clearer in my head. Obviously, unnecessary in the tests.

time_axis = MapAxis(nodes=np.linspace(0.0, 1000.0, 4), unit='s', name='time', node_type=node_type, interp=interp)
geom = WcsGeom.create(skydir=(0, 0), binsz=1, npix=10, axes=[energy_axis, time_axis])
m = Map.from_geom(geom, unit='m2')
m.write('test.fits', overwrite=True)

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

When writing to files from tests, it should always be done in a tempdir to avoid creating file in the current working directory. If you search Gammapy for tmpdir you will find many examples.

However, usually it's better to just serialise / deserialise in-memory to an HDUList, no need to create temp dirs and files which is slow. So suggest to do that.

I would also suggest to add asserts on the header keys you are adding here after serialisation. At the moment you just assert that you tround-trip the geom. You could change the serialisation format however you like and the test would still pass. If you add a few asserts on header keys or column names you can directly show in the test that the serialisation is correct, it will establish behaviour better. The round-trip assert for geom is good, keep that.

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Added assert for the INTERP keywords.
I have still kept the fits r/w (now in tmpdir) because it is a few lines of code to check the geom otherwise.

m.write('test.fits', overwrite=True)
m2 = Map.read('test.fits')
assert (m2.geom == m.geom)

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

Remove parentheses.

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Done

def test_read_write(node_type, interp):
energy_axis = MapAxis(nodes=np.logspace(-2., 2, 5), unit='TeV', name='energy', node_type=node_type, interp=interp)
time_axis = MapAxis(nodes=np.linspace(0.0, 1000.0, 4), unit='s', name='time', node_type=node_type, interp=interp)
geom = WcsGeom.create(skydir=(0, 0), binsz=1, npix=10, axes=[energy_axis, time_axis])

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

skydir=(0, 0) is default and irrelevant here. Remove.

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Removed

@pytest.mark.parametrize("node_type", ["edges", "center"])
@pytest.mark.parametrize("interp", ["log", "lin", "sqrt"])
def test_read_write(node_type, interp):
energy_axis = MapAxis(nodes=np.logspace(-2., 2, 5), unit='TeV', name='energy', node_type=node_type, interp=interp)

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

This line is very long. Did you run make black? Presumably it will break into many lines?
Maybe shorten it by passing args by position, or by changing np.logspace(-2., 2, 5) to the slightly shorter [1, 2]?

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

ran black. Breaks it into many lines.

@@ -122,7 +120,10 @@ def find_and_read_bands(hdu, header=None):
name = re.search("(.+)_MIN", cols[0]).group(1)
else:
name = cols[0]
intp = "INTERP%i" % (i+1)

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

"intp" is a pretty cryptic name. Maybe call it "interp_key" or something a bit more descriptive?

This comment has been minimized.

@registerrier

registerrier Nov 15, 2018

Contributor

The keyword length is limited to 8 characters unfortunately.

This comment has been minimized.

@cdeil

cdeil Nov 15, 2018

Member

I just meant the Python variable name, not the FITS header key name.

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Changed to interp_key

intp = "INTERP%i" % (i+1)
if header is not None:
if intp in header:
interp = header[intp]

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

Add input validation, i.e. raise ValueError or warning or silently skip if an interp that Gammapy doesn't understand occurs? Otherwise you put a value in, and it will probably crash or silently misbehave later? Might be useful to have a test for that?

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Added a log.warning()`. Will take default vaules in such cases

@@ -1349,6 +1350,10 @@ def _fill_header_from_axes(self, header):
else:
raise ValueError("Invalid node type {!r}".format(ax.node_type))
key_interp = "INTERP%i" % idx

This comment has been minimized.

@cdeil

cdeil Nov 14, 2018

Member

If it's just used once you could avoid the extra local variable since this is still short and readable:

header["INTERP%i" % idx] = ax._interp

This comment has been minimized.

@AtreyeeS

AtreyeeS Nov 15, 2018

Contributor

Now it is used twice, so keeping a variable name

@cdeil cdeil changed the title from Improvement in fits serialisation. to Fix MapAxis interpolation FITS serialisation Nov 14, 2018

@cdeil cdeil self-assigned this Nov 14, 2018

@cdeil cdeil added the bug label Nov 14, 2018

@cdeil cdeil added this to To do in Map analysis via automation Nov 14, 2018

@cdeil cdeil added this to the 0.9 milestone Nov 14, 2018

@cdeil

This comment has been minimized.

Member

cdeil commented Nov 15, 2018

@AtreyeeS - Are you still working on this? Or should I merge it in?

@cdeil cdeil removed the request for review from registerrier Nov 15, 2018

@AtreyeeS

This comment has been minimized.

Contributor

AtreyeeS commented Nov 15, 2018

@AtreyeeS - Are you still working on this? Or should I merge it in?

Done from my side

@cdeil

cdeil approved these changes Nov 15, 2018

@cdeil

cdeil approved these changes Nov 15, 2018

@AtreyeeS - Thanks!

@cdeil cdeil merged commit e2c6ee3 into gammapy:master Nov 15, 2018

0 of 3 checks passed

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

Map analysis automation moved this from To do to Done Nov 15, 2018

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