Add XSPEC table absorption model to spectral table model #824

jjlk · 2016-12-18T15:02:53Z

Hi @cdeil and @joleroi ,
I added a class method to TableModel to import models from XSPEC (see flx2tab).

Maybe we could propose this format on https://gamma-astro-data-formats.readthedocs.io/en/latest/? @cdeil, @registerrier , some volunteers?

I also made optional the logarithm scale value on the y-axis and I added a plot function because, as I understood it, the SpectralModel plot method was really attached to have some physical units for the y-values.

Thanks!

cdeil · 2016-12-19T07:52:16Z

Please add $GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz so that tests run:
https://travis-ci.org/gammapy/gammapy/jobs/184942465#L997

You can directly add the file to master in gammapy-extra, unless you want it reviewed.
The usual things are: please add a README file what it is in the folder, and don't commit too large files without discussion (ask first if it's >> 1 MB).

Concerning format -- yes, I'd pretty much always would prefer formats are described over at gamma-astro-data-formats, instead of in the Gammapy docs. Especially for things like this, where the same data is of interest to several science tools.
I don't have time for this one. If you do, the start would be to either make a pull request with a very short and simple description of what you have here (as a starting point) or an issue proposing it, and CC a few people that might be able to give useful feedback.

cdeil

I left a few inline comments.

cdeil · 2016-12-19T07:54:18Z

gammapy/spectrum/models.py

+        filename = str(make_path(filename))
+
+        # Check if parameter value is in range
+        table_param = fits.open(filename)['PARAMETERS']


Can you try and use astropy.table.Table here?
We use it almost everywhere else instead of astropy.io.fits.BinTable directly.
To read a specific HDU:

Table.read(filename, hdu='MY_HDU_NAME')

cdeil · 2016-12-19T07:55:24Z

gammapy/spectrum/models.py

+        filename = '$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz'
+        table_model = TableModel.from_file(filename=filename, param=0.3)
+        """
+        energy_col = 'ENERG'


You're using this energy_col = 'ENERG' two times below, with string formatting.
IMO it would be simpler to remove it, and just use the column names directly below.

cdeil · 2016-12-19T07:56:25Z

gammapy/spectrum/models.py

+
+        Examples
+        --------
+        Fill table from a model


To get proper syntax highlighting, you have to indent the code in the docstring by four spaces, like this:

Here's some code: print('hello!')

cdeil · 2016-12-19T07:58:41Z

gammapy/spectrum/models.py

                                   fill_value=-np.Inf,
                                   bounds_error=False,
                                   kind='cubic')

+    @classmethod
+    def from_file(cls, filename, param):


What is param???

Mention it with one line (or add a Parameters section) in the docstring?

…eter y_label in plot function

jjlk · 2016-12-19T10:05:51Z

Hi @cdeil ,

Please add $GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz so that tests run:
https://travis-ci.org/gammapy/gammapy/jobs/184942465#L997

The files are already there. The test failed on Travis because of a precision problem. Locally it works fine for me. For now I just changed the parameter value for the test.

Concerning format -- yes, I'd pretty much always would prefer formats are described over at gamma-astro-data-formats, instead of in the Gammapy docs. Especially for things like this, where the same data is of interest to several science tools.
I don't have time for this one. If you do, the start would be to either make a pull request with a very short and simple description of what you have here (as a starting point) or an issue proposing it,
and CC a few people that might be able to give useful feedback.

Oké, I'll have a look at that later.

Thanks for your comments, I took all of them into account!

joleroi · 2016-12-19T11:02:51Z

I added a plot function because, as I understood it, the SpectralModel plot method was really attached to have some physical units for the y-values.

Doesn't it work if you pass flux_unit=''?

I had a discussion with @adonath about this and we decided it might be better do rename flux_unit -> y_unit everywhere. Its just not top priority I think, but if you want to do it, go ahead (would also need to be done in gammapy/spectrum/results.py)

joleroi

Minor stuff

joleroi · 2016-12-19T11:04:01Z

gammapy/spectrum/models.py

-                                   kind='cubic')
+    @classmethod
+    def from_file(cls, filename, param):
+        """A Table containing aborbed values from a XSPEC model


typo: absorbed (?)

joleroi · 2016-12-19T11:05:30Z

gammapy/spectrum/models.py

+        table_param = Table.read(filename, hdu='PARAMETERS')
+        param_min = table_param['MINIMUM']
+        param_max = table_param['MAXIMUM']
+        assert(param >= param_min and param <= param_max), "Parameter\


I've never seen a hard line break within a string literal coded like this, So I don't know if it's "allowed". What I've seen many time is

a = " Very long string, It's so long that we need a new " "line somewhere"

A coding rule for Gammapy (and Astropy / ... most projects): Never use assert to check inputs, raise an exception instead.

Here:

if param_min <= param <= param_max: raise ValueError(<message here>)

joleroi · 2016-12-19T11:07:19Z

gammapy/spectrum/models.py

+        energy = energy_bounds.log_centers
+
+        # Get spectrum values (no interpolation, take closest value for param)
+        table_spectra = Table.read(filename, hdu='SPECTRA')


Is it ok to hardcode the hdu names? Is this always the same for the format you're using, or is it better to pass the hdu names as kwargs?

cdeil · 2016-12-19T11:14:07Z

gammapy/spectrum/models.py

-                                   kind='cubic')
+    @classmethod
+    def from_file(cls, filename, param):
+        """A Table containing aborbed values from a XSPEC model


Maybe call the method read_<something> instead of from_file?
E.g. read_xspec_absorbed or read_ebl?

I think we always use the term "read" for classmethods reading from files, or at least I'm +1 to do that.

cdeil

Left two minor inline comments.

cdeil · 2016-12-19T11:16:30Z

gammapy/spectrum/models.py

+            File containing the XSPEC model
+        param : float
+            Model parameter value
+        Examples


I don't think this docstring will render correctly.
You need to put an empty line before the Examples heading, and I think also colon before the code block.

Can you please check locally with

python setup.py build_docs open docs/_build/html/index.html

and then check this docstring on the HGPS page?

We will add automated checks for docstrings in the future, there's a discussion for Astropy generally how to best do it. For now, checking docstrings is done manually.

cdeil · 2016-12-19T11:17:49Z

gammapy/spectrum/tests/test_models.py

+@requires_data('gammapy-extra')
+def test_table_model_from_file():
+    filename = '$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz'
+    energy_range = (0.03, 10)


No need for local variables here IMO, since only used once.
Just pass as keyword arguments below.

jjlk · 2016-12-19T12:58:19Z

Hi @joleroi and @cdeil ,

Doesn't it work if you pass flux_unit=''?
I had a discussion with @adonath about this and we decided it might be better do rename flux_unit -> y_unit everywhere. Its just not top priority I think, but if you want to do it, go ahead (would also need to be done in gammapy/spectrum/results.py)

I tried that and it didn't work. I didn't push a lot to investiguate where the problem was.

Maybe call the method read_ instead of from_file?
E.g. read_xspec_absorbed or read_ebl?

Done, read_xspec_model

I don't think this docstring will render correctly.

Fixed

A coding rule for Gammapy (and Astropy / ... most projects): Never use assert to check inputs, raise an exception instead.

Fixed

Is it ok to hardcode the hdu names? Is this always the same for the format you're using, or is it better to pass the hdu names as kwargs?

Yes, it is a standard

No need for local variables here IMO, since only used once.
Just pass as keyword arguments below.

Fixed

And multiple typos fixed.

Thanks!

cdeil

Looks good to me.

@joleroi - Do you want to do a final review, or should we merge?

(I'd suggest postponing plot y-axis clean-up to a separate issue / PR)

Modify table to handle absorption from XSPEC model and add a test

0ff1ad9

cdeil assigned joleroi Dec 19, 2016

cdeil added the feature label Dec 19, 2016

cdeil added this to the 0.6 milestone Dec 19, 2016

cdeil suggested changes Dec 19, 2016

View reviewed changes

Use Table instead of BinTable, improve doc, modify test, remove param…

2b6ce87

…eter y_label in plot function

joleroi reviewed Dec 19, 2016

View reviewed changes

cdeil reviewed Dec 19, 2016

View reviewed changes

Doc fixed, change frome_file to read and +

62af36e

cdeil approved these changes Dec 19, 2016

View reviewed changes

joleroi merged commit 1fec65b into gammapy:master Dec 19, 2016

jjlk deleted the xspec_model branch December 19, 2016 15:10

cdeil changed the title ~~Modify TableModel to handle absorption model from XSPEC table~~ Add XSPEC table absorption model to spectral table model Dec 20, 2016

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add XSPEC table absorption model to spectral table model #824

Add XSPEC table absorption model to spectral table model #824

jjlk commented Dec 18, 2016

cdeil commented Dec 19, 2016

cdeil left a comment

cdeil Dec 19, 2016

cdeil Dec 19, 2016

cdeil Dec 19, 2016

cdeil Dec 19, 2016

jjlk commented Dec 19, 2016

joleroi commented Dec 19, 2016

joleroi left a comment

joleroi Dec 19, 2016

joleroi Dec 19, 2016

cdeil Dec 19, 2016

joleroi Dec 19, 2016

cdeil Dec 19, 2016

cdeil left a comment

cdeil Dec 19, 2016

cdeil Dec 19, 2016

jjlk commented Dec 19, 2016

cdeil left a comment

Add XSPEC table absorption model to spectral table model #824

Add XSPEC table absorption model to spectral table model #824

Conversation

jjlk commented Dec 18, 2016

cdeil commented Dec 19, 2016

cdeil left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jjlk commented Dec 19, 2016

joleroi commented Dec 19, 2016

joleroi left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cdeil left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jjlk commented Dec 19, 2016

cdeil left a comment

Choose a reason for hiding this comment