Commit ce4d814d authored by Médéric Boquien's avatar Médéric Boquien
Browse files

When computing the flux in filters, np.trapz() become the bottleneck of the...

When computing the flux in filters, np.trapz() become the bottleneck of the code. A large part of the time is spent on safeguards and on operations for nD arrays. However in our specific case we only have 1D arrays. Also we can cache some computed variables, for instance dx which only depends on the wavelength sampling. We do that employing the same key that is used for caching the best sampling of each filter in compute_fnu(). That way the caches are consistent. The key is based on the size of the wavelength grid, the name of the filter, and the redshift.
parent 8bd942a9
......@@ -321,7 +321,7 @@ class SED(object):
dist = 10. * parsec
f_lambda = utils.luminosity_to_flux(
np.trapz(transmission_r * l_lambda_r, wavelength_r),
utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key),
dist)
# Return Fν in mJy. The 1e-9 factor is because λ is in nm and 1e29 for
......
......@@ -7,6 +7,9 @@ import numpy as np
from scipy.constants import c, pi
from scipy.interpolate import interp1d
# Cache of dx for integrate(y,dx) done by flux_trapz
dx_cache = {}
def lambda_to_nu(wavelength):
"""Convert wavelength (nm) to frequency (Hz)
......@@ -366,3 +369,40 @@ def interpolate_lumin(wl, lumin, wl_new, lumin_new):
lumin_out[-1, :] = np.interp(wl_best, wl_new, lumin_new, left=0., right=0.)
return (wl_best, lumin_out)
def flux_trapz(y, x, key):
"""
Light weight trapezoidal rule used to compute the flux in filters. It
assumes that all the x and y input arrays are already numpy arrays to save
time. Also the width between x elements are saved in cache using the "key"
parameter as for a given x sampling it will always be the same. We do not
compute they key ourselves because we do not have a proper way to hash it.
Using x[0], x[-1], and x.size is not sufficient as if the redshift is small
it is very likely that not of these elements will change. However the
calling function has this knowledge and actually uses them to get the
resampled filters from the cache. That way the two cache are sure to remain
consistent.
Parameters
----------
y: 1D array
The values to be integrated, typically fluxes
x: 1D array
Sampling of y, typically wavelengths
key: tuple
Contains the size of the x array, the name of the filter, and the
redshift. Those elements point to a unique x grid. This is used to
cache some computations are the x sampling will not change.
Returns
-------
flux_trapz: float
integral(y, dx)
"""
if key in dx_cache:
dx = dx_cache[key]
else:
dx = np.diff(x)
dx_cache[key] = dx
return (dx*(y[1:]+y[:-1])).sum()/2.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment