UnitConversionError in calculate_predicted_counts

[joleroi]

The following line in gammapy.spectrum.utils.calculate_predicted_counts is really strange

 true_energy = aeff.energy.data.to('TeV')

Obviously, the unit of an astropy quantity should not matter. However, if you change the line to

 true_energy = aeff.energy.data

the test

python setup.py test -t gammay/spectrum/tests/test_fit.py

fails with

UnitConversionError: 'TeV(13/10) / keV(13/10)' and '' (dimensionless) are not convertible

This looks like a bug in astropy (the unit should be convertible), but I cannot reproduce it with only astropy code.

[cdeil]

I get this, not sure if it's relevant or if the number at the end makes sense:

In [1]: import astropy.units as u

In [2]: u.Unit('TeV(13/10) / keV(13/10)')
WARNING: UnitsWarning: 'TeV(13/10) / keV(13/10)' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Out[2]: Unit("TeV(13/10) / keV(13/10)")

In [3]: u.Unit('TeV(13/10) / keV(13/10)').to('')
WARNING: UnitsWarning: 'TeV(13/10) / keV(13/10)' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Out[3]: 501187233627.27277

It should be possible to debug or reproduce this by printing what aeff.energy.data is.
Is it a Quantity object? What is it's unit exactly?

I don't think so, but it might be that there's some configuration in astropy.units that can result in an UnitConversionError when running tests, but otherwise is a warning when run from a normal script.

[joleroi]

Consider the following code

In [59]: energy1 = 5 * u.TeV

In [60]: energy2 = 6 * u.GeV

In [61]: index = -1.23125892847591837508761349875614 * u.Unit('')

In [62]: value = (energy1/energy2) ** index

In [65]: value.unit.__dict__
Out[65]: 
{'_bases': [Unit("GeV"), Unit("TeV")],
 '_decomposed_cache': None,
 '_hash': None,
 '_powers': [1.2312589284759183, -1.2312589284759183],
 '_scale': 1.0}

The issue could be cause be rounding differences in the _powers attribute. Possible solutions

• Always convert to dimensionless in PowerLaw like models before applying the index
• Always convert to fixed units in model evaluation

[cdeil]

I just now came across this oddity trying to help @vorugantia with an example:

>>> from gammapy.catalog import SourceCatalog3FGL
>>> cat = SourceCatalog3FGL()
>>> source = cat['3FGL J2021.1+3651']
>>> print(source.spectral_model)
ExponentialCutoffPowerLaw3FGL
ParameterList
Parameter(name='index', value=1.6358338594436646, unit=Unit(dimensionless), min=0, max=None, frozen=False)
Parameter(name='amplitude', value=1.2799539206298505e-10, unit=Unit("1 / (cm2 MeV s)"), min=0, max=None, frozen=False)
Parameter(name='reference', value=838.4884033203125, unit=Unit("MeV"), min=None, max=None, frozen=0)
Parameter(name='ecut', value=3042.514404296875, unit=Unit("MeV"), min=None, max=None, frozen=False)

Covariance: None
>>> import astropy.units as u
>>> source.spectral_model(1 * u.TeV)
<Quantity 1.8504322368224425e-148 MeV(5/8) / (cm2 s TeV(13/8))>
>>> source.spectral_model(1 * u.TeV).to('cm-2 s-1 TeV-1')
<Quantity 2.833153003692258e-152 1 / (cm2 s TeV)>

@joleroi @adonath - I think we have to discuss (offline) what to do about __call__ for spectral models, no?

[joleroi]

Postponing to 0.9

[adonath]

The calculate_predicted_counts does not exist anymore. @cdeil What do you think is odd in the example, that you showed above? I think it's exactly the way quantities are supposed to work...

[cdeil]

This is a very complicated unit that we return from Gammapy:

>>> source.spectral_model(1 * u.TeV)
<Quantity 1.8504322368224425e-148 MeV(5/8) / (cm2 s TeV(13/8))>

If possible, I think we should return simpler units here.
If not easily possible in the scheme how we handle units in Gammapy, I guess this might be a good candidate for a FAQ entry "Why does Gammapy return complicated units? How do I simplify them?"