Source code for gammapy.utils.integrate

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from astropy.units import Quantity

__all__ = ["integrate_spectrum"]


[docs]def integrate_spectrum(func, xmin, xmax, 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. xmin : `~astropy.units.Quantity` or array-like Integration range minimum xmax : `~astropy.units.Quantity` or array-like 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 """ is_quantity = False if isinstance(xmin, Quantity): unit = xmin.unit xmin = xmin.value xmax = xmax.to_value(unit) is_quantity = True if np.isscalar(xmin): logmin = np.log10(xmin) logmax = np.log10(xmax) n = int((logmax - logmin) * ndecade) x = np.logspace(logmin, logmax, n) else: x = np.append(xmin, xmax[-1]) if is_quantity: x = x * unit y = func(x) val = _trapz_loglog(y, x, intervals=intervals) return val
# This function is copied over from https://github.com/zblz/naima/blob/master/naima/utils.py#L261 # and slightly modified to allow use with the uncertainties package def _trapz_loglog(y, x, axis=-1, intervals=False): """Integrate using the composite trapezoidal rule in log-log space. Integrate `y` (`x`) along given axis in loglog space. Parameters ---------- y : array_like Input array to integrate. x : array_like, optional Independent variable to integrate over. axis : int, optional Specify the axis. intervals : bool, optional Return array of shape x not the total integral, default: False Returns ------- 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]), ) tozero = (y[slice1] == 0.0) + (y[slice2] == 0.0) + (x[slice1] == x[slice2]) trapzs[tozero] = 0.0 if intervals: return trapzs * x_unit * y_unit ret = np.add.reduce(trapzs, axis) * x_unit * y_unit return ret