Simpilify _trapz_loglog
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from astropy.units import Quantity
+from .interpolation import LogScale
__all__ = ["integrate_spectrum"]
@@ -77,52 +78,28 @@ def _trapz_loglog(y, x, axis=-1, intervals=False):
trapz : float
Definite integral as approximated by trapezoidal rule in loglog space.
"""
- log10 = np.log10
-
- try:
- y_unit = y.unit
- y = y.value
- except AttributeError:
- y_unit = 1.0
- try:
- x_unit = x.unit
- x = x.value
- except AttributeError:
- x_unit = 1.0
-
- y = np.asanyarray(y)
- x = np.asanyarray(x)
-
- slice1 = [slice(None)] * y.ndim
- slice2 = [slice(None)] * y.ndim
- slice1[axis] = slice(None, -1)
- slice2[axis] = slice(1, None)
- slice1, slice2 = tuple(slice1), tuple(slice2)
-
- if x.ndim == 1:
- shape = [1] * y.ndim
- shape[axis] = x.shape[0]
- x = x.reshape(shape)
-
- with np.errstate(invalid="ignore", divide="ignore"):
- # Compute the power law indices in each integration bin
- b = log10(y[slice2] / y[slice1]) / log10(x[slice2] / x[slice1])
-
- # if local powerlaw index is -1, use \int 1/x = log(x); otherwise use normal
- # powerlaw integration
- trapzs = np.where(
- np.abs(b + 1.0) > 1e-10,
- (y[slice1] * (x[slice2] * (x[slice2] / x[slice1]) ** b - x[slice1]))
- / (b + 1),
- x[slice1] * y[slice1] * np.log(x[slice2] / x[slice1]),
- )
+ from gammapy.modeling.models import PowerLawSpectralModel
- tozero = (y[slice1] == 0.0) + (y[slice2] == 0.0) + (x[slice1] == x[slice2])
- trapzs[tozero] = 0.0
+ # see https://stackoverflow.com/a/56840428
+ x, y = np.moveaxis(x, axis, 0), np.moveaxis(y, axis, 0)
- if intervals:
- return trapzs * x_unit * y_unit
+ emin, emax = x[:-1], x[1:]
+ vals_emin, vals_emax = y[:-1], y[1:]
+
+ # log scale has the build-in zero clipping
+ log = LogScale()
+ index = -log(vals_emin / vals_emax) / log(emin / emax)
+ index[np.isnan(index)] = np.inf
- ret = np.add.reduce(trapzs, axis) * x_unit * y_unit
+ integral = PowerLawSpectralModel.evaluate_integral(
+ emin=emin,
+ emax=emax,
+ index=index,
+ reference=emin,
+ amplitude=vals_emin
+ )
+
+ if intervals:
+ return integral
- return ret
+ return integral.sum(axis=axis)
