PIG 21 - Models improvements

  • Author: Axel Donath, Régis Terrier and Quentin Remy

  • Created: Jun 10, 2020

  • Accepted: July 31, 2020

  • Status: accepted

  • Discussion: GH 2944

Abstract

This PIG outlines further minor improvements to the modeling framework of Gammapy. This includes introduction of a spectra models with norm parameters, minimal support for energy dependent spatial models as well as simplifications of the YAML serialisation.

Proposal

Spectral Norm Models

For the purpose of adjusting template based models we propose to introduce a new class of spectral models. Those spectral models feature a norm parameter instead of amplitude and are named using the NormSpectralModel suffix:

from gammapy.modeling.models import (
    PowerLawNormSpectralModel,
    LogParabolaNormSpectralModel,
    PiecewiseBrokenPowerlawNormSpectralModel,
    ConstantNormSpectralModel,
)

pwl_norm = PowerLawNormSpectralModel(norm=, tilt=, reference=)
log_parabola_norm = LogParabolaNormSpectralModel(norm=, alpha=, beta=)
bpwl_norm = PiecewiseBrokenPowerlawNormSpectralModel(norms=, energies=)
const_norm = ConstantNormSpectralModel(norm=)

# typically used like
    template = TemplateSpectralModel()
model = template * pwl_norm

model = BackgroundModel(
            map=map,
            spectral_model=pwl_norm
    )

model = SkyDiffuseCube(
            map=map,
            spectral_model=pwl_norm
    )

This requires removing the hard-coded parameters of the TemplateSpectralModel and BackgroundModel.

Energy Dependent Spatial Models

A very common use-case in scientific analyses is to look for energy dependent morphology of extended sources. In the past this kind of analysis has been typically done by splitting the data into energy bands and doing individual fits of the morphology in these energy bands. In a combined spectral and spatial (“cube”) analysis this can be naturally achieved by allowing for an energy dependent spatial model, where the energy dependence is e.g. described by a parametric model expression with energy dependent parameters.

In the current model scheme we use a “factorised” representation of the source model, where the spatial, energy and time dependence are independent model components and not correlated:

\[f(l, b, E, t) = A \cdot F(l, b) \cdot G(E) \cdot H(t)\]

To represent energy dependent morphology we propose to introduce energy dependent spatial models:

\[f(l, b, E, t) = A \cdot F(l, b, E) \cdot G(E) \cdot H(t)\]

In general the energy dependence is optional. If the spatial model does not declare an energy dependence it assumes the same morphology for all energies. This also ensures backwards compatibility with the current behaviour.

To limit the implementation effort in this PIG we propose to only add an example energy dependent custom model to our documentation. We do not propose to introduce general dependence of arbitrary model parameters for any spatial model, such as GaussianSpatialModel or DiskSpatialModel. An example of how this can be achieved with a custom implemented model is given below:

from gammapy.modeling.models import SpatialModel
from astropy.coordinates.angle_utilities import angular_separation

class MyCustomGaussianModel(SpatialModel):
    """My custom 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"
    lon_0 = Parameter("lon_0", "0 deg")
    lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)

    sigma_1TeV = Parameter("sigma_1TeV", "1 deg", min=0)
    sigma_10TeV = Parameter("sigma_10TeV", "0.5 deg", min=0)

@staticmethod:
def evaluate(lon, lat, energy, lon_0, lat_0, sigma_1TeV, sigma_10TeV):
    """Evaluate custom Gaussian model"""
    sigmas = u.Quantity([sigma_1TeV, sigma_10TeV])
    energy_nodes = [1, 10] * u.TeV
    sigma = np.interp(energy, energy_nodes, sigmas)

    sep = angular_separation(lon, lat, lon_0, lat_0)

    exponent = -0.5 * (sep / sigma) ** 2
    norm = 1 / (2 * np.pi * sigma ** 2)
    return norm * np.exp(exponent)

@property
def evaluation_radius(self):
    """Evaluation radius (`~astropy.coordinates.Angle`)."""
    return 5 * self.sigma_1TeV.quantity

Spectral Absorption Model

In the current handling of absorbed spectral models we have a very special Absorption model, which is not a spectral model. To resolve this special case, we propose to refactor the existing code and handle the absorbed case using a CompoundSpectralModel. The new implementation is used as follows:

from gammapy.modeling.models import EBLAbsorptionNormSpectralModel, PowerLawSpectralModel

absorption = EBLAbsorptionNormSpectralModel.from_reference(
    redshift=0.1, alpha_norm=1, reference="dominguez"
)

pwl = PowerLawSpectralModel()

spectral_model = absorption * pwl

assert isinstance(spectral_model, CompoundSpectralModel)

In addition we propose to rename .table_model to .to_template_spectral_model(redshift, alpha_norm).

Additional Models

In addition we propose to implement the following models in gammapy.modeling.models:

  • SersicSpatialModel following the parametrisation of the Astropy Sersic2D model.

  • BrokenPowerLaw with the following parametrisation:

\[\begin{split}\begin{split}\phi(E) = \phi_0 \times \left \{ \begin{eqnarray} \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\ \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise} \end{eqnarray} \right .\end{split}\end{split}\]

  • We propose to re-introduce the PhaseCurveModel to compute mean fluxes over time

Simplify YAML Representation

Introduce Shorter YAML Alias Tags

To simplify the definition of models in YAML files as well as for interactive model creation we propose to introduce shorter YAML tags for all models in Gammapy. For backwards compatibility we propose to support the class name as well. We require that the short YAML tag is only unique within the model class, so that the same tag such as “gauss” can be re-used by spectral, spatial as well as temporal components. The uniqueness is guaranteed in the YAML file, because of the different sections for the model types. For writing the models the verbose class name is used by default.

We propose to introduce the following YAML tags:

Class Name

YAML Tag

ConstantSpectralModel

const

ConstantNormSpectralModel

const-norm

CompoundSpectralModel

compound

PowerLawSpectralModel

pl

PowerLawNormSpectralModel

pl-norm

PowerLaw2SpectralModel

pl-2

SmoothBrokenPowerLawSpectralModel

sbpl

BrokenPowerLawSpectralModel

bpl

PiecewiseBrokenPowerLawNormSpectraModel

pwbpl-norm

ExpCutoffPowerLawSpectralModel

ecpl

ExpCutoffPowerLaw3FGLSpectralModel

ecpl-3fgl

SuperExpCutoffPowerLaw3FGLSpectralModel

secpl-3fgl

SuperExpCutoffPowerLaw4FGLSpectralModel

secpl-4fgl

LogParabolaSpectralModel

lp

LogParabolaNormSpectralModel

lp-norm

TemplateSpectralModel

template

GaussianSpectralModel

gauss

EBLAbsorbtionNormSpectralModel

ebl-absorbtion-norm

NaimaSpectralModel

naima

ScaleSpectralModel

scale

Class Name

YAML Tag

ConstantSpatialModel

const

TemplateSpatialModel

template

GaussianSpatialModel

gauss

DiskSpatialModel

disk

PointSpatialModel

point

ShellSpatialModel

shell

SersicSpatialModel

sersic

Class Name

YAML Tag

ConstantTemporalModel

const

LightCurveTemplateTemporalModel

template

GaussianTemporalModel

gauss

ExpDecayTemporalModel

exp-decay

To simplify the interactive model creation we propose to introduce:

from gammapy.modeling.models import SkyModel

model = SkyModel.create(
    spectral_type="pl",
            spectral_pars={"index": 2},
            spatial_type="gauss",
            spatial_pars={"sigma": "0.1 deg"},
            temporal_type="const",
)

In addition the Model.create() factory function should be adapted to:

from gammapy.modeling.models import Model

pwl = Model.create(tag="gauss", model_type="spectral", **pars)

Simplify YAML Parameter Representation

To further simplify the structure of the YAML file we propose to remove parameter properties if they are equivalent to the default values. The current representation looks like:

spectral:
    type: PowerLawSpectralModel
    parameters:
    - name: index
      value: 2.0
      unit: ''
      min: .nan
      max: .nan
      frozen: false
      error: 0
    - name: amplitude
      value: 1.0e-12
      unit: cm-2 s-1 TeV-1
      min: .nan
      max: .nan
      frozen: false
      error: 0
    - name: reference
      value: 1.0
      unit: TeV
      min: .nan
      max: .nan
      frozen: true
      error: 0

After the simplification:

spectral:
    type: PowerLawSpectralModel
    parameters:
    - name: index
      value: 2.0
    - name: amplitude
      value: 1.0e-12
      unit: cm-2 s-1 TeV-1
    - name: reference
      value: 1.0
      unit: TeV
      frozen: true

Decision

This PIG received feedback in form of comments to the PR GH 2944 and via private Slack messages to the authors. Most concerns were raised against the introduction of NormSpectralModels. In a following discussion among the lead developers it was decided to stick with the proposal, but keep the abstraction of SkyDiffuseCube to avoid the case, where a spatial model carries the flux information. A final review announced on the Gammapy and CC mailing list did not yield any additional comments. Therefore the PIG was accepted on July 31, 2020.