files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Built-in models in Gammapy."""
  from gammapy.utils.registry import Registry
  from .core import *
  from .cube import *
  from .spatial import *
  from .spectral import *
  from .spectral_cosmic_ray import *
  from .spectral_crab import *
  from .temporal import *
  
  SPATIAL_MODEL_REGISTRY = Registry(
      [
          ConstantSpatialModel,
          TemplateSpatialModel,
          DiskSpatialModel,
          GaussianSpatialModel,
          GeneralizedGaussianSpatialModel,
          PointSpatialModel,
          ShellSpatialModel,
      ]
  )
  """Registry of spatial model classes."""
  
  SPECTRAL_MODEL_REGISTRY = Registry(
      [
          ConstantSpectralModel,
          CompoundSpectralModel,
          PowerLawSpectralModel,
          PowerLaw2SpectralModel,
          BrokenPowerLawSpectralModel,
          SmoothBrokenPowerLawSpectralModel,
          PiecewiseBrokenPowerLawSpectralModel,
          ExpCutoffPowerLawSpectralModel,
          ExpCutoffPowerLaw3FGLSpectralModel,
          SuperExpCutoffPowerLaw3FGLSpectralModel,
          SuperExpCutoffPowerLaw4FGLSpectralModel,
          LogParabolaSpectralModel,
          TemplateSpectralModel,
          GaussianSpectralModel,
          AbsorbedSpectralModel,
          NaimaSpectralModel,
          ScaleSpectralModel,
          PowerLawNormSpectralModel,
          LogParabolaNormSpectralModel,
          ExpCutoffPowerLawNormSpectralModel,
      ]
  )
  """Registry of spectral model classes."""
  
  TEMPORAL_MODEL_REGISTRY = Registry(
      [
          ConstantTemporalModel,
          LightCurveTemplateTemporalModel,
          ExpDecayTemporalModel,
          GaussianTemporalModel,
      ]
  )
  """Registry of temporal models classes."""
  
  MODEL_REGISTRY = Registry([SkyModel, BackgroundModel])
  """Registry of model classes"""
  
  
  __all__ = [
      "MODEL_REGISTRY",
      "SPATIAL_MODEL_REGISTRY",
      "TEMPORAL_MODEL_REGISTRY",
      "SPECTRAL_MODEL_REGISTRY",
      "SkyModelBase",
      "Models",
      "SkyModel",
      "BackgroundModel",
      "create_crab_spectral_model",
      "create_cosmic_ray_spectral_model",
      "Absorption",
      "SpatialModel",
      "SpectralModel",
      "TemporalModel",
  ]
  
  __all__.extend(cls.__name__ for cls in SPATIAL_MODEL_REGISTRY)
  __all__.extend(cls.__name__ for cls in SPECTRAL_MODEL_REGISTRY)
  __all__.extend(cls.__name__ for cls in TEMPORAL_MODEL_REGISTRY)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/__init__.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Built-in models in Gammapy."""
  from gammapy.utils.registry import Registry
  from .core import *
  from .cube import *
  from .spatial import *
  from .spectral import *
  from .spectral_cosmic_ray import *
  from .spectral_crab import *
  from .temporal import *
  
  SPATIAL_MODEL_REGISTRY = Registry(
      [
          ConstantSpatialModel,
          TemplateSpatialModel,
          DiskSpatialModel,
          GaussianSpatialModel,
          GeneralizedGaussianSpatialModel,
          PointSpatialModel,
          ShellSpatialModel,
      ]
  )
  """Registry of spatial model classes."""
  
  SPECTRAL_MODEL_REGISTRY = Registry(
      [
          ConstantSpectralModel,
          CompoundSpectralModel,
          PowerLawSpectralModel,
          PowerLaw2SpectralModel,
          BrokenPowerLawSpectralModel,
          SmoothBrokenPowerLawSpectralModel,
          ExpCutoffPowerLawSpectralModel,
          ExpCutoffPowerLaw3FGLSpectralModel,
          SuperExpCutoffPowerLaw3FGLSpectralModel,
          SuperExpCutoffPowerLaw4FGLSpectralModel,
          LogParabolaSpectralModel,
          TemplateSpectralModel,
          GaussianSpectralModel,
          AbsorbedSpectralModel,
          NaimaSpectralModel,
          ScaleSpectralModel,
          PowerLawNormSpectralModel,
          LogParabolaNormSpectralModel,
          ExpCutoffPowerLawNormSpectralModel,
      ]
  )
  """Registry of spectral model classes."""
  
  TEMPORAL_MODEL_REGISTRY = Registry(
      [
          ConstantTemporalModel,
          LightCurveTemplateTemporalModel,
          ExpDecayTemporalModel,
          GaussianTemporalModel,
      ]
  )
  """Registry of temporal models classes."""
  
  MODEL_REGISTRY = Registry([SkyModel, BackgroundModel])
  """Registry of model classes"""
  
  
  __all__ = [
      "MODEL_REGISTRY",
      "SPATIAL_MODEL_REGISTRY",
      "TEMPORAL_MODEL_REGISTRY",
      "SPECTRAL_MODEL_REGISTRY",
      "SkyModelBase",
      "Models",
      "SkyModel",
      "BackgroundModel",
      "create_crab_spectral_model",
      "create_cosmic_ray_spectral_model",
      "Absorption",
      "SpatialModel",
      "SpectralModel",
      "TemporalModel",
  ]
  
  __all__.extend(cls.__name__ for cls in SPATIAL_MODEL_REGISTRY)
  __all__.extend(cls.__name__ for cls in SPECTRAL_MODEL_REGISTRY)
  __all__.extend(cls.__name__ for cls in TEMPORAL_MODEL_REGISTRY)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/__init__.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
      _type = "spectral"
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically if no analytical solution defined.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          if hasattr(self, "evaluate_integral"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_integral(emin, emax, **kwargs)
          else:
              return integrate_spectrum(self, emin, emax, **kwargs)
  
      def integral_error(self, emin, emax):
          """Evaluate the error of the integral flux of a given spectrum in
          a given energy range.
  
          Parameters
          ----------
          emin, emax :  `~astropy.units.Quantity`
              Lower and upper bound of integration range.
  
          Returns
          -------
          flux, flux_err : tuple of `~astropy.units.Quantity`
              Integral flux and flux error betwen emin and emax.
          """
          energy = np.sqrt(emin * emax)
          flux = self.integral(emin, emax, intervals=True)
          dnde, dnde_err = self.evaluate_error(energy, epsilon=1e-4)
          flux_err = flux * dnde_err / dnde
          return u.Quantity([flux.value, flux_err.value], unit=flux.unit)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          if hasattr(self, "evaluate_energy_flux"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_energy_flux(emin, emax, **kwargs)
          else:
              return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit :
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
      _type = "spectral"
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically if no analytical solution defined.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          if hasattr(self, "evaluate_integral"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_integral(emin, emax, **kwargs)
          else:
              return integrate_spectrum(self, emin, emax, **kwargs)
  
      def integral_error(self, emin, emax):
          """Evaluate the error of the integral flux of a given spectrum in
          a given energy range.
  
          Parameters
          ----------
          emin, emax :  `~astropy.units.Quantity`
              Lower and upper bound of integration range.
  
          Returns
          -------
          flux, flux_err : tuple of `~astropy.units.Quantity`
              Integral flux and flux error betwen emin and emax.
          """
          energy = np.sqrt(emin * emax)
          flux = self.integral(emin, emax, intervals=True)
          dnde, dnde_err = self.evaluate_error(energy, epsilon=1e-4)
          flux_err = flux * dnde_err / dnde
          return u.Quantity([flux.value, flux_err.value], unit=flux.unit)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          if hasattr(self, "evaluate_energy_flux"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_energy_flux(emin, emax, **kwargs)
          else:
              return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit :
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
      PiecewiseBrokenPowerLawSpectralModel,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert "ExpCutoffPowerLawSpectralModel" in model0.spectral_model.tag
      assert "PointSpatialModel" in model0.spatial_model.tag
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert "pl" in model1.spectral_model.tag
      assert "PowerLawSpectralModel" in model1.spectral_model.tag
      assert "DiskSpatialModel" in model1.spatial_model.tag
      assert "disk" in model1.spatial_model.tag
      assert "LightCurveTemplateTemporalModel" in model1.temporal_model.tag
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s)"
  
      assert "TemplateSpectralModel" in model2.spectral_model.tag
      assert "TemplateSpatialModel" in model2.spatial_model.tag
  
      assert not model2.spatial_model.normalize
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  def test_PiecewiseBrokenPowerLawSpectralModel_io(tmp_path):
  
      energy = [1, 3, 7, 10] * u.TeV
      values = [1, 5, 3, 0.5] * u.Unit("cm-2 s-1 TeV-1")
      with pytest.raises(ValueError):
          PiecewiseBrokenPowerLawSpectralModel(
              energy=[1,] * u.TeV, values=[1, 5] * u.Unit("cm-2 s-1 TeV-1")
          )
      with pytest.raises(ValueError):
          PiecewiseBrokenPowerLawSpectralModel(
              energy=[1,] * u.TeV, values=[1,] * u.Unit("cm-2 s-1 TeV-1")
          )
      model = PiecewiseBrokenPowerLawSpectralModel(energy=energy, values=values)
      model.parameters[0].value = 2
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      for k in range(len(parnames)):
          assert parnames[k] == f"norm{k}"
  
      new_model = PiecewiseBrokenPowerLawSpectralModel.from_dict(model_dict)
      assert_allclose(new_model.parameters[0].value, 2)
      assert_allclose(new_model.energy, energy)
      assert_allclose(new_model.values, [2, 5, 3, 0.5] * values.unit)
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("pl", "spectral"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert "PowerLawSpectralModel" in new_model.spectral_model.tag
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel", "spatial")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", "spatial", map=map_constantmodel)
      yield Model.create(
          "DiskSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg", r_0="3 deg"
      )
      yield Model.create("gauss", "spatial", lon_0="1 deg", lat_0="2 deg
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_io.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert "ExpCutoffPowerLawSpectralModel" in model0.spectral_model.tag
      assert "PointSpatialModel" in model0.spatial_model.tag
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert "pl" in model1.spectral_model.tag
      assert "PowerLawSpectralModel" in model1.spectral_model.tag
      assert "DiskSpatialModel" in model1.spatial_model.tag
      assert "disk" in model1.spatial_model.tag
      assert "LightCurveTemplateTemporalModel" in model1.temporal_model.tag
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s)"
  
      assert "TemplateSpectralModel" in model2.spectral_model.tag
      assert "TemplateSpatialModel" in model2.spatial_model.tag
  
      assert not model2.spatial_model.normalize
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("pl", "spectral"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert "PowerLawSpectralModel" in new_model.spectral_model.tag
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel", "spatial")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", "spatial", map=map_constantmodel)
      yield Model.create(
          "DiskSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg", r_0="3 deg"
      )
      yield Model.create("gauss", "spatial", lon_0="1 deg", lat_0="2 deg", sigma="3 deg")
      yield Model.create("PointSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg")
      yield Model.create(
          "ShellSpatialModel",
          "spatial",
          lon_0="1 deg",
          lat_0="2 deg",
          radius="3 deg",
          width="4 deg",
      )
      yield Model.create("ConstantSpectralModel", "spectral", const="99 cm-2 s-1 TeV-1")
      yield Model.create(
          "CompoundSpectralModel",
          "spectral",
          model1=Model.create("PowerLawSpectralModel", "spectral"),
          model2=Model.create("PowerLawSpectralModel", "spectral"),
          operator=np.add,
      )
      yield Model.create("PowerLawSpectralModel", "spectral")
      yield Model.create("PowerLawNormSpectralModel", "spectral")
      yield Model.create("PowerLaw2SpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawNormSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw4FGLSpectralModel", "spectral")
      yield Model.create("LogParabolaSpectralModel", "spectral")
      yield Model.create("LogParabolaNormS
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_io.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
      PiecewiseBrokenPowerLawSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2SpectralModel(
              amplitude=u.Quantity(2.9227116204223784, "cm-2 s-1"),
              index=2.3 * u.Unit(""),
              emin=1 * u.TeV,
              emax=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.6 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.765838739678921, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV cm-2 s-1"),
          e_peak=4 * u.TeV,
      ),
      dict(
          name="norm-ecpl",
          model=ExpCutoffPowerLawNormSpectralModel(
              index=1.6 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, ""),
          integral_1_10TeV=u.Quantity(3.765838739678921, "TeV"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV2"),
      ),
      dict(
          name="ecpl_3fgl",
          model=ExpCutoffPowerLaw3FGLSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              ecut=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.7349563611124971, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.6034046173089, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(5.340285560055799, "TeV cm-2 s-1"),
      ),
      dict(
          name="plsec_4fgl",
          model=SuperExpCutoffPowerLaw4FGLSpectralModel(
              index_1=1.5,
              index_2=2,
              amplitude=1 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              expfactor=1e-2,
          ),
          val_at_2TeV=u.Quantity(0.3431043087721737, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(1.2125247, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.38072082, "TeV cm-2 s-1"),
      ),
      dict(
          name="logpar",
          model=LogParabolaSpectralModel(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="norm-logpar",
          model=LogParabolaNormSpectralModel(
              alpha=2.3 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, ""),
          integral_1_10TeV=u.Quantity(2.255689748270628, "TeV"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV2"),
      ),
      dict(
          name="logpar10",
          model=LogParabolaSpectralModel.from_log
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2SpectralModel(
              amplitude=u.Quantity(2.9227116204223784, "cm-2 s-1"),
              index=2.3 * u.Unit(""),
              emin=1 * u.TeV,
              emax=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.6 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.765838739678921, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV cm-2 s-1"),
          e_peak=4 * u.TeV,
      ),
      dict(
          name="norm-ecpl",
          model=ExpCutoffPowerLawNormSpectralModel(
              index=1.6 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, ""),
          integral_1_10TeV=u.Quantity(3.765838739678921, "TeV"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV2"),
      ),
      dict(
          name="ecpl_3fgl",
          model=ExpCutoffPowerLaw3FGLSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              ecut=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.7349563611124971, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.6034046173089, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(5.340285560055799, "TeV cm-2 s-1"),
      ),
      dict(
          name="plsec_4fgl",
          model=SuperExpCutoffPowerLaw4FGLSpectralModel(
              index_1=1.5,
              index_2=2,
              amplitude=1 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              expfactor=1e-2,
          ),
          val_at_2TeV=u.Quantity(0.3431043087721737, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(1.2125247, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.38072082, "TeV cm-2 s-1"),
      ),
      dict(
          name="logpar",
          model=LogParabolaSpectralModel(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="norm-logpar",
          model=LogParabolaNormSpectralModel(
              alpha=2.3 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, ""),
          integral_1_10TeV=u.Quantity(2.255689748270628, "TeV"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV2"),
      ),
      dict(
          name="logpar10",
          model=LogParabolaSpectralModel.from_log10(
              alpha=2.3 * u.
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_spectral.py