Improve error methods on spectral models

gammapy/catalog/gammacat.py

@@ -107,7 +107,7 @@ def spectral_model(self):
                 pars['emax'] = DEFAULT_EMAX
             else:
                 pars['emax'] = d['spec_erange_max']
-            np.isnan(d['spec_erange_max'])
+
             # TODO: remove this hack once this issue is resolved in gamma-cat
             # https://github.com/gammapy/gamma-cat/issues/101
             pars['amplitude'] = d['spec_norm'].value * u.Unit('cm-2 s-1')

gammapy/spectrum/models.py

@@ -103,10 +103,10 @@ def integral(self, emin, emax, **kwargs):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        emin, emax : `~astropy.units.Quantity`
+            Lower and upper bound of integration range.
+        **kwargs : dict
+            Keyword arguments passed to `integrate_spectrum`
         """
         return integrate_spectrum(self, emin, emax, **kwargs)
 
@@ -116,10 +116,10 @@ def integral_error(self, emin, emax, **kwargs):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        emin, emax : `~astropy.units.Quantity`
+            Lower adn upper  bound of integration range.
+        **kwargs : dict
+            Keyword arguments passed to `integrate_spectrum`
 
         Returns
         -------
@@ -147,10 +147,11 @@ def energy_flux(self, emin, emax, **kwargs):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        emin, emax : `~astropy.units.Quantity`
+            Lower and upper bound of integration range.
+        **kwargs : dict
+            Keyword arguments passed to `integrate_spectrum`
+
         """
 
         def f(x):
@@ -167,10 +168,10 @@ def energy_flux_error(self, emin, emax, **kwargs):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
+        emin, emax : `~astropy.units.Quantity`
             Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        **kwargs : dict
+            Keyword arguments passed to `integrate_spectrum`
 
         Returns
         -------
@@ -421,10 +422,8 @@ def integral(self, emin, emax, **kwargs):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range.
+        emin, emax : `~astropy.units.Quantity`
+            Lower and upper bound of integration range.
 
         """
         # kwargs are passed to this function but not used
@@ -445,6 +444,42 @@ def integral(self, emin, emax, **kwargs):
         integral = prefactor * (upper - lower)
         return integral
 
+
+    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
+
     def energy_flux(self, emin, emax):
         r"""
         Compute energy flux in given energy range analytically.
@@ -459,20 +494,14 @@ def energy_flux(self, emin, emax):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        emin, emax : `~astropy.units.Quantity`
+            Lower and upper bound of integration range.
         """
+
         pars = self.parameters
         val = -1 * pars['index'].value + 2
 
-        try:
-            val_zero = np.isclose(val.n, 0)
-        except AttributeError:
-            val_zero = np.isclose(val, 0)
-
-        if val_zero:
+        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
@@ -483,6 +512,41 @@ def energy_flux(self, emin, emax):
             lower = (emin / pars['reference'].quantity) ** val
             return prefactor * (upper - lower)
 
+    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
+
     def to_sherpa(self, name='default'):
         """Return Sherpa `~sherpa.models.PowLaw1d`
 
@@ -563,10 +627,8 @@ def integral(self, emin, emax):
 
         Parameters
         ----------
-        emin : `~astropy.units.Quantity`
-            Lower bound of integration range.
-        emax : `~astropy.units.Quantity`
-            Upper bound of integration range
+        emin, emax : `~astropy.units.Quantity`
+            Lower and upper bound of integration range.
 
         """
         pars = self.parameters
@@ -579,6 +641,35 @@ def integral(self, emin, emax):
 
         return pars['amplitude'].quantity * top / bottom
 
+    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
+
     def inverse(self, value):
         """
         Return energy for a given function value of the spectral model.

gammapy/spectrum/tests/test_models.py

@@ -198,3 +198,26 @@ def test_absorbed_spectral_model():
     desired = absorption.evaluate(energy=reference, amplitude=1.)
     actual = model_ref_energy/pwl_ref_energy
     assert_quantity_allclose(actual, desired)
+
+@requires_dependency('uncertainties')
+def test_pwl_index_2_error():
+    pars, errs = {}, {}
+    pars['amplitude'] = 1E-12 * u.Unit('TeV-1 cm-2 s-1')
+    pars['reference'] = 1 * u.Unit('TeV')
+    pars['index'] = 2 * u.Unit('')
+    errs['amplitude'] = 0.1E-12 * u.Unit('TeV-1 cm-2 s-1')
+
+    pwl = PowerLaw(**pars)
+    pwl.parameters.set_parameter_errors(errs)
+
+    val, val_err = pwl.evaluate_error(1 * u.TeV)
+    assert_quantity_allclose(val, 1E-12 * u.Unit('TeV-1 cm-2 s-1'))
+    assert_quantity_allclose(val_err, 0.1E-12 * u.Unit('TeV-1 cm-2 s-1'))
+
+    flux, flux_err = pwl.integral_error(1 * u.TeV, 10 * u.TeV)
+    assert_quantity_allclose(flux, 9E-13 * u.Unit('cm-2 s-1'))
+    assert_quantity_allclose(flux_err, 9E-14 * u.Unit('cm-2 s-1'))
+
+    eflux, eflux_err = pwl.energy_flux_error(1 * u.TeV, 10 * u.TeV)
+    assert_quantity_allclose(eflux, 2.302585E-12 * u.Unit('TeV cm-2 s-1'))
+    assert_quantity_allclose(eflux_err, 0.2302585E-12 * u.Unit('TeV cm-2 s-1'))