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

Implement Fermi 4FGL catalog spectral models and catalog #2297

Merged
merged 6 commits into from Jul 23, 2019

Conversation

@kaoriinakashima
Copy link
Contributor

@kaoriinakashima kaoriinakashima commented Jul 18, 2019

This PR adds a first version of the 4FGL catalog and source classes (see #2055), as well as a new spectral model. It only implements the spectral models (not sure if all of them), it's a first step, not a complete 4FGL implementation.

@cdeil cdeil self-assigned this Jul 18, 2019
@cdeil cdeil added the feature label Jul 18, 2019
@cdeil cdeil added this to the 0.13 milestone Jul 18, 2019
Copy link
Member

@cdeil cdeil left a comment

See inline comments.

Please also add one test for the new spectral model (copy & adapt the test for the similar existing model).

gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
gammapy/catalog/tests/test_fermi.py Outdated Show resolved Hide resolved
gammapy/spectrum/models.py Show resolved Hide resolved
@cdeil cdeil changed the title adding the Fermi 4FGL catalog and implementing the spectrum function Add Fermi 4FGL catalog spectral models Jul 19, 2019
@cdeil cdeil removed this from the 0.13 milestone Jul 22, 2019
@cdeil cdeil added this to the 0.14 milestone Jul 22, 2019
…t_fermi.py is ok, but the test_models.py is not completly ready yet.
gammapy/spectrum/models.py Outdated Show resolved Hide resolved
gammapy/spectrum/models.py Outdated Show resolved Hide resolved
gammapy/catalog/fermi.py Outdated Show resolved Hide resolved
@cdeil
Copy link
Member

@cdeil cdeil commented Jul 22, 2019

This PLSuperExpCutoff4FGL.evaluate implementation isn't correct yet, it's not working:

>>> from astropy import units as u 
>>> from gammapy.spectrum.models import PLSuperExpCutoff4FGL
>>> model = PLSuperExpCutoff4FGL()
>>> model(2 * u.TeV)
Traceback (most recent call last):
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 32, in get_converter
    scale = from_unit._to(to_unit)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 953, in _to
    "'{0!r}' is not a scaled version of '{1!r}'".format(self, other))
astropy.units.core.UnitConversionError: 'Unit("TeV2")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 147, in helper_dimensionless_to_dimensionless
    return ([get_converter(unit, dimensionless_unscaled)],
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 35, in get_converter
    from_unit, to_unit, get_current_unit_registry().equivalencies)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 890, in _apply_equivalencies
    unit_str, other_str))
astropy.units.core.UnitConversionError: 'TeV2' and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/work/code/gammapy2/gammapy/spectrum/models.py", line 51, in __call__
    return self.evaluate(energy, **kwargs)
  File "/Users/deil/work/code/gammapy2/gammapy/spectrum/models.py", line 1105, in evaluate
    cutoff = expfactor * np.exp(reference ** index_2 - energy ** index_2)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py", line 446, in __array_ufunc__
    converters, unit = converters_and_unit(function, method, *inputs)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/converters.py", line 166, in converters_and_unit
    converters, result_unit = ufunc_helper(function, *units)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 152, in helper_dimensionless_to_dimensionless
    .format(f.__name__))
astropy.units.core.UnitTypeError: Can only apply 'exp' function to dimensionless quantities

You currently have this code:

cutoff = expfactor * np.exp(reference ** index_2 - energy ** index_2)

So the error message makes sense, one can't call np.exp on energy to the power of something.

The model formula is in equation 3 in https://arxiv.org/pdf/1902.10045.pdf (please mention that in the model docstring as reference).

There you see that the factor a is inside the exp function and is supposed to have unit of MeV ** - index_2
You can try putting it there and see if it works then.

cutoff = np.exp(expfactor * (reference ** index_2 - energy ** index_2))

If not, e.g. because of float rounding errors in the unit computation, then we might have to do a more manual unit handling for this case.

gammapy/spectrum/tests/test_models.py Outdated Show resolved Hide resolved
Copy link
Contributor Author

@kaoriinakashima kaoriinakashima left a comment

This PLSuperExpCutoff4FGL.evaluate implementation isn't correct yet, it's not working:

>>> from astropy import units as u 
>>> from gammapy.spectrum.models import PLSuperExpCutoff4FGL
>>> model = PLSuperExpCutoff4FGL()
>>> model(2 * u.TeV)
Traceback (most recent call last):
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 32, in get_converter
    scale = from_unit._to(to_unit)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 953, in _to
    "'{0!r}' is not a scaled version of '{1!r}'".format(self, other))
astropy.units.core.UnitConversionError: 'Unit("TeV2")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 147, in helper_dimensionless_to_dimensionless
    return ([get_converter(unit, dimensionless_unscaled)],
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 35, in get_converter
    from_unit, to_unit, get_current_unit_registry().equivalencies)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 890, in _apply_equivalencies
    unit_str, other_str))
astropy.units.core.UnitConversionError: 'TeV2' and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/work/code/gammapy2/gammapy/spectrum/models.py", line 51, in __call__
    return self.evaluate(energy, **kwargs)
  File "/Users/deil/work/code/gammapy2/gammapy/spectrum/models.py", line 1105, in evaluate
    cutoff = expfactor * np.exp(reference ** index_2 - energy ** index_2)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py", line 446, in __array_ufunc__
    converters, unit = converters_and_unit(function, method, *inputs)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/converters.py", line 166, in converters_and_unit
    converters, result_unit = ufunc_helper(function, *units)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity_helper/helpers.py", line 152, in helper_dimensionless_to_dimensionless
    .format(f.__name__))
astropy.units.core.UnitTypeError: Can only apply 'exp' function to dimensionless quantities

You currently have this code:

cutoff = expfactor * np.exp(reference ** index_2 - energy ** index_2)

So the error message makes sense, one can't call np.exp on energy to the power of something.

The model formula is in equation 3 in https://arxiv.org/pdf/1902.10045.pdf (please mention that in the model docstring as reference).

There you see that the factor a is inside the exp function and is supposed to have unit of MeV ** - index_2
You can try putting it there and see if it works then.

cutoff = np.exp(expfactor * (reference ** index_2 - energy ** index_2))

If not, e.g. because of float rounding errors in the unit computation, then we might have to do a more manual unit handling for this case.

Done!

kaoriinakashima and others added 3 commits Jul 22, 2019
cdeil
cdeil approved these changes Jul 22, 2019
Copy link
Member

@cdeil cdeil left a comment

@kaoriinakashima - Awesome first Gammapy pull request. Thank you!

I think this is fine to be merged, but @adonath - could you please review this also and adjust if needed.
I'm not really sure if we absolutely can't avoid this extra spectral model (see description in 4FGL paper, they somehow also re-used a previous model), and I'm not sure if the uncertainties eval is tested - it looks like we don't have unit tests for evaluate_error in test_models.py.

Also, the 4FGL FITS file needs to be added in gammapy-extra and to the relevant data index files. @adonath - Would you be willing to do this also?
Alternatively - I have time again on Wednesday and could finish this up.

@cdeil cdeil assigned adonath and unassigned cdeil Jul 22, 2019
@cdeil cdeil added this to To do in gammapy.makers via automation Jul 22, 2019
@adonath
Copy link
Member

@adonath adonath commented Jul 23, 2019

Thanks a lot @kaoriinakashima, this is really nice work for a first PR. I don't have any further comments.

@cdeil I think its OK to keep the extra model and to have models in different parametrisations. I agree it's partly confusing to users to have this many models, but on the other hand it's also confusing to users to only provide models with a different parameterisation then a reference analysis. I'll go ahead now merge this PR and add the 4FGL to $GAMMAPY_DATA now.

@adonath adonath merged commit b127bc8 into gammapy:master Jul 23, 2019
9 checks passed
gammapy.makers automation moved this from To do to Done Jul 23, 2019
@adonath adonath removed this from the 0.14 milestone Jul 25, 2019
@adonath adonath added this to the 0.13 milestone Jul 25, 2019
@adonath adonath changed the title Add Fermi 4FGL catalog spectral models Implement Fermi 4FGL catalog spectral models Jul 25, 2019
@adonath adonath changed the title Implement Fermi 4FGL catalog spectral models Implement Fermi 4FGL catalog spectral models and catalog Jul 25, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants