.. include:: ../../references.txt .. _pig-007: ************** PIG 7 - Models ************** * Author: Axel Donath, Christoph Deil & Regis Terrier * Created: Dec 17, 2018 * Withdrawn: Oct 7, 2019 * Status: withdrawn * Discussion: `GH 1971`_ Introduction ============ One of the most important features of Gammapy is the modeling of gamma-ray data by parametric models. This includes the modeling of source spectra by spectral models, the modeling of source morphologies by spatial models, the modeling of light curves or phase curves with temporal models and combinations thereof. Instrument specific characteristics, such as PSF, energy resolution, effective area and background are also considered as part of the model framework. Models are used both in interactive analysis, as well as scripted or command line based analyses. For this reason they must implement an easy to use API as well as serialization to disk. Gammapy should provide a variety of built-in models, allow to combine existing models using limited arithmetics and offer the possibility for users to easily implement custom models. A lot of development work has already put into the modeling framework of Gammapy. This PIG outlines the missing development work on the model framework for Gammapy v1.0. Proposal ======== Introduce naming scheme for models ---------------------------------- Currently Gammapy implements mainly four categories of models: spectral models, morphology or spatial models, combined spectral and spatial models (``SkyModel``) and lightcurve or phase models. We propose a uniform naming scheme for these model classes, based on the prefixes ``Spectral...``, ``Spatial...``, ``Temporal``: +-----------------+----------------------------------------------+----------------------+----------------------+-----------------------------+ | Base class | Template class | IRF folded class |Background class | Model classes | +=================+==============================================+======================+======================+=============================+ |SpectralModel |SpectralTemplate |SpectralIRFModel |SpectralBackground | Powerlaw, LogParabola, ... | |(SpectralModel) |(TableModel) | | | | +-----------------+----------------------------------------------+----------------------+----------------------+-----------------------------+ |SpatialModel |SpatialTemplate |SpatialIRFModel? |SpatialBackground? | Shell, SpatialGaussian, ... | |(SkySpatialModel)|(SkyDiffuseMap) | | | | +-----------------+----------------------------------------------+----------------------+----------------------+-----------------------------+ |TemporalModel |PhaseCurveTemplate, LightCurveTemplate | |TemporalBackground? | Sine, Exponential, ... | | |(PhaseCurveTableModel), (LightCurveTableModel)| | | | +-----------------+----------------------------------------------+----------------------+----------------------+-----------------------------+ |SourceModel |SourceTemplate |SourceIRFModel |SourceBackground |Arithmetic spectral/spatial | |(SkyModel) |(SkyDiffuseCube) | | |models? | +-----------------+----------------------------------------------+----------------------+----------------------+-----------------------------+ The current class names are given in parentheses. Optional model classes, where the need is unclear are marked with a question mark. Parametric models are named after their parametrization, such as ``Powerlaw`` or ``Gaussian``. In case the name is not unique (e.g. for a gaussian spectral and spatial model, the corresponding prefix should be used, such as ``SpectralGaussian`` or ``SpatialGaussian`` or ``SpatialConstant`` and ``SpectralConstant``). In addition there is the ``SourceModels`` class to represent a sum of ``SourceModel`` objects. Unify calling interface for models ---------------------------------- For model evaluation we would like to have a nice user interface, that supports ``Quantity`` and ``SkyCoord`` objects on the one hand and a performance oriented and flexible interface for efficient evaluation for model fitting. The latter We propose a unified implementation scheme for ``__call__`` and ``.evaluate()`` for all parametric models. The actual arthmetics of the model is defined in ``.evaluate()`` and should be implemented using numpy expressions only. This ensure compatibility with all the possible higher level callers that use ``Quantity``, ``SkyCoords`` or autograd tools. :: # for "normal" models it can be a static method @staticmethod def evaluate(energy, amplitude, index, reference): value = np.power(energy / reference, -index) return amplitude * value # "template" models, that need to access ``self.interpolate()`` def evaluate(self, energy, amplitude): value = self.interpolate(energy) return amplitude * value # "IRF" models, that need to access ``self.geom``, ``self.psf`` etc. def evaluate(self, **kwargs): value = self.model(self.energy, **kwargs) npred = self.apply_psf(value) return npred A nice user interface as well as parameter dispatching is implemented in the model base classe's ``__call__`` method: :: # spectral model def __call__(self, energy): pars = self.parameters.quantities value = self.evaluate(energy, **pars) return value def __call__(self, skycoord): pars = self.parameters.quantities lon = lon.wrap_at("") value = self.evaluate(lon, lat, **pars) return value Model evaluation for fitting will be first implement using the ``__call__`` interface. In case we find performance issue related to ``Quantity`` and ``SkyCoord`` handling or need auto differentiation tools such as ``autograd``, we can later implement a more efficient interface relying directly on ``Model.evaluate()``. Introduction of background models --------------------------------- For any kind of IACTs analysis hadronic background models are required, which model the residual hadronic (gamma-like) background emission from templates provided by DL3 or by generic parametric models. The background template differ from source templates, because they are defined in reconstructed spatial and energy coordinates. We propose the introduction of the following two model classes: BackgroundModel ~~~~~~~~~~~~~~~ Base class for all background models. BackgroundIRFModel ~~~~~~~~~~~~~~~~~~ Class to represent a template background model. It is initialized with a background map template and introduces additional background parameters to adjust the spectral shape and amplitude. This model can only be evaluated on a fixed grid, defined by the input map. The model is interpolated in the preparation step. No IRFs are applied. :: from gammapy.cube import BackgroundTemplate background_map = Map() background = BackgroundTemplate(template=background_map, norm=1, tilt=0.1) npred = background.evaluate() Introduction of "forward folded" models --------------------------------------- To combine parametric model components with their correspoding instrument response functions, such as PSF, energy dispersion and exposure, we propose the introduction of "forward folded" models. These take a fixed data grid, typically defined by the exposure map, on which the model is evaluated. After application of instrument reponse functions, the number of prediced counts is returned by the ``.evaluate()`` method. We propose to introduce the following classes: SpectralIRFModel ~~~~~~~~~~~~~~~~ A "forward folded" spectral model, that applies energy dispersion and effective area to a ``SpectralModel``. It corresponds in functionality to the current ``CountsPredictor``. :: from gammapy.spectrum import Powerlaw, SpectralIRFModel pwl = Powerlaw() spectral_irf_model = SpectralIRFModel(model=pwl, exposure=exposure, edisp=edisp) npred = spectral_irf_model.evaluate() SpatialIRFModel ~~~~~~~~~~~~~~~ A convolved model, that applies the PSF and exposure to a ``SpatialModel``. The need for this model is not clear, as we might solve the use case of morphology fitting by a combined spatial and spectral analysis with one energy bin. SourceIRFModel ~~~~~~~~~~~~~~ A "forward folded" source model, that applies IRFs to a ``SourceModel`` instance and returns an integrated quantity corresponding to predicted counts. It can only be evaluated on a fixed grid, which is defined on initialization by the exposure map. In addition it takes reduced PSF and EDISP objects. In functionality it corresponds to the current ``MapEvaluator``, but with the model parameters attached and the background handling removed. Additional "hidden" IRF parameters e.g. to propagate systematics could be optionally introduced later. This class will replace the current ``MapEvaluator``. :: from gammapy.cube import SourceIRFModel, SourceModel model = SourceModel() source_irf_model = SourceIRFModel(model, exposure=exposure, edisp=edisp, psf=psf) npred = source_irf_model.evaluate() Improve SourceModels class -------------------------- The existing ``CompoundSkyModel`` is a nice very generic abstraction to support any kind of arithmetics between ``SkyModel`` objects, but the number of use cases for other operators, except for "+" is very limited and can always be achieved by implementing a custom model. Serilization and component handling of this hierarchical model component structure is intrinsically difficult. For this reason we propose to first support and improve the existing ``SourceModels`` that implements an easier to handle flat hierarchy for model components. Support for arbitrary model arithmetics can be introduced, if needed, after Gammapy v1.0. We propose to remove the ``CompoundSourceModel`` and reimplement the ``+`` operator using the ``SourceModels`` class. :: from gammapy.cube import SourceModel, SourceModels component_1 = SourceModel() component_2 = SourceModel() total_model = component_1 + component_2 assert isinstance(total_model, SourceModels) # becomes quivalent to total_model = SourceModels([component_1, component_2]) Possibly remove the ``SpectralCompoundModel`` for consistency. Introduction of model name attributes ------------------------------------- The ``SourceModel`` class should be improved to support model component names. This way the ``SourceModels`` object can be improved to access model components by name, as well as the XML serialization (which supports names) becomes complete. The following example should work: :: from gammapy.cube import SourceModel, SourceModels component_1 = SourceModel(name="source_1") component_2 = SourceModel(name="source_2") total_model = SourceModels([component_1, component_2]) total_model["source_1"].parameters["index"].value = 2.3 # or alternatively?? total_model.parameters["source_1.index"].value = 2.3 # delete a model component in place del model["source_2"] This is also simplifies the parameter access in the fitting back-end, because parameter names become unique. E.g. no need for cryptic ``par_00X_`` parameter names in the minuit backend, which simplifies debugging and interaction by the user with methods such as ``Fit.likelihood_profile()`` or ``Fit.confidence()``, where parameter identifiers must be given. Improve and implement model serilization ---------------------------------------- Serialization of models is an important feature for standard analyses. In ``gammapy.utils.serialization`` there is a prototype implementation of XML serialization for ``SourceModels``, We propose to fix this existing XML serilization prototype so that the following easily works: :: from gammapy.cube import SourceModels model = SourceModels.read("model.xml") print(model) model.write("model.xml") In addition we propose to implement an additional YAML serialization format, which results in human-readable model configuration files. Once available the YAML format should be the preferred serialization: :: model = SourceModels.read("model.yaml") print(model) model.write("model.yaml") The two file formats should be handled automatically in ``SourceModel.read()`` and ``SourceModel.write()``. Improve spatial models ---------------------- Implement sky coordinate handling ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Implement support for coordinate systems. Add an interface for ``SkyCoord``, by introducing a ``SpatialModel.skycood`` attribute. :: from astropy.coordinates import SkyCoord from gammapy.image.models import PointSource # instanciation using a SkyCoord object point_source = PointSource(skycood=SkyCoord(0, 0, frame="galactic", unit="deg")) # or point_source = PointSource(lon="0 deg", lat="0 deg", frame="galactic") skycoord = point_source.skycoord assert isinstance(skycoord, SkyCoord) # evaluation of models using sky coordinates coords = SkyCoord([0, 1], [0.5, 0.5], frame="galactic", unit="deg") values = point_source(coords) # maybe override evaluation for efficient evaluation for fitting values = point_source(lon, lat) Implement default parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For spectral models Gammapy has defined default parameters, which simplifes the interactive usage of the models. From this experience we propose to introduce default parameters for ``SpatialModel`` as well as ``SourceModel``. :: from gammapy.image.models import Shell from gammapy.cube.models import SourceModels shell = Shell() print(shell) Implement evaluation region specifications ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For efficient evaluation of model components, we propose that ``SpatialModel`` classes declare a circular region ``.evaluation_radius``, they are evaluated on. :: from gammapy.image.models import Shell shell = Shell() print(shell.evaluation_radius) The ``.evaluation_radius`` should be a property of the model, that is computed depending of the model parameter, that determines the spatial extension of the model. The handling of ``PointSource`` and additional evaluation margins, which depend on the pixel size, is done in the step, where the model component is evaluated. Expose model parameters as attributes ------------------------------------- To simplify the interactive handling of model instances, model parameters should be exposed as attribues on the model, so that the following works: :: from gammapy.spectrum.models import Powerlaw pwl = Powerlaw() pwl.amplitude.value = "1e-12 cm-2 s-1 TeV-1" # instead of pwl.parameters["amplitude"].value = "1e-12 cm-2 s-1 TeV-1" The attributes should be auto-generated in the ``Model`` base class, based on the list of parameters defined by the model. Add new parametric models ------------------------- We propose to implement the following parametric models: - SkyGaussianElongated: Non radial-symmetric Gaussian. In planar approximation... - SkyDiskElongated: Non radial-symmetric disk. In planar approximation... - SpectralGaussian: Spectral Gaussian model. ``LogGaussian``? Task list ========= This is a proposal for a list of pull requests implementing the proposed changes, ordered by priority: 1. Implement a ``BackgroundModel`` base class and a ``BackgroundIRFModel`` class. Use the ``BackgroundIRFModel`` object in the current ``MapEvaluator``. Join the background parameters with the source model parameters and expose the joined parameter list at the ``MapEvaluator`` object, so that it can be used be used in current ``MapFit``. (**v0.10**) 2. Reimplement the "+" operator for ``SkyModel`` using the ``SkyModels`` class. Remove the ``CompoundSkyModel`` class. (**v0.XX**) 3. Implement support for model component names in ``SkyModel``. Implement ``SkyModels.__getitem__`` to allow for access of model components by name. (**v0.XX**) 4. Fix existing XML serilization of ``SkyModels`` and add serialization of the model component name. (**v0.XX**) 5. Implement a prototype YAML serialization for ``SkyModels``. Add ``SkyModels.from_yaml()`` and ``SkyModels.to_yaml()`` methods. Extend ``SkyModels.read()`` to recognise the file type. (**v0.XX**) 6. Implement support for coordinate systems in the ``SpatialModel`` base class. Add a ``SpatialModel.skycoord`` attribute to return a ``SkyCoord`` object and allow initialization of the spatial model positions with ``SkyCoord``. (**v0.XX**) 7. Implement the ``evaluation_radius`` property for all spatial models. (**v0.XX**) 8. Implement default parameter values for spatial models and ``SourceModel``. (**v0.XX**) 9. Expose model parameters as attributes on the models. (**v0.XX**) 10. Rename all models according to the proposed naming scheme. (**v0.XX**) 11. Rename existing ``MapEvaluator`` to ``SourceIRFModel`` and remove background handling. Temporarily move the background handling to the ``MapFit`` class. (**v0.XX**) 12. Rename existing ``CountsPredictor`` to ``SpectralIRFModel`` and unify API with ``SourceIRFModel`` as much as possible. (**v0.XX**) Alternatives ============ The proposed model names are open for discussion. Instead of ``Spectral...``, ``Spatial...`` and ``Source...`` one could also use the suffixes ``...1D``, ``...2D`` or ``...3D``. Instead of ``...Template`` one could use ``...TableModel``. The ``...IRFModel`` class are basically evaluating models on given coordinate grids and correspond in functionality to the current ``MapEvaluator``. A naming scheme based on ``ModelEvaluator`` seems a good alternative proposal. This would include a ``SpectralEvaluator`` and a ``SourceEvaluator`` and ``BackgroundEvaluator`` class. Decision ======== The authors decided to withdraw the PIG. Most of the proposed changes have been discussed and implemented independently. .. _gammapy: https://github.com/gammapy/gammapy .. _gammapy-web: https://github.com/gammapy/gammapy-webpage .. _GH 1971: https://github.com/gammapy/gammapy/pull/1971