<jupyter_start><jupyter_text>Gammapy ModelsThis is an introduction and overview on how to work with models in Gammapy. The sub-package `~gammapy.modeling` contains all the functionality related to modeling and fittingdata. This includes spectral, spatial and temporal model classes, as well as the fitand parameter API. We will cover the follwing topics in order:1. [Spectral Models](Spectral-Models)1. [Spatial Models](Spatial-Models)1. [SkyModel and SkyDiffuseCube](SkyModel-and-SkyDiffuseCube)1. [Model Lists and Serialisation](Model-Lists-and-Serialisation)1. [Implementing as Custom Model](Implementing-a-Custom-Model)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](https://docs.gammapy.org/dev/modeling/gallery/index.htmlspectral-models).Note that there is a separate tutorial [modeling](modeling.ipynb) that explains about `~gammapy.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<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from gammapy.maps import Map<jupyter_output><empty_output><jupyter_text>Spectral ModelsAll models are imported from the `~gammapy.modeling.models` namespace. Let's start with a `PowerLawSpectralModel`:<jupyter_code>from gammapy.modeling.models import PowerLawSpectralModel
pwl = PowerLawSpectralModel()
print(pwl)<jupyter_output><empty_output><jupyter_text>To get a list of all available spectral models you can import and print the spectral model registry or take a look at the [model gallery](https://docs.gammapy.org/dev/modeling/gallery/index.htmlspectral-models):<jupyter_code>from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY

print(SPECTRAL_MODEL_REGISTRY)<jupyter_output><empty_output><jupyter_text>Spectral models all come with default parameters. Different parametervalues can be passed on creation of the model, either as a string definingthe value and unit or as an `astropy.units.Quantity` object directly:<jupyter_code>amplitude = 1e-12 * u.Unit("TeV-1 cm-2 s-1")
pwl = PowerLawSpectralModel(amplitude=amplitude, index=2.2)<jupyter_output><empty_output><jupyter_text>For convenience a `str` specifying the value and unit can be passed as well:<jupyter_code>pwl = PowerLawSpectralModel(amplitude="2.7e-12 TeV-1 cm-2 s-1", index=2.2)
print(pwl)<jupyter_output><empty_output><jupyter_text>The model can be evaluated at given energies by calling the model instance:<jupyter_code>energy = [1, 3, 10, 30] * u.TeV
dnde = pwl(energy)
print(dnde)<jupyter_output><empty_output><jupyter_text>The returned quantity is a differential photon flux. For spectral models you can computed in addition the integrated and energy fluxin a given energy range:<jupyter_code>flux = pwl.integral(emin=1 * u.TeV, emax=10 * u.TeV)
print(flux)

eflux = pwl.energy_flux(emin=1 * u.TeV, emax=10 * u.TeV)
print(eflux)<jupyter_output><empty_output><jupyter_text>This also works for a list or an array of integration boundaries:<jupyter_code>energy = [1, 3, 10, 30] * u.TeV
flux = pwl.integral(emin=energy[:-1], emax=energy[1:])
print(flux)<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>dnde = 2.7e-12 * u.Unit("TeV-1 cm-2 s-1")
energy = pwl.inverse(dnde)
print(energy)<jupyter_output><empty_output><jupyter_text>As a convenience you can also plot any spectral model in a given energy range:<jupyter_code>pwl.plot(energy_range=[1, 100] * u.TeV)<jupyter_output><empty_output><jupyter_text>Spatial Models Spatial models are imported from the same `~gammapy.modeling.models` namespace, let's start with a `GaussianSpatialModel`:<jupyter_code>from gammapy.modeling.models import GaussianSpatialModel
gauss = GaussianSpatialModel(lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg")
print(gauss)<jupyter_output><empty_output><jupyter_text>Again you can check the `SPATIAL_MODELS` registry to see which models are available or take a look at the [model gallery](https://docs.gammapy.org/dev/modeling/gallery/index.htmlspatial-models).<jupyter_code>from gammapy.modeling.models import SPATIAL_MODEL_REGISTRY

print(SPATIAL_MODEL_REGISTRY)<jupyter_output><empty_output><jupyter_text>The default coordinate frame for all spatial models is ``"icrs"``, but the frame can be modified using the``frame`` argument:<jupyter_code>gauss = GaussianSpatialModel(
    lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg", frame="galactic"
)<jupyter_output><empty_output><jupyter_text>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`:<jupyter_code>print(gauss.position)<jupyter_output><empty_output><jupyter_text>Spatial models can be evaluated again by calling the instance:<jupyter_code>lon = [0, 0.1] * u.deg
lat = [0, 0.1] * u.deg

flux_per_omega = gauss(lon, lat)
print(flux_per_omega)<jupyter_output><empty_output><jupyter_text>The returned quantity corresponds to a surface brightness. Spatial modelcan be also evaluated using `~gammapy.maps.Map` and `~gammapy.maps.Geom` objects:<jupyter_code>m = Map.create(skydir=(0, 0), width=(1, 1), binsz=0.02, frame="galactic")
m.quantity = gauss.evaluate_geom(m.geom)
m.plot(add_cbar=True);<jupyter_output><empty_output><jupyter_text>Again for convenience the model can be plotted directly:<jupyter_code>gauss.plot(add_cbar=True);<jupyter_output><empty_output><jupyter_text>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 `regions.SkyRegion` object:<jupyter_code>print(gauss.to_region())<jupyter_output><empty_output><jupyter_text>Now we can plot the region on an sky image:<jupyter_code># create and plot the model
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());<jupyter_output><empty_output><jupyter_text>The `.to_region()` method can also be useful to write e.g. ds9 region files using `write_ds9` from the `regions` package:<jupyter_code>from regions import write_ds9

regions = [gauss.to_region(), gauss_elongated.to_region()]

filename = "regions.reg"
write_ds9(regions, filename, coordsys="galactic", fmt=".4f", radunit="deg")
!cat regions.reg<jupyter_output><empty_output><jupyter_text>SkyModel and SkyDiffuseCube The `~gammapy.modeling.models.SkyModel` class combines a spectral and a spatial model. It can be createdfrom existing spatial and spectral model components:<jupyter_code>from gammapy.modeling.models import SkyModel

model = SkyModel(spectral_model=pwl, spatial_model=gauss, name="my-source")
print(model)<jupyter_output><empty_output><jupyter_text>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 serilisation. If you don't define a name, a unique random name is generated:<jupyter_code>model_without_name = SkyModel(spectral_model=pwl, spatial_model=gauss)
print(model_without_name.name)<jupyter_output><empty_output><jupyter_text>The spectral and spatial component of the source model can be accessed using `.spectral_model` and `.spatial_model`:<jupyter_code>model.spectral_model
model.spatial_model<jupyter_output><empty_output><jupyter_text>And can be used as you have seen already seen above:<jupyter_code>model.spectral_model.plot(energy_range=[1, 10] * u.TeV);<jupyter_output><empty_output><jupyter_text>In some cases (e.g. when doing a spectral analysis) there is only a spectral model associated with the source. So the spatial model is optional:<jupyter_code>model_spectrum = SkyModel(spectral_model=pwl, name="source-spectrum")
print(model_spectrum)<jupyter_output><empty_output><jupyter_text>Additionally the `~gammapy.modeling.models.SkyDiffuseCube` 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:<jupyter_code>from gammapy.modeling.models import SkyDiffuseCube
diffuse = SkyDiffuseCube.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz"
)
print(diffuse)<jupyter_output><empty_output><jupyter_text>Model Lists and SerialisationIn a typical analysis scenario a model consists of mutiple model components, or a "catalog" or "source library". To handle this list of multiple model components, Gammapy has a `Models` class:<jupyter_code>from gammapy.modeling.models import Models
models = Models([model, diffuse])
print(models)<jupyter_output><empty_output><jupyter_text>Individual model components in the list can be accessed by their name:<jupyter_code>print(models["my-source"])<jupyter_output><empty_output><jupyter_text>**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:<jupyter_code>print(models.names)<jupyter_output><empty_output><jupyter_text>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 as explained in other analysis tutorials (see for example the [modeling](modeling.ipynb) notebook).The `Models` class also has in place `.append()` and `.extend()` methods:<jupyter_code>model_copy = model.copy(name="my-source-copy")
models.append(model_copy)<jupyter_output><empty_output><jupyter_text>This list of models can be also serialised toa custom YAML based format:<jupyter_code>models_yaml = models.to_yaml()
print(models_yaml)<jupyter_output><empty_output><jupyter_text>The structure of the yaml files follows the structure of the python objects.The `components` listed correspond to the `SkyModel` and `SkyDiffuseCube` components of the `Models`. For each `SkyModel` we have  informations 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 defiend 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:<jupyter_code>models.write("models.yaml", overwrite=True)
models_read = Models.read("models.yaml")<jupyter_output><empty_output><jupyter_text>Additionally the models can exported and imported togeter with the data using the `Datasets.read()` and `Datasets.write()` methods as shown in the [analysis_mwl](analysis_mwl.ipynb) notebook.  Implementing a Custom ModelIn 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`:<jupyter_code>from gammapy.modeling.models import SpectralModel, Parameter


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)
    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<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>my_custom_model = MyCustomSpectralModel(mean="3 TeV")
print(my_custom_model)
my_custom_model.integral(1 * u.TeV, 10 * u.TeV)
my_custom_model.plot(energy_range=[1, 10] * u.TeV)<jupyter_output><empty_output><jupyter_text>As a next step we can also register the custom model in the `SPECTRAL_MODELS` registry, so that it becomes available for serilisation:<jupyter_code>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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09434770196974022, 'fraction_numerical': 0.009063659720239795, 'mean_word_length': 6.171443193449335, 'pattern_counts': {'":': 0, '<': 141, '<?xml version=': 0, '>': 141, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '35', 'dir': 'jupyter-structured-clean-dedup', 'id': '353008', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4590, 'n_tokens_neox': 4390, 'n_words': 1367, 'path': 'docs/tutorials/models.ipynb', 'repo_name': 'vikasj78/gammapy', 'size': '14012'}