Source code for beast.fitting.trim_grid

import numpy as np
import tables

from beast.physicsmodel.grid import SEDGrid
from astropy.table 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.filter_aliases[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]) # link to the noisemodel values model_bias = sedgrid_noisemodel["bias"] model_unc = 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"] # has to be complete in all filters - otherwise observation model not defined # toothpick model means that if compl = 0, then bias = 0, and sigma = 0 from ASTs above_ast = model_compl > 0 sum_above_ast = np.sum(above_ast, axis=1) (indxs,) = np.where(sum_above_ast >= n_filters) n_ast_indxs = len(indxs) print("number of original models = ", len(sedgrid.seds[:, 0])) print("number of ast trimmed models = ", n_ast_indxs) if n_ast_indxs <= 0: raise ValueError("no models are brighter than the minimum ASTs run") # 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] # print(k, min(model_val), max(model_val), min(model_bias[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 = SEDGrid( sedgrid.lamb, seds=sedgrid.seds[indxs], grid=Table(cols), backend="memory" ) filternames = obsdata.filters g.header["filters"] = " ".join(filternames) # trimmed grid name g.write(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] )