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)
  
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)
  
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 : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          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 error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623).
  
              When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your
              plotting code and the error band extends to negative values, it is not
              shown correctly. To circumvent this issue also use
              ``plt.loglog(nonposx='clip', nonposy='clip')``
              or ``plt.semilogy(nonposy='clip')``.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
          **kwargs : dict
              Keyword arguments forwarded to `matplotlib.pyplot.fill_between`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("facecolor", "black")
          kwargs.setdefault("alpha", 0.2)
          kwargs.setdefault("linewidth", 0)
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          flux, flux_err = self.evaluate_error(energy).to(flux_unit)
  
          y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power)
          y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power)
  
          where = (energy >= energy_range[0]) & (energy <= energy_range[1])
          ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs)
  
          self._plot_format_ax(ax, energy, y_lo, energy_power)
          return ax
  
      @staticmethod
      def _plot_format_ax(ax, energy, y, energy_power):
          ax.set_xlabel(f"Energy [{energy.unit}]")
          if energy_power > 0:
              ax.set_ylabel(f"E{energy_power} * Flux [{y.unit}]")
          else:
              ax.set_ylabel(f"Flux [{y.unit}]")
  
          ax.set_xscale("log", nonposx="clip")
          ax.set_yscale("log", nonposy="clip")
  
      @staticmethod
      def _plot_scale_flux(energy, flux, energy_power):
          try:
              eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0]
          except IndexError:
              eunit = energy.unit
          y = flux * np.power(energy, energy_power)
          return y.to(flux.unit * eunit ** energy_power)
  
      def spectral_index(self, energy, epsilon=1e-5):
          """Compute spectral index at given energy.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to estimate the index
          epsilon : float
              Fractional energy increment to use for determining the spectral index.
  
          Returns
          -------
          index : float
              Estimated spectral index.
          """
          f1 = self(energy)
          f2 = self(energy * (1 + epsilon))
          return np.log(f1 / f2) / np.log(1 + epsilon)
  
      def inverse(self, value, emin=0.1 * u.TeV, emax=100 * u.TeV):
          """Return energy for a given function value of the spectral model.
  
          Calls the `scipy.optimize.brentq` numerical root finding method.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          emin : `~astropy.units.Quantity`
              Lower bracket value in case solution is not unique.
          emax : `~astropy.units.Quantity`
              Upper bracket value in case solution is not unique.
  
          Returns
          -------
          energy : `~astropy.units.Quantity`
              Energies at which the model has the given ``value``.
          """
          eunit = "TeV"
  
          energies = []
          for val in np.atleast_1d(value):
  
              def f(x):
                  # scale by 1e12 to achieve better precision
                  energy = u.Quantity(x, eunit, copy=False)
                  y = self(energy).to_value(value.unit)
                  return 1e12 * (y - val.value)
  
              energy = scipy.optimize.brentq(
                  f, emin.to_value(eunit), emax.to_value(eunit)
              )
              energies.append(energy)
  
          return u.Quantity(energies, eunit, copy=False)
  
  
  class ConstantSpectralModel(SpectralModel):
      r"""Constant model.
  
      For more information see :ref:`constant-spectral-model`.
  
      Parameters
      ----------
      const : `~astropy.units.Quantity`
          :math:`k`
      """
  
      tag = ["ConstantSpectralModel", "const"]
      const = Parameter("const", "1e-12 cm-2 s-1 TeV-1")
  
      @staticmethod
      def evaluate(energy, const):
          """Evaluate the model (static function)."""
          return np.ones(np.atleast_1d(energy).shape) * const
  
  
  class CompoundSpectralModel(SpectralModel):
      """Arithmetic combination of two spectral models.
  
      For more information see :ref:`compound-spectral-model`.
      """
  
      tag = ["CompoundSpectralModel", "compound"]
  
      def __init__(self, model1, model2, operator):
          self.model1 = model1
          self.model2 = model2
          self.operator = operator
          super().__init__()
  
      @property
      def parameters(self):
          return self.model1.parameters + self.model2.parameters
  
      def __str__(self):
          return (
              f"{self.__class__.__name__}\n"
              f"    Component 1 : {self.model1}\n"
              f"    Component 2 : {self.model2}\n"
              f"    Operator : {self.operator.__name__}\n"
          )
  
      def __call__(self, energy):
          val1 = self.model1(energy)
          val2 = self.model2(energy)
          return self.operator(val1, val2)
  
      def to_dict(self):
          return {
              "type": self.tag[0],
              "model1": self.model1.to_dict(),
              "model2": self.model2.to_dict(),
              "operator": self.operator.__name__,
          }
  
      @classmethod
      def from_dict(cls, data):
          from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY
  
          model1_cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["model1"]["type"])
          model1 = model1_cls.from_dict(data["model1"])
          model2_cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["model2"]["type"])
          model2 = model2_cls.from_dict(data["model2"])
          op = getattr(operator, data["operator"])
          return cls(model1, model2, op)
  
  
  class PowerLawSpectralModel(SpectralModel):
      r"""Spectral power-law model.
  
      For more information see :ref:`powerlaw-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
          
      See Also
      --------
      PowerLaw2SpectralModel, PowerLawNormSpectralModel
      """
  
      tag = ["PowerLawSpectralModel", "pl"]
      index = Parameter("index", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference):
          """Evaluate the model (static function)."""
          return amplitude * np.power((energy / reference), -index)
  
      @staticmethod
      def evaluate_integral(emin, emax, index, amplitude, reference):
          r"""Integrate power law analytically (static function).
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE = \left.
              \phi_0 \frac{E_0}{-\Gamma + 1} \left( \frac{E}{E_0} \right)^{-\Gamma + 1}
              \right \vert _{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range
          """
          val = -1 * index + 1
  
          prefactor = amplitude * reference / val
          upper = np.power((emax / reference), val)
          lower = np.power((emin / reference), val)
          integral = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              integral[mask] = (amplitude * reference * np.log(emax / emin))[mask]
  
          return integral
  
      @staticmethod
      def evaluate_energy_flux(emin, emax, index, amplitude, reference):
          r"""Compute energy flux in given energy range analytically (static function).
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE = \left.
              \phi_0 \frac{E_0^2}{-\Gamma + 2} \left( \frac{E}{E_0} \right)^{-\Gamma + 2}
              \right \vert _{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          val = -1 * index + 2
  
          prefactor = amplitude * reference ** 2 / val
          upper = (emax / reference) ** val
          lower = (emin / reference) ** val
          energy_flux = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)
              # for reference
              energy_flux[mask] = amplitude * reference ** 2 * np.log(emax / emin)[mask]
  
          return energy_flux
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          base = value / self.amplitude.quantity
          return self.reference.quantity * np.power(base, -1.0 / self.index.value)
  
      @property
      def pivot_energy(self):
          r"""The decorrelation energy is defined as:
  
          .. math::
  
              E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
  
          Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
          """
          index_err = self.index.error
          reference = self.reference.quantity
          amplitude = self.amplitude.quantity
          cov_index_ampl = self.covariance.data[0, 1] * amplitude.unit
          return reference * np.exp(cov_index_ampl / (amplitude * index_err ** 2))
  
  
  class PowerLawNormSpectralModel(SpectralModel):
      r"""Spectral power-law model with normalized amplitude parameter.
  
      Parameters
      ----------
      tilt : `~astropy.units.Quantity`
          :math:`\Gamma`
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
          
      See Also
      --------
      PowerLawSpectralModel, PowerLaw2SpectralModel
      """
  
      tag = ["PowerLawNormSpectralModel", "pl-norm"]
      norm = Parameter("norm", 1, unit="")
      tilt = Parameter("tilt", 0)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, tilt, norm, reference):
          """Evaluate the model (static function)."""
          return norm * np.power((energy / reference), -tilt)
  
      @staticmethod
      def evaluate_integral(emin, emax, tilt, norm, reference):
          """Evaluate pwl integral."""
          val = -1 * tilt + 1
  
          prefactor = norm * reference / val
          upper = np.power((emax / reference), val)
          lower = np.power((emin / reference), val)
          integral = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              integral[mask] = (norm * reference * np.log(emax / emin))[mask]
  
          return integral
  
      @staticmethod
      def evaluate_energy_flux(emin, emax, tilt, norm, reference):
          """Evaluate the energy flux (static function)"""
          val = -1 * tilt + 2
  
          prefactor = norm * reference ** 2 / val
          upper = (emax / reference) ** val
          lower = (emin / reference) ** val
          energy_flux = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)
              # for reference
              energy_flux[mask] = norm * reference ** 2 * np.log(emax / emin)[mask]
  
          return energy_flux
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          base = value / self.norm.quantity
          return self.reference.quantity * np.power(base, -1.0 / self.tilt.value)
  
      @property
      def pivot_energy(self):
          r"""The decorrelation energy is defined as:
  
          .. math::
  
              E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
  
          Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
          """
          tilt_err = self.tilt.error
          reference = self.reference.quantity
          norm = self.norm.quantity
          cov_tilt_norm = self.covariance.data[0, 1] * norm.unit
          return reference * np.exp(cov_tilt_norm / (norm * tilt_err ** 2))
  
  
  class PowerLaw2SpectralModel(SpectralModel):
      r"""Spectral power-law model with integral as amplitude parameter.
  
      For more information see :ref:`powerlaw2-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          Spectral index :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          Integral flux :math:`F_0`.
      emin : `~astropy.units.Quantity`
          Lower energy limit :math:`E_{0, min}`.
      emax : `~astropy.units.Quantity`
          Upper energy limit :math:`E_{0, max}`.
      
      See Also
      --------
      PowerLawSpectralModel, PowerLawNormSpectralModel
      """
      tag = ["PowerLaw2SpectralModel", "pl-2"]
  
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1")
      index = Parameter("index", 2)
      emin = Parameter("emin", "0.1 TeV", frozen=True)
      emax = Parameter("emax", "100 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, amplitude, index, emin, emax):
          """Evaluate the model (static function)."""
          top = -index + 1
  
          # to get the energies dimensionless we use a modified formula
          bottom = emax - emin * (emin / emax) ** (-index)
          return amplitude * (top / bottom) * np.power(energy / emax, -index)
  
      @staticmethod
      def evaluate_integral(emin_lim, emax_lim, amplitude, index, emin, emax):
          r"""Integrate power law analytically.
  
          .. math::
              F(E_{min}, E_{max}) = F_0 \cdot \frac{E_{max}^{\Gamma + 1} \
                                  - E_{min}^{\Gamma + 1}}{E_{0, max}^{\Gamma + 1} \
                                  - E_{0, min}^{\Gamma + 1}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          temp1 = np.power(emax_lim, -index.value + 1)
          temp2 = np.power(emin_lim, -index.value + 1)
          top = temp1 - temp2
  
          temp1 = np.power(emax, -index.value + 1)
          temp2 = np.power(emin, -index.value + 1)
          bottom = temp1 - temp2
  
          return amplitude * top / bottom
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          amplitude = self.amplitude.quantity
          index = self.index.value
          emin = self.emin.quantity
          emax = self.emax.quantity
  
          # to get the energies dimensionless we use a modified formula
          top = -index + 1
          bottom = emax - emin * (emin / emax) ** (-index)
          term = (bottom / top) * (value / amplitude)
          return np.power(term.to_value(""), -1.0 / index) * emax
  
  
  class BrokenPowerLawSpectralModel(SpectralModel):
      r"""Spectral broken power-law model.
  
      For more information see :ref:`broken-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index1 : `~astropy.units.Quantity`
          :math:`\Gamma1`
      index2 : `~astropy.units.Quantity`
          :math:`\Gamma2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      ebreak : `~astropy.units.Quantity`
          :math:`E_{break}`
          
      See Also
      --------
      SmoothBrokenPowerLawSpectralModel
      """
  
      tag = ["BrokenPowerLawSpectralModel", "bpl"]
      index1 = Parameter("index1", 2.0)
      index2 = Parameter("index2", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      ebreak = Parameter("ebreak", "1 TeV")
  
      @staticmethod
      def evaluate(energy, index1, index2, amplitude, ebreak):
          """Evaluate the model (static function)."""
          energy = np.atleast_1d(energy)
          cond = energy < ebreak
          bpwl = amplitude * np.ones(energy.shape)
          bpwl[cond] *= (energy[cond] / ebreak) ** (-index1)
          bpwl[~cond] *= (energy[~cond] / ebreak) ** (-index2)
          return bpwl
  
  
  class SmoothBrokenPowerLawSpectralModel(SpectralModel):
      r"""Spectral smooth broken power-law model.
  
      For more information see :ref:`smooth-broken-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index1 : `~astropy.units.Quantity`
          :math:`\Gamma1`
      index2 : `~astropy.units.Quantity`
          :math:`\Gamma2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ebreak : `~astropy.units.Quantity`
          :math:`E_{break}`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
          
      See Also
      --------
      BrokenPowerLawSpectralModel
      """
  
      tag = ["SmoothBrokenPowerLawSpectralModel", "sbpl"]
      index1 = Parameter("index1", 2.0)
      index2 = Parameter("index2", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      ebreak = Parameter("ebreak", "1 TeV")
      reference = Parameter("reference", "1 TeV", frozen=True)
      beta = Parameter("beta", 1, frozen=True)
  
      @staticmethod
      def evaluate(energy, index1, index2, amplitude, ebreak, reference, beta):
          """Evaluate the model (static function)."""
          beta *= np.sign(index2 - index1)
          pwl = amplitude * (energy / reference) ** (-index1)
          brk = (1 + (energy / ebreak) ** ((index2 - index1) / beta)) ** (-beta)
          return pwl * brk
  
  
  class PiecewiseBrokenPowerLawSpectralModel(SpectralModel):
      """Piecewise broken power-law at fixed energy nodes.
  
      Parameters
      ----------
      energy : `~astropy.units.Quantity`
          Array of energies at which the model values are given (nodes).
      values : array
          Array with the initial values of the model at energies ``energy``.
          A normalisation parameters is created for each value.
      """
  
      tag = "PiecewiseBrokenPowerLawSpectralModel"
  
      def __init__(self, energy, values, parameters=None):
          self._energy = energy
          self.init_values = values
          if len(values) != len(energy):
              raise ValueError("dimension mismatch")
          if len(values) < 2:
              raise ValueError("Input arrays must contians at least 2 elements")
          if parameters is None:
              parameters = Parameters(
                  [Parameter(f"norm{k}", 1.0) for k in range(len(values))]
              )
          for parameter in parameters:
              setattr(self, parameter.name, parameter)
          self.default_parameters = parameters
  
      @classmethod
      def from_parameters(cls, parameters, energy, values):
          init = cls(energy, values, parameters=parameters)
          return init
  
      @property
      def values(self):
          return np.array([p.value for p in self.parameters]) * self.init_values
  
      @property
      def energy(self):
          return self._energy
  
      def __call__(self, energy):
          return self.evaluate(energy)
  
      def evaluate(self, energy):
          logedata = np.log10(np.atleast_1d(energy.value))
          loge = np.log10(self.energy.to(energy.unit).value)
          logv = np.log10(self.values.value)
          ne = len(loge)
          conds = (
              [(logedata < loge[1])]
              + [
                  (logedata >= loge[k]) & (logedata < loge[k + 1])
                  for k in range(1, ne - 2)
              ]
              + [(logedata >= loge[-2])]
          )
          a = (logv[1:] - logv[:-1]) / (loge[1:] - loge[:-1])
          b = logv[1:] - a * loge[1:]
  
          output = np.zeros(logedata.shape)
          for k in range(ne - 1):
              output[conds[k]] = 10 ** (a[k] * logedata[conds[k]] + b[k])
          return output * self.values.unit
  
      def to_dict(self):
          return {
              "type": self.tag,
              "parameters": self.parameters.to_dict(),
              "energy": {
                  "data": self.energy.data.tolist(),
                  "unit": str(self.energy.unit),
              },
              "values": {
                  "data": self.init_values.data.tolist(),
                  "unit": str(self.values.unit),
              },
          }
  
      @classmethod
      def from_dict(cls, data):
          energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
          values = u.Quantity(data["values"]["data"], data["values"]["unit"])
          if "parameters" in data:
              parameters = Parameters.from_dict(data["parameters"])
              return cls.from_parameters(parameters, energy=energy, values=values)
          else:
              return cls(energy=energy, values=values)
  
  
  class ExpCutoffPowerLawSpectralModel(SpectralModel):
      r"""Spectral exponential cutoff power-law model.
  
      For more information see :ref:`exp-cutoff-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      lambda_ : `~astropy.units.Quantity`
          :math:`\lambda`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
  
      See Also
      --------
      ExpCutoffPowerLawNormSpectralModel
      """
  
      tag = ["ExpCutoffPowerLawSpectralModel", "ecpl"]
  
      index = Parameter("index", 1.5)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      lambda_ = Parameter("lambda_", "0.1 TeV-1")
      alpha = Parameter("alpha", "1.0", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference, lambda_, alpha):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index)
          cutoff = np.exp(-np.power(energy * lambda_, alpha))
  
          return pwl * cutoff
  
      @property
      def e_peak(self):
          r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`).
  
          This is the peak in E^2 x dN/dE and is given by:
  
          .. math::
              E_{Peak} =  \left(\frac{2 - \Gamma}{\alpha}\right)^{1/\alpha} / \lambda
          """
          reference = self.reference.quantity
          index = self.index.quantity
          lambda_ = self.lambda_.quantity
          alpha = self.alpha.quantity
  
          if index >= 2 or lambda_ == 0.0 or alpha == 0.0:
              return np.nan * reference.unit
          else:
              return np.power((2 - index) / alpha, 1 / alpha) / lambda_
  
  
  class ExpCutoffPowerLawNormSpectralModel(SpectralModel):
      r"""Norm spectral exponential cutoff power-law model.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      lambda_ : `~astropy.units.Quantity`
          :math:`\lambda`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
  
      See Also
      --------
      ExpCutoffPowerLawSpectralModel
      """
      tag = ["ExpCutoffPowerLawNormSpectralModel", "ecpl-norm"]
  
      index = Parameter("index", 1.5)
      norm = Parameter("norm", 1, unit="")
      reference = Parameter("reference", "1 TeV", frozen=True)
      lambda_ = Parameter("lambda_", "0.1 TeV-1")
      alpha = Parameter("alpha", "1.0", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, norm, reference, lambda_, alpha):
          """Evaluate the model (static function)."""
          pwl = norm * (energy / reference) ** (-index)
          cutoff = np.exp(-np.power(energy * lambda_, alpha))
  
          return pwl * cutoff
  
  
  class ExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
      r"""Spectral exponential cutoff power-law model used for 3FGL.
  
      For more information see :ref:`exp-cutoff-powerlaw-3fgl-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ecut : `~astropy.units.Quantity`
          :math:`E_{C}`
      """
  
      tag = ["ExpCutoffPowerLaw3FGLSpectralModel", "ecpl-3fgl"]
      index = Parameter("index", 1.5)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      ecut = Parameter("ecut", "10 TeV")
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference, ecut):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index)
          cutoff = np.exp((reference - energy) / ecut)
          return pwl * cutoff
  
  
  class SuperExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
      r"""Spectral super exponential cutoff power-law model used for 3FGL.
  
      For more information see :ref:`super-exp-cutoff-powerlaw-3fgl-spectral-model`.
  
      .. math::
          \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma_1}
                    \exp \left( \left(\frac{E_0}{E_{C}} \right)^{\Gamma_2} -
                                \left(\frac{E}{E_{C}} \right)^{\Gamma_2}
                                \right)
  
      Parameters
      ----------
      index_1 : `~astropy.units.Quantity`
          :math:`\Gamma_1`
      index_2 : `~astropy.units.Quantity`
          :math:`\Gamma_2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ecut : `~astropy.units.Quantity`
          :math:`E_{C}`
      """
  
      tag = ["SuperExpCutoffPowerLaw3FGLSpectralModel", "secpl-3fgl"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      ecut = Parameter("ecut", "10 TeV")
      index_1 = Parameter("index_1", 1.5)
      index_2 = Parameter("index_2", 2)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, ecut, index_1, index_2):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index_1)
          cutoff = np.exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2)
          return pwl * cutoff
  
  
  class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
      r"""Spectral super exponential cutoff power-law model used for 4FGL.
  
      For more information see :ref:`super-exp-cutoff-powerlaw-4fgl-spectral-model`.
  
      Parameters
      ----------
      index_1 : `~astropy.units.Quantity`
          :math:`\Gamma_1`
      index_2 : `~astropy.units.Quantity`
          :math:`\Gamma_2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      expfactor : `~astropy.units.Quantity`
          :math:`a`, given as dimensionless value but
          internally assumes unit of :math:`[E_0]` power :math:`-\Gamma_2`
      """
  
      tag = ["SuperExpCutoffPowerLaw4FGLSpectralModel", "secpl-4fgl"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      expfactor = Parameter("expfactor", "1e-2")
      index_1 = Parameter("index_1", 1.5)
      index_2 = Parameter("index_2", 2)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, expfactor, index_1, index_2):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index_1)
          cutoff = np.exp(
              expfactor
              / reference.unit ** index_2
              * (reference ** index_2 - energy ** index_2)
          )
          return pwl * cutoff
  
  
  class LogParabolaSpectralModel(SpectralModel):
      r"""Spectral log parabola model.
  
      For more information see :ref:`logparabola-spectral-model`.
  
      Parameters
      ----------
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
  
      See Also
      --------
      LogParabolaNormSpectralModel
  
      """
      tag = ["LogParabolaSpectralModel", "lp"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "10 TeV", frozen=True)
      alpha = Parameter("alpha", 2)
      beta = Parameter("beta", 1)
  
      @classmethod
      def from_log10(cls, amplitude, reference, alpha, beta):
          """Construct from :math:`log_{10}` parametrization."""
          beta_ = beta / np.log(10)
          return cls(amplitude=amplitude, reference=reference, alpha=alpha, beta=beta_)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, alpha, beta):
          """Evaluate the model (static function)."""
          xx = energy / reference
          exponent = -alpha - beta * np.log(xx)
          return amplitude * np.power(xx, exponent)
  
      @property
      def e_peak(self):
          r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`).
  
          This is the peak in E^2 x dN/dE and is given by:
  
          .. math::
              E_{Peak} = E_{0} \exp{ (2 - \alpha) / (2 * \beta)}
          """
          reference = self.reference.quantity
          alpha = self.alpha.quantity
          beta = self.beta.quantity
          return reference * np.exp((2 - alpha) / (2 * beta))
  
  
  class LogParabolaNormSpectralModel(SpectralModel):
      r"""Norm spectral log parabola model.
  
      Parameters
      ----------
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
  
      See Also
      --------
      LogParabolaSpectralModel
      """
      tag = ["LogParabolaNormSpectralModel", "lp-norm"]
      norm = Parameter("norm", 1, unit="")
      reference = Parameter("reference", "10 TeV", frozen=True)
      alpha = Parameter("alpha", 2)
      beta = Parameter("beta", 1)
  
      @classmethod
      def from_log10(cls, norm, reference, alpha, beta):
          """Construct from :math:`log_{10}` parametrization."""
          beta_ = beta / np.log(10)
          return cls(norm=norm, reference=reference, alpha=alpha, beta=beta_)
  
      @staticmethod
      def evaluate(energy, norm, reference, alpha, beta):
          """Evaluate the model (static function)."""
          xx = energy / reference
          exponent = -alpha - beta * np.log(xx)
          return norm * np.power(xx, exponent)
  
  
  class TemplateSpectralModel(SpectralModel):
      """A model generated from a table of energy and value arrays.
  
      For more information see :ref:`template-spectral-model`.
  
      Parameters
      ----------
      energy : `~astropy.units.Quantity`
          Array of energies at which the model values are given
      values : array
          Array with the values of the model at energies ``energy``.
      interp_kwargs : dict
          Interpolation keyword arguments pass to `scipy.interpolate.interp1d`.
          By default all values outside the interpolation range are set to zero.
          If you want to apply linear extrapolation you can pass `interp_kwargs={'fill_value':
          'extrapolate', 'kind': 'linear'}`. If you want to choose the interpolation
          scaling applied to values, you can use `interp_kwargs={"values_scale": "log"}`.
      meta : dict, optional
          Meta information, meta['filename'] will be used for serialization
      """
  
      tag = ["TemplateSpectralModel", "template"]
  
      def __init__(
          self, energy, values, interp_kwargs=None, meta=None,
      ):
          self.energy = energy
          self.values = u.Quantity(values, copy=False)
          self.meta = dict() if meta is None else meta
          interp_kwargs = interp_kwargs or {}
          interp_kwargs.setdefault("values_scale", "log")
          interp_kwargs.setdefault("points_scale", ("log",))
  
          self._evaluate = ScaledRegularGridInterpolator(
              points=(energy,), values=values, **interp_kwargs
          )
  
          super().__init__()
  
      @classmethod
      def read_xspec_model(cls, filename, param, **kwargs):
          """Read XSPEC table model.
  
          The input is a table containing absorbed values from a XSPEC model as a
          function of energy.
  
          TODO: Format of the file should be described and discussed in
          https://gamma-astro-data-formats.readthedocs.io/en/latest/index.html
  
          Parameters
          ----------
          filename : str
              File containing the XSPEC model
          param : float
              Model parameter value
  
          Examples
          --------
          Fill table from an EBL model (Franceschini, 2008)
  
          >>> from gammapy.modeling.models import TemplateSpectralModel
          >>> filename = '$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz'
          >>> table_model = TemplateSpectralModel.read_xspec_model(filename=filename, param=0.3)
          """
          filename = make_path(filename)
  
          # Check if parameter value is in range
          table_param = Table.read(filename, hdu="PARAMETERS")
          pmin = table_param["MINIMUM"]
          pmax = table_param["MAXIMUM"]
          if param < pmin or param > pmax:
              raise ValueError(f"Out of range: param={param}, min={pmin}, max={pmax}")
  
          # Get energy values
          table_energy = Table.read(filename, hdu="ENERGIES")
          energy_lo = table_energy["ENERG_LO"]
          energy_hi = table_energy["ENERG_HI"]
  
          # set energy to log-centers
          energy = np.sqrt(energy_lo * energy_hi)
  
          # Get spectrum values (no interpolation, take closest value for param)
          table_spectra = Table.read(filename, hdu="SPECTRA")
          idx = np.abs(table_spectra["PARAMVAL"] - param).argmin()
          values = u.Quantity(table_spectra[idx][1], "", copy=False)  # no dimension
  
          kwargs.setdefault("interp_kwargs", {"values_scale": "lin"})
          return cls(energy=energy, values=values, **kwargs)
  
      def evaluate(self, energy):
          """Evaluate the model (static function)."""
          return self._evaluate((energy,), clip=True)
  
      def to_dict(self):
          return {
              "type": self.tag[0],
              "energy": {
                  "data": self.energy.data.tolist(),
                  "unit": str(self.energy.unit),
              },
              "values": {
                  "data": self.values.data.tolist(),
                  "unit": str(self.values.unit),
              },
          }
  
      @classmethod
      def from_dict(cls, data):
          energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
          values = u.Quantity(data["values"]["data"], data["values"]["unit"])
          return cls(energy=energy, values=values)
  
  
  class ScaleSpectralModel(SpectralModel):
      """Wrapper to scale another spectral model by a norm factor.
  
      Parameters
      ----------
      model : `SpectralModel`
          Spectral model to wrap.
      norm : float
          Multiplicative norm factor for the model value.
      """
  
      tag = ["ScaleSpectralModel", "scale"]
      norm = Parameter("norm", 1, unit="")
  
      def __init__(self, model, norm=norm.quantity):
          self.model = model
          self._covariance = None
          super().__init__(norm=norm)
  
      def evaluate(self, energy, norm):
          return norm * self.model(energy)
  
      def integral(self, emin, emax, **kwargs):
          return self.norm.value * self.model.integral(emin, emax, **kwargs)
  
  
  class Absorption:
      r"""Gamma-ray absorption models.
  
      For more information see :ref:`absorption-spectral-model`.
  
      Parameters
      ----------
      energy : `~astropy.units.Quantity`
          Energy node values
      param : `~astropy.units.Quantity`
          Parameter node values
      data : `~astropy.units.Quantity`
          Model value
      filename : str
          Filename of the absorption model used for serialisation.
      interp_kwargs : dict
          Interpolation option passed to `ScaledRegularGridInterpolator`.
          By default the models are extrapolated outside the range. To prevent
          this and raise an error instead use interp_kwargs = {"extrapolate": False}
      """
  
      tag = "Absorption"
  
      def __init__(self, energy, param, data, filename=None, interp_kwargs=None):
          self.filename = filename
          # set values log centers
          self.param = param
          self.energy = energy
          self.energy = energy
          self.data = u.Quantity(data, copy=False)
  
          interp_kwargs = interp_kwargs or {}
          interp_kwargs.setdefault("points_scale", ("log", "lin"))
          interp_kwargs.setdefault("extrapolate", True)
  
          self._evaluate = ScaledRegularGridInterpolator(
              points=(self.param, self.energy), values=self.data, **interp_kwargs
          )
  
      def to_dict(self):
          if self.filename is None:
              return {
                  "type": self.tag,
                  "energy": {
                      "data": self.energy.data.tolist(),
                      "unit": str(self.energy.unit),
                  },
                  "param": {
                      "data": self.param.data.tolist(),
                      "unit": str(self.param.unit),
                  },
                  "values": {
                      "data": self.data.data.tolist(),
                      "unit": str(self.data.unit),
                  },
              }
          else:
              return {"type": self.tag, "filename": self.filename}
  
      @classmethod
      def from_dict(cls, data):
  
          if "filename" in data:
              return cls.read(data["filename"])
          else:
              energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
              param = u.Quantity(data["param"]["data"], data["param"]["unit"])
              values = u.Quantity(data["values"]["data"], data["values"]["unit"])
              return cls(energy=energy, param=param, data=values)
  
      @classmethod
      def read(cls, filename, interp_kwargs=None):
          """Build object from an XSPEC model.
  
          Todo: Format of XSPEC binary files should be referenced at https://gamma-astro-data-formats.readthedocs.io/en/latest/
  
          Parameters
          ----------
          filename : str
              File containing the model.
          interp_kwargs : dict
              Interpolation option passed to `ScaledRegularGridInterpolator`.
  
          Returns
          -------
          absorption : `Absorption`
              Absorption model.
          """
          # Create EBL data array
          filename = make_path(filename)
          table_param = Table.read(filename, hdu="PARAMETERS")
  
          # TODO: for some reason the table contain duplicated values
          param, idx = np.unique(table_param[0]["VALUE"], return_index=True)
  
          # Get energy values
          table_energy = Table.read(filename, hdu="ENERGIES")
          energy_lo = u.Quantity(
              table_energy["ENERG_LO"], "keV", copy=False
          )  # unit not stored in file
          energy_hi = u.Quantity(
              table_energy["ENERG_HI"], "keV", copy=False
          )  # unit not stored in file
          energy = np.sqrt(energy_lo * energy_hi)
  
          # Get spectrum values
          table_spectra = Table.read(filename, hdu="SPECTRA")
          data = table_spectra["INTPSPEC"].data[idx, :]
          return cls(
              energy=energy,
              param=param,
              data=data,
              filename=filename,
              interp_kwargs=interp_kwargs,
          )
  
      @classmethod
      def read_builtin(cls, name, interp_kwargs=None):
          """Read one of the built-in absorption models.
  
          Parameters
          ----------
          name : {'franceschini', 'dominguez', 'finke'}
              name of one of the available model in gammapy-data
  
          References
          ----------
          .. [1] Franceschini et al., "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity",
              `Link <https://ui.adsabs.harvard.edu/abs/2008A%26A...487..837F>`__
          .. [2] Dominguez et al., " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions"
              `Link <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__
          .. [3] Finke et al., "Modeling the Extragalactic Background Light from Stars and Dust"
              `Link <https://ui.adsabs.harvard.edu/abs/2010ApJ...712..238F>`__
  
          Returns
          -------
          absorption : `Absorption`
              Absorption model.
  
          """
          models = dict()
          models["franceschini"] = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz"
          models["dominguez"] = "$GAMMAPY_DATA/ebl/ebl_dominguez11.fits.gz"
          models["finke"] = "$GAMMAPY_DATA/ebl/frd_abs.fits.gz"
  
          return cls.read(models[name], interp_kwargs=interp_kwargs)
  
      def table_model(self, parameter):
          """Table model for a given parameter value.
  
          Parameters
          ----------
          parameter : float
              Parameter value.
  
          Returns
          -------
          template_model : `TemplateSpectralModel`
              Template spectral model.
          """
          energy = self.energy
          values = self.evaluate(energy=energy, parameter=parameter)
          return TemplateSpectralModel(
              energy=energy, values=values, interp_kwargs={"values_scale": "log"}
          )
  
      def evaluate(self, energy, parameter):
          """Evaluate model for energy and parameter value."""
          return np.clip(self._evaluate((parameter, energy)), 0, 1)
  
  
  class AbsorbedSpectralModel(SpectralModel):
      r"""Spectral model with EBL absorption.
  
      For more information see :ref:`absorbed-spectral-model`.
  
      Parameters
      ----------
      spectral_model : `SpectralModel`
          Spectral model.
      absorption : `Absorption`
          Absorption model.
      redshift : float
          Redshift of the absorption model
      alpha_norm: float
          Norm of the EBL model
      """
  
      tag = ["AbsorbedSpectralModel", "ebl-absorption"]
      alpha_norm = Parameter("alpha_norm", 1.0, frozen=True)
      redshift = Parameter("redshift", 0.1, frozen=True)
  
      def __init__(
          self, spectral_model, absorption, redshift, alpha_norm=alpha_norm.quantity,
      ):
          self.spectral_model = spectral_model
          self.absorption = absorption
  
          min_ = self.absorption.param.min()
          max_ = self.absorption.param.max()
  
          redshift = Parameter("redshift", redshift, frozen=True, min=min_, max=max_)
          super().__init__(redshift=redshift, alpha_norm=alpha_norm)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          self._covariance.set_subcovariance(self.spectral_model.covariance)
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          subcovar = self._covariance.get_subcovariance(
              self.spectral_model.covariance.parameters
          )
          self.spectral_model.covariance = subcovar
  
      @property
      def parameters(self):
          return (
              Parameters([self.redshift, self.alpha_norm])
              + self.spectral_model.parameters
          )
  
      def evaluate(self, energy, **kwargs):
          """Evaluate the model at a given energy."""
          # assign redshift value and remove it from dictionary
          # since it does not belong to the spectral model
          parameter = kwargs.pop("redshift")
          alpha_norm = kwargs.pop("alpha_norm")
  
          dnde = self.spectral_model.evaluate(energy=energy, **kwargs)
          absorption = self.absorption.evaluate(energy=energy, parameter=parameter)
          # Power rule: (e ^ a) ^ b = e ^ (a * b)
          absorption = np.power(absorption, alpha_norm)
          return dnde * absorption
  
      def to_dict(self):
          return {
              "type": self.tag,
              "base_model": self.spectral_model.to_dict(),
              "absorption": self.absorption.to_dict(),
              "absorption_parameter": {"name": "redshift", "value": self.redshift.value,},
              "parameters": Parameters([self.redshift, self.alpha_norm]).to_dict(),
          }
  
      @classmethod
      def from_dict(cls, data):
          from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY
  
          model_class = SPECTRAL_MODEL_REGISTRY.get_cls(data["base_model"]["type"])
  
          model = cls(
              spectral_model=model_class.from_dict(data["base_model"]),
              absorption=Absorption.from_dict(data["absorption"]),
              redshift=data["absorption_parameter"]["value"],
          )
          return model
  
  
  class NaimaSpectralModel(SpectralModel):
      r"""A wrapper for Naima models.
  
      For more information see :ref:`naima-spectral-model`.
  
      Parameters
      ----------
      radiative_model : `~naima.models.BaseRadiative`
          An instance of a radiative model defined in `~naima.models`
      distance : `~astropy.units.Quantity`, optional
          Distance to the source. If set to 0, the intrinsic differential
          luminosity will be returned. Default is 1 kpc
      seed : str or list of str, optional
          Seed photon field(s) to be considered for the `radiative_model` flux computation,
          in case of a `~naima.models.InverseCompton` model. It can be a subset of the
          `seed_photon_fields` list defining the `radiative_model`. Default is the whole list
          of photon fields
      nested_models : dict
          Additionnal parameters for nested models not supplied by the radiative model,
          for now this is used  only for synchrotron self-compton model
      """
  
      tag = ["NaimaSpectralModel", "naima"]
  
      def __init__(
          self, radiative_model, distance=1.0 * u.kpc, seed=None, nested_models=None
      ):
          import naima
  
          self.radiative_model = radiative_model
          self._particle_distribution = self.radiative_model.particle_distribution
          self.distance = u.Quantity(distance)
          self.seed = seed
  
          if nested_models is None:
              nested_models = {}
  
          self.nested_models = nested_models
  
          if isinstance(self._particle_distribution, naima.models.TableModel):
              param_names = ["amplitude"]
          else:
              param_names = self._particle_distribution.param_names
  
          parameters = []
  
          for name in param_names:
              value = getattr(self._particle_distribution, name)
              parameter = Parameter(name, value)
              parameters.append(parameter)
  
          # In case of a synchrotron radiative model, append B to the fittable parameters
          if "B" in self.radiative_model.param_names:
              value = getattr(self.radiative_model, "B")
              parameter = Parameter("B", value)
              parameters.append(parameter)
  
          # In case of a synchrotron self compton model, append B and Rpwn to the fittable parameters
          if (
              isinstance(self.radiative_model, naima.models.InverseCompton)
              and "SSC" in self.nested_models
          ):
              B = self.nested_models["SSC"]["B"]
              radius = self.nested_models["SSC"]["radius"]
              parameters.append(Parameter("B", B))
              parameters.append(Parameter("radius", radius, frozen=True))
  
          self.default_parameters = Parameters(parameters)
          super().__init__()
  
      def _evaluate_ssc(
          self, energy,
      ):
          """
          Compute photon density spectrum from synchrotron emission for synchrotron self-compton model,
          assuming uniform synchrotron emissivity inside a sphere of radius R
          (see Section 4.1 of Atoyan & Aharonian 1996)
  
          based on :
          "https://naima.readthedocs.io/en/latest/examples.html#crab-nebula-ssc-model"
  
          """
          import naima
  
          SYN = naima.models.Synchrotron(
              self._particle_distribution,
              B=self.B.quantity,
              Eemax=self.radiative_model.Eemax,
              Eemin=self.radiative_model.Eemin,
          )
  
          Esy = np.logspace(-7, 9, 100) * u.eV
          Lsy = SYN.flux(Esy, distance=0 * u.cm)  # use distance 0 to get luminosity
          phn_sy = Lsy / (4 * np.pi * self.radius.quantity ** 2 * const.c) * 2.24
          # The factor 2.24 comes from the assumption on uniform synchrotron
          # emissivity inside a sphere
  
          if "SSC" not in self.radiative_model.seed_photon_fields:
              self.radiative_model.seed_photon_fields["SSC"] = {
                  "isotropic": True,
                  "type": "array",
                  "energy": Esy,
                  "photon_density": phn_sy,
              }
          else:
              self.radiative_model.seed_photon_fields["SSC"]["photon_density"] = phn_sy
  
          dnde = self.radiative_model.flux(
              energy, seed=self.seed, distance=self.distance
          ) + SYN.flux(energy, distance=self.distance)
          return dnde
  
      def evaluate(self, energy, **kwargs):
          """Evaluate the model."""
          import naima
  
          for name, value in kwargs.items():
              setattr(self._particle_distribution, name, value)
  
          if "B" in self.radiative_model.param_names:
              self.radiative_model.B = self.B.quantity
  
          if (
              isinstance(self.radiative_model, naima.models.InverseCompton)
              and "SSC" in self.nested_models
          ):
              dnde = self._evaluate_ssc(energy.flatten())
          elif self.seed is not None:
              dnde = self.radiative_model.flux(
                  energy.flatten(), seed=self.seed, distance=self.distance
              )
          else:
              dnde = self.radiative_model.flux(energy.flatten(), distance=self.distance)
  
          dnde = dnde.reshape(energy.shape)
          unit = 1 / (energy.unit * u.cm ** 2 * u.s)
          return dnde.to(unit)
  
      @classmethod
      def from_dict(cls, data):
          raise NotImplementedError(
              "Currently the NaimaSpectralModel cannot be read from YAML"
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          raise NotImplementedError(
              "Currently the NaimaSpectralModel cannot be built from a list of parameters."
          )
  
  
  class GaussianSpectralModel(SpectralModel):
      r"""Gaussian spectral model.
  
      For more information see :ref:`gaussian-spectral-model`.
  
      Parameters
      ----------
      norm : `~astropy.units.Quantity`
          :math:`N_0`
      mean : `~astropy.units.Quantity`
          :math:`\bar{E}`
      sigma : `~astropy.units.Quantity`
          :math:`\sigma`
      """
  
      tag = ["GaussianSpectralModel", "gauss"]
      norm = Parameter("norm", 1e-12 * u.Unit("cm-2 s-1"))
      mean = Parameter("mean", 1 * u.TeV)
      sigma = Parameter("sigma", 2 * u.TeV)
  
      @staticmethod
      def evaluate(energy, norm, mean, sigma):
          return (
              norm
              / (sigma * np.sqrt(2 * np.pi))
              * np.exp(-((energy - mean) ** 2) / (2 * sigma ** 2))
          )
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate Gaussian analytically.
  
          .. math::
              F(E_{min}, E_{max}) = \frac{N_0}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma})\right]_{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range
          """
          # kwargs are passed to this function but not used
          # this is to get a consistent API with SpectralModel.integral()
          u_min = (
              (emin - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          u_max = (
              (emax - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
  
          return (
              self.norm.quantity
              / 2
              * (scipy.special.erf(u_max) - scipy.special.erf(u_min))
          )
  
      def energy_flux(self, emin, emax):
          r"""Compute energy flux in given energy range analytically.
  
          .. math::
              G(E_{min}, E_{max}) =  \frac{N_0 \sigma}{\sqrt{2*\pi}}* \left[ - \exp(\frac{E_{min}-\bar{E}}{\sqrt{2} \sigma})
              \right]_{E_{min}}^{E_{max}} + \frac{N_0 * \bar{E}}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma})
               \right]_{E_{min}}^{E_{max}}
  
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          u_min = (
              (emin - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          u_max = (
              (emax - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          a = self.norm.quantity * self.sigma.quantity / np.sqrt(2 * np.pi)
          b = self.norm.quantity * self.mean.quantity / 2
          return a * (np.exp(-(u_min ** 2)) - np.exp(-(u_max ** 2))) + b * (
              scipy.special.erf(u_max) - scipy.special.erf(u_min)
          )
  
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 : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          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 error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623).
  
              When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your
              plotting code and the error band extends to negative values, it is not
              shown correctly. To circumvent this issue also use
              ``plt.loglog(nonposx='clip', nonposy='clip')``
              or ``plt.semilogy(nonposy='clip')``.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
          **kwargs : dict
              Keyword arguments forwarded to `matplotlib.pyplot.fill_between`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("facecolor", "black")
          kwargs.setdefault("alpha", 0.2)
          kwargs.setdefault("linewidth", 0)
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          flux, flux_err = self.evaluate_error(energy).to(flux_unit)
  
          y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power)
          y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power)
  
          where = (energy >= energy_range[0]) & (energy <= energy_range[1])
          ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs)
  
          self._plot_format_ax(ax, energy, y_lo, energy_power)
          return ax
  
      @staticmethod
      def _plot_format_ax(ax, energy, y, energy_power):
          ax.set_xlabel(f"Energy [{energy.unit}]")
          if energy_power > 0:
              ax.set_ylabel(f"E{energy_power} * Flux [{y.unit}]")
          else:
              ax.set_ylabel(f"Flux [{y.unit}]")
  
          ax.set_xscale("log", nonposx="clip")
          ax.set_yscale("log", nonposy="clip")
  
      @staticmethod
      def _plot_scale_flux(energy, flux, energy_power):
          try:
              eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0]
          except IndexError:
              eunit = energy.unit
          y = flux * np.power(energy, energy_power)
          return y.to(flux.unit * eunit ** energy_power)
  
      def spectral_index(self, energy, epsilon=1e-5):
          """Compute spectral index at given energy.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to estimate the index
          epsilon : float
              Fractional energy increment to use for determining the spectral index.
  
          Returns
          -------
          index : float
              Estimated spectral index.
          """
          f1 = self(energy)
          f2 = self(energy * (1 + epsilon))
          return np.log(f1 / f2) / np.log(1 + epsilon)
  
      def inverse(self, value, emin=0.1 * u.TeV, emax=100 * u.TeV):
          """Return energy for a given function value of the spectral model.
  
          Calls the `scipy.optimize.brentq` numerical root finding method.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          emin : `~astropy.units.Quantity`
              Lower bracket value in case solution is not unique.
          emax : `~astropy.units.Quantity`
              Upper bracket value in case solution is not unique.
  
          Returns
          -------
          energy : `~astropy.units.Quantity`
              Energies at which the model has the given ``value``.
          """
          eunit = "TeV"
  
          energies = []
          for val in np.atleast_1d(value):
  
              def f(x):
                  # scale by 1e12 to achieve better precision
                  energy = u.Quantity(x, eunit, copy=False)
                  y = self(energy).to_value(value.unit)
                  return 1e12 * (y - val.value)
  
              energy = scipy.optimize.brentq(
                  f, emin.to_value(eunit), emax.to_value(eunit)
              )
              energies.append(energy)
  
          return u.Quantity(energies, eunit, copy=False)
  
  
  class ConstantSpectralModel(SpectralModel):
      r"""Constant model.
  
      For more information see :ref:`constant-spectral-model`.
  
      Parameters
      ----------
      const : `~astropy.units.Quantity`
          :math:`k`
      """
  
      tag = ["ConstantSpectralModel", "const"]
      const = Parameter("const", "1e-12 cm-2 s-1 TeV-1")
  
      @staticmethod
      def evaluate(energy, const):
          """Evaluate the model (static function)."""
          return np.ones(np.atleast_1d(energy).shape) * const
  
  
  class CompoundSpectralModel(SpectralModel):
      """Arithmetic combination of two spectral models.
  
      For more information see :ref:`compound-spectral-model`.
      """
  
      tag = ["CompoundSpectralModel", "compound"]
  
      def __init__(self, model1, model2, operator):
          self.model1 = model1
          self.model2 = model2
          self.operator = operator
          super().__init__()
  
      @property
      def parameters(self):
          return self.model1.parameters + self.model2.parameters
  
      def __str__(self):
          return (
              f"{self.__class__.__name__}\n"
              f"    Component 1 : {self.model1}\n"
              f"    Component 2 : {self.model2}\n"
              f"    Operator : {self.operator.__name__}\n"
          )
  
      def __call__(self, energy):
          val1 = self.model1(energy)
          val2 = self.model2(energy)
          return self.operator(val1, val2)
  
      def to_dict(self):
          return {
              "type": self.tag[0],
              "model1": self.model1.to_dict(),
              "model2": self.model2.to_dict(),
              "operator": self.operator.__name__,
          }
  
      @classmethod
      def from_dict(cls, data):
          from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY
  
          model1_cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["model1"]["type"])
          model1 = model1_cls.from_dict(data["model1"])
          model2_cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["model2"]["type"])
          model2 = model2_cls.from_dict(data["model2"])
          op = getattr(operator, data["operator"])
          return cls(model1, model2, op)
  
  
  class PowerLawSpectralModel(SpectralModel):
      r"""Spectral power-law model.
  
      For more information see :ref:`powerlaw-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
          
      See Also
      --------
      PowerLaw2SpectralModel, PowerLawNormSpectralModel
      """
  
      tag = ["PowerLawSpectralModel", "pl"]
      index = Parameter("index", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference):
          """Evaluate the model (static function)."""
          return amplitude * np.power((energy / reference), -index)
  
      @staticmethod
      def evaluate_integral(emin, emax, index, amplitude, reference):
          r"""Integrate power law analytically (static function).
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE = \left.
              \phi_0 \frac{E_0}{-\Gamma + 1} \left( \frac{E}{E_0} \right)^{-\Gamma + 1}
              \right \vert _{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range
          """
          val = -1 * index + 1
  
          prefactor = amplitude * reference / val
          upper = np.power((emax / reference), val)
          lower = np.power((emin / reference), val)
          integral = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              integral[mask] = (amplitude * reference * np.log(emax / emin))[mask]
  
          return integral
  
      @staticmethod
      def evaluate_energy_flux(emin, emax, index, amplitude, reference):
          r"""Compute energy flux in given energy range analytically (static function).
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE = \left.
              \phi_0 \frac{E_0^2}{-\Gamma + 2} \left( \frac{E}{E_0} \right)^{-\Gamma + 2}
              \right \vert _{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          val = -1 * index + 2
  
          prefactor = amplitude * reference ** 2 / val
          upper = (emax / reference) ** val
          lower = (emin / reference) ** val
          energy_flux = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)
              # for reference
              energy_flux[mask] = amplitude * reference ** 2 * np.log(emax / emin)[mask]
  
          return energy_flux
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          base = value / self.amplitude.quantity
          return self.reference.quantity * np.power(base, -1.0 / self.index.value)
  
      @property
      def pivot_energy(self):
          r"""The decorrelation energy is defined as:
  
          .. math::
  
              E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
  
          Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
          """
          index_err = self.index.error
          reference = self.reference.quantity
          amplitude = self.amplitude.quantity
          cov_index_ampl = self.covariance.data[0, 1] * amplitude.unit
          return reference * np.exp(cov_index_ampl / (amplitude * index_err ** 2))
  
  
  class PowerLawNormSpectralModel(SpectralModel):
      r"""Spectral power-law model with normalized amplitude parameter.
  
      Parameters
      ----------
      tilt : `~astropy.units.Quantity`
          :math:`\Gamma`
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
          
      See Also
      --------
      PowerLawSpectralModel, PowerLaw2SpectralModel
      """
  
      tag = ["PowerLawNormSpectralModel", "pl-norm"]
      norm = Parameter("norm", 1, unit="")
      tilt = Parameter("tilt", 0)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, tilt, norm, reference):
          """Evaluate the model (static function)."""
          return norm * np.power((energy / reference), -tilt)
  
      @staticmethod
      def evaluate_integral(emin, emax, tilt, norm, reference):
          """Evaluate pwl integral."""
          val = -1 * tilt + 1
  
          prefactor = norm * reference / val
          upper = np.power((emax / reference), val)
          lower = np.power((emin / reference), val)
          integral = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              integral[mask] = (norm * reference * np.log(emax / emin))[mask]
  
          return integral
  
      @staticmethod
      def evaluate_energy_flux(emin, emax, tilt, norm, reference):
          """Evaluate the energy flux (static function)"""
          val = -1 * tilt + 2
  
          prefactor = norm * reference ** 2 / val
          upper = (emax / reference) ** val
          lower = (emin / reference) ** val
          energy_flux = prefactor * (upper - lower)
  
          mask = np.isclose(val, 0)
  
          if mask.any():
              # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)
              # for reference
              energy_flux[mask] = norm * reference ** 2 * np.log(emax / emin)[mask]
  
          return energy_flux
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          base = value / self.norm.quantity
          return self.reference.quantity * np.power(base, -1.0 / self.tilt.value)
  
      @property
      def pivot_energy(self):
          r"""The decorrelation energy is defined as:
  
          .. math::
  
              E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
  
          Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
          """
          tilt_err = self.tilt.error
          reference = self.reference.quantity
          norm = self.norm.quantity
          cov_tilt_norm = self.covariance.data[0, 1] * norm.unit
          return reference * np.exp(cov_tilt_norm / (norm * tilt_err ** 2))
  
  
  class PowerLaw2SpectralModel(SpectralModel):
      r"""Spectral power-law model with integral as amplitude parameter.
  
      For more information see :ref:`powerlaw2-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          Spectral index :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          Integral flux :math:`F_0`.
      emin : `~astropy.units.Quantity`
          Lower energy limit :math:`E_{0, min}`.
      emax : `~astropy.units.Quantity`
          Upper energy limit :math:`E_{0, max}`.
      
      See Also
      --------
      PowerLawSpectralModel, PowerLawNormSpectralModel
      """
      tag = ["PowerLaw2SpectralModel", "pl-2"]
  
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1")
      index = Parameter("index", 2)
      emin = Parameter("emin", "0.1 TeV", frozen=True)
      emax = Parameter("emax", "100 TeV", frozen=True)
  
      @staticmethod
      def evaluate(energy, amplitude, index, emin, emax):
          """Evaluate the model (static function)."""
          top = -index + 1
  
          # to get the energies dimensionless we use a modified formula
          bottom = emax - emin * (emin / emax) ** (-index)
          return amplitude * (top / bottom) * np.power(energy / emax, -index)
  
      @staticmethod
      def evaluate_integral(emin_lim, emax_lim, amplitude, index, emin, emax):
          r"""Integrate power law analytically.
  
          .. math::
              F(E_{min}, E_{max}) = F_0 \cdot \frac{E_{max}^{\Gamma + 1} \
                                  - E_{min}^{\Gamma + 1}}{E_{0, max}^{\Gamma + 1} \
                                  - E_{0, min}^{\Gamma + 1}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          temp1 = np.power(emax_lim, -index.value + 1)
          temp2 = np.power(emin_lim, -index.value + 1)
          top = temp1 - temp2
  
          temp1 = np.power(emax, -index.value + 1)
          temp2 = np.power(emin, -index.value + 1)
          bottom = temp1 - temp2
  
          return amplitude * top / bottom
  
      def inverse(self, value):
          """Return energy for a given function value of the spectral model.
  
          Parameters
          ----------
          value : `~astropy.units.Quantity`
              Function value of the spectral model.
          """
          amplitude = self.amplitude.quantity
          index = self.index.value
          emin = self.emin.quantity
          emax = self.emax.quantity
  
          # to get the energies dimensionless we use a modified formula
          top = -index + 1
          bottom = emax - emin * (emin / emax) ** (-index)
          term = (bottom / top) * (value / amplitude)
          return np.power(term.to_value(""), -1.0 / index) * emax
  
  
  class BrokenPowerLawSpectralModel(SpectralModel):
      r"""Spectral broken power-law model.
  
      For more information see :ref:`broken-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index1 : `~astropy.units.Quantity`
          :math:`\Gamma1`
      index2 : `~astropy.units.Quantity`
          :math:`\Gamma2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      ebreak : `~astropy.units.Quantity`
          :math:`E_{break}`
          
      See Also
      --------
      SmoothBrokenPowerLawSpectralModel
      """
  
      tag = ["BrokenPowerLawSpectralModel", "bpl"]
      index1 = Parameter("index1", 2.0)
      index2 = Parameter("index2", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      ebreak = Parameter("ebreak", "1 TeV")
  
      @staticmethod
      def evaluate(energy, index1, index2, amplitude, ebreak):
          """Evaluate the model (static function)."""
          energy = np.atleast_1d(energy)
          cond = energy < ebreak
          bpwl = amplitude * np.ones(energy.shape)
          bpwl[cond] *= (energy[cond] / ebreak) ** (-index1)
          bpwl[~cond] *= (energy[~cond] / ebreak) ** (-index2)
          return bpwl
  
  
  class SmoothBrokenPowerLawSpectralModel(SpectralModel):
      r"""Spectral smooth broken power-law model.
  
      For more information see :ref:`smooth-broken-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index1 : `~astropy.units.Quantity`
          :math:`\Gamma1`
      index2 : `~astropy.units.Quantity`
          :math:`\Gamma2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ebreak : `~astropy.units.Quantity`
          :math:`E_{break}`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
          
      See Also
      --------
      BrokenPowerLawSpectralModel
      """
  
      tag = ["SmoothBrokenPowerLawSpectralModel", "sbpl"]
      index1 = Parameter("index1", 2.0)
      index2 = Parameter("index2", 2.0)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      ebreak = Parameter("ebreak", "1 TeV")
      reference = Parameter("reference", "1 TeV", frozen=True)
      beta = Parameter("beta", 1, frozen=True)
  
      @staticmethod
      def evaluate(energy, index1, index2, amplitude, ebreak, reference, beta):
          """Evaluate the model (static function)."""
          beta *= np.sign(index2 - index1)
          pwl = amplitude * (energy / reference) ** (-index1)
          brk = (1 + (energy / ebreak) ** ((index2 - index1) / beta)) ** (-beta)
          return pwl * brk
  
  
  class ExpCutoffPowerLawSpectralModel(SpectralModel):
      r"""Spectral exponential cutoff power-law model.
  
      For more information see :ref:`exp-cutoff-powerlaw-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      lambda_ : `~astropy.units.Quantity`
          :math:`\lambda`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
  
      See Also
      --------
      ExpCutoffPowerLawNormSpectralModel
      """
  
      tag = ["ExpCutoffPowerLawSpectralModel", "ecpl"]
  
      index = Parameter("index", 1.5)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      lambda_ = Parameter("lambda_", "0.1 TeV-1")
      alpha = Parameter("alpha", "1.0", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference, lambda_, alpha):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index)
          cutoff = np.exp(-np.power(energy * lambda_, alpha))
  
          return pwl * cutoff
  
      @property
      def e_peak(self):
          r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`).
  
          This is the peak in E^2 x dN/dE and is given by:
  
          .. math::
              E_{Peak} =  \left(\frac{2 - \Gamma}{\alpha}\right)^{1/\alpha} / \lambda
          """
          reference = self.reference.quantity
          index = self.index.quantity
          lambda_ = self.lambda_.quantity
          alpha = self.alpha.quantity
  
          if index >= 2 or lambda_ == 0.0 or alpha == 0.0:
              return np.nan * reference.unit
          else:
              return np.power((2 - index) / alpha, 1 / alpha) / lambda_
  
  
  class ExpCutoffPowerLawNormSpectralModel(SpectralModel):
      r"""Norm spectral exponential cutoff power-law model.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      lambda_ : `~astropy.units.Quantity`
          :math:`\lambda`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
  
      See Also
      --------
      ExpCutoffPowerLawSpectralModel
      """
      tag = ["ExpCutoffPowerLawNormSpectralModel", "ecpl-norm"]
  
      index = Parameter("index", 1.5)
      norm = Parameter("norm", 1, unit="")
      reference = Parameter("reference", "1 TeV", frozen=True)
      lambda_ = Parameter("lambda_", "0.1 TeV-1")
      alpha = Parameter("alpha", "1.0", frozen=True)
  
      @staticmethod
      def evaluate(energy, index, norm, reference, lambda_, alpha):
          """Evaluate the model (static function)."""
          pwl = norm * (energy / reference) ** (-index)
          cutoff = np.exp(-np.power(energy * lambda_, alpha))
  
          return pwl * cutoff
  
  
  class ExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
      r"""Spectral exponential cutoff power-law model used for 3FGL.
  
      For more information see :ref:`exp-cutoff-powerlaw-3fgl-spectral-model`.
  
      Parameters
      ----------
      index : `~astropy.units.Quantity`
          :math:`\Gamma`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ecut : `~astropy.units.Quantity`
          :math:`E_{C}`
      """
  
      tag = ["ExpCutoffPowerLaw3FGLSpectralModel", "ecpl-3fgl"]
      index = Parameter("index", 1.5)
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      ecut = Parameter("ecut", "10 TeV")
  
      @staticmethod
      def evaluate(energy, index, amplitude, reference, ecut):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index)
          cutoff = np.exp((reference - energy) / ecut)
          return pwl * cutoff
  
  
  class SuperExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):
      r"""Spectral super exponential cutoff power-law model used for 3FGL.
  
      For more information see :ref:`super-exp-cutoff-powerlaw-3fgl-spectral-model`.
  
      .. math::
          \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma_1}
                    \exp \left( \left(\frac{E_0}{E_{C}} \right)^{\Gamma_2} -
                                \left(\frac{E}{E_{C}} \right)^{\Gamma_2}
                                \right)
  
      Parameters
      ----------
      index_1 : `~astropy.units.Quantity`
          :math:`\Gamma_1`
      index_2 : `~astropy.units.Quantity`
          :math:`\Gamma_2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      ecut : `~astropy.units.Quantity`
          :math:`E_{C}`
      """
  
      tag = ["SuperExpCutoffPowerLaw3FGLSpectralModel", "secpl-3fgl"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      ecut = Parameter("ecut", "10 TeV")
      index_1 = Parameter("index_1", 1.5)
      index_2 = Parameter("index_2", 2)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, ecut, index_1, index_2):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index_1)
          cutoff = np.exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2)
          return pwl * cutoff
  
  
  class SuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):
      r"""Spectral super exponential cutoff power-law model used for 4FGL.
  
      For more information see :ref:`super-exp-cutoff-powerlaw-4fgl-spectral-model`.
  
      Parameters
      ----------
      index_1 : `~astropy.units.Quantity`
          :math:`\Gamma_1`
      index_2 : `~astropy.units.Quantity`
          :math:`\Gamma_2`
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      expfactor : `~astropy.units.Quantity`
          :math:`a`, given as dimensionless value but
          internally assumes unit of :math:`[E_0]` power :math:`-\Gamma_2`
      """
  
      tag = ["SuperExpCutoffPowerLaw4FGLSpectralModel", "secpl-4fgl"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "1 TeV", frozen=True)
      expfactor = Parameter("expfactor", "1e-2")
      index_1 = Parameter("index_1", 1.5)
      index_2 = Parameter("index_2", 2)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, expfactor, index_1, index_2):
          """Evaluate the model (static function)."""
          pwl = amplitude * (energy / reference) ** (-index_1)
          cutoff = np.exp(
              expfactor
              / reference.unit ** index_2
              * (reference ** index_2 - energy ** index_2)
          )
          return pwl * cutoff
  
  
  class LogParabolaSpectralModel(SpectralModel):
      r"""Spectral log parabola model.
  
      For more information see :ref:`logparabola-spectral-model`.
  
      Parameters
      ----------
      amplitude : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
  
      See Also
      --------
      LogParabolaNormSpectralModel
  
      """
      tag = ["LogParabolaSpectralModel", "lp"]
      amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
      reference = Parameter("reference", "10 TeV", frozen=True)
      alpha = Parameter("alpha", 2)
      beta = Parameter("beta", 1)
  
      @classmethod
      def from_log10(cls, amplitude, reference, alpha, beta):
          """Construct from :math:`log_{10}` parametrization."""
          beta_ = beta / np.log(10)
          return cls(amplitude=amplitude, reference=reference, alpha=alpha, beta=beta_)
  
      @staticmethod
      def evaluate(energy, amplitude, reference, alpha, beta):
          """Evaluate the model (static function)."""
          xx = energy / reference
          exponent = -alpha - beta * np.log(xx)
          return amplitude * np.power(xx, exponent)
  
      @property
      def e_peak(self):
          r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`).
  
          This is the peak in E^2 x dN/dE and is given by:
  
          .. math::
              E_{Peak} = E_{0} \exp{ (2 - \alpha) / (2 * \beta)}
          """
          reference = self.reference.quantity
          alpha = self.alpha.quantity
          beta = self.beta.quantity
          return reference * np.exp((2 - alpha) / (2 * beta))
  
  
  class LogParabolaNormSpectralModel(SpectralModel):
      r"""Norm spectral log parabola model.
  
      Parameters
      ----------
      norm : `~astropy.units.Quantity`
          :math:`\phi_0`
      reference : `~astropy.units.Quantity`
          :math:`E_0`
      alpha : `~astropy.units.Quantity`
          :math:`\alpha`
      beta : `~astropy.units.Quantity`
          :math:`\beta`
  
      See Also
      --------
      LogParabolaSpectralModel
      """
      tag = ["LogParabolaNormSpectralModel", "lp-norm"]
      norm = Parameter("norm", 1, unit="")
      reference = Parameter("reference", "10 TeV", frozen=True)
      alpha = Parameter("alpha", 2)
      beta = Parameter("beta", 1)
  
      @classmethod
      def from_log10(cls, norm, reference, alpha, beta):
          """Construct from :math:`log_{10}` parametrization."""
          beta_ = beta / np.log(10)
          return cls(norm=norm, reference=reference, alpha=alpha, beta=beta_)
  
      @staticmethod
      def evaluate(energy, norm, reference, alpha, beta):
          """Evaluate the model (static function)."""
          xx = energy / reference
          exponent = -alpha - beta * np.log(xx)
          return norm * np.power(xx, exponent)
  
  
  class TemplateSpectralModel(SpectralModel):
      """A model generated from a table of energy and value arrays.
  
      For more information see :ref:`template-spectral-model`.
  
      Parameters
      ----------
      energy : `~astropy.units.Quantity`
          Array of energies at which the model values are given
      values : array
          Array with the values of the model at energies ``energy``.
      interp_kwargs : dict
          Interpolation keyword arguments pass to `scipy.interpolate.interp1d`.
          By default all values outside the interpolation range are set to zero.
          If you want to apply linear extrapolation you can pass `interp_kwargs={'fill_value':
          'extrapolate', 'kind': 'linear'}`. If you want to choose the interpolation
          scaling applied to values, you can use `interp_kwargs={"values_scale": "log"}`.
      meta : dict, optional
          Meta information, meta['filename'] will be used for serialization
      """
  
      tag = ["TemplateSpectralModel", "template"]
  
      def __init__(
          self, energy, values, interp_kwargs=None, meta=None,
      ):
          self.energy = energy
          self.values = u.Quantity(values, copy=False)
          self.meta = dict() if meta is None else meta
          interp_kwargs = interp_kwargs or {}
          interp_kwargs.setdefault("values_scale", "log")
          interp_kwargs.setdefault("points_scale", ("log",))
  
          self._evaluate = ScaledRegularGridInterpolator(
              points=(energy,), values=values, **interp_kwargs
          )
  
          super().__init__()
  
      @classmethod
      def read_xspec_model(cls, filename, param, **kwargs):
          """Read XSPEC table model.
  
          The input is a table containing absorbed values from a XSPEC model as a
          function of energy.
  
          TODO: Format of the file should be described and discussed in
          https://gamma-astro-data-formats.readthedocs.io/en/latest/index.html
  
          Parameters
          ----------
          filename : str
              File containing the XSPEC model
          param : float
              Model parameter value
  
          Examples
          --------
          Fill table from an EBL model (Franceschini, 2008)
  
          >>> from gammapy.modeling.models import TemplateSpectralModel
          >>> filename = '$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz'
          >>> table_model = TemplateSpectralModel.read_xspec_model(filename=filename, param=0.3)
          """
          filename = make_path(filename)
  
          # Check if parameter value is in range
          table_param = Table.read(filename, hdu="PARAMETERS")
          pmin = table_param["MINIMUM"]
          pmax = table_param["MAXIMUM"]
          if param < pmin or param > pmax:
              raise ValueError(f"Out of range: param={param}, min={pmin}, max={pmax}")
  
          # Get energy values
          table_energy = Table.read(filename, hdu="ENERGIES")
          energy_lo = table_energy["ENERG_LO"]
          energy_hi = table_energy["ENERG_HI"]
  
          # set energy to log-centers
          energy = np.sqrt(energy_lo * energy_hi)
  
          # Get spectrum values (no interpolation, take closest value for param)
          table_spectra = Table.read(filename, hdu="SPECTRA")
          idx = np.abs(table_spectra["PARAMVAL"] - param).argmin()
          values = u.Quantity(table_spectra[idx][1], "", copy=False)  # no dimension
  
          kwargs.setdefault("interp_kwargs", {"values_scale": "lin"})
          return cls(energy=energy, values=values, **kwargs)
  
      def evaluate(self, energy):
          """Evaluate the model (static function)."""
          return self._evaluate((energy,), clip=True)
  
      def to_dict(self):
          return {
              "type": self.tag[0],
              "energy": {
                  "data": self.energy.data.tolist(),
                  "unit": str(self.energy.unit),
              },
              "values": {
                  "data": self.values.data.tolist(),
                  "unit": str(self.values.unit),
              },
          }
  
      @classmethod
      def from_dict(cls, data):
          energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
          values = u.Quantity(data["values"]["data"], data["values"]["unit"])
          return cls(energy=energy, values=values)
  
  
  class ScaleSpectralModel(SpectralModel):
      """Wrapper to scale another spectral model by a norm factor.
  
      Parameters
      ----------
      model : `SpectralModel`
          Spectral model to wrap.
      norm : float
          Multiplicative norm factor for the model value.
      """
  
      tag = ["ScaleSpectralModel", "scale"]
      norm = Parameter("norm", 1, unit="")
  
      def __init__(self, model, norm=norm.quantity):
          self.model = model
          self._covariance = None
          super().__init__(norm=norm)
  
      def evaluate(self, energy, norm):
          return norm * self.model(energy)
  
      def integral(self, emin, emax, **kwargs):
          return self.norm.value * self.model.integral(emin, emax, **kwargs)
  
  
  class Absorption:
      r"""Gamma-ray absorption models.
  
      For more information see :ref:`absorption-spectral-model`.
  
      Parameters
      ----------
      energy : `~astropy.units.Quantity`
          Energy node values
      param : `~astropy.units.Quantity`
          Parameter node values
      data : `~astropy.units.Quantity`
          Model value
      filename : str
          Filename of the absorption model used for serialisation.
      interp_kwargs : dict
          Interpolation option passed to `ScaledRegularGridInterpolator`.
          By default the models are extrapolated outside the range. To prevent
          this and raise an error instead use interp_kwargs = {"extrapolate": False}
      """
  
      tag = "Absorption"
  
      def __init__(self, energy, param, data, filename=None, interp_kwargs=None):
          self.filename = filename
          # set values log centers
          self.param = param
          self.energy = energy
          self.energy = energy
          self.data = u.Quantity(data, copy=False)
  
          interp_kwargs = interp_kwargs or {}
          interp_kwargs.setdefault("points_scale", ("log", "lin"))
          interp_kwargs.setdefault("extrapolate", True)
  
          self._evaluate = ScaledRegularGridInterpolator(
              points=(self.param, self.energy), values=self.data, **interp_kwargs
          )
  
      def to_dict(self):
          if self.filename is None:
              return {
                  "type": self.tag,
                  "energy": {
                      "data": self.energy.data.tolist(),
                      "unit": str(self.energy.unit),
                  },
                  "param": {
                      "data": self.param.data.tolist(),
                      "unit": str(self.param.unit),
                  },
                  "values": {
                      "data": self.data.data.tolist(),
                      "unit": str(self.data.unit),
                  },
              }
          else:
              return {"type": self.tag, "filename": self.filename}
  
      @classmethod
      def from_dict(cls, data):
  
          if "filename" in data:
              return cls.read(data["filename"])
          else:
              energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
              param = u.Quantity(data["param"]["data"], data["param"]["unit"])
              values = u.Quantity(data["values"]["data"], data["values"]["unit"])
              return cls(energy=energy, param=param, data=values)
  
      @classmethod
      def read(cls, filename, interp_kwargs=None):
          """Build object from an XSPEC model.
  
          Todo: Format of XSPEC binary files should be referenced at https://gamma-astro-data-formats.readthedocs.io/en/latest/
  
          Parameters
          ----------
          filename : str
              File containing the model.
          interp_kwargs : dict
              Interpolation option passed to `ScaledRegularGridInterpolator`.
  
          Returns
          -------
          absorption : `Absorption`
              Absorption model.
          """
          # Create EBL data array
          filename = make_path(filename)
          table_param = Table.read(filename, hdu="PARAMETERS")
  
          # TODO: for some reason the table contain duplicated values
          param, idx = np.unique(table_param[0]["VALUE"], return_index=True)
  
          # Get energy values
          table_energy = Table.read(filename, hdu="ENERGIES")
          energy_lo = u.Quantity(
              table_energy["ENERG_LO"], "keV", copy=False
          )  # unit not stored in file
          energy_hi = u.Quantity(
              table_energy["ENERG_HI"], "keV", copy=False
          )  # unit not stored in file
          energy = np.sqrt(energy_lo * energy_hi)
  
          # Get spectrum values
          table_spectra = Table.read(filename, hdu="SPECTRA")
          data = table_spectra["INTPSPEC"].data[idx, :]
          return cls(
              energy=energy,
              param=param,
              data=data,
              filename=filename,
              interp_kwargs=interp_kwargs,
          )
  
      @classmethod
      def read_builtin(cls, name, interp_kwargs=None):
          """Read one of the built-in absorption models.
  
          Parameters
          ----------
          name : {'franceschini', 'dominguez', 'finke'}
              name of one of the available model in gammapy-data
  
          References
          ----------
          .. [1] Franceschini et al., "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity",
              `Link <https://ui.adsabs.harvard.edu/abs/2008A%26A...487..837F>`__
          .. [2] Dominguez et al., " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions"
              `Link <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__
          .. [3] Finke et al., "Modeling the Extragalactic Background Light from Stars and Dust"
              `Link <https://ui.adsabs.harvard.edu/abs/2010ApJ...712..238F>`__
  
          Returns
          -------
          absorption : `Absorption`
              Absorption model.
  
          """
          models = dict()
          models["franceschini"] = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz"
          models["dominguez"] = "$GAMMAPY_DATA/ebl/ebl_dominguez11.fits.gz"
          models["finke"] = "$GAMMAPY_DATA/ebl/frd_abs.fits.gz"
  
          return cls.read(models[name], interp_kwargs=interp_kwargs)
  
      def table_model(self, parameter):
          """Table model for a given parameter value.
  
          Parameters
          ----------
          parameter : float
              Parameter value.
  
          Returns
          -------
          template_model : `TemplateSpectralModel`
              Template spectral model.
          """
          energy = self.energy
          values = self.evaluate(energy=energy, parameter=parameter)
          return TemplateSpectralModel(
              energy=energy, values=values, interp_kwargs={"values_scale": "log"}
          )
  
      def evaluate(self, energy, parameter):
          """Evaluate model for energy and parameter value."""
          return np.clip(self._evaluate((parameter, energy)), 0, 1)
  
  
  class AbsorbedSpectralModel(SpectralModel):
      r"""Spectral model with EBL absorption.
  
      For more information see :ref:`absorbed-spectral-model`.
  
      Parameters
      ----------
      spectral_model : `SpectralModel`
          Spectral model.
      absorption : `Absorption`
          Absorption model.
      redshift : float
          Redshift of the absorption model
      alpha_norm: float
          Norm of the EBL model
      """
  
      tag = ["AbsorbedSpectralModel", "ebl-absorption"]
      alpha_norm = Parameter("alpha_norm", 1.0, frozen=True)
      redshift = Parameter("redshift", 0.1, frozen=True)
  
      def __init__(
          self, spectral_model, absorption, redshift, alpha_norm=alpha_norm.quantity,
      ):
          self.spectral_model = spectral_model
          self.absorption = absorption
  
          min_ = self.absorption.param.min()
          max_ = self.absorption.param.max()
  
          redshift = Parameter("redshift", redshift, frozen=True, min=min_, max=max_)
          super().__init__(redshift=redshift, alpha_norm=alpha_norm)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          self._covariance.set_subcovariance(self.spectral_model.covariance)
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          subcovar = self._covariance.get_subcovariance(
              self.spectral_model.covariance.parameters
          )
          self.spectral_model.covariance = subcovar
  
      @property
      def parameters(self):
          return (
              Parameters([self.redshift, self.alpha_norm])
              + self.spectral_model.parameters
          )
  
      def evaluate(self, energy, **kwargs):
          """Evaluate the model at a given energy."""
          # assign redshift value and remove it from dictionary
          # since it does not belong to the spectral model
          parameter = kwargs.pop("redshift")
          alpha_norm = kwargs.pop("alpha_norm")
  
          dnde = self.spectral_model.evaluate(energy=energy, **kwargs)
          absorption = self.absorption.evaluate(energy=energy, parameter=parameter)
          # Power rule: (e ^ a) ^ b = e ^ (a * b)
          absorption = np.power(absorption, alpha_norm)
          return dnde * absorption
  
      def to_dict(self):
          return {
              "type": self.tag,
              "base_model": self.spectral_model.to_dict(),
              "absorption": self.absorption.to_dict(),
              "absorption_parameter": {"name": "redshift", "value": self.redshift.value,},
              "parameters": Parameters([self.redshift, self.alpha_norm]).to_dict(),
          }
  
      @classmethod
      def from_dict(cls, data):
          from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY
  
          model_class = SPECTRAL_MODEL_REGISTRY.get_cls(data["base_model"]["type"])
  
          model = cls(
              spectral_model=model_class.from_dict(data["base_model"]),
              absorption=Absorption.from_dict(data["absorption"]),
              redshift=data["absorption_parameter"]["value"],
          )
          return model
  
  
  class NaimaSpectralModel(SpectralModel):
      r"""A wrapper for Naima models.
  
      For more information see :ref:`naima-spectral-model`.
  
      Parameters
      ----------
      radiative_model : `~naima.models.BaseRadiative`
          An instance of a radiative model defined in `~naima.models`
      distance : `~astropy.units.Quantity`, optional
          Distance to the source. If set to 0, the intrinsic differential
          luminosity will be returned. Default is 1 kpc
      seed : str or list of str, optional
          Seed photon field(s) to be considered for the `radiative_model` flux computation,
          in case of a `~naima.models.InverseCompton` model. It can be a subset of the
          `seed_photon_fields` list defining the `radiative_model`. Default is the whole list
          of photon fields
      nested_models : dict
          Additionnal parameters for nested models not supplied by the radiative model,
          for now this is used  only for synchrotron self-compton model
      """
  
      tag = ["NaimaSpectralModel", "naima"]
  
      def __init__(
          self, radiative_model, distance=1.0 * u.kpc, seed=None, nested_models=None
      ):
          import naima
  
          self.radiative_model = radiative_model
          self._particle_distribution = self.radiative_model.particle_distribution
          self.distance = u.Quantity(distance)
          self.seed = seed
  
          if nested_models is None:
              nested_models = {}
  
          self.nested_models = nested_models
  
          if isinstance(self._particle_distribution, naima.models.TableModel):
              param_names = ["amplitude"]
          else:
              param_names = self._particle_distribution.param_names
  
          parameters = []
  
          for name in param_names:
              value = getattr(self._particle_distribution, name)
              parameter = Parameter(name, value)
              parameters.append(parameter)
  
          # In case of a synchrotron radiative model, append B to the fittable parameters
          if "B" in self.radiative_model.param_names:
              value = getattr(self.radiative_model, "B")
              parameter = Parameter("B", value)
              parameters.append(parameter)
  
          # In case of a synchrotron self compton model, append B and Rpwn to the fittable parameters
          if (
              isinstance(self.radiative_model, naima.models.InverseCompton)
              and "SSC" in self.nested_models
          ):
              B = self.nested_models["SSC"]["B"]
              radius = self.nested_models["SSC"]["radius"]
              parameters.append(Parameter("B", B))
              parameters.append(Parameter("radius", radius, frozen=True))
  
          self.default_parameters = Parameters(parameters)
          super().__init__()
  
      def _evaluate_ssc(
          self, energy,
      ):
          """
          Compute photon density spectrum from synchrotron emission for synchrotron self-compton model,
          assuming uniform synchrotron emissivity inside a sphere of radius R
          (see Section 4.1 of Atoyan & Aharonian 1996)
  
          based on :
          "https://naima.readthedocs.io/en/latest/examples.html#crab-nebula-ssc-model"
  
          """
          import naima
  
          SYN = naima.models.Synchrotron(
              self._particle_distribution,
              B=self.B.quantity,
              Eemax=self.radiative_model.Eemax,
              Eemin=self.radiative_model.Eemin,
          )
  
          Esy = np.logspace(-7, 9, 100) * u.eV
          Lsy = SYN.flux(Esy, distance=0 * u.cm)  # use distance 0 to get luminosity
          phn_sy = Lsy / (4 * np.pi * self.radius.quantity ** 2 * const.c) * 2.24
          # The factor 2.24 comes from the assumption on uniform synchrotron
          # emissivity inside a sphere
  
          if "SSC" not in self.radiative_model.seed_photon_fields:
              self.radiative_model.seed_photon_fields["SSC"] = {
                  "isotropic": True,
                  "type": "array",
                  "energy": Esy,
                  "photon_density": phn_sy,
              }
          else:
              self.radiative_model.seed_photon_fields["SSC"]["photon_density"] = phn_sy
  
          dnde = self.radiative_model.flux(
              energy, seed=self.seed, distance=self.distance
          ) + SYN.flux(energy, distance=self.distance)
          return dnde
  
      def evaluate(self, energy, **kwargs):
          """Evaluate the model."""
          import naima
  
          for name, value in kwargs.items():
              setattr(self._particle_distribution, name, value)
  
          if "B" in self.radiative_model.param_names:
              self.radiative_model.B = self.B.quantity
  
          if (
              isinstance(self.radiative_model, naima.models.InverseCompton)
              and "SSC" in self.nested_models
          ):
              dnde = self._evaluate_ssc(energy.flatten())
          elif self.seed is not None:
              dnde = self.radiative_model.flux(
                  energy.flatten(), seed=self.seed, distance=self.distance
              )
          else:
              dnde = self.radiative_model.flux(energy.flatten(), distance=self.distance)
  
          dnde = dnde.reshape(energy.shape)
          unit = 1 / (energy.unit * u.cm ** 2 * u.s)
          return dnde.to(unit)
  
      @classmethod
      def from_dict(cls, data):
          raise NotImplementedError(
              "Currently the NaimaSpectralModel cannot be read from YAML"
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          raise NotImplementedError(
              "Currently the NaimaSpectralModel cannot be built from a list of parameters."
          )
  
  
  class GaussianSpectralModel(SpectralModel):
      r"""Gaussian spectral model.
  
      For more information see :ref:`gaussian-spectral-model`.
  
      Parameters
      ----------
      norm : `~astropy.units.Quantity`
          :math:`N_0`
      mean : `~astropy.units.Quantity`
          :math:`\bar{E}`
      sigma : `~astropy.units.Quantity`
          :math:`\sigma`
      """
  
      tag = ["GaussianSpectralModel", "gauss"]
      norm = Parameter("norm", 1e-12 * u.Unit("cm-2 s-1"))
      mean = Parameter("mean", 1 * u.TeV)
      sigma = Parameter("sigma", 2 * u.TeV)
  
      @staticmethod
      def evaluate(energy, norm, mean, sigma):
          return (
              norm
              / (sigma * np.sqrt(2 * np.pi))
              * np.exp(-((energy - mean) ** 2) / (2 * sigma ** 2))
          )
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate Gaussian analytically.
  
          .. math::
              F(E_{min}, E_{max}) = \frac{N_0}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma})\right]_{E_{min}}^{E_{max}}
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range
          """
          # kwargs are passed to this function but not used
          # this is to get a consistent API with SpectralModel.integral()
          u_min = (
              (emin - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          u_max = (
              (emax - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
  
          return (
              self.norm.quantity
              / 2
              * (scipy.special.erf(u_max) - scipy.special.erf(u_min))
          )
  
      def energy_flux(self, emin, emax):
          r"""Compute energy flux in given energy range analytically.
  
          .. math::
              G(E_{min}, E_{max}) =  \frac{N_0 \sigma}{\sqrt{2*\pi}}* \left[ - \exp(\frac{E_{min}-\bar{E}}{\sqrt{2} \sigma})
              \right]_{E_{min}}^{E_{max}} + \frac{N_0 * \bar{E}}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma})
               \right]_{E_{min}}^{E_{max}}
  
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          """
          u_min = (
              (emin - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          u_max = (
              (emax - self.mean.quantity) / (np.sqrt(2) * self.sigma.quantity)
          ).to_value("")
          a = self.norm.quantity * self.sigma.quantity / np.sqrt(2 * np.pi)
          b = self.norm.quantity * self.mean.quantity / 2
          return a * (np.exp(-(u_min ** 2)) - np.exp(-(u_max ** 2))) + b * (
              scipy.special.erf(u_max) - scipy.special.erf(u_min)
          )
  
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", 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("LogParabolaNormSpectralModel", "spectral")
      yield Model.create(
          "TemplateSpectralModel", "spectral", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create(
          "PiecewiseBrokenPowerLawSpectralModel", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )
      yield Model.create("GaussianSpectralModel")
      # TODO: yield Model.create("AbsorbedSpectralModel")
      # TODO: yield Model.create("NaimaSpectralModel")
      # TODO: yield Model.create("ScaleSpectralModel")
      yield Model.create("ConstantTemporalModel", "temporal")
      yield Model.create("LightCurveTemplateTemporalModel", "temporal", Table())
      yield Model.create(
          "SkyModel",
          spatial_model=Model.create("ConstantSpatialModel", "spatial"),
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
      )
      m1 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_nodes([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("TemplateSpatialModel", "spatial", map=m1)
      m2 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_edges([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("BackgroundModel", map=m2)
  
  
  @pytest.mark.parametrize("model_class", MODEL_REGISTRY)
  def test_all_model_classes(model_class):
      if isinstance(model_class.tag, list):
          assert model_class.tag[0] == model_class.__name__
      else:
          assert model_class.tag == model_class.__name__
  
  
  @pytest.mark.parametrize("model", make_all_models())
  def test_all_model_instances(model):
      tag = model.tag[0] if isinstance(model.tag, list) else model.tag
      assert tag == model.__class__.__name__
  
  
  @requires_data()
  def test_missing_parameters():
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      assert models["source1"].spatial_model.e in models.parameters
      assert len(models["source1"].spatial_model.parameters) == 6
  
  
  def test_registries_print():
      print(MODEL_REGISTRY)
  
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("LogParabolaNormSpectralModel", "spectral")
      yield Model.create(
          "TemplateSpectralModel", "spectral", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create("GaussianSpectralModel", "spectral")
      # TODO: yield Model.create("AbsorbedSpectralModel")
      # TODO: yield Model.create("NaimaSpectralModel")
      # TODO: yield Model.create("ScaleSpectralModel")
      yield Model.create("ConstantTemporalModel", "temporal")
      yield Model.create("LightCurveTemplateTemporalModel", "temporal", Table())
      yield Model.create(
          "SkyModel",
          spatial_model=Model.create("ConstantSpatialModel", "spatial"),
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
      )
      m1 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_nodes([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("TemplateSpatialModel", "spatial", map=m1)
      m2 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_edges([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("BackgroundModel", map=m2)
  
  
  @pytest.mark.parametrize("model_class", MODEL_REGISTRY)
  def test_all_model_classes(model_class):
      if isinstance(model_class.tag, list):
          assert model_class.tag[0] == model_class.__name__
      else:
          assert model_class.tag == model_class.__name__
  
  
  @pytest.mark.parametrize("model", make_all_models())
  def test_all_model_instances(model):
      tag = model.tag[0] if isinstance(model.tag, list) else model.tag
      assert tag == model.__class__.__name__
  
  
  @requires_data()
  def test_missing_parameters():
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      assert models["source1"].spatial_model.e in models.parameters
      assert len(models["source1"].spatial_model.parameters) == 6
  
  
  def test_registries_print():
      print(MODEL_REGISTRY)
  
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_log10(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=1.151292546497023 * 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="constant",
          model=ConstantSpectralModel(const=4 / u.cm ** 2 / u.s / u.TeV),
          val_at_2TeV=u.Quantity(4, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(35.9999999999999, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(198.00000000000006, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw_index1",
          model=PowerLawSpectralModel(
              index=1 * u.Unit(""),
              amplitude=2 / 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(4.605170185, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(18.0, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl_2",
          model=ExpCutoffPowerLawSpectralModel(
              index=2.0 * 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(0.81873075, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.83075297, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.41406327, "TeV cm-2 s-1"),
          e_peak=np.nan * u.TeV,
      ),
      dict(
          name="GaussianSpectralModel",
          model=GaussianSpectralModel(
              norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2 * u.TeV
          ),
          val_at_2TeV=u.Quantity(7.978845608028654, "cm-2 s-1 TeV-1"),
          val_at_3TeV=u.Quantity(2.973439029468601e-05, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.9999988533937123, "cm-2 s-1"),
          integral_infinity=u.Quantity(4, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.999998896163037, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.8 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
              alpha=0.8,
          ),
          val_at_2TeV=u.Quantity(0.871694294554192, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.026342, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.38652453, "TeV cm-2 s-1"),
          e_peak=1.7677669529663684 * u.TeV,
      ),
      dict(
          name="bpl",
          model=BrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.125, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.45649740094103286, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(0.9669999668731384, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
              reference=1 * u.TeV,
              beta=1,
          ),
          val_at_2TeV=u.Quantity(0.28284271247461906, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.9956923907948155, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(2.2372256145972207, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl-hard",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=2.5 * u.Unit(""),
              index2=1.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
              reference=1 * u.TeV,
              beta=1,
          ),
          val_at_2TeV=u.Quantity(3.5355339059327378, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(13.522782989735022, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(40.06681812966845, "TeV cm-2 s-1"),
      ),
      dict(
          name="pbpl",
          model=PiecewiseBrokenPowerLawSpectralModel(
              energy=[1, 3, 7, 10] * u.TeV,
              values=[1, 5, 3, 0.5] * u.Unit("cm-2 s-1 TeV-1"),
          ),
          val_at_2TeV=u.Quantity(2.76058404, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(24.757573885411876, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(117.74087966682515, "TeV cm-2 s-1"),
      ),
  ]
  
  # Add compound models
  
  TEST_MODELS.append(
      dict(
          name="compound3",
          model=TEST_MODELS[0]["model"] + TEST_MODELS[0]["model"],
          val_at_2TeV=TEST_MODELS[0]["val_at_2TeV"] * 2,
          integral_1_10TeV=TEST_MODELS[0]["integral_1_10TeV"] * 2,
          eflux_1_10TeV=TEST_MODELS[0]["eflux_1_10TeV"] * 2,
      )
  )
  
  TEST_MODELS.append(
      dict(
          name="compound6",
          model=TEST_MODELS[11]["model"] + u.Quantity(4, "cm-2 s-1 TeV-1"),
          val_at_2TeV=TEST_MODELS[11]["val_at_2TeV"] * 2,
          integral_1_10TeV=TEST_MODELS[11]["integral_1_10TeV"] * 2,
          eflux_1_10TeV=TEST_MODELS[11]["eflux_1_10TeV"] * 2,
      )
  )
  
  TEST_MODELS.append(
      dict(
          name="table_model",
          model=table_model(),
          # Values took from power law expectation
          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"),
      )
  )
  
  
  @requires_dependency("scipy")
  @pytest.mark.parametrize("spectrum", TEST_MODELS, ids=lambda _: _["name"])
  def test_models(spectrum):
      model = spectrum["model"]
      energy = 2 * u.TeV
      value = model(energy)
      assert_quantity_allclose(value, spectrum["val_at_2TeV"], rtol=1e-7)
      if "val_at_3TeV" in spectrum:
          energy = 3 * u.TeV
          value = model(energy)
          assert_quantity_allclose(value, spectrum["val_at_3TeV"], rtol=1e-7)
  
      emin = 1 * u.TeV
      emax = 10 * u.TeV
      assert_quantity_allclose(
          model.integral(emin=emin, emax=emax), spectrum["integral_1_10TeV"], rtol=1e-5
      )
      assert_quantity_allclose(
          model.energy_flux(emin=emin, emax=emax), spectrum["eflux_1_10TeV"], rtol=1e-5
      )
  
      if "e_peak" in spectrum:
          assert_quantity_allclose(model.e_peak, spectrum["e_peak"], rtol=1e-2)
  
      # inverse for ConstantSpectralModel is irrelevant.
      # inverse for Gaussian and PiecewiseBrokenPowerLawSpectralModel have a degeneracy
      if not (
          isinstance(model, ConstantSpectralModel)
          or spectrum["name"] == "compound6"
          or spectrum["name"] == "GaussianSpectralModel"
          or spectrum["name"] == "pbpl"
      ):
          assert_quantity_allclose(model.inverse(value), 2 * u.TeV, rtol=0.01)
  
      if "integral_infinity" in spectrum:
          emin = 0 * u.TeV
          emax = 10000 * u.TeV
          assert_quantity_allclose(
              model.integral(emin=emin, emax=emax), spectrum["integral_infinity"]
          )
  
      model.to_dict()
  
      assert "" in str(model)
  
      # check that an array evaluation works (otherwise e.g. plotting raises an error)
      e_array = [2, 10, 20] * u.TeV
      e_array = e_array[:, np.newaxis, np.newaxis]
      val = model(e_array)
      assert val.shape == e_array.shape
      assert_quantity_allclose(val[0], spectrum["val_at_2TeV"])
  
  
  def test_model_unit():
      pwl = PowerLawSpectralModel()
      value = pwl(500 * u.MeV)
      assert value.unit == "cm-2 s-1 TeV-1"
  
  
  @requires_dependency("matplotlib")
  def test_model_plot():
      pwl = PowerLawSpectralModel(
          amplitude=1e-12 * u.Unit("TeV-1 cm-2 s-1"), reference=1 * u.Unit("TeV"), index=2
      )
      pwl.amplitude.error = 0.1e-12 * u.Unit("TeV-1 cm-2 s-1")
  
      with mpl_plot_check():
          pwl.plot((1 * u.TeV, 10 * u.TeV))
  
      with mpl_plot_check():
          pwl.plot_error((1 * u.TeV, 10 * u.TeV))
  
  
  def test_to_from_dict():
      spectrum = TEST_MODELS[0]
      model = spectrum["model"]
  
      model_dict = model.to_dict()
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
  
      assert isinstance(new_model, PowerLawSpectralModel)
  
      actual = [par.value for par in new_model.parameters]
      desired = [par.value for par in model.parameters]
      assert_quantity_allclose(actual, desired)
  
  
  def test_to_from_dict_compound():
      spectrum = TEST_MODELS[-2]
      model = spectrum["model"]
      assert spectrum["name"] == "compound6"
      model_dict = model.to_dict()
      assert model_dict["operator"] == "add"
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
  
      assert isinstance(new_model, CompoundSpectralModel)
  
      actual = [par.value for par in new_model.parameters]
      desired = [par.value for par in model.parameters]
      assert_quantity_allclose(actual, desired)
  
  
  @requires_dependency("matplotlib")
  @requires_data()
  def test_table_model_from_file():
      filename = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz"
      absorption_z03 = TemplateSpectralModel.read_xspec_model(
          filename=filename, param=0.3
      )
      with mpl_plot_check():
          absorption_z03.plot(energy_range=(0.03, 10), energy_unit=u.TeV, flux_unit="")
  
  
  @requires_data()
  def test_absorption():
      # absorption values for given redshift
      redshift = 0.117
      absorption = Absorption.read_builtin("dominguez")
  
      # Spectral model corresponding to PKS 2155-304 (quiescent state)
      index = 3.53
      amplitude = 1.81 * 1e-12 * u.Unit("cm-2 s-1 TeV-1")
      reference = 1 * u.TeV
      pwl = PowerLawSpectralModel(index=index, amplitude=amplitude, reference=reference)
  
      # EBL + PWL model
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, redshift=redshift
      )
  
      desired = u.Quantity(5.140765e-13, "TeV-1 s-1 cm-2")
      assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
      assert model.alpha_norm.value == 1.0
  
      # EBL + PWL model: test if norm of EBL=0: it mean model =pwl
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, alpha_norm=0, redshift=redshift
      )
      assert_quantity_allclose(model(1 * u.TeV), pwl(1 * u.TeV), rtol=1e-3)
  
      # EBL + PWL model: Test with a norm different of 1
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, alpha_norm=1.5, redshift=redshift
      )
      desired = u.Quantity(2.739695e-13, "TeV-1 s-1 cm-2")
      assert model.alpha_norm.value == 1.5
      assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
  
      # Test error propagation
      model.spectral_model.amplitude.error = 0.1 * model.spectral_model.amplitude.value
      dnde, dnde_err = model.evaluate_error(1 * u.TeV)
      assert_allclose(dnde_err / dnde, 0.1)
  
  
  @requires_data()
  def test_absorbed_extrapolate():
      ebl_model = "dominguez"
      z = 0.001
      absorption = Absorption.read_builtin(ebl_model)
  
      model = absorption.table_model(z)
      assert_allclose(model(1 * u.TeV), 1)
  
  
  def test_ecpl_integrate():
      # regression test to check the numerical integration for small energy bins
      ecpl = ExpCutoffPowerLawSpectralModel()
      value = ecpl.integral(1 * u.TeV, 1.1 * u.TeV)
      assert_quantity_allclose(value, 8.380761e-14 * u.Unit("s-1 cm-2"))
  
  
  def test_pwl_pivot_energy():
      pwl = PowerLawSpectralModel(amplitude="5.35510540e-11 cm-2 s-1 TeV-1")
  
      pwl.covariance = [
          [0.0318377 ** 2, 6.56889442e-14, 0],
          [6.56889442e-14, 0, 0],
          [0, 0, 0],
      ]
  
      assert_quantity_allclose(pwl.pivot_energy, 3.3540034240210987 * u.TeV)
  
  
  def test_TemplateSpectralModel_evaluate_tiny():
      energy = np.array([1.00000000e06, 1.25892541e06, 1.58489319e06, 1.99526231e06])
      values = np.array([4.39150790e-38, 1.96639562e-38, 8.80497507e-39, 3.94262401e-39])
  
      model = TemplateSpectralModel(
          energy=energy, values=values * u.Unit("MeV-1 s-1 sr-1")
      )
      result = model.evaluate(energy)
      tiny = np.finfo(np.float32).tiny
      mask = abs(values) - tiny > tiny
      np.testing.assert_allclose(
          values[mask] / values.max(), result[mask].value / values.max()
      )
      mask = abs(result.value) - tiny <= tiny
      assert np.all(result[mask] == 0.0)
  
  
  def test_TemplateSpectralModel_compound():
      energy = [1.00e06, 1.25e06, 1.58e06, 1.99e06] * u.MeV
      values = [4.39e-7, 1.96e-7, 8.80e-7, 3.94e-7] * u.Unit("MeV-1 s-1 sr-1")
  
      template = TemplateSpectralModel(energy=energy, values=values)
      correction = PowerLawNormSpectralModel(norm=2)
      model = CompoundSpectralModel(template, correction, operator=operator.mul)
      assert np.allclose(model(energy), 2 * values)
  
      model_mul = template * correction
      assert isinstance(model_mul, CompoundSpectralModel)
      assert np.allclose(model_mul(energy), 2 * values)
  
      model_dict = model.to_dict()
      assert model_dict["operator"] == "mul"
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
      assert isinstance(new_model, CompoundSpectralModel)
      assert np.allclose(new_model(energy), 2 * values)
  
  
  @requires_dependency("naima")
  class TestNaimaModel:
      # Used to test model value at 2 TeV
      energy = 2 * u.TeV
  
      # Used to test model integral and energy flux
      emin = 1 * u.TeV
      emax = 10 * u.TeV
  
      # Used to that if array evaluation works
      e_array = [2, 10, 20] * u.TeV
      e_array = e_array[:, np.newaxis, np.newaxis]
  
      def test_pion_decay(self):
          import naima
  
          particle_distribution = naima.models.PowerLaw(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=2.5
          )
          radiative_model = naima.radiative.PionDecay(
              particle_distribution, nh=1 * u.cm ** -3
          )
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 9.725347355450884e-14 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 3.530537143620737e-13 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 7.643559573105779e-13 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
          model.amplitude.error = 0.1 * model.amplitude.value
  
          out = model.evaluate_error(1 * u.TeV)
          assert_allclose(out.data, [5.266068e-13, 5.266068e-14], rtol=1e-3)
  
      def test_ic(self):
          import naima
  
          particle_distribution = naima.models.ExponentialCutoffBrokenPowerLaw(
              amplitude=2e33 / u.eV,
              e_0=10 * u.TeV,
              alpha_1=2.5,
              alpha_2=2.7,
              e_break=900 * u.GeV,
              e_cutoff=10 * u.TeV,
          )
          radiative_model = naima.radiative.InverseCompton(
              particle_distribution, seed_photon_fields=["CMB"]
          )
  
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 4.347836316893546e-12 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 1.595813e-11 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 2.851283e-11 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV, rtol=1e-5
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV, rtol=1e-5
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
      def test_synchrotron(self):
          import naima
  
          particle_distribution = naima.models.LogParabola(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=1.3, beta=0.5
          )
          radiative_model = naima.radiative.Synchrotron(particle_distribution, B=2 * u.G)
  
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 1.0565840392550432e-24 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 4.449303e-13 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 4.594242e-13 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV, rtol=1e-5
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV, rtol=1e-5
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
          model.B.value = 3  # update B
          val_at_2TeV = 5.1985064062296e-16 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
  
      def test_ssc(self):
          import naima
  
          ECBPL = naima.models.ExponentialCutoffBrokenPowerLaw(
              amplitude=3.699e36 / u.eV,
              e_0=1 * u.TeV,
              e_break=0.265 * u.TeV,
              alpha_1=1.5,
              alpha_2=3.233,
              e_cutoff=1863 * u.TeV,
              beta=2.0,
          )
  
          radiative_model = naima.radiative.InverseCompton(
              ECBPL,
              seed_photon_fields=[
                  "CMB",
                  ["FIR", 70 * u.K, 0.5 * u.eV / u.cm ** 3],
                  ["NIR", 5000 * u.K, 1 * u.eV / u.cm ** 3],
              ],
              Eemax=50 * u.PeV,
              Eemin=0.1 * u.GeV,
          )
          B = 125 * u.uG
          radius = 2.1 * u.pc
          nested_models = {"SSC": {"B": B, "radius": radius}}
          model = NaimaSpectralModel(radiative_model, nested_models=nested_models)
          assert_quantity_allclose(model.B.quantity, B)
          assert_quantity_allclose(model.radius.quantity, radius)
          val_at_2TeV = 1.6703761561806372e-11 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV, rtol=1e-5)
  
          model.parameters["B"].value = 100
          val_at_2TeV = 1.441331153167876e-11 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV, rtol=1e-5)
  
      def test_bad_init(self):
          import naima
  
          particle_distribution = naima.models.PowerLaw(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=2.5
          )
          radiative_model = naima.radiative.PionDecay(
              particle_distribution, nh=1 * u.cm ** -3
          )
          model = NaimaSpectralModel(radiative_model)
  
          with pytest.raises(NotImplementedError):
              NaimaSpectralModel.from_dict(model.to_dict())
          with pytest.raises(NotImplementedError):
              NaimaSpectralModel.from_parameters(model.parameters)
  
  
  class TestSpectralModelErrorPropagation:
      """Test spectral model error propagation.
  
      https://github.com/gammapy/gammapy/blob/master/docs/development/pigs/pig-014.rst#proposal
      https://nbviewer.jupyter.org/github/gammapy/gammapy-extra/blob/master/experiments/uncertainty_estimation_prototype.ipynb
      """
  
      def setup(self):
          self.model = LogParabolaSpectralModel(
              amplitude=3.76e-11 * u.Unit("cm-2 s-1 TeV-1"),
              reference=1 * u.TeV,
              alpha=2.44,
              beta=0.25,
          )
          self.model.covariance = [
              [1.31e-23, 0, -6.80e-14, 3.04e-13],
              [0, 0, 0, 0],
              [-6.80e-14, 0, 0.00899, 0.00904],
              [3.04e-13, 0, 0.00904, 0.0284],
          ]
  
      def test_evaluate_error_scalar(self):
          # evaluate_error on scalar
          out = self.model.evaluate_error(1 * u.TeV)
          assert isinstance(out, u.Quantity)
          assert out.unit == "cm-2 s-1 TeV-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [3.7600e-11, 3.6193e-12], rtol=1e-3)
  
      def test_evaluate_error_array(self):
          out = self.model.evaluate_error([1, 100] * u.TeV)
          assert out.shape == (2, 2)
          expected = [[3.76e-11, 2.469e-18], [3.619e-12, 9.375e-18]]
          assert_allclose(out.data, expected, rtol=1e-3)
  
      def test_evaluate_error_unit(self):
          out = self.model.evaluate_error(1e6 * u.MeV)
          assert out.unit == "cm-2 s-1 TeV-1"
          assert_allclose(out.data, [3.760e-11, 3.6193e-12], rtol=1e-3)
  
      @pytest.mark.xfail(reason="FIXME, do we need this method?")
      def test_integral_error(self):
          out = self.model.integral_error(1 * u.TeV, 10 * u.TeV)
          assert out.unit == "cm-2 s-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [2.197e-11, 2.796e-12], rtol=1e-3)
  
      @pytest.mark.xfail(reason="FIXME, do we need this method?")
      def test_energy_flux_error(self):
          out = self.model.energy_flux_error(1 * u.TeV, 10 * u.TeV)
          assert out.unit == "TeV cm-2 s-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [4.119e-11, 8.157e-12], rtol=1e-3)
  
      def test_ecpl_model(self):
          # Regression test for ECPL model
          # https://github.com/gammapy/gammapy/issues/2007
          model = ExpCutoffPowerLawSpectralModel(
              amplitude=2.076183759227292e-12 * u.Unit("cm-2 s-1 TeV-1"),
              index=1.8763343736076483,
              lambda_=0.08703226432146616 * u.Unit("TeV-1"),
              reference=1 * u.TeV,
          )
          model.covariance = [
              [0.00204191498, -1.507724e-14, 0.0, -0.001834819, 0.0],
              [-1.507724e-14, 1.6864740e-25, 0.0, 1.854251e-14, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0],
              [-0.001834819175, 1.8542517e-14, 0.0, 0.0032559101, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0],
          ]
  
          out = model.evaluate_error(1 * u.TeV)
          assert_allclose(out.data, [1.903129e-12, 2.979976e-13], rtol=1e-3)
  
          out = model.evaluate_error(0.1 * u.TeV)
          assert_allclose(out.data, [1.548176e-10, 1.933612e-11], rtol=1e-3)
  
  
  def test_integral_error_PowerLaw():
      energy = np.linspace(1 * u.TeV, 10 * u.TeV, 10)
      emin = energy[:-1]
      emax = energy[1:]
  
      powerlaw = PowerLawSpectralModel()
      powerlaw.parameters['index'].error = 0.4
      powerlaw.parameters['amplitude'].error = 1e-13
  
      flux, flux_error = powerlaw.integral_error(emin,emax)
  
      assert_allclose(flux.value[0]/1e-13, 5.0, rtol=0.1)
      assert_allclose(flux_error.value[0]/1e-14, 8.546615432273905, rtol=0.01)
  
  
  def test_integral_error_ExpCutOffPowerLaw():
      energy = np.linspace(1 * u.TeV, 10 * u.TeV, 10)
      emin = energy[:-1]
      emax = energy[1:]
      
      exppowerlaw = ExpCutoffPowerLawSpectralModel()
      exppowerlaw.parameters['index'].error = 0.4
      exppowerlaw.parameters['amplitude'].error = 1e-13
      exppowerlaw.parameters['lambda_'].error = 0.03
      
      flux, flux_error = exppowerlaw.integral_error(emin, emax)
      
      assert_allclose(flux.value[0]/1e-13, 5.05855622, rtol=0.01)
      assert_allclose(flux_error.value[0]/1e-14, 8.90907063, rtol=0.01)
  
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.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=1.151292546497023 * 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="constant",
          model=ConstantSpectralModel(const=4 / u.cm ** 2 / u.s / u.TeV),
          val_at_2TeV=u.Quantity(4, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(35.9999999999999, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(198.00000000000006, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw_index1",
          model=PowerLawSpectralModel(
              index=1 * u.Unit(""),
              amplitude=2 / 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(4.605170185, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(18.0, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl_2",
          model=ExpCutoffPowerLawSpectralModel(
              index=2.0 * 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(0.81873075, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.83075297, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.41406327, "TeV cm-2 s-1"),
          e_peak=np.nan * u.TeV,
      ),
      dict(
          name="GaussianSpectralModel",
          model=GaussianSpectralModel(
              norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2 * u.TeV
          ),
          val_at_2TeV=u.Quantity(7.978845608028654, "cm-2 s-1 TeV-1"),
          val_at_3TeV=u.Quantity(2.973439029468601e-05, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.9999988533937123, "cm-2 s-1"),
          integral_infinity=u.Quantity(4, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.999998896163037, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.8 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
              alpha=0.8,
          ),
          val_at_2TeV=u.Quantity(0.871694294554192, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.026342, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.38652453, "TeV cm-2 s-1"),
          e_peak=1.7677669529663684 * u.TeV,
      ),
      dict(
          name="bpl",
          model=BrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.125, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.45649740094103286, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(0.9669999668731384, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
              reference=1 * u.TeV,
              beta=1,
          ),
          val_at_2TeV=u.Quantity(0.28284271247461906, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.9956923907948155, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(2.2372256145972207, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl-hard",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=2.5 * u.Unit(""),
              index2=1.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
              reference=1 * u.TeV,
              beta=1,
          ),
          val_at_2TeV=u.Quantity(3.5355339059327378, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(13.522782989735022, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(40.06681812966845, "TeV cm-2 s-1"),
      ),
  ]
  
  # Add compound models
  
  TEST_MODELS.append(
      dict(
          name="compound3",
          model=TEST_MODELS[0]["model"] + TEST_MODELS[0]["model"],
          val_at_2TeV=TEST_MODELS[0]["val_at_2TeV"] * 2,
          integral_1_10TeV=TEST_MODELS[0]["integral_1_10TeV"] * 2,
          eflux_1_10TeV=TEST_MODELS[0]["eflux_1_10TeV"] * 2,
      )
  )
  
  TEST_MODELS.append(
      dict(
          name="compound6",
          model=TEST_MODELS[11]["model"] + u.Quantity(4, "cm-2 s-1 TeV-1"),
          val_at_2TeV=TEST_MODELS[11]["val_at_2TeV"] * 2,
          integral_1_10TeV=TEST_MODELS[11]["integral_1_10TeV"] * 2,
          eflux_1_10TeV=TEST_MODELS[11]["eflux_1_10TeV"] * 2,
      )
  )
  
  TEST_MODELS.append(
      dict(
          name="table_model",
          model=table_model(),
          # Values took from power law expectation
          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"),
      )
  )
  
  
  @requires_dependency("scipy")
  @pytest.mark.parametrize("spectrum", TEST_MODELS, ids=lambda _: _["name"])
  def test_models(spectrum):
      model = spectrum["model"]
      energy = 2 * u.TeV
      value = model(energy)
      assert_quantity_allclose(value, spectrum["val_at_2TeV"], rtol=1e-7)
      if "val_at_3TeV" in spectrum:
          energy = 3 * u.TeV
          value = model(energy)
          assert_quantity_allclose(value, spectrum["val_at_3TeV"], rtol=1e-7)
  
      emin = 1 * u.TeV
      emax = 10 * u.TeV
      assert_quantity_allclose(
          model.integral(emin=emin, emax=emax), spectrum["integral_1_10TeV"], rtol=1e-5
      )
      assert_quantity_allclose(
          model.energy_flux(emin=emin, emax=emax), spectrum["eflux_1_10TeV"], rtol=1e-5
      )
  
      if "e_peak" in spectrum:
          assert_quantity_allclose(model.e_peak, spectrum["e_peak"], rtol=1e-2)
  
      # inverse for ConstantSpectralModel is irrelevant.
      # inverse for Gaussian has a degeneracy
      if not (
          isinstance(model, ConstantSpectralModel)
          or spectrum["name"] == "compound6"
          or spectrum["name"] == "GaussianSpectralModel"
      ):
          assert_quantity_allclose(model.inverse(value), 2 * u.TeV, rtol=0.01)
  
      if "integral_infinity" in spectrum:
          emin = 0 * u.TeV
          emax = 10000 * u.TeV
          assert_quantity_allclose(
              model.integral(emin=emin, emax=emax), spectrum["integral_infinity"]
          )
  
      model.to_dict()
  
      assert "" in str(model)
  
      # check that an array evaluation works (otherwise e.g. plotting raises an error)
      e_array = [2, 10, 20] * u.TeV
      e_array = e_array[:, np.newaxis, np.newaxis]
      val = model(e_array)
      assert val.shape == e_array.shape
      assert_quantity_allclose(val[0], spectrum["val_at_2TeV"])
  
  
  def test_model_unit():
      pwl = PowerLawSpectralModel()
      value = pwl(500 * u.MeV)
      assert value.unit == "cm-2 s-1 TeV-1"
  
  
  @requires_dependency("matplotlib")
  def test_model_plot():
      pwl = PowerLawSpectralModel(
          amplitude=1e-12 * u.Unit("TeV-1 cm-2 s-1"), reference=1 * u.Unit("TeV"), index=2
      )
      pwl.amplitude.error = 0.1e-12 * u.Unit("TeV-1 cm-2 s-1")
  
      with mpl_plot_check():
          pwl.plot((1 * u.TeV, 10 * u.TeV))
  
      with mpl_plot_check():
          pwl.plot_error((1 * u.TeV, 10 * u.TeV))
  
  
  def test_to_from_dict():
      spectrum = TEST_MODELS[0]
      model = spectrum["model"]
  
      model_dict = model.to_dict()
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
  
      assert isinstance(new_model, PowerLawSpectralModel)
  
      actual = [par.value for par in new_model.parameters]
      desired = [par.value for par in model.parameters]
      assert_quantity_allclose(actual, desired)
  
  
  def test_to_from_dict_compound():
      spectrum = TEST_MODELS[-2]
      model = spectrum["model"]
      assert spectrum["name"] == "compound6"
      model_dict = model.to_dict()
      assert model_dict["operator"] == "add"
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
  
      assert isinstance(new_model, CompoundSpectralModel)
  
      actual = [par.value for par in new_model.parameters]
      desired = [par.value for par in model.parameters]
      assert_quantity_allclose(actual, desired)
  
  
  @requires_dependency("matplotlib")
  @requires_data()
  def test_table_model_from_file():
      filename = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz"
      absorption_z03 = TemplateSpectralModel.read_xspec_model(
          filename=filename, param=0.3
      )
      with mpl_plot_check():
          absorption_z03.plot(energy_range=(0.03, 10), energy_unit=u.TeV, flux_unit="")
  
  
  @requires_data()
  def test_absorption():
      # absorption values for given redshift
      redshift = 0.117
      absorption = Absorption.read_builtin("dominguez")
  
      # Spectral model corresponding to PKS 2155-304 (quiescent state)
      index = 3.53
      amplitude = 1.81 * 1e-12 * u.Unit("cm-2 s-1 TeV-1")
      reference = 1 * u.TeV
      pwl = PowerLawSpectralModel(index=index, amplitude=amplitude, reference=reference)
  
      # EBL + PWL model
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, redshift=redshift
      )
  
      desired = u.Quantity(5.140765e-13, "TeV-1 s-1 cm-2")
      assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
      assert model.alpha_norm.value == 1.0
  
      # EBL + PWL model: test if norm of EBL=0: it mean model =pwl
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, alpha_norm=0, redshift=redshift
      )
      assert_quantity_allclose(model(1 * u.TeV), pwl(1 * u.TeV), rtol=1e-3)
  
      # EBL + PWL model: Test with a norm different of 1
      model = AbsorbedSpectralModel(
          spectral_model=pwl, absorption=absorption, alpha_norm=1.5, redshift=redshift
      )
      desired = u.Quantity(2.739695e-13, "TeV-1 s-1 cm-2")
      assert model.alpha_norm.value == 1.5
      assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
  
      # Test error propagation
      model.spectral_model.amplitude.error = 0.1 * model.spectral_model.amplitude.value
      dnde, dnde_err = model.evaluate_error(1 * u.TeV)
      assert_allclose(dnde_err / dnde, 0.1)
  
  
  @requires_data()
  def test_absorbed_extrapolate():
      ebl_model = "dominguez"
      z = 0.001
      absorption = Absorption.read_builtin(ebl_model)
  
      model = absorption.table_model(z)
      assert_allclose(model(1 * u.TeV), 1)
  
  
  def test_ecpl_integrate():
      # regression test to check the numerical integration for small energy bins
      ecpl = ExpCutoffPowerLawSpectralModel()
      value = ecpl.integral(1 * u.TeV, 1.1 * u.TeV)
      assert_quantity_allclose(value, 8.380761e-14 * u.Unit("s-1 cm-2"))
  
  
  def test_pwl_pivot_energy():
      pwl = PowerLawSpectralModel(amplitude="5.35510540e-11 cm-2 s-1 TeV-1")
  
      pwl.covariance = [
          [0.0318377 ** 2, 6.56889442e-14, 0],
          [6.56889442e-14, 0, 0],
          [0, 0, 0],
      ]
  
      assert_quantity_allclose(pwl.pivot_energy, 3.3540034240210987 * u.TeV)
  
  
  def test_TemplateSpectralModel_evaluate_tiny():
      energy = np.array([1.00000000e06, 1.25892541e06, 1.58489319e06, 1.99526231e06])
      values = np.array([4.39150790e-38, 1.96639562e-38, 8.80497507e-39, 3.94262401e-39])
  
      model = TemplateSpectralModel(
          energy=energy, values=values * u.Unit("MeV-1 s-1 sr-1")
      )
      result = model.evaluate(energy)
      tiny = np.finfo(np.float32).tiny
      mask = abs(values) - tiny > tiny
      np.testing.assert_allclose(
          values[mask] / values.max(), result[mask].value / values.max()
      )
      mask = abs(result.value) - tiny <= tiny
      assert np.all(result[mask] == 0.0)
  
  
  def test_TemplateSpectralModel_compound():
      energy = [1.00e06, 1.25e06, 1.58e06, 1.99e06] * u.MeV
      values = [4.39e-7, 1.96e-7, 8.80e-7, 3.94e-7] * u.Unit("MeV-1 s-1 sr-1")
  
      template = TemplateSpectralModel(energy=energy, values=values)
      correction = PowerLawNormSpectralModel(norm=2)
      model = CompoundSpectralModel(template, correction, operator=operator.mul)
      assert np.allclose(model(energy), 2 * values)
  
      model_mul = template * correction
      assert isinstance(model_mul, CompoundSpectralModel)
      assert np.allclose(model_mul(energy), 2 * values)
  
      model_dict = model.to_dict()
      assert model_dict["operator"] == "mul"
      model_class = SPECTRAL_MODEL_REGISTRY.get_cls(model_dict["type"])
      new_model = model_class.from_dict(model_dict)
      assert isinstance(new_model, CompoundSpectralModel)
      assert np.allclose(new_model(energy), 2 * values)
  
  
  @requires_dependency("naima")
  class TestNaimaModel:
      # Used to test model value at 2 TeV
      energy = 2 * u.TeV
  
      # Used to test model integral and energy flux
      emin = 1 * u.TeV
      emax = 10 * u.TeV
  
      # Used to that if array evaluation works
      e_array = [2, 10, 20] * u.TeV
      e_array = e_array[:, np.newaxis, np.newaxis]
  
      def test_pion_decay(self):
          import naima
  
          particle_distribution = naima.models.PowerLaw(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=2.5
          )
          radiative_model = naima.radiative.PionDecay(
              particle_distribution, nh=1 * u.cm ** -3
          )
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 9.725347355450884e-14 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 3.530537143620737e-13 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 7.643559573105779e-13 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
          model.amplitude.error = 0.1 * model.amplitude.value
  
          out = model.evaluate_error(1 * u.TeV)
          assert_allclose(out.data, [5.266068e-13, 5.266068e-14], rtol=1e-3)
  
      def test_ic(self):
          import naima
  
          particle_distribution = naima.models.ExponentialCutoffBrokenPowerLaw(
              amplitude=2e33 / u.eV,
              e_0=10 * u.TeV,
              alpha_1=2.5,
              alpha_2=2.7,
              e_break=900 * u.GeV,
              e_cutoff=10 * u.TeV,
          )
          radiative_model = naima.radiative.InverseCompton(
              particle_distribution, seed_photon_fields=["CMB"]
          )
  
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 4.347836316893546e-12 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 1.595813e-11 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 2.851283e-11 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV, rtol=1e-5
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV, rtol=1e-5
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
      def test_synchrotron(self):
          import naima
  
          particle_distribution = naima.models.LogParabola(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=1.3, beta=0.5
          )
          radiative_model = naima.radiative.Synchrotron(particle_distribution, B=2 * u.G)
  
          model = NaimaSpectralModel(radiative_model)
  
          val_at_2TeV = 1.0565840392550432e-24 * u.Unit("cm-2 s-1 TeV-1")
          integral_1_10TeV = 4.449303e-13 * u.Unit("cm-2 s-1")
          eflux_1_10TeV = 4.594242e-13 * u.Unit("TeV cm-2 s-1")
  
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
          assert_quantity_allclose(
              model.integral(emin=self.emin, emax=self.emax), integral_1_10TeV, rtol=1e-5
          )
          assert_quantity_allclose(
              model.energy_flux(emin=self.emin, emax=self.emax), eflux_1_10TeV, rtol=1e-5
          )
          val = model(self.e_array)
          assert val.shape == self.e_array.shape
  
          model.B.value = 3  # update B
          val_at_2TeV = 5.1985064062296e-16 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV)
  
      def test_ssc(self):
          import naima
  
          ECBPL = naima.models.ExponentialCutoffBrokenPowerLaw(
              amplitude=3.699e36 / u.eV,
              e_0=1 * u.TeV,
              e_break=0.265 * u.TeV,
              alpha_1=1.5,
              alpha_2=3.233,
              e_cutoff=1863 * u.TeV,
              beta=2.0,
          )
  
          radiative_model = naima.radiative.InverseCompton(
              ECBPL,
              seed_photon_fields=[
                  "CMB",
                  ["FIR", 70 * u.K, 0.5 * u.eV / u.cm ** 3],
                  ["NIR", 5000 * u.K, 1 * u.eV / u.cm ** 3],
              ],
              Eemax=50 * u.PeV,
              Eemin=0.1 * u.GeV,
          )
          B = 125 * u.uG
          radius = 2.1 * u.pc
          nested_models = {"SSC": {"B": B, "radius": radius}}
          model = NaimaSpectralModel(radiative_model, nested_models=nested_models)
          assert_quantity_allclose(model.B.quantity, B)
          assert_quantity_allclose(model.radius.quantity, radius)
          val_at_2TeV = 1.6703761561806372e-11 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV, rtol=1e-5)
  
          model.parameters["B"].value = 100
          val_at_2TeV = 1.441331153167876e-11 * u.Unit("cm-2 s-1 TeV-1")
          value = model(self.energy)
          assert_quantity_allclose(value, val_at_2TeV, rtol=1e-5)
  
      def test_bad_init(self):
          import naima
  
          particle_distribution = naima.models.PowerLaw(
              amplitude=2e33 / u.eV, e_0=10 * u.TeV, alpha=2.5
          )
          radiative_model = naima.radiative.PionDecay(
              particle_distribution, nh=1 * u.cm ** -3
          )
          model = NaimaSpectralModel(radiative_model)
  
          with pytest.raises(NotImplementedError):
              NaimaSpectralModel.from_dict(model.to_dict())
          with pytest.raises(NotImplementedError):
              NaimaSpectralModel.from_parameters(model.parameters)
  
  
  class TestSpectralModelErrorPropagation:
      """Test spectral model error propagation.
  
      https://github.com/gammapy/gammapy/blob/master/docs/development/pigs/pig-014.rst#proposal
      https://nbviewer.jupyter.org/github/gammapy/gammapy-extra/blob/master/experiments/uncertainty_estimation_prototype.ipynb
      """
  
      def setup(self):
          self.model = LogParabolaSpectralModel(
              amplitude=3.76e-11 * u.Unit("cm-2 s-1 TeV-1"),
              reference=1 * u.TeV,
              alpha=2.44,
              beta=0.25,
          )
          self.model.covariance = [
              [1.31e-23, 0, -6.80e-14, 3.04e-13],
              [0, 0, 0, 0],
              [-6.80e-14, 0, 0.00899, 0.00904],
              [3.04e-13, 0, 0.00904, 0.0284],
          ]
  
      def test_evaluate_error_scalar(self):
          # evaluate_error on scalar
          out = self.model.evaluate_error(1 * u.TeV)
          assert isinstance(out, u.Quantity)
          assert out.unit == "cm-2 s-1 TeV-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [3.7600e-11, 3.6193e-12], rtol=1e-3)
  
      def test_evaluate_error_array(self):
          out = self.model.evaluate_error([1, 100] * u.TeV)
          assert out.shape == (2, 2)
          expected = [[3.76e-11, 2.469e-18], [3.619e-12, 9.375e-18]]
          assert_allclose(out.data, expected, rtol=1e-3)
  
      def test_evaluate_error_unit(self):
          out = self.model.evaluate_error(1e6 * u.MeV)
          assert out.unit == "cm-2 s-1 TeV-1"
          assert_allclose(out.data, [3.760e-11, 3.6193e-12], rtol=1e-3)
  
      @pytest.mark.xfail(reason="FIXME, do we need this method?")
      def test_integral_error(self):
          out = self.model.integral_error(1 * u.TeV, 10 * u.TeV)
          assert out.unit == "cm-2 s-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [2.197e-11, 2.796e-12], rtol=1e-3)
  
      @pytest.mark.xfail(reason="FIXME, do we need this method?")
      def test_energy_flux_error(self):
          out = self.model.energy_flux_error(1 * u.TeV, 10 * u.TeV)
          assert out.unit == "TeV cm-2 s-1"
          assert out.shape == (2,)
          assert_allclose(out.data, [4.119e-11, 8.157e-12], rtol=1e-3)
  
      def test_ecpl_model(self):
          # Regression test for ECPL model
          # https://github.com/gammapy/gammapy/issues/2007
          model = ExpCutoffPowerLawSpectralModel(
              amplitude=2.076183759227292e-12 * u.Unit("cm-2 s-1 TeV-1"),
              index=1.8763343736076483,
              lambda_=0.08703226432146616 * u.Unit("TeV-1"),
              reference=1 * u.TeV,
          )
          model.covariance = [
              [0.00204191498, -1.507724e-14, 0.0, -0.001834819, 0.0],
              [-1.507724e-14, 1.6864740e-25, 0.0, 1.854251e-14, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0],
              [-0.001834819175, 1.8542517e-14, 0.0, 0.0032559101, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0],
          ]
  
          out = model.evaluate_error(1 * u.TeV)
          assert_allclose(out.data, [1.903129e-12, 2.979976e-13], rtol=1e-3)
  
          out = model.evaluate_error(0.1 * u.TeV)
          assert_allclose(out.data, [1.548176e-10, 1.933612e-11], rtol=1e-3)
  
  
  def test_integral_error_PowerLaw():
      energy = np.linspace(1 * u.TeV, 10 * u.TeV, 10)
      emin = energy[:-1]
      emax = energy[1:]
  
      powerlaw = PowerLawSpectralModel()
      powerlaw.parameters['index'].error = 0.4
      powerlaw.parameters['amplitude'].error = 1e-13
  
      flux, flux_error = powerlaw.integral_error(emin,emax)
  
      assert_allclose(flux.value[0]/1e-13, 5.0, rtol=0.1)
      assert_allclose(flux_error.value[0]/1e-14, 8.546615432273905, rtol=0.01)
  
  
  def test_integral_error_ExpCutOffPowerLaw():
      energy = np.linspace(1 * u.TeV, 10 * u.TeV, 10)
      emin = energy[:-1]
      emax = energy[1:]
      
      exppowerlaw = ExpCutoffPowerLawSpectralModel()
      exppowerlaw.parameters['index'].error = 0.4
      exppowerlaw.parameters['amplitude'].error = 1e-13
      exppowerlaw.parameters['lambda_'].error = 0.03
      
      flux, flux_error = exppowerlaw.integral_error(emin, emax)
      
      assert_allclose(flux.value[0]/1e-13, 5.05855622, rtol=0.01)
      assert_allclose(flux_error.value[0]/1e-14, 8.90907063, rtol=0.01)
  
files: gammapy/modeling/models/tests/test_spectral.py