# 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 evaluate_integral_pwl, 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()
[docs]class SpectralModel(Model):
"""Spectral model base class."""
[docs] 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)
@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 __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
[docs] 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)
[docs] def integral(self, emin, emax, **kwargs):
r"""Integrate spectral model numerically.
.. 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`
"""
return integrate_spectrum(self, emin, emax, **kwargs)
[docs] 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)
return integrate_spectrum(f, emin, emax, **kwargs)
[docs] 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
[docs] 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)
[docs] 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)
[docs] 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)
[docs]class ConstantSpectralModel(SpectralModel):
r"""Constant model.
For more information see :ref:`constant-spectral-model`.
Parameters
----------
const : `~astropy.units.Quantity`
:math:`k`
"""
tag = "ConstantSpectralModel"
const = Parameter("const", "1e-12 cm-2 s-1 TeV-1")
[docs] @staticmethod
def evaluate(energy, const):
"""Evaluate the model (static function)."""
return np.ones(np.atleast_1d(energy).shape) * const
[docs]class CompoundSpectralModel(SpectralModel):
"""Arithmetic combination of two spectral models.
For more information see :ref:`compound-spectral-model`.
"""
tag = "CompoundSpectralModel"
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}\n"
)
[docs] def __call__(self, energy):
val1 = self.model1(energy)
val2 = self.model2(energy)
return self.operator(val1, val2)
[docs] def to_dict(self):
return {
"model1": self.model1.to_dict(),
"model2": self.model2.to_dict(),
"operator": self.operator,
}
[docs]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`
"""
tag = "PowerLawSpectralModel"
index = Parameter("index", 2.0)
amplitude = Parameter("amplitude", "1e-12 cm-2 s-1 TeV-1")
reference = Parameter("reference", "1 TeV", frozen=True)
evaluate_integral = staticmethod(evaluate_integral_pwl)
[docs] @staticmethod
def evaluate(energy, index, amplitude, reference):
"""Evaluate the model (static function)."""
return amplitude * np.power((energy / reference), -index)
[docs] @staticmethod
def evaluate_energy_flux(emin, emax, index, amplitude, reference):
"""Evaluate the energy flux (static function)"""
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
[docs] def integral(self, emin, emax, **kwargs):
r"""Integrate power law analytically.
.. 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
"""
kwargs = {par.name: par.quantity for par in self.parameters}
kwargs = self._convert_evaluate_unit(kwargs, emin)
return self.evaluate_integral(emin=emin, emax=emax, **kwargs)
[docs] def energy_flux(self, emin, emax):
r"""Compute energy flux in given energy range analytically.
.. 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.
"""
kwargs = {par.name: par.quantity for par in self.parameters}
kwargs = self._convert_evaluate_unit(kwargs, emin)
return self.evaluate_energy_flux(emin=emin, emax=emax, **kwargs)
[docs] 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))
[docs]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}`.
"""
tag = "PowerLaw2SpectralModel"
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)
[docs] @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)
[docs] def integral(self, emin, emax, **kwargs):
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, -self.index.value + 1)
temp2 = np.power(emin, -self.index.value + 1)
top = temp1 - temp2
temp1 = np.power(self.emax.quantity, -self.index.value + 1)
temp2 = np.power(self.emin.quantity, -self.index.value + 1)
bottom = temp1 - temp2
return self.amplitude.quantity * top / bottom
[docs] 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
[docs]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`
"""
tag = "SmoothBrokenPowerLawSpectralModel"
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)
[docs] @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
[docs]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`
"""
tag = "ExpCutoffPowerLawSpectralModel"
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)
[docs] @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_
[docs]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"
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")
[docs] @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
[docs]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"
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)
[docs] @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
[docs]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"
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)
[docs] @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
[docs]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`
"""
tag = "LogParabolaSpectralModel"
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)
[docs] @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_)
[docs] @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))
[docs]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``.
norm : float
Model scale that is multiplied to the supplied arrays. Defaults to 1.
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"
norm = Parameter("norm", 1, unit="")
tilt = Parameter("tilt", 0, unit="", frozen=True)
reference = Parameter("reference", "1 TeV", frozen=True)
def __init__(
self,
energy,
values,
norm=norm.quantity,
tilt=tilt.quantity,
reference=reference.quantity,
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__(norm=norm, tilt=tilt, reference=reference)
[docs] @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)
[docs] def evaluate(self, energy, norm, tilt, reference):
"""Evaluate the model (static function)."""
values = self._evaluate((energy,), clip=True)
tilt_factor = np.power(energy / reference, -tilt)
return norm * values * tilt_factor
[docs] 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.values.data.tolist(),
"unit": str(self.values.unit),
},
}
[docs] @classmethod
def from_dict(cls, data):
parameters = Parameters.from_dict(data["parameters"])
energy = u.Quantity(data["energy"]["data"], data["energy"]["unit"])
values = u.Quantity(data["values"]["data"], data["values"]["unit"])
return cls.from_parameters(parameters, energy=energy, values=values)
[docs]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"
norm = Parameter("norm", 1, unit="")
def __init__(self, model, norm=norm.quantity):
self.model = model
self._covariance = None
super().__init__(norm=norm)
[docs] def evaluate(self, energy, norm):
return norm * self.model(energy)
[docs] def integral(self, emin, emax, **kwargs):
return self.norm.value * self.model.integral(emin, emax, **kwargs)
[docs]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.data = data
self.filename = filename
# set values log centers
self.param = param
self.energy = energy
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=data, **interp_kwargs
)
[docs] 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}
[docs] @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)
[docs] @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,
)
[docs] @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)
[docs] 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"}
)
[docs] def evaluate(self, energy, parameter):
"""Evaluate model for energy and parameter value."""
return np.clip(self._evaluate((parameter, energy)), 0, 1)
[docs]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"
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
)
[docs] 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
[docs] 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(),
}
[docs] @classmethod
def from_dict(cls, data):
from gammapy.modeling.models import SPECTRAL_MODELS
model_class = SPECTRAL_MODELS.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
[docs]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"
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
[docs] 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)
[docs] @classmethod
def from_dict(cls, data):
raise NotImplementedError(
"Currently the NaimaSpectralModel cannot be read from YAML"
)
[docs]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"
norm = Parameter("norm", 1e-12 * u.Unit("cm-2 s-1"))
mean = Parameter("mean", 1 * u.TeV)
sigma = Parameter("sigma", 2 * u.TeV)
[docs] @staticmethod
def evaluate(energy, norm, mean, sigma):
return (
norm
/ (sigma * np.sqrt(2 * np.pi))
* np.exp(-((energy - mean) ** 2) / (2 * sigma ** 2))
)
[docs] 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))
)
[docs] 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)
)