import numpy as np
import tables
from beast.physicsmodel.grid import SpectralGrid
from beast.external.eztables import Table
__all__ = ["trim_models"]
[docs]def trim_models(
sedgrid,
sedgrid_noisemodel,
obsdata,
sed_outname,
noisemodel_outname,
sigma_fac=3.0,
n_detected=4,
inFlux=True,
trunchen=False,
):
"""
For a given set of observations, there will be models that are so
bright or faint that they will always have ~0 probability of fitting
the data. This program trims those models out of the SED grid
so that time is not spent calculating model points that are always
zero probability.
Parameters
----------
sedgrid : grid.SEDgrid instance
model grid
sedgrid_noisemodel : beast noisemodel instance
noise model data
obsdata : Observation object instance
observation catalog
sed_outname : str
name for output sed file
noisemodel_outname : str
name for output noisemodel file
sigma_fac : float, optional
factor for trimming the upper and lower range of grid so that
the model range cuts off sigma_fac above and below the brightest
and faintest models, respectively (default: 3.)
n_detected : int, optional
minimum number of bands where ASTs yielded a detection for
a given model, if fewer detections than n_detected this model
gets eliminated (default: 4)
inFlux : bool, optional
if true data are in fluxes (default: True)
trunchen : bool, optional
if true use the trunchen noise model (default: False)
"""
# Store the brigtest and faintest fluxes in each band (for data and asts)
n_filters = len(obsdata.filters)
min_data = np.zeros(n_filters)
max_data = np.zeros(n_filters)
min_models = np.zeros(n_filters)
max_models = np.zeros(n_filters)
for k, filtername in enumerate(obsdata.filters):
sfiltname = obsdata.data.resolve_alias(filtername)
if inFlux:
min_data[k] = np.amin(obsdata.data[sfiltname] * obsdata.vega_flux[k])
max_data[k] = np.amax(obsdata.data[sfiltname] * obsdata.vega_flux[k])
else:
min_data[k] = np.amin(
10 ** (-0.4 * obsdata.data[sfiltname]) * obsdata.vega_flux[k]
)
max_data[k] = np.amax(
10 ** (-0.4 * obsdata.data[sfiltname]) * obsdata.vega_flux[k]
)
min_models[k] = np.amin(sedgrid.seds[:, k])
max_models[k] = np.amax(sedgrid.seds[:, k])
# first remove all models that have any band with fluxes below the
# faintest ASTs run
# when the noisemodel was computed, models with fluxes below the
# faintest ASTs were tagged with a negative error/uncertainty
# identify the models that have been detected in enough bands
# the idea here is that if the ASTs are not measured that means
# that *none* were recovered and this implies
# that no model with these values would be recovered and thus the
# probability should always be zero
model_unc = sedgrid_noisemodel["error"]
above_ast = model_unc > 0
sum_above_ast = np.sum(above_ast, axis=1)
indxs, = np.where(sum_above_ast >= n_detected)
# cache the noisemodel values
model_bias = sedgrid_noisemodel["bias"]
model_unc = np.fabs(sedgrid_noisemodel["error"])
model_compl = sedgrid_noisemodel["completeness"]
if trunchen:
model_q_norm = sedgrid_noisemodel["q_norm"]
model_icov_diag = sedgrid_noisemodel["icov_diag"]
model_icov_offdiag = sedgrid_noisemodel["icov_offdiag"]
if len(indxs) <= 0:
raise ValueError("no models are brighter than the minimum ASTs run")
n_ast_indxs = len(indxs)
# Find models with fluxes (with margin) between faintest and brightest data
for k in range(n_filters):
print("working on filter # = ", k)
# Get upper and lower values for the models given the noise model
# sigma_fac defaults to 3.
model_val = sedgrid.seds[indxs, k] + model_bias[indxs, k]
model_down = model_val - sigma_fac * model_unc[indxs, k]
model_up = model_val + sigma_fac * model_unc[indxs, k]
nindxs, = np.where((model_up >= min_data[k]) & (model_down <= max_data[k]))
if len(nindxs) > 0:
indxs = indxs[nindxs]
if len(indxs) == 0:
raise ValueError("no models that are within the data range")
print("number of original models = ", len(sedgrid.seds[:, 0]))
print("number of ast trimmed models = ", n_ast_indxs)
print("number of trimmed models = ", len(indxs))
# Save the grid
print("Writing trimmed sedgrid to disk into {0:s}".format(sed_outname))
cols = {}
for key in list(sedgrid.grid.keys()):
cols[key] = sedgrid.grid[key][indxs]
# New column to save the index of the model in the full grid
cols["fullgrid_idx"] = indxs.astype(int)
g = SpectralGrid(
sedgrid.lamb, seds=sedgrid.seds[indxs], grid=Table(cols), backend="memory"
)
filternames = obsdata.filters
g.grid.header["filters"] = " ".join(filternames)
# trimmed grid name
g.writeHDF(sed_outname)
# save the trimmed noise model
print("Writing trimmed noisemodel to disk into {0:s}".format(noisemodel_outname))
with tables.open_file(noisemodel_outname, "w") as outfile:
outfile.create_array(outfile.root, "bias", model_bias[indxs])
outfile.create_array(outfile.root, "error", model_unc[indxs])
outfile.create_array(outfile.root, "completeness", model_compl[indxs])
if trunchen:
outfile.create_array(outfile.root, "q_norm", model_q_norm[indxs])
outfile.create_array(outfile.root, "icov_diag", model_icov_diag[indxs])
outfile.create_array(
outfile.root, "icov_offdiag", model_icov_offdiag[indxs]
)