Commit 38e7f35b authored by Médéric Boquien's avatar Médéric Boquien

Implement a cache for the filters to compute the fluxes. Rather than passing...

Implement a cache for the filters to compute the fluxes. Rather than passing the transmission table as an argument to compute_fnu(), which makes caching difficult and/or slow to compute the hash, we rather pass the filter name. Then the SED object fetches the filter from the database and resamples it on the optimal wavelength grid. The result is stored in cache to avoid carrying out this operation repeatedly.
parent 899be64f
......@@ -37,7 +37,6 @@ from ...utils import read_table
from .. import AnalysisModule, complete_obs_table
from .utils import save_table_analysis, save_table_best, analyse_chi2
from ...warehouse import SedWarehouse
from ...data import Database
from .workers import sed as worker_sed
from .workers import init_sed as init_worker_sed
from .workers import init_analysis as init_worker_analysis
......@@ -134,14 +133,7 @@ class PdfAnalysis(AnalysisModule):
lim_flag = config["lim_flag"].lower() == "true"
mock_flag = config["mock_flag"].lower() == "true"
# Get the needed filters in the pcigale database. We use an ordered
# dictionary because we need the keys to always be returned in the
# same order. We also put the filters in the shared modules as they
# are needed to compute the fluxes during the models generation.
with Database() as base:
filters = OrderedDict([(name, base.get_filter(name))
for name in column_list
if not name.endswith('_err')])
filters = [name for name in column_list if not name.endswith('_err')]
n_filters = len(filters)
# Read the observation table and complete it by adding error where
......
......@@ -30,8 +30,8 @@ def init_sed(params, filters, analysed, redshifts, fluxes, variables,
----------
params: ParametersHandler
Handles the parameters from a 1D index.
filters: OrderedDict
Contains filters to compute the fluxes.
filters: List
Contains the names of the filters to compute the fluxes.
analysed: list
Variable names to be analysed.
redshifts: RawArray and tuple containing the shape
......@@ -173,9 +173,8 @@ def sed(idx):
model_fluxes = -99. * np.ones(len(gbl_filters))
model_variables = -99. * np.ones(len(gbl_analysed_variables))
else:
model_fluxes = np.array([sed.compute_fnu(filter_.trans_table,
filter_.effective_wavelength)
for filter_ in gbl_filters.values()])
model_fluxes = np.array([sed.compute_fnu(filter_) for filter_ in
gbl_filters])
model_variables = np.array([sed.info[name]
for name in gbl_analysed_variables])
......
......@@ -26,7 +26,6 @@ import time
import numpy as np
from .. import AnalysisModule
from ...data import Database
from ..utils import ParametersHandler, backup_dir, save_fluxes
from ...utils import read_table
from ...warehouse import SedWarehouse
......@@ -98,14 +97,7 @@ class SaveFluxes(AnalysisModule):
out_format = parameters["output_format"]
save_sed = parameters["save_sed"].lower() == "true"
# Get the needed filters in the pcigale database. We use an ordered
# dictionary because we need the keys to always be returned in the
# same order. We also put the filters in the shared modules as they
# are needed to compute the fluxes during the models generation.
with Database() as base:
filters = OrderedDict([(name, base.get_filter(name))
for name in column_list
if not name.endswith('_err')])
filters = [name for name in column_list if not name.endswith('_err')]
n_filters = len(filters)
w_redshifting = creation_modules.index('redshifting')
......
......@@ -22,8 +22,8 @@ def init_fluxes(params, filters, save_sed, fluxes, info, t_begin, n_computed):
----------
params: ParametersHandler
Handles the parameters from a 1D index.
filters: OrderedDict
Contains filters to compute the fluxes.
filters: List
Contains the names of the filters to compute the fluxes.
save_sed: boolean
Indicates whether the SED should be saved.
fluxes: RawArray and tuple containing the shape
......@@ -84,9 +84,8 @@ def fluxes(idx):
if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
model_fluxes = -99. * np.ones(len(gbl_filters))
else:
model_fluxes = np.array([sed.compute_fnu(filter_.trans_table,
filter_.effective_wavelength)
for filter_ in gbl_filters.values()])
model_fluxes = np.array([sed.compute_fnu(filter_) for filter_ in
gbl_filters])
gbl_model_fluxes[idx, :] = model_fluxes
gbl_model_info[idx, :] = list(sed.info.values())
......
......@@ -36,6 +36,8 @@ from . import utils
from .io.vo import save_sed_to_vo
from scipy.constants import c, parsec
from scipy.interpolate import interp1d
from ..data import Database
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
......@@ -66,6 +68,7 @@ class SED(object):
self.luminosities = None
self.info = OrderedDict()
self.mass_proportional_info = []
self.cache_filters = {}
@property
def sfh(self):
......@@ -242,18 +245,14 @@ class SED(object):
- self.contribution_names[::-1].index(name))
return self.luminosities[idx]
def compute_fnu(self, transmission, lambda_eff):
def compute_fnu(self, filter_name):
"""
Compute the Fν flux density corresponding the filter which
transmission is given.
Compute the Fν flux density in a given filter
As the SED stores the Lλ luminosity density, we first compute the Fλ
flux density. Fλ is the integration of the Lλ luminosity multiplied by
the filter transmission, normalised to this transmission and corrected
by the luminosity distance of the source. This is done by the
pcigale.sed.utils.luminosity_to_flux function.
Fλ = luminosity_to_flux( integ( LλT(λ)dλ ) / integ( T(λ)dλ ) )
by the luminosity distance of the source.
Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
to compute Fν, we make the approximation:
......@@ -267,12 +266,9 @@ class SED(object):
Parameters
----------
transmission: 2D array of floats
A numpy 2D array containing the filter response profile
wavelength[nm] vs transmission).
lambda_eff: float
Effective wavelength of the filter in nm.
filter_name: string
Name of the filter to integrate into. It must be presnt in the
database.
Return
------
......@@ -280,47 +276,55 @@ class SED(object):
The integrated Fν density in mJy.
"""
# Filter limits
lambda_min = transmission[0][0]
lambda_max = transmission[0][-1]
wavelength = self.wavelength_grid
l_lambda = self.luminosity
# Test if the spectrum cover all the filter extend
if ((wavelength[0] > lambda_min) or
(wavelength[-1] < lambda_max)):
f_nu = -99.
# First we try to fetch the filter's wavelength, transmission and
# effective wavelength from the cache. The two keys are the size of the
# spectrum wavelength grid and the name of the filter. The first key is
# necessary because different spectra may have different sampling. To
# avoid having the resample the filter every time on the optimal grid
# (spectrum+filter), we store the resampled filter. That way we only
# have to resample to spectrum.
key = (wavelength.size, filter_name)
if key in list(self.cache_filters.keys()):
wavelength_r, transmission_r, lambda_eff = self.cache_filters[key]
else:
with Database() as db:
filter_ = db.get_filter(filter_name)
trans_table = filter_.trans_table
lambda_eff = filter_.effective_wavelength
lambda_min = filter_.trans_table[0][0]
lambda_max = filter_.trans_table[0][-1]
# Test if the filter covers all the spectrum extent. If not then
# the flux is not defined
if ((wavelength[0] > lambda_min) or (wavelength[-1] < lambda_max)):
return -99.
# We regrid both spectrum and filter to the best wavelength grid
# to avoid interpolating a high wavelength density curve to a low
# density one. Also, we limit the work wavelength domain to the
# filter one.
w = np.where((wavelength >= lambda_min)&(wavelength <= lambda_max))
wavelength_r = utils.best_grid(wavelength[w], transmission[0])
l_lambda_r = np.interp(wavelength_r, wavelength, l_lambda)
transmission_r = np.interp(wavelength_r, transmission[0],
transmission[1])
wavelength_r = utils.best_grid(wavelength[w], trans_table[0])
transmission_r = np.interp(wavelength_r, trans_table[0], trans_table[1])
if 'universe.luminosity_distance' in self.info.keys():
dist = self.info['universe.luminosity_distance']
else:
dist = 10. * parsec
self.cache_filters[key] = (wavelength_r, transmission_r, lambda_eff)
f_lambda = utils.luminosity_to_flux(
np.trapz(transmission_r * l_lambda_r, wavelength_r),
dist)
l_lambda_r = np.interp(wavelength_r, wavelength, self.luminosity)
# Fν in W/m²/Hz. The 1e-9 factor is because λ is in nm.
f_nu = lambda_eff * f_lambda * lambda_eff * 1e-9 / c
if 'universe.luminosity_distance' in self.info.keys():
dist = self.info['universe.luminosity_distance']
else:
dist = 10. * parsec
# Conversion from W/m²/Hz to mJy
f_nu *= 1e+29
f_lambda = utils.luminosity_to_flux(
np.trapz(transmission_r * l_lambda_r, wavelength_r),
dist)
return f_nu
# Return Fν in mJy. The 1e-9 factor is because λ is in nm and 1e29 for
# convert from W/m²/Hz to mJy.
return lambda_eff * lambda_eff* f_lambda * 1e-9 / c * 1e29
def to_votable(self, filename, mass=1.):
"""
......@@ -352,5 +356,6 @@ class SED(object):
sed.contribution_names = self.contribution_names[:]
sed.info = self.info.copy()
sed.mass_proportional_info = self.mass_proportional_info[:]
sed.cache_filters = self.cache_filters
return sed
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