Source code for gammapy.spectrum.models

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Spectral models for Gammapy.
"""
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
import copy
import operator
import astropy.units as u
from astropy.table import Table
from ..utils.energy import EnergyBounds
from ..utils.nddata import NDDataArray, BinnedDataAxis
from .utils import integrate_spectrum
from ..utils.scripts import make_path
from ..utils.modeling import Parameter, ParameterList

__all__ = [
    'SpectralModel',
    'PowerLaw',
    'PowerLaw2',
    'ExponentialCutoffPowerLaw',
    'ExponentialCutoffPowerLaw3FGL',
    'PLSuperExpCutoff3FGL',
    'LogParabola',
    'TableModel',
    'AbsorbedSpectralModel',
    'Absorption',
]


[docs]class SpectralModel(object): """Spectral model base class. Derived classes should store their parameters as `~gammapy.utils.modeling.ParameterList` See for example return pardict of `~gammapy.spectrum.models.PowerLaw`. """ def __repr__(self): fmt = '{}()' return fmt.format(self.__class__.__name__) def __str__(self): ss = self.__class__.__name__ ss += '\n\nParameters: \n\n\t' table = self.parameters.to_table() ss += '\n\t'.join(table.pformat()) if self.parameters.covariance is not None: ss += '\n\nCovariance: \n\n\t' covar = self.parameters.covariance_to_table() ss += '\n\t'.join(covar.pformat()) return ss
[docs] def __call__(self, energy): """Call evaluate method of derived classes""" kwargs = dict() for par in self.parameters.parameters: kwargs[par.name] = par.quantity return self.evaluate(energy, **kwargs)
def __mul__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.mul) def __rmul__(self, model): # This is needed to support e.g. 5 * model return self.__mul__(model) def __add__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(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 = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.sub) def __rsub__(self, model): return self.__sub__(model) def __truediv__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.truediv) def __rtruediv__(self, model): return self.__div__(model) def _parse_uarray(self, uarray): from uncertainties import unumpy values = unumpy.nominal_values(uarray) errors = unumpy.std_devs(uarray) return values, errors def _convert_energy(self, energy): try: energy = energy.to(self.parameters['reference'].unit) except IndexError: energy = energy.to(self.parameters['emin'].unit) return energy
[docs] def evaluate_error(self, energy): """Evaluate spectral model with error propagation. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to evaluate Returns ------- flux, flux_error : tuple of `~astropy.units.Quantity` Tuple of flux and flux error. """ energy = self._convert_energy(energy) unit = self(energy).unit upars = self.parameters._ufloats uarray = self.evaluate(energy.value, **upars) return self._parse_uarray(uarray) * unit
[docs] def integral(self, emin, emax, **kwargs): """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.spectrum.integrate_spectrum` """ return integrate_spectrum(self, emin, emax, **kwargs)
[docs] def integral_error(self, emin, emax, **kwargs): """Integrate spectral model numerically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower adn upper bound of integration range. **kwargs : dict Keyword arguments passed to func:`~gammapy.spectrum.integrate_spectrum` Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats def f(x): return self.evaluate(x, **upars) uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) return self._parse_uarray(uarray) * unit
[docs] def energy_flux(self, emin, emax, **kwargs): """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.spectrum.integrate_spectrum` """ def f(x): return x * self(x) return integrate_spectrum(f, emin, emax, **kwargs)
[docs] def energy_flux_error(self, emin, emax, **kwargs): """Compute energy flux in given energy range with error propagation. .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower bound of integration range. **kwargs : dict Keyword arguments passed to `func:`~gammapy.spectrum.integrate_spectrum` Returns ------- energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` Tuple of energy flux and energy flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.energy_flux(emin, emax, **kwargs).unit upars = self.parameters._ufloats def f(x): return x * self.evaluate(x, **upars) uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) return self._parse_uarray(uarray) * unit
[docs] def to_dict(self): """Convert to dict.""" retval = self.parameters.to_dict() retval['name'] = self.__class__.__name__ return retval
[docs] @classmethod def from_dict(cls, val): """Create from dict.""" classname = val.pop('name') parameters = ParameterList.from_dict(val) model = globals()[classname]() model.parameters = parameters model.parameters.covariance = parameters.covariance return model
[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` 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 = EnergyBounds.equal_log_spacing( emin, emax, n_points, energy_unit) # 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 = EnergyBounds.equal_log_spacing( emin, emax, n_points, energy_unit) 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
def _plot_format_ax(self, ax, energy, y, energy_power): ax.set_xlabel('Energy [{}]'.format(energy.unit)) if energy_power > 0: ax.set_ylabel('E{} * Flux [{}]'.format(energy_power, y.unit)) else: ax.set_ylabel('Flux [{}]'.format(y.unit)) ax.set_xscale("log", nonposx='clip') ax.set_yscale("log", nonposy='clip') def _plot_scale_flux(self, energy, flux, energy_power): eunit = [_ for _ in flux.unit.bases if _.physical_type == 'energy'][0] y = flux * np.power(energy, energy_power) return y.to(flux.unit * eunit ** energy_power)
[docs] def to_sherpa(self, name='default'): """Convert to Sherpa model. To be implemented by subclasses """ raise NotImplementedError('{}'.format(self.__class__.__name__))
[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``. """ from scipy.optimize import brentq energies = [] for val in np.atleast_1d(value): def f(x): # scale by 1e12 to achieve better precision y = self(x * u.TeV).to(value.unit).value return 1e12 * (y - val.value) energy = brentq(f, emin.to('TeV').value, emax.to('TeV').value) energies.append(energy) return energies * u.TeV
[docs] def copy(self): """A deep copy.""" return copy.deepcopy(self)
class ConstantModel(SpectralModel): r"""Constant model .. math:: \phi(E) = k Parameters ---------- const : `~astropy.units.Quantity` :math:`k` """ def __init__(self, const): self.parameters = ParameterList([ Parameter('const', const) ]) @staticmethod def evaluate(energy, const): """Evaluate the model (static function).""" return const class CompoundSpectralModel(SpectralModel): """Represents the algebraic combination of two `~gammapy.spectrum.models.SpectralModel` """ def __init__(self, model1, model2, operator): self.model1 = model1 self.model2 = model2 self.operator = operator @property def parameters(self): val = self.model1.parameters.parameters + self.model2.parameters.parameters return ParameterList(val) def __str__(self): ss = self.__class__.__name__ ss += '\n Component 1 : {}'.format(self.model1) ss += '\n Component 2 : {}'.format(self.model2) ss += '\n Operator : {}'.format(self.operator) return ss def __call__(self, energy): val1 = self.model1(energy) val2 = self.model2(energy) return self.operator(val1, val2) def to_dict(self): retval = dict() retval['model1'] = self.model1.to_dict() retval['model2'] = self.model2.to_dict() retval['operator'] = self.operator
[docs]class PowerLaw(SpectralModel): r"""Spectral power-law model. .. math:: \phi(E) = \phi_0 \cdot \left( \frac{E}{E_0} \right)^{-\Gamma} Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma` amplitude : `~astropy.units.Quantity` :math:`Phi_0` reference : `~astropy.units.Quantity` :math:`E_0` Examples -------- This is how to plot the default `PowerLaw` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import PowerLaw pwl = PowerLaw() pwl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, index=2., amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=1 * u.TeV): self.parameters = ParameterList([ Parameter('index', index), Parameter('amplitude', amplitude), Parameter('reference', reference, parmin=0, frozen=True) ])
[docs] @staticmethod def evaluate(energy, index, amplitude, reference): """Evaluate the model (static function).""" return amplitude * np.power((energy / reference), -index)
[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 are passed to this function but not used # this is to get a consistent API with SpectralModel.integral() pars = self.parameters if np.isclose(pars['index'].value, 1): e_unit = emin.unit prefactor = pars['amplitude'].quantity * pars['reference'].quantity.to(e_unit) upper = np.log(emax.to(e_unit).value) lower = np.log(emin.value) else: val = -1 * pars['index'].value + 1 prefactor = pars['amplitude'].quantity * pars['reference'].quantity / val upper = np.power((emax / pars['reference'].quantity), val) lower = np.power((emin / pars['reference'].quantity), val) return prefactor * (upper - lower)
[docs] def integral_error(self, emin, emax, **kwargs): r"""Integrate power law analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ # kwargs are passed to this function but not used # this is to get a consistent API with SpectralModel.integral() emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats if np.isclose(upars['index'].nominal_value, 1): prefactor = upars['amplitude'] * upars['reference'] upper = np.log(emax.value) lower = np.log(emin.value) else: val = -1 * upars['index'] + 1 prefactor = upars['amplitude'] * upars['reference'] / val upper = np.power((emax.value / upars['reference']), val) lower = np.power((emin.value / upars['reference']), val) uarray = prefactor * (upper - lower) return self._parse_uarray(uarray) * unit
[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. """ pars = self.parameters val = -1 * pars['index'].value + 2 if np.isclose(val, 0): # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) # for reference temp = pars['amplitude'].quantity * pars['reference'].quantity ** 2 return temp * np.log(emax / emin) else: prefactor = pars['amplitude'].quantity * pars['reference'].quantity ** 2 / val upper = (emax / pars['reference'].quantity) ** val lower = (emin / pars['reference'].quantity) ** val return prefactor * (upper - lower)
[docs] def energy_flux_error(self, emin, emax, **kwargs): r"""Compute energy flux in given energy range analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` Tuple of energy flux and energy flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.energy_flux(emin, emax, **kwargs).unit upars = self.parameters._ufloats val = -1 * upars['index'] + 2 if np.isclose(val.nominal_value, 0): # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) # for reference temp = upars['amplitude'] * upars['reference'] ** 2 uarray = temp * np.log(emax.value / emin.value) else: prefactor = upars['amplitude'] * upars['reference'] ** 2 / val upper = (emax.value / upars['reference']) ** val lower = (emin.value / upars['reference']) ** val uarray = prefactor * (upper - lower) return self._parse_uarray(uarray) * unit
[docs] def to_sherpa(self, name='default'): """Convert to `sherpa.models.PowLaw1D`. Parameters ---------- name : str, optional Name of the sherpa model instance """ import sherpa.models as m model = m.PowLaw1D('powlaw1d.' + name) model.gamma = self.parameters['index'].value model.ref = self.parameters['reference'].quantity.to('keV').value model.ampl = self.parameters['amplitude'].quantity.to('cm-2 s-1 keV-1').value return model
[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. """ p = self.parameters base = value / p['amplitude'].quantity return p['reference'].quantity * np.power(base, - 1. / p['index'].value)
[docs]class PowerLaw2(SpectralModel): r"""Spectral power-law model with integral as amplitude parameter. See also: https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html .. math:: \phi(E) = F_0 \cdot \frac{\Gamma + 1}{E_{0, max}^{\Gamma + 1} - E_{0, min}^{\Gamma + 1}} \cdot E^{-\Gamma} 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}`. Examples -------- This is how to plot the default `PowerLaw2` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import PowerLaw2 pwl2 = PowerLaw2() pwl2.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, amplitude=1E-12 * u.Unit('cm-2 s-1'), index=2, emin=0.1 * u.TeV, emax=100 * u.TeV): self.parameters = ParameterList([ Parameter('amplitude', amplitude), Parameter('index', index), Parameter('emin', emin, frozen=True), Parameter('emax', emax, frozen=True) ])
[docs] @staticmethod def evaluate(energy, amplitude, index, emin, emax): """Evaluate the model (static function).""" top = -index + 1 bottom = emax ** (-index + 1) - emin ** (-index + 1) return amplitude * (top / bottom) * np.power(energy, -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. """ pars = self.parameters temp1 = np.power(emax, -pars['index'].value + 1) temp2 = np.power(emin, -pars['index'].value + 1) top = temp1 - temp2 temp1 = np.power(pars['emax'].quantity, -pars['index'].value + 1) temp2 = np.power(pars['emin'].quantity, -pars['index'].value + 1) bottom = temp1 - temp2 return pars['amplitude'].quantity * top / bottom
[docs] def integral_error(self, emin, emax, **kwargs): r"""Integrate power law analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats temp1 = np.power(emax.value, -upars['index'] + 1) temp2 = np.power(emin.value, -upars['index'] + 1) top = temp1 - temp2 temp1 = np.power(upars['emax'], -upars['index'] + 1) temp2 = np.power(upars['emin'], -upars['index'] + 1) bottom = temp1 - temp2 uarray = upars['amplitude'] * top / bottom return self._parse_uarray(uarray) * unit
[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. """ p = self.parameters index = p['index'].value top = -index + 1 bottom = (p['emax'].quantity.to('TeV').value ** (-index + 1) - p['emin'].quantity.to('TeV').value ** (-index + 1)) term = (bottom / top) * (value / p['amplitude'].quantity).to('1 / TeV') return np.power(term.value, -1. / index) * u.TeV
[docs]class ExponentialCutoffPowerLaw(SpectralModel): r"""Spectral exponential cutoff power-law model. .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp(-\lambda E) 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` Examples -------- This is how to plot the default `ExponentialCutoffPowerLaw` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import ExponentialCutoffPowerLaw ecpl = ExponentialCutoffPowerLaw() ecpl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, index=1.5, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=1 * u.TeV, lambda_=0.1 / u.TeV): self.parameters = ParameterList([ Parameter('index', index), Parameter('amplitude', amplitude), Parameter('reference', reference, frozen=True), Parameter('lambda_', lambda_) ])
[docs] @staticmethod def evaluate(energy, index, amplitude, reference, lambda_): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index) try: cutoff = np.exp(-energy * lambda_) except AttributeError: from uncertainties.unumpy import exp cutoff = exp(-energy * lambda_) return pwl * cutoff
[docs] def to_sherpa(self, name='default'): """Convert to a `~sherpa.models.ArithmeticModel`. Parameters ---------- name : str, optional Name of the sherpa model instance """ from .sherpa_utils import SherpaExponentialCutoffPowerLaw model = SherpaExponentialCutoffPowerLaw(name='ecpl.' + name) pars = self.parameters model.gamma = pars['index'].value model.ref = pars['reference'].quantity.to('keV').value model.ampl = pars['amplitude'].quantity.to('cm-2 s-1 keV-1').value # Sherpa ExponentialCutoffPowerLaw expects cutoff in 1/TeV model.cutoff = pars['lambda_'].quantity.to('TeV-1').value return model
@property def e_peak(self): r"""Spectral energy distribution peak energy (`~astropy.utils.Quantity`). This is the peak in E^2 x dN/dE and is given by: .. math:: E_{Peak} = (2 - \Gamma) / \lambda """ p = self.parameters reference = p['reference'].quantity index = p['index'].quantity lambda_ = p['lambda_'].quantity if index >= 2: return np.nan * reference.unit else: return (2 - index) / lambda_
[docs]class ExponentialCutoffPowerLaw3FGL(SpectralModel): r"""Spectral exponential cutoff power-law model used for 3FGL. Note that the parametrization is different from `ExponentialCutoffPowerLaw`: .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp \left( \frac{E_0 - E}{E_{C}} \right) 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}` Examples -------- This is how to plot the default `ExponentialCutoffPowerLaw3FGL` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import ExponentialCutoffPowerLaw3FGL ecpl_3fgl = ExponentialCutoffPowerLaw3FGL() ecpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, index=1.5, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=1 * u.TeV, ecut=10 * u.TeV): self.parameters = ParameterList([ Parameter('index', index), Parameter('amplitude', amplitude), Parameter('reference', reference, frozen=True), Parameter('ecut', ecut) ])
[docs] @staticmethod def evaluate(energy, index, amplitude, reference, ecut): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index) try: cutoff = np.exp((reference - energy) / ecut) except AttributeError: from uncertainties.unumpy import exp cutoff = exp((reference - energy) / ecut) return pwl * cutoff
[docs]class PLSuperExpCutoff3FGL(SpectralModel): r"""Spectral super exponential cutoff power-law model used for 3FGL. .. 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}` Examples -------- This is how to plot the default `PLSuperExpCutoff3FGL` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import PLSuperExpCutoff3FGL secpl_3fgl = PLSuperExpCutoff3FGL() secpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, index_1=1.5, index_2=2, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=1 * u.TeV, ecut=10 * u.TeV): # TODO: order or parameters is different from argument list / docstring. Make uniform! self.parameters = ParameterList([ Parameter('amplitude', amplitude), Parameter('reference', reference, frozen=True), Parameter('ecut', ecut), Parameter('index_1', index_1), Parameter('index_2', index_2), ])
[docs] @staticmethod def evaluate(energy, amplitude, reference, ecut, index_1, index_2): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index_1) try: cutoff = np.exp((reference / ecut) ** (index_2) - (energy / ecut) ** (index_2)) except AttributeError: from uncertainties.unumpy import exp cutoff = exp((reference / ecut) ** (index_2) - (energy / ecut) ** (index_2)) return pwl * cutoff
[docs]class LogParabola(SpectralModel): r"""Spectral log parabola model. .. math:: \phi(E) = \phi_0 \left( \frac{E}{E_0} \right) ^ { - \alpha - \beta \log{ \left( \frac{E}{E_0} \right) } } Note that :math:`log` refers to the natural logarithm. This is consistent with the `Fermi Science Tools <https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html>`_ and `ctools <http://cta.irap.omp.eu/ctools-devel/users/user_manual/getting_started/models.html#log-parabola>`_. The `Sherpa <http://cxc.harvard.edu/sherpa/ahelp/logparabola.html_ package>`_ package, however, uses :math:`log_{10}`. If you have parametrization based on :math:`log_{10}` you can use the :func:`~gammapy.spectrum.models.LogParabola.from_log10` method. 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` Examples -------- This is how to plot the default `LogParabola` model: .. code:: python from astropy import units as u from gammapy.spectrum.models import LogParabola log_parabola = LogParabola() log_parabola.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, amplitude=1E-12 * u.Unit('cm-2 s-1 TeV-1'), reference=10 * u.TeV, alpha=2, beta=1): self.parameters = ParameterList([ Parameter('amplitude', amplitude), Parameter('reference', reference, frozen=True), Parameter('alpha', alpha), Parameter('beta', beta) ])
[docs] @classmethod def from_log10(cls, amplitude, reference, alpha, beta): """Construct LogParabola 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).""" # TODO: can this comment be removed? # cast dimensionless values as np.array, because of bug in Astropy < v1.2 # https://github.com/astropy/astropy/issues/4764 try: xx = (energy / reference).to('') exponent = -alpha - beta * np.log(xx) except AttributeError: from uncertainties.unumpy import log xx = energy / reference exponent = -alpha - beta * log(xx) return amplitude * np.power(xx, exponent)
@property def e_peak(self): r"""Spectral energy distribution peak energy (`~astropy.utils.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)} """ p = self.parameters reference = p['reference'].quantity alpha = p['alpha'].quantity beta = p['beta'].quantity return reference * np.exp((2 - alpha) / (2 * beta))
[docs]class TableModel(SpectralModel): """A model generated from a table of energy and value arrays. the units returned will be the units of the values array provided at initialization. The model will return values interpolated in log-space, returning 0 for energies outside of the limits of the provided energy array. Class implementation follows closely what has been done in `naima.models.TableModel` Parameters ---------- energy : `~astropy.units.Quantity` array Array of energies at which the model values are given values : array Array with the values of the model at energies ``energy``. scale : float Model scale that is multiplied to the supplied arrays. Defaults to 1. scale_logy : boolean interpolation can be done linearly or in logarithm """ def __init__(self, energy, values, scale=1, scale_logy=True): from scipy.interpolate import interp1d self.parameters = ParameterList([ Parameter('scale', scale, parmin=0, unit='') ]) self.energy = energy self.values = values self.scale_logy = scale_logy self.lo_threshold = energy[0] self.hi_threshold = energy[-1] loge = np.log10(self.energy.to('eV').value) try: self.unit = self.values.unit if scale_logy is True: y = np.log10(self.values.value) else: y = self.values.value except AttributeError: self.unit = u.Unit('') if scale_logy is True: y = np.log10(self.values) else: y = self.values # The type conversion is a fix for: # https://travis-ci.org/gammapy/gammapy/jobs/210576260 self.interpy = interp1d(loge.astype(float), y.astype(float), fill_value=-np.Inf, bounds_error=False, kind='cubic')
[docs] @classmethod def read_xspec_model(cls, filename, param): """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.spectrum.models import TableModel >>> filename = '$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz' >>> table_model = TableModel.read_xspec_model(filename=filename, param=0.3) """ filename = str(make_path(filename)) # Check if parameter value is in range table_param = Table.read(filename, hdu='PARAMETERS') param_min = table_param['MINIMUM'] param_max = table_param['MAXIMUM'] if param < param_min or param > param_max: err = 'Parameter out of range, param={}, param_min={}, param_max={}'.format( param, param_min, param_max) raise ValueError(err) # Get energy values table_energy = Table.read(filename, hdu='ENERGIES') energy_lo = table_energy['ENERG_LO'] energy_hi = table_energy['ENERG_HI'] # Hack while format is not fixed, energy values are in keV energy_bounds = EnergyBounds.from_lower_and_upper_bounds(lower=energy_lo, upper=energy_hi, unit=u.keV) energy = energy_bounds.log_centers # 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 = table_spectra[idx][1] * u.Unit('') # no dimension return cls(energy=energy, values=values, scale_logy=False)
[docs] def evaluate(self, energy, scale): """Evaluate the model (static function).""" # What's with all this checking? # TODO: Try `np.asanyarray` and always return an array (even for scalar input)? is_array = True try: len(energy) except: is_array = False # Not working for astropy.units.quantity.Quantity # if isinstance(energy, (np.ndarray, np.generic)): if is_array: # Test if array # initialise array value to zero (dim energy) values = np.zeros(len(energy), dtype=float) # mask for energy range mask = (energy >= self.lo_threshold) & ( energy <= self.hi_threshold) # apply interpolation for masked values values[mask] = self.interpy(np.log10(energy[mask].to('eV').value)) # Get rid of negative values (due to interpolation) # Needed because of the rand.poisson used in SpectrumSimulation class # Should be fixed in the class itself ? if self.scale_logy is False: values[values < 0] = 0. else: # if not array # test if energy is in range if (energy >= self.lo_threshold or energy <= self.hi_threshold): values = self.interpy(np.log10(energy.to('eV').value)) else: values = 0 if self.scale_logy: values = np.power(10, values) return scale * values * self.unit
[docs] def plot(self, energy_range, ax=None, energy_unit='TeV', n_points=100, **kwargs): """Plot spectral model curve. kwargs are forwarded to :func:`~matplotlib.pyplot.plot` Parameters ---------- energy_range : `~astropy.units.Quantity` Plot range ax : `~matplotlib.axes.Axes`, optional Axis energy_unit : str, `~astropy.units.Unit`, optional Unit of the energy axis 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 = EnergyBounds.equal_log_spacing(emin, emax, n_points, energy_unit) y = self.interpy(np.log10(energy.to('eV').value)) * self.parameters['scale'].quantity if self.scale_logy: y = np.power(10, y) ax.plot(energy.value, y, **kwargs) ax.set_xlabel('Energy [{}]'.format(energy.unit)) ax.set_ylabel('Table model') ax.set_xscale("log", nonposx='clip') if self.scale_logy: ax.set_yscale("log", nonposy='clip') return ax
[docs]class Absorption(object): """Gamma-ray absorption models. Parameters ---------- energy_lo, energy_hi : `~astropy.units.Quantity` Lower and upper bin edges of energy axis param_lo, param_hi : `~astropy.units.Quantity` Lower and upper bin edges of parameter axis data : `~astropy.units.Quantity` Model value Examples -------- Create and plot EBL absorption models for a redshift of 0.5: .. plot:: :include-source: import matplotlib.pyplot as plt import astropy.units as u from gammapy.spectrum.models import Absorption # Load tables for z=0.5 redshift = 0.5 dominguez = Absorption.read_builtin('dominguez').table_model(redshift) franceschini = Absorption.read_builtin('franceschini').table_model(redshift) finke = Absorption.read_builtin('finke').table_model(redshift) # start customised plot energy_range = [0.08, 3] * u.TeV ax = plt.gca() opts = dict(energy_range=energy_range, energy_unit='TeV', ax=ax) franceschini.plot(label='Franceschini 2008', **opts) finke.plot(label='Finke 2010', **opts) dominguez.plot(label='Dominguez 2011', **opts) # tune plot ax.set_ylabel(r'Absorption coefficient [$\exp{(-\tau(E))}$]') ax.set_xlim(energy_range.value) # we get ride of units ax.set_ylim([1.e-4, 2.]) ax.set_yscale('log') ax.set_title('EBL models (z=' + str(redshift) + ')') plt.grid(which='both') plt.legend(loc='best') # legend # show plot plt.show() """ def __init__(self, energy_lo, energy_hi, param_lo, param_hi, data): axes = [ BinnedDataAxis(param_lo, param_hi, interpolation_mode='linear', name='parameter'), BinnedDataAxis(energy_lo, energy_hi, interpolation_mode='log', name='energy'), ] self.data = NDDataArray(axes=axes, data=data) self.data.default_interp_kwargs['fill_value'] = None
[docs] @classmethod def read(cls, filename): """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. """ # Create EBL data array filename = str(make_path(filename)) table_param = Table.read(filename, hdu='PARAMETERS') par_min = table_param['MINIMUM'] par_max = table_param['MAXIMUM'] par_array = table_param[0]['VALUE'] par_delta = np.diff(par_array) * 0.5 param_lo, param_hi = par_array, par_array # initialisation param_lo[0] = par_min - par_delta[0] param_lo[1:] -= par_delta param_hi[:-1] += par_delta param_hi[-1] = par_max # Get energy values table_energy = Table.read(filename, hdu='ENERGIES') energy_lo = table_energy['ENERG_LO'] * u.keV # unit not stored in file energy_hi = table_energy['ENERG_HI'] * u.keV # unit not stored in file # Energies are in keV energy_bounds = EnergyBounds.from_lower_and_upper_bounds(lower=energy_lo, upper=energy_hi, unit=u.keV) # Get spectrum values table_spectra = Table.read(filename, hdu='SPECTRA') data = table_spectra['INTPSPEC'].data return cls( energy_lo=energy_bounds.lower_bounds, energy_hi=energy_bounds.upper_bounds, param_lo=param_lo, param_hi=param_hi, data=data, )
[docs] @classmethod def read_builtin(cls, name): """Read one of the built-in absorption models. Parameters ---------- name : {'franceschini', 'dominguez', 'finke'} name of one of the available model in gammapy-extra References ---------- .. [1] Franceschini et al., "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity", `Link <http://adsabs.harvard.edu/abs/2008A%26A...487..837F>`__ .. [2] Dominguez et al., " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions" `Link <http://adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__ .. [3] Finke et al., "Modeling the Extragalactic Background Light from Stars and Dust" `Link <http://adsabs.harvard.edu/abs/2010ApJ...712..238F>`__ """ models = dict() models['franceschini'] = '$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz' models['dominguez'] = '$GAMMAPY_EXTRA/datasets/ebl/ebl_dominguez11.fits.gz' models['finke'] = '$GAMMAPY_EXTRA/datasets/ebl/frd_abs.fits.gz' return cls.read(models[name])
[docs] def table_model(self, parameter, unit='TeV'): """Table model for a given parameter (`~gammapy.spectrum.models.TableModel`). Parameters ---------- parameter : `float` Parameter value. unit : `str`, (optional) desired value for energy axis """ energy_axis = self.data.axes[1] energy = (energy_axis.log_center()).to(unit) values = self.evaluate(energy=energy, parameter=parameter) return TableModel(energy=energy, values=values, scale_logy=False)
[docs] def evaluate(self, energy, parameter): """Evaluate model for energy and parameter value.""" return self.data.evaluate(energy=energy, parameter=parameter)
[docs]class AbsorbedSpectralModel(SpectralModel): """Spectral model with EBL absorption. Parameters ---------- spectral_model : `~gammapy.spectrum.models.SpectralModel` Spectral model. absorption : `~gammapy.spectrum.models.Absorption` Absorption model. parameter : `float` parameter value for absorption model parameter_name : `str`, optional parameter name """ def __init__(self, spectral_model, absorption, parameter, parameter_name='redshift'): self.spectral_model = spectral_model self.absorption = absorption self.parameter = parameter self.parameter_name = parameter_name # initialise list parameters from spectral model param_list = [] for param in spectral_model.parameters.parameters: param_list.append(param) # Add parameter to the list param_min = self.absorption.data.axes[0].lo[0] param_max = self.absorption.data.axes[0].lo[-1] par = Parameter(parameter_name, parameter, parmin=param_min, parmax=param_max, frozen=True) param_list.append(par) self.parameters = ParameterList(param_list)
[docs] def evaluate(self, energy, **kwargs): """Evaluate the model at a given energy.""" # assign redshift value and remove it from dictionnary # since it does not belong to the spectral model parameter = kwargs[self.parameter_name] del kwargs[self.parameter_name] flux = self.spectral_model.evaluate(energy=energy, **kwargs) absorption = self.absorption.evaluate(energy=energy, parameter=parameter) return flux * absorption