Converting a astromodels Model¶
The core concept behind that plugin is the model converter that maps a astromodels Model to a gammapy SkyModel, so we can perform the forward folding of said model in gammapy. Let’s have a look!
We start by importing everything for a simple astromodels Model
[1]:
from astromodels.core.model import Model
from astromodels.functions import Powerlaw
from astromodels.functions import Disk_on_sphere
from astromodels.sources import ExtendedSource, PointSource
from astromodels.utils.io import display
import astropy.units as u
13:25:22 WARNING The naima package is not available. Models that depend on it will not be functions.py:43 available
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:65 will not be available.
WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:33 available
which we then define
[2]:
pl = Powerlaw()
disk = Disk_on_sphere()
source = ExtendedSource("test_source", spectral_shape=pl, spatial_shape=disk)
disk.lon0 = 1 * u.deg
disk.lat0 = 10 * u.deg
model = Model(source)
and take a look at what we created
[3]:
display(model)
| N | |
|---|---|
| Point sources | 0 |
| Extended sources | 1 |
| Particle sources | 0 |
Free parameters (5):
| value | min_value | max_value | unit | |
|---|---|---|---|---|
| test_source.Disk_on_sphere.lon0 | 1 | 0 | 360 | deg |
| test_source.Disk_on_sphere.lat0 | 10 | -90 | 90 | deg |
| test_source.Disk_on_sphere.radius | 15 | 0 | 20 | deg |
| test_source.spectrum.main.Powerlaw.K | 1 | 1e-30 | 1000 | keV-1 s-1 cm-2 |
| test_source.spectrum.main.Powerlaw.index | -2.01 | -10 | 10 |
Fixed parameters (1):
(abridged. Use complete=True to see all fixed parameters)
Properties (0):
(none)
Linked parameters (0):
(none)
Independent variables:
(none)
Linked functions (0):
(none)
B.E.A.U.T.I.F.U.L.
Ok, now for the conversion - first we need to import stuff again:
[4]:
from gammapy_plugin.converter import AstromodelConverter
13:25:23 INFO Starting 3ML! __init__.py:44
WARNING WARNINGs here are NOT errors __init__.py:45
WARNING but are inform you about optional packages that can be installed __init__.py:46
WARNING to disable these messages, turn off start_warning in your config file __init__.py:47
WARNING no display variable set. using backend for graphics without display (agg) __init__.py:53
13:25:24 WARNING ROOT minimizer not available minimization.py:1208
WARNING Multinest minimizer not available minimization.py:1218
WARNING PyGMO is not available minimization.py:1228
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:126 software installed and configured?
WARNING No fermitools installed lat_transient_builder.py:44
WARNING Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:335 performances in 3ML
WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:335 performances in 3ML
WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:335 performances in 3ML
and it is indeed as simple as running
[5]:
conv = AstromodelConverter(model=model)
We now have an instance of the converter which has already perfomed the conversion up on initialization. We only created one single source so we only have one single model which we can index
[6]:
conv.gammapy_models[0]
[6]:
SkyModel
Name : test_source
Datasets names : None
Spectral model type : SpectralModelConverted
Spatial model type : SpatialModelConverted
Temporal model type :
Parameters:
test_source.spectrum.main.Powerlaw.K : 1.000 +/- 0.00 1 / (keV s cm2)
test_source.spectrum.main.Powerlaw.piv (frozen): 1.000 keV
test_source.spectrum.main.Powerlaw.index : -2.010 +/- 0.00
lon0 : 1.000 +/- 0.00 deg
lat0 : 10.000 +/- 0.00 deg
radius : 15.000 +/- 0.00 deg
Nice! We see that both our lon0and lat0 values are the ones we assigned before :)
Now let’s update them. If we simply run
[7]:
pl.K.value = 10
conv.gammapy_models[0]
[7]:
SkyModel
Name : test_source
Datasets names : None
Spectral model type : SpectralModelConverted
Spatial model type : SpatialModelConverted
Temporal model type :
Parameters:
test_source.spectrum.main.Powerlaw.K : 1.000 +/- 0.00 1 / (keV s cm2)
test_source.spectrum.main.Powerlaw.piv (frozen): 1.000 keV
test_source.spectrum.main.Powerlaw.index : -2.010 +/- 0.00
lon0 : 1.000 +/- 0.00 deg
lat0 : 10.000 +/- 0.00 deg
radius : 15.000 +/- 0.00 deg
We see that nothing has changed :/ We have to invoke the update() function:
[8]:
conv.update()
conv.gammapy_models[0]
[8]:
SkyModel
Name : test_source
Datasets names : None
Spectral model type : SpectralModelConverted
Spatial model type : SpatialModelConverted
Temporal model type :
Parameters:
test_source.spectrum.main.Powerlaw.K : 10.000 +/- 0.00 1 / (keV s cm2)
test_source.spectrum.main.Powerlaw.piv (frozen): 1.000 keV
test_source.spectrum.main.Powerlaw.index : -2.010 +/- 0.00
lon0 : 1.000 +/- 0.00 deg
lat0 : 10.000 +/- 0.00 deg
radius : 15.000 +/- 0.00 deg
Perfect!
This calls the indiviudal update methods of the sources which then update all the models.
Point Sources¶
In case you are only analysing a spectrum dataset and do not care about pixelatiion on a map or the PSF at all you can set the convert_ps flag to False.
By default this is set to true, meaning when you transform an astromodels PointSource you will get a PointSpatialModel (see e.g. here):
[9]:
pl = Powerlaw()
ps = PointSource("test_ps", spectral_shape=pl, ra=0, dec=12.3)
model_ps = Model(ps)
conv_ps = AstromodelConverter(model_ps)
conv_ps.gammapy_models[0]
[9]:
SkyModel
Name : test_ps
Datasets names : None
Spectral model type : SpectralModelConverted
Spatial model type : PointSourceModelConverted
Temporal model type :
Parameters:
test_ps.spectrum.main.Powerlaw.K : 1.000 +/- 0.00 1 / (keV s cm2)
test_ps.spectrum.main.Powerlaw.piv (frozen): 1.000 keV
test_ps.spectrum.main.Powerlaw.index : -2.010 +/- 0.00
lon_0 (frozen): 0.000 deg
lat_0 (frozen): 12.300 deg
If you set it to False you will not get any spatial component:
[10]:
pl = Powerlaw()
ps = PointSource("test_ps", spectral_shape=pl, ra=0, dec=12.3)
model_ps = Model(ps)
conv_ps = AstromodelConverter(model_ps, convert_ps=False)
conv_ps.gammapy_models[0]
[10]:
SkyModel
Name : test_ps
Datasets names : None
Spectral model type : SpectralModelConverted
Spatial model type :
Temporal model type :
Parameters:
test_ps.spectrum.main.Powerlaw.K : 1.000 +/- 0.00 1 / (keV s cm2)
test_ps.spectrum.main.Powerlaw.piv (frozen): 1.000 keV
test_ps.spectrum.main.Powerlaw.index : -2.010 +/- 0.00
Plotting¶
Stealing from the gammapy docs we can now also use the same plotting routines for the models:
But first let’s update our parameters to something in the “gammapy energy range”
[11]:
pl.K = 1 * u.Unit("TeV-1 cm-2 s-1")
pl.piv = 1 * u.TeV
and please be aware that astromodels uses \(K\times\left(\frac{E}{piv}\right)^{i}\) not \(K\times\left(\frac{E}{piv}\right)^{-i}\) for a powerlaw definition.
[12]:
from astropy import units as u
import matplotlib.pyplot as plt
energy_bounds = [0.1, 100] * u.TeV
conv_ps.gammapy_models[0].spectral_model.plot(energy_bounds)
plt.grid(which="both")
We can also do the same for the Disk spatial model before
[13]:
from gammapy.maps import Map, WcsGeom
pl = Powerlaw()
disk = Disk_on_sphere()
source = ExtendedSource("test_source", spectral_shape=pl, spatial_shape=disk)
disk.lon0 = 20 * u.deg
disk.lat0 = 0 * u.deg
disk.radius = 1 * u.deg
model = Model(source)
conv = AstromodelConverter(model=model)
lon_0 = disk.lon0.value
lat_0 = disk.lat0.value
reval = 2 * disk.radius.value
dr = 0.02
geom = WcsGeom.create(
skydir=(lon_0, lat_0),
binsz=dr,
width=(2 * reval, 2 * reval),
frame="icrs",
)
fig, ax = plt.subplots(1, figsize=(9, 6))
meval = conv.gammapy_models[0].spatial_model.evaluate_geom(geom)
Map.from_geom(geom=geom, data=meval.value, unit=meval.unit).plot(ax=ax)
plt.tight_layout()