Note
Click here to download the full example code or to run this example in your browser via Binder
Models#
This is an introduction and overview on how to work with models in Gammapy.
The sub-package modeling contains all the functionality
related to modeling and fitting data. This includes spectral, spatial
and temporal model classes, as well as the fit and parameter API.The
models follow a naming scheme which contains the category as a suffix to
the class name. An overview of all the available models can be found in
the model gallery.
Note that there are separate tutorials,
Modelling and
Fitting that explains about
modeling, the Gammapy modeling and fitting framework. You
have to read that to learn how to work with models in order to analyse
data.
Setup#
# %matplotlib inline
import numpy as np
from astropy import units as u
import matplotlib.pyplot as plt
from IPython.display import display
from gammapy.maps import Map, MapAxis, WcsGeom
Check setup#
from gammapy.utils.check import check_tutorials_setup
check_tutorials_setup()
System:
python_executable : /Users/terrier/Code/anaconda3/envs/gammapy-dev/bin/python
python_version : 3.8.13
machine : x86_64
system : Darwin
Gammapy package:
version : 1.0rc2
path : /Users/terrier/Code/gammapy-dev/gammapy/gammapy
Other packages:
numpy : 1.22.4
scipy : 1.9.3
astropy : 5.1
regions : 0.6
click : 8.1.3
yaml : 6.0
IPython : 8.4.0
jupyterlab : 3.4.8
matplotlib : 3.5.3
pandas : 1.5.0
healpy : 1.16.1
iminuit : 2.17.0
sherpa : 4.15.0
naima : 0.10.0
emcee : 3.1.3
corner : 2.2.1
Gammapy environment variables:
GAMMAPY_DATA : /Users/terrier/Code/gammapy-dev/gammapy-data
Spectral models#
All models are imported from the models namespace.
Let’s start with a PowerLawSpectralModel:
from gammapy.modeling.models import PowerLawSpectralModel
pwl = PowerLawSpectralModel()
print(pwl)
PowerLawSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral index 2.0000e+00 ... nan False False
spectral amplitude 1.0000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral reference 1.0000e+00 TeV ... nan True False
To get a list of all available spectral models you can import and print the spectral model registry or take a look at the spectral-models-gallery
from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY
print(SPECTRAL_MODEL_REGISTRY)
Registry
--------
ConstantSpectralModel : ['ConstantSpectralModel', 'const']
CompoundSpectralModel : ['CompoundSpectralModel', 'compound']
PowerLawSpectralModel : ['PowerLawSpectralModel', 'pl']
PowerLaw2SpectralModel : ['PowerLaw2SpectralModel', 'pl-2']
BrokenPowerLawSpectralModel : ['BrokenPowerLawSpectralModel', 'bpl']
SmoothBrokenPowerLawSpectralModel : ['SmoothBrokenPowerLawSpectralModel', 'sbpl']
PiecewiseNormSpectralModel : ['PiecewiseNormSpectralModel', 'piecewise-norm']
ExpCutoffPowerLawSpectralModel : ['ExpCutoffPowerLawSpectralModel', 'ecpl']
ExpCutoffPowerLaw3FGLSpectralModel : ['ExpCutoffPowerLaw3FGLSpectralModel', 'ecpl-3fgl']
SuperExpCutoffPowerLaw3FGLSpectralModel : ['SuperExpCutoffPowerLaw3FGLSpectralModel', 'secpl-3fgl']
SuperExpCutoffPowerLaw4FGLDR3SpectralModel: ['SuperExpCutoffPowerLaw4FGLDR3SpectralModel', 'secpl-4fgl-dr3']
SuperExpCutoffPowerLaw4FGLSpectralModel : ['SuperExpCutoffPowerLaw4FGLSpectralModel', 'secpl-4fgl']
LogParabolaSpectralModel : ['LogParabolaSpectralModel', 'lp']
TemplateSpectralModel : ['TemplateSpectralModel', 'template']
TemplateNDSpectralModel : ['TemplateNDSpectralModel', 'templateND']
GaussianSpectralModel : ['GaussianSpectralModel', 'gauss']
EBLAbsorptionNormSpectralModel : ['EBLAbsorptionNormSpectralModel', 'ebl-norm']
NaimaSpectralModel : ['NaimaSpectralModel', 'naima']
ScaleSpectralModel : ['ScaleSpectralModel', 'scale']
PowerLawNormSpectralModel : ['PowerLawNormSpectralModel', 'pl-norm']
LogParabolaNormSpectralModel : ['LogParabolaNormSpectralModel', 'lp-norm']
ExpCutoffPowerLawNormSpectralModel : ['ExpCutoffPowerLawNormSpectralModel', 'ecpl-norm']
DarkMatterAnnihilationSpectralModel : ['DarkMatterAnnihilationSpectralModel', 'dm-annihilation']
Spectral models all come with default parameters. Different parameter
values can be passed on creation of the model, either as a string
defining the value and unit or as an astropy.units.Quantity object
directly:
For convenience a str specifying the value and unit can be passed as
well:
pwl = PowerLawSpectralModel(amplitude="2.7e-12 TeV-1 cm-2 s-1", index=2.2)
print(pwl)
PowerLawSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral index 2.2000e+00 ... nan False False
spectral amplitude 2.7000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral reference 1.0000e+00 TeV ... nan True False
The model can be evaluated at given energies by calling the model instance:
[2.70000000e-12 2.40822469e-13 1.70358483e-14 1.51948705e-15] 1 / (cm2 s TeV)
The returned quantity is a differential photon flux.
For spectral models you can additionally compute the integrated and energy flux in a given energy range:
flux = pwl.integral(energy_min=1 * u.TeV, energy_max=10 * u.TeV)
print(flux)
eflux = pwl.energy_flux(energy_min=1 * u.TeV, energy_max=10 * u.TeV)
print(eflux)
2.108034597491956e-12 1 / (cm2 s)
4.982075849517389e-12 TeV / (cm2 s)
This also works for a list or an array of integration boundaries:
[1.64794383e-12 4.60090769e-13 1.03978226e-13] 1 / (cm2 s)
In some cases it can be useful to find use the inverse of a spectral model, to find the energy at which a given flux is reached:
1.0 TeV
As a convenience you can also plot any spectral model in a given energy range:
plt.figure()
pwl.plot(energy_bounds=[1, 100] * u.TeV)
<AxesSubplot:xlabel='Energy [TeV]', ylabel='dnde [1 / (cm2 s TeV)]'>
Norm Spectral Models#
Normed spectral models are a special class of Spectral Models, which
have a dimension-less normalisation. These spectral models feature a
norm parameter instead of amplitude and are named using the
NormSpectralModel suffix. They must be used along with another
spectral model, as a multiplicative correction factor according to their
spectral shape. They can be typically used for adjusting template based
models, or adding a EBL correction to some analytic model.
To check if a given SpectralModel is a norm model, you can simply
look at the is_norm_spectral_model property
# To see the available norm models shipped with gammapy:
for model in SPECTRAL_MODEL_REGISTRY:
if model.is_norm_spectral_model:
print(model)
<class 'gammapy.modeling.models.spectral.PiecewiseNormSpectralModel'>
<class 'gammapy.modeling.models.spectral.EBLAbsorptionNormSpectralModel'>
<class 'gammapy.modeling.models.spectral.PowerLawNormSpectralModel'>
<class 'gammapy.modeling.models.spectral.LogParabolaNormSpectralModel'>
<class 'gammapy.modeling.models.spectral.ExpCutoffPowerLawNormSpectralModel'>
As an example, we see the PowerLawNormSpectralModel
from gammapy.modeling.models import PowerLawNormSpectralModel
pwl_norm = PowerLawNormSpectralModel(tilt=0.1)
print(pwl_norm)
PowerLawNormSpectralModel
type name value unit error min max frozen is_norm link
-------- --------- ---------- ---- --------- --- --- ------ ------- ----
spectral norm 1.0000e+00 0.000e+00 nan nan False True
spectral tilt 1.0000e-01 0.000e+00 nan nan True False
spectral reference 1.0000e+00 TeV 0.000e+00 nan nan True False
We can check the correction introduced at each energy
[1.12794487 1. 0.89595846 0.79432823 0.7116851 ]
A typical use case of a norm model would be in applying spectral
correction to a TemplateSpectralModel. A template model is defined
by custom tabular values provided at initialization.
from gammapy.modeling.models import TemplateSpectralModel
plt.figure()
energy = [0.3, 1, 3, 10, 30] * u.TeV
values = [40, 30, 20, 10, 1] * u.Unit("TeV-1 s-1 cm-2")
template = TemplateSpectralModel(energy, values)
template.plot(energy_bounds=[0.2, 50] * u.TeV, label="template model")
normed_template = template * pwl_norm
normed_template.plot(energy_bounds=[0.2, 50] * u.TeV, label="normed_template model")
plt.legend()
<matplotlib.legend.Legend object at 0x14940d580>
Compound Spectral Model#
A CompoundSpectralModel is an arithmetic combination of two spectral
models. The model normed_template created in the preceding example
is an example of a CompoundSpectralModel
print(normed_template)
CompoundSpectralModel
Component 1 : TemplateSpectralModel
type name value unit error min max frozen is_norm link
-------- ---- ---------- ---- --------- --- --- ------ ------- ----
spectral norm 1.0000e+00 0.000e+00 nan nan True True
Component 2 : PowerLawNormSpectralModel
type name value unit error min max frozen is_norm link
-------- --------- ---------- ---- --------- --- --- ------ ------- ----
spectral norm 1.0000e+00 0.000e+00 nan nan False True
spectral tilt 1.0000e-01 0.000e+00 nan nan True False
spectral reference 1.0000e+00 TeV 0.000e+00 nan nan True False
Operator : mul
To create an additive model, you can do simply:
CompoundSpectralModel
Component 1 : PowerLawSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral index 2.2000e+00 ... nan False False
spectral amplitude 2.7000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral reference 1.0000e+00 TeV ... nan True False
Component 2 : TemplateSpectralModel
type name value unit error min max frozen is_norm link
-------- ---- ---------- ---- --------- --- --- ------ ------- ----
spectral norm 1.0000e+00 0.000e+00 nan nan True True
Operator : add
Spatial models#
Spatial models are imported from the same models
namespace, let’s start with a GaussianSpatialModel:
from gammapy.modeling.models import GaussianSpatialModel
gauss = GaussianSpatialModel(lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg")
print(gauss)
GaussianSpatialModel
type name value unit error min max frozen is_norm link
------- ----- ---------- ---- --------- ---------- --------- ------ ------- ----
spatial lon_0 0.0000e+00 deg 0.000e+00 nan nan False False
spatial lat_0 0.0000e+00 deg 0.000e+00 -9.000e+01 9.000e+01 False False
spatial sigma 2.0000e-01 deg 0.000e+00 0.000e+00 nan False False
spatial e 0.0000e+00 0.000e+00 0.000e+00 1.000e+00 True False
spatial phi 0.0000e+00 deg 0.000e+00 nan nan True False
Again you can check the SPATIAL_MODELS registry to see which models
are available or take a look at the spatial-models-gallery
from gammapy.modeling.models import SPATIAL_MODEL_REGISTRY
print(SPATIAL_MODEL_REGISTRY)
Registry
--------
ConstantSpatialModel : ['ConstantSpatialModel', 'const']
TemplateSpatialModel : ['TemplateSpatialModel', 'template']
DiskSpatialModel : ['DiskSpatialModel', 'disk']
GaussianSpatialModel : ['GaussianSpatialModel', 'gauss']
GeneralizedGaussianSpatialModel: ['GeneralizedGaussianSpatialModel', 'gauss-general']
PointSpatialModel : ['PointSpatialModel', 'point']
ShellSpatialModel : ['ShellSpatialModel', 'shell']
Shell2SpatialModel : ['Shell2SpatialModel', 'shell2']
The default coordinate frame for all spatial models is "icrs", but
the frame can be modified using the frame argument:
gauss = GaussianSpatialModel(
lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg", frame="galactic"
)
You can specify any valid astropy.coordinates frame. The center
position of the model can be retrieved as a
astropy.coordinates.SkyCoord object using SpatialModel.position:
print(gauss.position)
<SkyCoord (Galactic): (l, b) in deg
(0., 0.)>
Spatial models can be evaluated again by calling the instance:
lon = [0, 0.1] * u.deg
lat = [0, 0.1] * u.deg
flux_per_omega = gauss(lon, lat)
print(flux_per_omega)
[13061.88470839 10172.60603928] 1 / sr
The returned quantity corresponds to a surface brightness. Spatial model
can be also evaluated using Map and
Geom objects:
m = Map.create(skydir=(0, 0), width=(1, 1), binsz=0.02, frame="galactic")
m.quantity = gauss.evaluate_geom(m.geom)
plt.figure()
m.plot(add_cbar=True)
<WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>
Again for convenience the model can be plotted directly:
plt.figure()
gauss.plot(add_cbar=True)
<WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>
All spatial models have an associated sky region to it e.g. to
illustrate the extend of the model on a sky image. The returned object
is an SkyRegion object:
print(gauss.to_region())
Region: EllipseSkyRegion
center: <SkyCoord (Galactic): (l, b) in deg
(0., 0.)>
width: 0.6000000000000001 deg
height: 0.6000000000000001 deg
angle: 0.0 deg
Now we can plot the region on an sky image:
# create and plot the model
plt.figure()
gauss_elongated = GaussianSpatialModel(
lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg", e=0.7, phi="45 deg"
)
ax = gauss_elongated.plot(add_cbar=True)
# add region illustration
region = gauss_elongated.to_region()
region_pix = region.to_pixel(ax.wcs)
ax.add_artist(region_pix.as_artist(ec="w", fc="None"))
<matplotlib.patches.Ellipse object at 0x147bd19a0>
The to_region() method can also be useful to write e.g. ds9 region
files using write_ds9 from the regions package:
from regions import Regions
regions = Regions([gauss.to_region(), gauss_elongated.to_region()])
filename = "regions.reg"
regions.write(
filename,
format="ds9",
overwrite=True,
)
# !cat regions.reg
Temporal models#
Temporal models are imported from the same models
namespace, let’s start with a GaussianTemporalModel:
from gammapy.modeling.models import GaussianTemporalModel
gauss_temp = GaussianTemporalModel(t_ref=59240.0 * u.d, sigma=2.0 * u.d)
print(gauss_temp)
GaussianTemporalModel
type name value unit error min max frozen is_norm link
-------- ----- ---------- ---- --------- --- --- ------ ------- ----
temporal t_ref 5.9240e+04 d 0.000e+00 nan nan False False
temporal sigma 2.0000e+00 d 0.000e+00 nan nan False False
To check the TEMPORAL_MODELS registry to see which models are
available:
from gammapy.modeling.models import TEMPORAL_MODEL_REGISTRY
print(TEMPORAL_MODEL_REGISTRY)
Registry
--------
ConstantTemporalModel : ['ConstantTemporalModel', 'const']
LinearTemporalModel : ['LinearTemporalModel', 'linear']
LightCurveTemplateTemporalModel : ['LightCurveTemplateTemporalModel', 'template']
ExpDecayTemporalModel : ['ExpDecayTemporalModel', 'exp-decay']
GaussianTemporalModel : ['GaussianTemporalModel', 'gauss']
GeneralizedGaussianTemporalModel: ['GeneralizedGaussianTemporalModel', 'gengauss']
PowerLawTemporalModel : ['PowerLawTemporalModel', 'powerlaw']
SineTemporalModel : ['SineTemporalModel', 'sinus']
TemplatePhaseCurveTemporalModel : ['TemplatePhaseCurveTemporalModel', 'template-phase']
Temporal models can be evaluated on astropy.time.Time objects. The
returned quantity is a dimensionless number
from astropy.time import Time
time = Time("2021-01-29 00:00:00.000")
gauss_temp(time)
<Quantity 0.32465247>
As for other models, they can be plotted in a given time range
time = Time([59233.0, 59250], format="mjd")
gauss_temp.plot(time)
<AxesSubplot:xlabel='Time [iso]', ylabel='Norm / A.U.'>
SkyModel#
The SkyModel class combines a spectral, and
optionally, a spatial model and a temporal. It can be created from
existing spectral, spatial and temporal model components:
from gammapy.modeling.models import SkyModel
model = SkyModel(
spectral_model=pwl,
spatial_model=gauss,
temporal_model=gauss_temp,
name="my-source",
)
print(model)
SkyModel
Name : my-source
Datasets names : None
Spectral model type : PowerLawSpectralModel
Spatial model type : GaussianSpatialModel
Temporal model type : GaussianTemporalModel
Parameters:
index : 2.200 +/- 0.00
amplitude : 2.70e-12 +/- 0.0e+00 1 / (cm2 s TeV)
reference (frozen): 1.000 TeV
lon_0 : 0.000 +/- 0.00 deg
lat_0 : 0.000 +/- 0.00 deg
sigma : 0.200 +/- 0.00 deg
e (frozen): 0.000
phi (frozen): 0.000 deg
t_ref : 59240.000 +/- 0.00 d
sigma : 2.000 +/- 0.00 d
It is good practice to specify a name for your sky model, so that you can access it later by name and have meaningful identifier you serialisation. If you don’t define a name, a unique random name is generated:
model_without_name = SkyModel(spectral_model=pwl, spatial_model=gauss)
print(model_without_name.name)
nuXM_oVi
The individual components of the source model can be accessed using
.spectral_model, .spatial_model and .temporal_model:
print(model.spectral_model)
PowerLawSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral index 2.2000e+00 ... nan False False
spectral amplitude 2.7000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral reference 1.0000e+00 TeV ... nan True False
print(model.spatial_model)
GaussianSpatialModel
type name value unit error min max frozen is_norm link
------- ----- ---------- ---- --------- ---------- --------- ------ ------- ----
spatial lon_0 0.0000e+00 deg 0.000e+00 nan nan False False
spatial lat_0 0.0000e+00 deg 0.000e+00 -9.000e+01 9.000e+01 False False
spatial sigma 2.0000e-01 deg 0.000e+00 0.000e+00 nan False False
spatial e 0.0000e+00 0.000e+00 0.000e+00 1.000e+00 True False
spatial phi 0.0000e+00 deg 0.000e+00 nan nan True False
print(model.temporal_model)
GaussianTemporalModel
type name value unit error min max frozen is_norm link
-------- ----- ---------- ---- --------- --- --- ------ ------- ----
temporal t_ref 5.9240e+04 d 0.000e+00 nan nan False False
temporal sigma 2.0000e+00 d 0.000e+00 nan nan False False
And can be used as you have seen already seen above:
plt.figure()
model.spectral_model.plot(energy_bounds=[1, 10] * u.TeV)
<AxesSubplot:xlabel='Energy [TeV]', ylabel='dnde [1 / (cm2 s TeV)]'>
Note that the gammapy fitting can interface only with a SkyModel and
not its individual components. So, it is customary to work with
SkyModel even if you are not doing a 3D fit. Since the amplitude
parameter resides on the SpectralModel, specifying a spectral
component is compulsory. The temporal and spatial components are
optional. The temporal model needs to be specified only for timing
analysis. In some cases (e.g. when doing a spectral analysis) there is
no need for a spatial component either, and only a spectral model is
associated with the source.
model_spectrum = SkyModel(spectral_model=pwl, name="source-spectrum")
print(model_spectrum)
SkyModel
Name : source-spectrum
Datasets names : None
Spectral model type : PowerLawSpectralModel
Spatial model type :
Temporal model type :
Parameters:
index : 2.200 +/- 0.00
amplitude : 2.70e-12 +/- 0.0e+00 1 / (cm2 s TeV)
reference (frozen): 1.000 TeV
Additionally the spatial model of SkyModel
can be used to represent source models based on templates, where the
spatial and energy axes are correlated. It can be created e.g. from an
existing FITS file:
from gammapy.modeling.models import PowerLawNormSpectralModel, TemplateSpatialModel
diffuse_cube = TemplateSpatialModel.read(
"$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz", normalize=False
)
diffuse = SkyModel(PowerLawNormSpectralModel(), diffuse_cube)
print(diffuse)
SkyModel
Name : VLGilO4j
Datasets names : None
Spectral model type : PowerLawNormSpectralModel
Spatial model type : TemplateSpatialModel
Temporal model type :
Parameters:
norm : 1.000 +/- 0.00
tilt (frozen): 0.000
reference (frozen): 1.000 TeV
Note that if the spatial model is not normalized over the sky it has to
be combined with a normalized spectral model, for example
PowerLawNormSpectralModel. This is the only
case in gammapy.models.SkyModel where the unit is fully attached to
the spatial model.
Modifying model parameters#
Model parameters can be modified (eg: frozen, values changed, etc at any point), eg:
# Freezing a parameter
model.spectral_model.index.frozen = True
# Making a parameter free
model.spectral_model.index.frozen = False
# Changing a value
model.spectral_model.index.value = 3
# Setting min and max ranges on parameters
model.spectral_model.index.min = 1.0
model.spectral_model.index.max = 5.0
# Visualise the model as a table
display(model.parameters.to_table())
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --------- ------ ------- ----
spectral index 3.0000e+00 ... 5.000e+00 False False
spectral amplitude 2.7000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral reference 1.0000e+00 TeV ... nan True False
spatial lon_0 0.0000e+00 deg ... nan False False
spatial lat_0 0.0000e+00 deg ... 9.000e+01 False False
spatial sigma 2.0000e-01 deg ... nan False False
spatial e 0.0000e+00 ... 1.000e+00 True False
spatial phi 0.0000e+00 deg ... nan True False
temporal t_ref 5.9240e+04 d ... nan False False
temporal sigma 2.0000e+00 d ... nan False False
You can use the interactive boxes to choose model parameters by name, type or other attrributes mentioned in the column names.
Model lists and serialisation#
In a typical analysis scenario a model consists of multiple model
components, or a “catalog” or “source library”. To handle this list of
multiple model components, Gammapy has a Models class:
Models
Component 0: SkyModel
Name : my-source
Datasets names : None
Spectral model type : PowerLawSpectralModel
Spatial model type : GaussianSpatialModel
Temporal model type : GaussianTemporalModel
Parameters:
index : 3.000 +/- 0.00
amplitude : 2.70e-12 +/- 0.0e+00 1 / (cm2 s TeV)
reference (frozen): 1.000 TeV
lon_0 : 0.000 +/- 0.00 deg
lat_0 : 0.000 +/- 0.00 deg
sigma : 0.200 +/- 0.00 deg
e (frozen): 0.000
phi (frozen): 0.000 deg
t_ref : 59240.000 +/- 0.00 d
sigma : 2.000 +/- 0.00 d
Component 1: SkyModel
Name : VLGilO4j
Datasets names : None
Spectral model type : PowerLawNormSpectralModel
Spatial model type : TemplateSpatialModel
Temporal model type :
Parameters:
norm : 1.000 +/- 0.00
tilt (frozen): 0.000
reference (frozen): 1.000 TeV
Individual model components in the list can be accessed by their name:
print(models["my-source"])
SkyModel
Name : my-source
Datasets names : None
Spectral model type : PowerLawSpectralModel
Spatial model type : GaussianSpatialModel
Temporal model type : GaussianTemporalModel
Parameters:
index : 3.000 +/- 0.00
amplitude : 2.70e-12 +/- 0.0e+00 1 / (cm2 s TeV)
reference (frozen): 1.000 TeV
lon_0 : 0.000 +/- 0.00 deg
lat_0 : 0.000 +/- 0.00 deg
sigma : 0.200 +/- 0.00 deg
e (frozen): 0.000
phi (frozen): 0.000 deg
t_ref : 59240.000 +/- 0.00 d
sigma : 2.000 +/- 0.00 d
Note: To make the access by name unambiguous, models are required to have a unique name, otherwise an error will be thrown.
To see which models are available you can use the .names attribute:
print(models.names)
['my-source', 'VLGilO4j']
Note that a SkyModel object can be evaluated for a given longitude,
latitude, and energy, but the Models object cannot. This Models
container object will be assigned to Dataset or Datasets
together with the data to be fitted. Checkout e.g. the
Modelling tutorial for details.
The Models class also has in place .append() and .extend()
methods:
model_copy = model.copy(name="my-source-copy")
models.append(model_copy)
This list of models can be also serialised to a custom YAML based format:
models_yaml = models.to_yaml()
print(models_yaml)
components:
- name: my-source
type: SkyModel
spectral:
type: PowerLawSpectralModel
parameters:
- name: index
value: 3.0
min: 1.0
max: 5.0
- name: amplitude
value: 2.7e-12
unit: cm-2 s-1 TeV-1
- name: reference
value: 1.0
unit: TeV
spatial:
type: GaussianSpatialModel
frame: galactic
parameters:
- name: lon_0
value: 0.0
unit: deg
- name: lat_0
value: 0.0
unit: deg
- name: sigma
value: 0.2
unit: deg
- name: e
value: 0.0
- name: phi
value: 0.0
unit: deg
temporal:
type: GaussianTemporalModel
parameters:
- name: t_ref
value: 59240.0
unit: d
- name: sigma
value: 2.0
unit: d
- name: VLGilO4j
type: SkyModel
spectral:
type: PowerLawNormSpectralModel
parameters:
- name: norm
value: 1.0
- name: tilt
value: 0.0
- name: reference
value: 1.0
unit: TeV
spatial:
type: TemplateSpatialModel
frame: galactic
parameters: []
filename: /Users/terrier/Code/gammapy-dev/gammapy-data/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz
normalize: false
unit: 1 / (cm2 MeV s sr)
- name: my-source-copy
type: SkyModel
spectral:
type: PowerLawSpectralModel
parameters:
- name: index
value: 3.0
min: 1.0
max: 5.0
- name: amplitude
value: 2.7e-12
unit: cm-2 s-1 TeV-1
- name: reference
value: 1.0
unit: TeV
spatial:
type: GaussianSpatialModel
frame: galactic
parameters:
- name: lon_0
value: 0.0
unit: deg
- name: lat_0
value: 0.0
unit: deg
- name: sigma
value: 0.2
unit: deg
- name: e
value: 0.0
- name: phi
value: 0.0
unit: deg
temporal:
type: GaussianTemporalModel
parameters:
- name: t_ref
value: 59240.0
unit: d
- name: sigma
value: 2.0
unit: d
The structure of the yaml files follows the structure of the python
objects. The components listed correspond to the SkyModel and
components of the Models. For each SkyModel
we have information about its name, type (corresponding to the
tag attribute) and sub-mobels (i.e spectral model and eventually
spatial model). Then the spatial and spectral models are defined by
their type and parameters. The parameters keys name/value/unit are
mandatory, while the keys min/max/frozen are optionnals (so you can
prepare shorter files).
If you want to write this list of models to disk and read it back later you can use:
models.write("models.yaml", overwrite=True)
models_read = Models.read("models.yaml")
Additionally the models can exported and imported togeter with the data
using the Datasets.read() and Datasets.write() methods as shown
in the Multi instrument joint 3D and 1D analysis
notebook.
Implementing a custom model#
In order to add a user defined spectral model you have to create a SpectralModel subclass. This new model class should include:
a tag used for serialization (it can be the same as the class name)
an instantiation of each Parameter with their unit, default values and frozen status
the evaluate function where the mathematical expression for the model is defined.
As an example we will use a PowerLawSpectralModel plus a Gaussian (with
fixed width). First we define the new custom model class that we name
MyCustomSpectralModel:
from gammapy.modeling import Parameter
from gammapy.modeling.models import SpectralModel
class MyCustomSpectralModel(SpectralModel):
"""My custom spectral model, parametrising a power law plus a Gaussian spectral line.
Parameters
----------
amplitude : `astropy.units.Quantity`
Amplitude of the spectra model.
index : `astropy.units.Quantity`
Spectral index of the model.
reference : `astropy.units.Quantity`
Reference energy of the power law.
mean : `astropy.units.Quantity`
Mean value of the Gaussian.
width : `astropy.units.Quantity`
Sigma width of the Gaussian line.
"""
tag = "MyCustomSpectralModel"
amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1", min=0, is_norm=True)
index = Parameter("index", 2, min=0)
reference = Parameter("reference", "1 TeV", frozen=True)
mean = Parameter("mean", "1 TeV", min=0)
width = Parameter("width", "0.1 TeV", min=0, frozen=True)
@staticmethod
def evaluate(energy, index, amplitude, reference, mean, width):
pwl = PowerLawSpectralModel.evaluate(
energy=energy,
index=index,
amplitude=amplitude,
reference=reference,
)
gauss = amplitude * np.exp(-((energy - mean) ** 2) / (2 * width**2))
return pwl + gauss
It is good practice to also implement a docstring for the model,
defining the parameters and also definig a .tag, which specifies the
name of the model for serialisation. Also note that gammapy assumes that
all SpectralModel evaluate functions return a flux in unit of
"cm-2 s-1 TeV-1" (or equivalent dimensions).
This model can now be used as any other spectral model in Gammapy:
my_custom_model = MyCustomSpectralModel(mean="3 TeV")
print(my_custom_model)
print(my_custom_model.integral(1 * u.TeV, 10 * u.TeV))
plt.figure()
my_custom_model.plot(energy_bounds=[1, 10] * u.TeV)
MyCustomSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral amplitude 1.0000e-12 cm-2 s-1 TeV-1 ... nan False True
spectral index 2.0000e+00 ... nan False False
spectral reference 1.0000e+00 TeV ... nan True False
spectral mean 3.0000e+00 TeV ... nan False False
spectral width 1.0000e-01 TeV ... nan True False
1.1442739329466746e-12 1 / (cm2 s)
<AxesSubplot:xlabel='Energy [TeV]', ylabel='dnde [1 / (cm2 s TeV)]'>
As a next step we can also register the custom model in the
SPECTRAL_MODELS registry, so that it becomes available for
serilisation:
SPECTRAL_MODEL_REGISTRY.append(MyCustomSpectralModel)
model = SkyModel(spectral_model=my_custom_model, name="my-source")
models = Models([model])
models.write("my-custom-models.yaml", overwrite=True)
# !cat my-custom-models.yaml
Similarly you can also create custom spatial models and add them to the
SPATIAL_MODELS registry. In that case gammapy assumes that the
evaluate function return a normalized quantity in “sr-1” such as the
model integral over the whole sky is one.
Models with energy dependent morphology#
A common science case in the study of extended sources is to probe for energy dependent morphology, eg: in Supernova Remnants or Pulsar Wind Nebulae. Traditionally, this has been done by splitting the data into energy bands and doing individual fits of the morphology in these energy bands.
SkyModel offers a natural framework to simultaneously model the
energy and morphology, e.g. spatial extent described by a parametric
model expression with energy dependent parameters.
The models shipped within gammapy use a “factorised” representation of the source model, where the spatial (\(l,b\)), energy (\(E\)) and time (\(t\)) dependence are independent model components and not correlated:
:raw-latex:`\begin{align}f(l, b, E, t) = F(l, b) \cdot G(E) \cdot H(t)\end{align}`
To use full 3D models, ie $f(l, b, E) = F(l, b, E)
:raw-latex:`\cdot`G(E) $, you have to implement your own custom
SpatialModel. Note that it is still necessary to multiply by a
SpectralModel, \(G(E)\) to be dimensionally consistent.
In this example, we create Gaussian Spatial Model with the extension varying with energy. For simplicity, we assume a linear dependence on energy and parameterize this by specifying the extension at 2 energies. You can add more complex dependences, probably motivated by physical models.
from astropy.coordinates.angle_utilities import angular_separation
from gammapy.modeling.models import SpatialModel
class MyCustomGaussianModel(SpatialModel):
"""My custom Energy Dependent Gaussian model.
Parameters
----------
lon_0, lat_0 : `~astropy.coordinates.Angle`
Center position
sigma_1TeV : `~astropy.coordinates.Angle`
Width of the Gaussian at 1 TeV
sigma_10TeV : `~astropy.coordinates.Angle`
Width of the Gaussian at 10 TeV
"""
tag = "MyCustomGaussianModel"
is_energy_dependent = True
lon_0 = Parameter("lon_0", "0 deg")
lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
sigma_1TeV = Parameter("sigma_1TeV", "2.0 deg", min=0)
sigma_10TeV = Parameter("sigma_10TeV", "0.2 deg", min=0)
@staticmethod
def evaluate(lon, lat, energy, lon_0, lat_0, sigma_1TeV, sigma_10TeV):
sep = angular_separation(lon, lat, lon_0, lat_0)
# Compute sigma for the given energy using linear interpolation in log energy
sigma_nodes = u.Quantity([sigma_1TeV, sigma_10TeV])
energy_nodes = [1, 10] * u.TeV
log_s = np.log(sigma_nodes.to("deg").value)
log_en = np.log(energy_nodes.to("TeV").value)
log_e = np.log(energy.to("TeV").value)
sigma = np.exp(np.interp(log_e, log_en, log_s)) * u.deg
exponent = -0.5 * (sep / sigma) ** 2
norm = 1 / (2 * np.pi * sigma**2)
return norm * np.exp(exponent)
Serialisation of this model can be achieved as explained in the previous
section. You can now use it as standard SpatialModel in your
analysis. Note that this is still a SpatialModel and not a
SkyModel, so it needs to be multiplied by a SpectralModel as
before.
spatial_model = MyCustomGaussianModel()
spectral_model = PowerLawSpectralModel()
sky_model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
print(spatial_model.evaluation_radius)
None
To visualise it, we evaluate it on a 3D geom.
energy_axis = MapAxis.from_energy_bounds(
energy_min=0.1 * u.TeV, energy_max=10.0 * u.TeV, nbin=3, name="energy_true"
)
geom = WcsGeom.create(skydir=(0, 0), width=5.0 * u.deg, binsz=0.1, axes=[energy_axis])
spatial_model.plot_grid(geom=geom, add_cbar=True, figsize=(14, 3))
plt.show()
For computational purposes, it is useful to specify a
evaluation_radius for SpatialModels - this gives a size on which
to compute the model. Though optional, it is highly recommended for
Custom Spatial Models. This can be done, for ex, by defining the
following function inside the above class:
@property
def evaluation_radius(self):
"""Evaluation radius (`~astropy.coordinates.Angle`)."""
return 5 * np.max([self.sigma_1TeV.value, self.sigma_10TeV.value]) * u.deg
