Improved integration for common Astromodels Functions, and custom function cleanup#466
Conversation
…including Band functions * minor cleanups to custom functions and flux integration
including Band spectra * assorted small cleanups related to flux integrals
all the divides in question have floating-point denominators and so will return the same result as regular /
attribute tests in eval. * don't reload file and rebuild interpolator on every call to evaluate in SpecFromDat; do so only if the filename changes * retire Wide_Asymmetric_Gaussian function, because we can get exactly the same behavior from Astromodels' Asymmetric_Gaussian (which has identical code) by changing the upper limit on one of its parameters * change notebooks that used Wide_Asymmetric_Gaussian to reflect the previous change
it is read from a file or created directly via __init__() * simplify sanity check in ExtendedResponse's get_expectation() (now also enforces unit agreement between first two axes of response and input) * remove some unneeded imports
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
|
Thanks @jdbuhler . This look good to me. Just one thing: can you please commit the tutorial with notebook executed --i.e. with the output cells? While outputs are not ideal for version control, they are used directly for the website with the documentation. As a side note, the integration functionality is a candidate to be added directly into astromodels. This is something to think a out for later. Maybe we can discuss it in the 3ML meeting at the end of the month. |
|
Will do. And yeah, I would much prefer to have the integral as a method of the Function class, with a fallback to adaptive quadrature, inside astromodels. At least one of their Band implementations has an integration function that is currently broken (only works for integer alpha), so maybe they can use my code. |
…l_fit notebook * produce output for notebooks that were changed for Asymm_Gaussian fix
|
Thanks @jdbuhler. As an additional test I check the output of the Crab spectral fit, and it looks as expected. I'm merging. |
2 participants
This PR was motivated by an investigation into the cost of spectral fitting for short-duration GRBs. Fitting requires repeatedly integrating the spectrum for every Ei bin whenever its parameters change. Currently, this integration is done with adaptive quadrature (scipy's integrate.quad), and it cannot be vectorized over bins.
We introduce a new function, get_integral_values(), that takes an Astromodels function and a vector of x's and computes the integral between each successive pair of x's. For commonly used functions, we implement the integrals using analytic expressions or special functions, thereby avoiding the need for adaptive quadrature. This is much faster to evaluate for most of the functions mentioned in response/functions.py. Notably, we provide integral functions for a Gaussian and for a Band spectrum; the latter is 20-25x faster than quadrature for typical values of alpha. The computed integrals are also somewhat more accurate for step functions (which are poorly sampled by integrate.quad) and Band (for which quadrature loses accuracy around the cutoff energy).
This work also prompted a review and cleanup of the custom functions in the threeml subdirectory and of response-related code that makes use of spectral integration. COSILike has been left alone for now because it is due to be substantially refactored with the introduction of the interfaces redesign, but it does use the new integral code via the PointSourceResponse class.
While this work benefits users of common Astromodels spectral functions, it doesn't accelerate all such functions. We don't accelerate integrals over functions evaluated by 1D linear interpolation (e.g., SpecFromDat), though that could be done in the future using trapezoidal rebinning to replace quadrature. And we don't try to decompose Composite functions representing linear combinations of simple functions -- that should be done once, rather than on every call to integration, and so would require some additional infrastructure.
There are a few uses of integrate.quad in cosipy that don't involve Astromodels spectra, notably in LineBackground (involving sums of powerlaws and Gaussians) and 3D functions (involving integrals of a 1D linearly interpolated function). These could, with a bit of work, be converted to use the new infrastructure, but there is no point in doing so until we address the two extensions suggested above.