Skip to content
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 HAWC 2HWC catalog #1108

Merged
merged 4 commits into from Aug 31, 2017
Merged

Add HAWC 2HWC catalog #1108

merged 4 commits into from Aug 31, 2017

Conversation

@pdeiml
Copy link
Contributor

@pdeiml pdeiml commented Aug 24, 2017

I was geared to the existing files for e.g. gamma cat, fermi and hess.
But some things where not clear to me how to transform them to hwc2, e.g. the description of the catalog. There, I added some TODOs.
@cdeil Please have a look at them and tell me what to do at these points. (We can do it via Skype if you want)


EDIT by @cdeil - this is a PR implementing #887

@cdeil cdeil self-requested a review Aug 24, 2017
@cdeil cdeil self-assigned this Aug 24, 2017
@cdeil cdeil added the feature label Aug 24, 2017
@cdeil cdeil added this to the 0.7 milestone Aug 24, 2017
Copy link
Member

@cdeil cdeil left a comment

@pdeiml - This is on a good way, welcome to Gammapy and thanks for the contribution!

I left a bunch of inline comments. Please address them and add the spectral and morphology model properties and test for those, and more questions if you have them, then ping me for a second review.

"""
HAWC 2FHL catalog
TODO: Add reference_id to the paper?

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

For the docstring at the top of gammapy/catalog/hawc.py I would suggest:

HAWC catalogs (https://www.hawc-observatory.org/)

Note that it doesn't appear in the HTML docs, so it's of limited use to put more info here.

You can check the new docs locally via

python setup.py buid_docs
make doc-show

I would suggest to put a link to ADS and paper reference for 2HWC in the class-level docstring of SourceCatalog2HWC. If you don't mind doing it, it should really be added for the Fermi catalogs as well, and that would be a small and simple addition that you could just do in this PR now that we're noticing / discussing it.

Here's a good example of how to add a reference to a docstring:

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed this.
Unfortunately, I can not build the documentation at my Laptop at home. A long error occurs.
On the office PC it worked yesterday. I don't know why.

ss = self.__str__(info=info)
print(ss)

def __str__(self, info='all'):

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

There is no way to pass parameters to __str__. So it doesn't make sense to have an info='all' argument here.
Look at the other catalogs how it's done, and if the same mistake is there as well, fix it or ask how to fix it if unclear.

This comment has been minimized.

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

I'll fix the __str__ to how I prefer it in a separate PR now. Coming in 15 min ..

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

See #1109 for how to write the __str__ and info methods.
(merged in master just now; but no need to rebase this PR)

ss += '{:20s} : {:1}'.format('Test radius', d['spec1_radius'])
return ss

# TODO: Add other information to def ___str___(self, info)???

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Can you show the printout for one source in a comment in the PR?
That's much easier to review than the code.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Currently it looks like this for the second source of the catalog:

*** Basic info ***

Catalog row index (zero-based) : 1
HWC2_FHL_name: : 2HWC J0631+169

*** Position info ***

Measurement:
RA : 97.998 deg
DEC : 16.997 deg
GLON : 195.614 deg
GLAT : 3.507 deg
Position error : 0.114 deg

*** Spectral info ***

Spectrum 1:
Flux at 7TeV : 6.71e-15 +- 1.47e-15 1 / (cm2 s TeV)
index : -2.570 +- 0.150
Test radius : 0.0 deg

Spectrum 2:
Flux at 7TeV : 4.87e-14 +- 6.85e-15 1 / (cm2 s TeV)
index : -2.230 +- 0.080
Test radius : 2.0 deg

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Looks mostly good. Some small suggestions:

  • 'HWC2_FHL_name' -> 'Source name'
  • '7TeV' -> '7 TeV' (insert space)
  • 'index' -> 'Spectral index'
  • '6.71e-15 +- 1.47e-15 1 / (cm2 s TeV)' -> '6.71e-15 +- 1.47e-15 cm-2 s-1 TeV-1' (write unit by hand; I find the number followed by '1 / ' hard to read and have seen others confused by it.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed.

The output is now:

*** Basic info ***

Catalog row index (zero-based) : 1
Source name: : 2HWC J0631+169

*** Position info ***

Measurement:
RA : 97.998 deg
DEC : 16.997 deg
GLON : 195.614 deg
GLAT : 3.507 deg
Position error : 0.114 deg

*** Spectral info ***

Spectrum 1:
Flux at 7 TeV : 6.71e-15 +- 1.47e-15 cm-2 s-1 TeV-1
Spectral index : -2.570 +- 0.150
Test radius : 0.0 deg

Spectrum 2:
Flux at 7 TeV : 4.87e-14 +- 6.85e-15 cm-2 s-1 TeV-1
Spectral index : -2.230 +- 0.080
Test radius : 2.0 deg

TODO: Add examples
"""

name = '2hwc' #TODO: Different name, e.g. '2hwc fhl' (if changed, needs to be changed in test_hawc.py as well)

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Remove TODO comment `name='2hwc' is good.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed


#TODO: initialise self.filename? In gammacat it is, but Hess and fermi it is not

table = Table.read(filename, format='ascii.ecsv', delimiter=' ')

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Remove delimiter argument. Not needed.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed.

def __init__(self, filename='$GAMMAPY_EXTRA/datasets/catalogs/2HWC.ecsv'):
filename = str(make_path(filename))

#TODO: initialise self.filename? In gammacat it is, but Hess and fermi it is not

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Remove TODO. not needed.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed


table = Table.read(filename, format='ascii.ecsv', delimiter=' ')

source_name_key = 'source_name' #TODO: Call it e.g. 2HWC_FHL_name?

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Keep column names as-is in the ECSV file. Remove TODO here.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed

import pytest
from ..hawc import SourceCatalog2HWC

@pytest.fixture #TODO: add here a scope like in test_gammacat.py?

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Yes, add session scope here. No need to read the files multiple times during a test run.

def hwc2_fhl():
filename='$GAMMAPY_EXTRA/datasets/catalogs/2HWC.ecsv' #TODO: Not sure whether this makes sence, __init__ of SourceCatalog2HWC
# has this string as default argument
return SourceCatalog2HWC(filename = filename)

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Run pep8 on your files and fix up style violations like whitespace on this line.
See http://pep8.org/ if you're not familar.
I've been using PyCharm (free edition: https://www.jetbrains.com/pycharm-edu/download/ ) for a few years, and never had to run pep8 from the command line. It's awesome, highly recommend it.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

The white spaces are only due to the comment at the end.
Will delete it, when it is clear what to do here.

# has this string as default argument
return SourceCatalog2HWC(filename = filename)

#TODO: Add here @requires_data?

This comment has been minimized.

@cdeil

cdeil Aug 24, 2017
Member

Yes, you need a @requires_data('gammapy-extra') here.
To see it fail locally without adding that, you can do

export GAMMAPY_EXTRA=""

and then this test will fail when trying to read the ECSV file.
If you put the requires_data decorator the test will be skipped.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

It works, but if GAMMAPY_EXTRA=""
then the tests are only skipped, they do not return an error?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Yip. That's what the requires_data decorator is for. Skip test if required test data is not available, don't error on that test.

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 25, 2017

There are still some TODOs that you missed.

test_hawc.py: Line 9
hawc.py: Line 50, Line 98 (ok, related to your comment above) and especially line 121 which is the description which is shown by executing 'source_catalogs.info()'


@pytest.fixture(scope='session')
def hwc2_fhl():
filename='$GAMMAPY_EXTRA/datasets/catalogs/2HWC.ecsv'#TODO: Not sure whether this makes sence, __init__ of SourceCatalog2HWC

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Indeed, you can remove the filename and just have a one-line return SourceCatalog2HWC()

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed.

d = self.data
ss = '\n*** Position info ***\n\n'

# TODO: Add here the SIMBAD position???

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

No. Remove todo.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed.

One source is represented by `~gammapy.catalog.SourceCatalogObjectGammaCat`
See: TODO: Add here the reference_id?

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Should I add the webadress or the reference_id?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Please add the link to ADS.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Fixed

Examples
--------
TODO: Add examples

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Remove it?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

As you like. Some catalog docstrings have an example: http://docs.gammapy.org/en/latest/api/gammapy.catalog.SourceCatalogGammaCat.html
Others (e.g. 3FGL don't).

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

I removed it.

"""

name = '2hwc'
description = ' ' # TODO: Add description

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Description for source_catalogs.info() necessary? I think yes

This comment has been minimized.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

pls look whether the description is ok

Parameters
----------
info : {'all', 'basic', 'position', 'spectra} #TODO: call spectra different

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Is the word 'spectra' as argument ok?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

I see that HGPS has "spec" and Fermi has "spectral".
I would suggest "spectrum" everywhere.
If you can make it consistent here, please go a head.
If no, just leave as-is and I'll do some cleanup in a follow-up commit when merging this PR.
(I almost always do this nowadays to avoid going back and forth on details when I could just as well make the edits directly when I review the code of the PR after opening it up in PyCharm.)

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Changed it.

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 25, 2017

I have no TODOs left, hence - if your agree with descrition etc. - this PR is

RTM

@cdeil
cdeil approved these changes Aug 25, 2017
Copy link
Member

@cdeil cdeil left a comment

The code that's there now looks good.

@pdeiml - A major point of having these classes in Gammapy at all is that they translate the source parameter infos to Gammapy source models. Do you have time to add this via a spectral_model and morphology_model property and corresponding tests (again following other catalogs as example)? (If no, then I would just add it by adding a commit to your PR in the coming days)

ss += '{:20s} : {:.3f} +- {:.3f}\n'.format('Spectral index', d['spec0_index'], d['spec0_index_err'])
ss += '{:20s} : {:1}\n\n'.format('Test radius', d['spec0_radius'])

if(np.isnan(d['spec1_dnde'])):

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Does this work? Shouldn't it be "spec2" here in the if instead of "spec1"?
Can you please, in a comment, share the printout for a source that has only 1 spec, and for a sources that has 2 specs?

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

It works, in the ecsv-file the spectra are called 'spec0' and 'spec1'.
I can change the info-output to
Spectrum 0
Spectrum 1
but I don't like counting from 0 for non-binary/ digital things.

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

As I said: please past the two examples here so that I can review without having to check locally.

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

This here is a source with only one spectrum:

print(cat[0])

*** Basic info ***

Catalog row index (zero-based) : 0
Source name: : 2HWC J0534+220

*** Position info ***

Measurement:
RA : 83.628 deg
DEC : 22.024 deg
GLON : 184.547 deg
GLAT : -5.783 deg
Position error : 0.057 deg

*** Spectral info ***

Spectrum 1:
Flux at 7 TeV : 1.85e-13 +- 2.38e-15 cm-2 s-1 TeV-1
Spectral index : -2.580 +- 0.010
Test radius : 0.0 deg

No second spectrum available for this source

and this here a source with two spectra:

*** Basic info ***

Catalog row index (zero-based) : 1
Source name: : 2HWC J0631+169

*** Position info ***

Measurement:
RA : 97.998 deg
DEC : 16.997 deg
GLON : 195.614 deg
GLAT : 3.507 deg
Position error : 0.114 deg

*** Spectral info ***

Spectrum 1:
Flux at 7 TeV : 6.71e-15 +- 1.47e-15 cm-2 s-1 TeV-1
Spectral index : -2.570 +- 0.150
Test radius : 0.0 deg

Spectrum 2:
Flux at 7 TeV : 4.87e-14 +- 6.85e-15 cm-2 s-1 TeV-1
Spectral index : -2.230 +- 0.080
Test radius : 2.0 deg

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Perfect!
(I was confused that "1" is the second, you could add a code comment saying that, or leave as-is)

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 25, 2017

I can add a spectral_model function. That is not difficult.
But how should I handle the fact that there may be two models?
Should I write two spectral_model method which return the first and the second model seperately?
Or one method which returns both spectra as e.g. a list of SpectralModels?

What should be written in the morphology method?
We have only the position as information but nothing further, haven't we?

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 25, 2017

For the spectra, I see a few options how to expose it, and don't have a preference on API yet (but might later after using it for a while), choose the one you prefer.

Either a function with a parameter to select, or two separate properties. I guess returnning a length 1 or 2 list with 1 or 2 spectral models is also fine.
It might also be useful to add a bool property has_two_spec (or int n_spec or whatever you like), no?


I thought they measured source extensions. But looking at the paper I see that they didn't really. So forget about the morphology model, the info doesn't exist.


return ss

# TODO: Add an @property?

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

I don't know what this @Property thing ist, but should I add it here?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Put this:

@property
def n_spectra(self):
    """Number of measured spectra (1 or 2)."""
    return 1 if np.isnan(self.data['spec1_dnde']) else 2

And Google "Python property" if you want to learn about that "property thing" (but not really needed, the other source classes contain working examples in the code / tests that show how it works)

errs['amplitude'] = data['spec0_dnde_err']
pars['index'] = data['spec0_index']*u.Unit('')
errs['index'] = data['spec0_index_err']*u.Unit('')
pars['reference'] = 7*u.TeV #TODO: Check whether 7Tev is correct by definition of

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

I think it has to be 7TeV because all fluxes in the catalog are at 7TeV. But I don't know exactly how this pars['reference'] is handled in line 120 and whether the model is correct afterwards.

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Yes, 7 TeV is correct here. Remove the TODO.


assert source.data['source_name'] == '2HWC J0534+220'

#TODO: Add this @pytest.mark.parametrize thing?

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

I will write a test method here but can you tell me what I have to write into this @pytest.mark.parametrize thing.
Thx

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

To learn how pytest works, please read this: https://github.com/jiffyclub/pytest-features#parametrized-tests or copy / execute the other working examples in gammapy.catalog. You can use it if you like, but I'm not really sure if it's useful here to parametrize tests. I would probably not do it and just have one test for a single-spec source and one separate test for a source with 2 spectra.

Copy link
Member

@cdeil cdeil left a comment

See inline comments.


return ss

# TODO: Add an @property?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Put this:

@property
def n_spectra(self):
    """Number of measured spectra (1 or 2)."""
    return 1 if np.isnan(self.data['spec1_dnde']) else 2

And Google "Python property" if you want to learn about that "property thing" (but not really needed, the other source classes contain working examples in the code / tests that show how it works)

errs['amplitude'] = data['spec0_dnde_err']
pars['index'] = data['spec0_index']*u.Unit('')
errs['index'] = data['spec0_index_err']*u.Unit('')
pars['reference'] = 7*u.TeV #TODO: Check whether 7Tev is correct by definition of

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Yes, 7 TeV is correct here. Remove the TODO.

return N

@property
def spectral_model(self):

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

It's usually considered bad (confusing / hard to deal with return value / error prone) to have functions / methods / properties that sometimes return one thing, and sometimes something else.
At the moment you sometimes return a SpectralModel object, and sometimes a list of two spectral models.

I left a few suggestions how to do it before, but maybe I think maybe always returning a list of spectral models would be good and you like that solution? So make an empty list at the top, and then append one or two SpectralModel objects and return that.

  • The implementation should be improved to avoid code duplication (currently you have the same code to copy over parameters to a SpectralModel 3 times -> should be there only once). See my make_2hwc.py script for an example how to do that.
  • The power-law amplitude should be assigned a quantity wiht the correct units. In the tests you can evaluate the model e.g. at 1 TEV or 7 TeV to confirm that it works and gives the expected result.
  • The docstring needs some more info, i.e. explain that it's a list of 1 or 2 spectral models, and give a sentence what they are (consult the paper description if needed).

assert source.data['source_name'] == '2HWC J0534+220'

#TODO: Add this @pytest.mark.parametrize thing?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

To learn how pytest works, please read this: https://github.com/jiffyclub/pytest-features#parametrized-tests or copy / execute the other working examples in gammapy.catalog. You can use it if you like, but I'm not really sure if it's useful here to parametrize tests. I would probably not do it and just have one test for a single-spec source and one separate test for a source with 2 spectra.

ss += '{:20s} : {:.3f} +- {:.3f}\n'.format('Spectral index', d['spec0_index'], d['spec0_index_err'])
ss += '{:20s} : {:1}\n\n'.format('Test radius', d['spec0_radius'])

if(np.isnan(d['spec1_dnde'])):

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Perfect!
(I was confused that "1" is the second, you could add a code comment saying that, or leave as-is)


return 1 if np.isnan(self.data['spec1_dnde']) else 2

#TODO: Maybe as public method?

This comment has been minimized.

@pdeiml

pdeiml Aug 25, 2017
Author Contributor

Shall this be a public method?

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Keep it private. Remove TODO.

source in the catalog.
"""

model_class = PowerLaw

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

No need for model_class = PowerLaw local variable. Please remove and use PowerLaw directly below.

model_class = PowerLaw
models = []

models.append(model_class(**(self._get_spec_pars(0)[0])))

This comment has been minimized.

@cdeil

cdeil Aug 25, 2017
Member

Suggest to first make the model and then put it in a list. Easier to read.

model = PowerLaw(**(self._get_spec_pars(0)[0]))
model.parameters.set_parameter_errors(self._get_spec_pars(0)[1])
models = [model]
@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 25, 2017

Should be fine now. By the way: The tests are working without errors.

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 28, 2017

Running pep8 on hawc.py still returns:

gammapy/catalog/hawc.py:78:13: E128 continuation line under-indented for visual indent
gammapy/catalog/hawc.py:87:17: E128 continuation line under-indented for visual indent

which is because of the line break in these two lines.
Is it ok like it is? Avoiding the two errors would lead to a indent of the newline which is very long.

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 29, 2017

@cdeil Please review

@pdeiml
Copy link
Contributor Author

@pdeiml pdeiml commented Aug 30, 2017

@cdeil Please review and tell me what the problem of travis CI is. It tells me that SourceCatalogObject2HWC has no getitem? But I don't know why this is needed and all other test_*.py scripts are written analogeously to test_hawc.py.
Thanks.

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 30, 2017

@pdeiml - I did some changes in 2ea3dc3 and pushed the commit to this PR.

I still want to improve the spectral parameter printout to avoid the code duplication and improve the tests, but I have to go now. Will finish this and merge tomorrow.

@cdeil
Copy link
Member

@cdeil cdeil commented Aug 31, 2017

Some more fixes and improvements in 99161ff . Merging now.

@pdeiml - Thank you for making your first contribution to Gammapy; it's a very nice addition!

@cdeil cdeil merged commit 52e46a4 into gammapy:master Aug 31, 2017
0 of 2 checks passed
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 changed the title Add HAWC 2FHL catalog Add HAWC 2HWC catalog Aug 31, 2017
@pdeiml pdeiml deleted the pdeiml:add_hawc_2fhl_catalog branch Aug 31, 2017
@cdeil
Copy link
Member

@cdeil cdeil commented Aug 31, 2017

Docs for this are now online: http://docs.gammapy.org/en/latest/api/gammapy.catalog.SourceCatalog2HWC.html

@pdeiml - I've added you to the Gammapy contributor list: 34f2e93 (will appear in the docs after the next docs build)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants