diff --git a/rubin_sim/maf/maf_contrib/lss_metrics.py b/rubin_sim/maf/maf_contrib/lss_metrics.py index 287b1e83..ebb31786 100644 --- a/rubin_sim/maf/maf_contrib/lss_metrics.py +++ b/rubin_sim/maf/maf_contrib/lss_metrics.py @@ -4,7 +4,9 @@ import numpy as np import scipy -from rubin_sim.maf.metrics import BaseMetric, ExgalM5 +from rubin_sim.maf.metrics.base_metric import BaseMetric +from rubin_sim.maf.metrics.simple_metrics import Coaddm5Metric +from rubin_sim.maf.metrics.exgal_m5 import ExgalM5 class GalaxyCountsMetric(BaseMetric): diff --git a/rubin_sim/maf/maf_contrib/lss_obs_strategy/galaxy_counts_metric_extended.py b/rubin_sim/maf/maf_contrib/lss_obs_strategy/galaxy_counts_metric_extended.py index 894ec137..ba57149f 100644 --- a/rubin_sim/maf/maf_contrib/lss_obs_strategy/galaxy_counts_metric_extended.py +++ b/rubin_sim/maf/maf_contrib/lss_obs_strategy/galaxy_counts_metric_extended.py @@ -29,7 +29,9 @@ power_law_const_a, power_law_const_b, ) -from rubin_sim.maf.metrics import BaseMetric, Coaddm5Metric, ExgalM5 +from rubin_sim.maf.metrics.base_metric import BaseMetric +from rubin_sim.maf.metrics.simple_metrics import Coaddm5Metric +from rubin_sim.maf.metrics.exgal_m5 import ExgalM5 class GalaxyCountsMetricExtended(BaseMetric): diff --git a/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py b/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py index 8cf8f3e4..b0de8383 100644 --- a/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py +++ b/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py @@ -14,7 +14,11 @@ from rubin_scheduler.data import get_data_dir from rubin_sim.maf.maf_contrib.lss_obs_strategy import GalaxyCountsMetricExtended -from rubin_sim.maf.metrics import BaseMetric, ExgalM5, StarDensityMetric + +from rubin_sim.maf.metrics.base_metric import BaseMetric +from rubin_sim.maf.metrics.simple_metrics import Coaddm5Metric +from rubin_sim.maf.metrics.exgal_m5 import ExgalM5 +from rubin_sim.maf.metrics.star_density import StarDensityMetric from rubin_sim.maf.slicers import UserPointsSlicer diff --git a/rubin_sim/maf/maf_contrib/static_probes_fom_summary_metric.py b/rubin_sim/maf/maf_contrib/static_probes_fom_summary_metric.py index 83793b90..e215b641 100644 --- a/rubin_sim/maf/maf_contrib/static_probes_fom_summary_metric.py +++ b/rubin_sim/maf/maf_contrib/static_probes_fom_summary_metric.py @@ -34,6 +34,8 @@ def __init__( super().__init__(col=col, **kwargs) if col is None: self.col = "metricdata" + else: + self.col = col self.shear_m = shear_m self.sigma_z = sigma_z self.sig_delta_z = sig_delta_z diff --git a/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py b/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py index de6ad6bb..fe55df64 100644 --- a/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py +++ b/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py @@ -9,7 +9,8 @@ import scipy.integrate as integrate from rubin_sim.maf.maps import DustMap, DustMap3D, StellarDensityMap -from rubin_sim.maf.metrics import BaseMetric, CrowdingM5Metric +from rubin_sim.maf.metrics.base_metric import BaseMetric +from rubin_sim.maf.metrics.crowding_metric import CrowdingM5Metric from rubin_sim.phot_utils import DustValues diff --git a/rubin_sim/maf/metrics/__init__.py b/rubin_sim/maf/metrics/__init__.py index 3c39eb15..f0683ff4 100644 --- a/rubin_sim/maf/metrics/__init__.py +++ b/rubin_sim/maf/metrics/__init__.py @@ -6,6 +6,7 @@ from .cadence_metrics import * from .calibration_metrics import * from .chip_vendor_metric import * +from .cosmology_summary_metrics import * from .coverage_metric import * from .crowding_metric import * from .cumulative_metric import * @@ -40,8 +41,10 @@ from .surfb_metric import * from .technical_metrics import * from .tgaps import * +from .tomography_models import * from .transient_metrics import * from .use_metrics import * +from .uniformity_metrics import * from .vector_metrics import * from .visit_groups_metric import * from .weak_lensing_systematics_metric import * diff --git a/rubin_sim/maf/metrics/area_summary_metrics.py b/rubin_sim/maf/metrics/area_summary_metrics.py index 1890953c..12e4a643 100644 --- a/rubin_sim/maf/metrics/area_summary_metrics.py +++ b/rubin_sim/maf/metrics/area_summary_metrics.py @@ -96,7 +96,7 @@ def __init__( self.lower_threshold = lower_threshold self.mask_val = np.nan # Include so all values get passed self.col = col - self.units = "degrees" + self.units = "square degrees" def run(self, data_slice, slice_point=None): # find out what nside we have @@ -113,7 +113,7 @@ def run(self, data_slice, slice_point=None): npix = len( np.where( (data_slice[self.col] > self.lower_threshold) - and (data_slice[self.col] < self.upper_threshold) + & (data_slice[self.col] < self.upper_threshold) )[0] ) area = pix_area * npix diff --git a/rubin_sim/maf/metrics/cosmology_summary_metrics.py b/rubin_sim/maf/metrics/cosmology_summary_metrics.py new file mode 100644 index 00000000..709f125c --- /dev/null +++ b/rubin_sim/maf/metrics/cosmology_summary_metrics.py @@ -0,0 +1,1140 @@ +__all__ = ( + "TotalPowerMetric", + "StaticProbesFoMEmulatorMetricSimple", + "TomographicClusteringSigma8biasMetric", + "MultibandMeanzBiasMetric", +) + +import warnings + +import healpy as hp +import numpy as np +from scipy import interpolate +from scipy.stats import median_abs_deviation +from astropy import units as u +from astropy.coordinates import SkyCoord +from sklearn.cluster import KMeans +from copy import deepcopy + +from ..maf_contrib.static_probes_fom_summary_metric import StaticProbesFoMEmulatorMetric +from .area_summary_metrics import AreaThresholdMetric +from .base_metric import BaseMetric +from .simple_metrics import RmsMetric + + +# Cosmology-related summary metrics. +# These generally calculate a FoM for various DESC metrics. + + +class TotalPowerMetric(BaseMetric): + """Calculate the total power in the angular power spectrum, + between lmin/lmax. + + Parameters + ---------- + lmin : `float`, optional + Minimum ell value to include when calculating total power. + lmax : `float`, optional + Maximum ell value to include when calculating total power. + remove_monopole : `bool`, optional + Flag to remove monopole when calculating total power. + remove_dipole : `bool`, optional + Flag to remove dipole when calculating total power. + col : `str`, optional + The column name to operate on. + For summary metrics, this is almost always `metricdata`. + mask_val : `float` or np.nan, optional + The mask value to apply to the metric values when passed. + If this attribute exists, the metric_values will be passed + using metric_values.filled(mask_val). + If mask_val is `None` for a metric, metric_values will be passed + using metric_values.compressed(). + """ + + def __init__( + self, + lmin=100.0, + lmax=300.0, + remove_monopole=True, + remove_dipole=True, + col="metricdata", + mask_val=np.nan, + **kwargs, + ): + self.lmin = lmin + self.lmax = lmax + self.remove_monopole = remove_monopole + self.remove_dipole = remove_dipole + super().__init__(col=col, mask_val=mask_val, **kwargs) + + def run(self, data_slice, slice_point=None): + # Calculate the power spectrum. + data = data_slice[self.colname] + if self.remove_monopole: + data = hp.remove_monopole(data, verbose=False, bad=self.mask_val) + if self.remove_dipole: + data = hp.remove_dipole(data, verbose=False, bad=self.mask_val) + cl = hp.anafast(data) + ell = np.arange(np.size(cl)) + condition = np.where((ell <= self.lmax) & (ell >= self.lmin))[0] + totalpower = np.sum(cl[condition] * (2 * ell[condition] + 1)) + return totalpower + + +class StaticProbesFoMEmulatorMetricSimple(BaseMetric): + """Calculate the FoM for the combined static probes + (3x2pt, i.e. Weak Lensing, LSS, Clustering). + + Parameters + ---------- + year : `int`, optional + The year of the survey to calculate FoM. + This calibrates expected depth and area. + + Returns + ------- + result : `float` + The simple 3x2pt FoM emulator value, for the + years where the correlation between area/depth and value is defined. + + Notes + ----- + This FoM is purely statistical and does not factor in systematics. + The implementation here is simpler than in + `rubin_sim.maf.mafContrib.StaticProbesFoMEmulatorMetric`, and that + more sophisticated version should replace this metric. + + This version of the emulator was used to generate the results in + https://ui.adsabs.harvard.edu/abs/2018arXiv181200515L/abstract + + Note that this is truly a summary metric and should be run on the + output of Exgalm5_with_cuts. + """ + + def __init__(self, year=10, **kwargs): + self.year = year + super().__init__(col="metricdata", mask_val=-666, **kwargs) + + def run(self, data_slice, slice_point=None): + # derive nside from length of data slice + nside = hp.npix2nside(len(data_slice)) + pix_area = hp.nside2pixarea(nside, degrees=True) + + # Chop off any outliers (and also the masked value) + good_pix = np.where(data_slice[self.col] > 0)[0] + + # Calculate area and med depth from + area = pix_area * np.size(good_pix) + median_depth = np.median(data_slice[self.col][good_pix]) + + # FoM is calculated at the following values + if self.year == 1: + areas = [7500, 13000, 16000] + depths = [24.9, 25.2, 25.5] + fom_arr = [ + [1.212257e02, 1.462689e02, 1.744913e02], + [1.930906e02, 2.365094e02, 2.849131e02], + [2.316956e02, 2.851547e02, 3.445717e02], + ] + elif self.year == 3: + areas = [10000, 15000, 20000] + depths = [25.5, 25.8, 26.1] + fom_arr = [ + [1.710645e02, 2.246047e02, 2.431472e02], + [2.445209e02, 3.250737e02, 3.516395e02], + [3.173144e02, 4.249317e02, 4.595133e02], + ] + + elif self.year == 6: + areas = [10000, 15000, 2000] + depths = [25.9, 26.1, 26.3] + fom_arr = [ + [2.346060e02, 2.414678e02, 2.852043e02], + [3.402318e02, 3.493120e02, 4.148814e02], + [4.452766e02, 4.565497e02, 5.436992e02], + ] + + elif self.year == 10: + areas = [10000, 15000, 20000] + depths = [26.3, 26.5, 26.7] + fom_arr = [ + [2.887266e02, 2.953230e02, 3.361616e02], + [4.200093e02, 4.292111e02, 4.905306e02], + [5.504419e02, 5.624697e02, 6.441837e02], + ] + else: + warnings.warn("FoMEmulator is not defined for this year") + return self.badval + + # Interpolate FoM to the actual values for this sim + areas = [[i] * 3 for i in areas] + depths = [depths] * 3 + f = interpolate.interp2d(areas, depths, fom_arr, bounds_error=False) + fom = f(area, median_depth)[0] + return fom + + +class TomographicClusteringSigma8biasMetric(BaseMetric): + """Compute bias on sigma8 due to spurious contamination of density maps. + Run as summary metric on NestedLinearMultibandModelMetric. + + Points of contact / contributors: Boris Leistedt + + Parameters + ---------- + density_tomograph_model : `dict` + dictionary containing models calculated for fiducial N(z)s and Cells: + lmax : numpy.array of int, of shape (Nbins, ) + lmax corresponding to kmax of 0.05 + poly1d_coefs_loglog : numpy.array of float, of shape (Nbins, ) + polynomial fits to log(C_ell) vs log(ell) computed for CCL + sigma8square_model (float) + value of sigma8^2 used as fiducal model for CCL + power_multiplier : `float`, optional + fraction of power (variance) which is uncorrected + and thus biases sigma8 + lmin : `int`, optional + lmin for the analysis + convert_to_sigma8 : `str`, optional + Convert the bias to sigma8 instead of sigma8^2 + (via change of variables for the uncertainty) + + Returns + ------- + result : `float` + Value of sigma8 bias calculated from this model: + (sigma8^2_obs - sigma^2_model) / error on sigma8^2_obs + if `convert_to_sigma8` is True, + then it is about sigma8 instead of sigma8^2. + + Notes + ----- + This is a summary metric to be run on the results + of the NestedLinearMultibandModelMetric. + + NestedLinearMultibandModelMetric converts 6-band depth maps into + a set of maps (e.g tomographic redshift bins) which describe + spurious density fluctuations in each bin. + + This summary metric multiplies the maps by the parameter power_multiplier, + which can be used to describe the fraction of power uncorrected by + systematics mitigation schemes, computes the total power + (via angular power spectra with lmin-lmax limits) + and then infers sigma8^2 via a model of the angular power spectra. + By taking sigma8_obs minus sigma8_model divided by the uncertainty, + one derives a bias. + """ + + def __init__( + self, + density_tomography_model, + power_multiplier=0.1, + lmin=10, + convert_to_sigma8=True, + **kwargs, + ): + super().__init__(col="metricdata", **kwargs) + # Set mask_val, so that we receive metric_values.filled(mask_val) + self.mask_val = hp.UNSEEN + self.badval = hp.UNSEEN + + self.convert_to_sigma8 = convert_to_sigma8 + + self.power_multiplier = power_multiplier + self.lmin = lmin + self.density_tomography_model = density_tomography_model + # to compute angular power spectra and total power, + # initialize an array of metrics, with the right lmin and lmax. + self.totalPowerMetrics = [ + TotalPowerMetric(lmin=lmin, lmax=lmax, mask_val=self.mask_val) + for lmax in density_tomography_model["lmax"] + ] + self.areaThresholdMetric = AreaThresholdMetric( + lower_threshold=hp.UNSEEN, + upper_threshold=np.inf, + mask_val=self.mask_val, + ) + + def run(self, data_slice, slice_point=None): + # need to define an array of bad values for the masked pixels + badval_arr = np.repeat(self.badval, len(self.density_tomography_model["lmax"])) + # converts the input recarray to an array + data_slice_list = [ + badval_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + # should be (nbins, npix) + data_slice_arr = np.asarray(data_slice_list, dtype=float).T + # hp.mollview(data_slice_arr[0]) + ### data_slice_arr[data_slice_arr == -666] = hp.UNSEEN ############ + data_slice_arr[~np.isfinite(data_slice_arr)] = ( + hp.UNSEEN + ) # need to work with TotalPowerMetric and healpix + + # measure valid sky fractions and total power + # (via angular power spectra) in each bin. + # The original metric returns an array at each slice_point (of the + # original slicer) -- so there is a bit of "rearrangement" that + # has to happen to be able to pass a np.array with right dtype + # (i.e. dtype = [("metricdata", float)]) to each call to + # the AreaThresholdMetric and TotalPowerMetric `run` methods. + totalsky = 42000 + fskys = np.array( + [ + self.areaThresholdMetric.run(np.core.records.fromrecords(x, dtype=[("metricdata", float)])) + / totalsky + for x in data_slice_arr + ] + ) # sky fraction + spuriousdensitypowers = ( + np.array( + [ + self.totalPowerMetrics[i].run( + np.core.records.fromrecords(x, dtype=[("metricdata", float)]) + ) + for i, x in enumerate(data_slice_arr) + ] + ) + / fskys + ) + print("spuriousdensitypowers:", spuriousdensitypowers) + print("fskys:", fskys) + + def solve_for_multiplicative_factor(spurious_powers, model_cells, fskys, lmin, power_multiplier): + """ + Infer multiplicative factor sigma8^2 (and uncertainty) + from the model Cells and observed total powers + since it os a Gaussian posterior distribution. + """ + # solve for multiplicative sigma8^2 term between + # measured angular power spectra + # (spurious measured Cells times power_multiplier) + # and model ones (polynomial model from CCL). + n_bins = model_cells["lmax"].size + assert len(spurious_powers) == n_bins + assert len(fskys) == n_bins + assert model_cells["poly1d_coefs_loglog"].shape[0] == n_bins + totalvar_mod = np.zeros((n_bins, 1)) + totalvar_obs = np.zeros((n_bins, 1)) + totalvar_var = np.zeros((n_bins, 1)) + # loop over tomographic bins + # hardcoded; assumed CCL cosmology + sigma8square_model = model_cells["sigma8square_model"] + for i in range(n_bins): + # get model Cells from polynomial model (in log log space) + ells = np.arange(lmin, model_cells["lmax"][i]) + polynomial_model = np.poly1d(model_cells["poly1d_coefs_loglog"][i, :]) + cells_model = np.exp(polynomial_model(np.log(ells))) + + # model variance is sum of cells x (2l+1) + totalvar_mod[i, 0] = np.sum(cells_model * (2 * ells + 1)) + + # observations is spurious power noiseless model + totalvar_obs[i, 0] = totalvar_mod[i, 0] + spurious_powers[i] * power_multiplier + + # simple model variance of cell based on Gaussian covariance + cells_var = 2 * cells_model ** 2 / (2 * ells + 1) / fskys[i] + totalvar_var[i, 0] = np.sum(cells_var * (2 * ells + 1) ** 2) + + # model assumed sigma8 = 0.8 + # (add CCL cosmology here? or how I obtained them + documentation) + # results_fractional_spurious_power = + # totalvar_obs / totalvar_mod - 1.0 + + # model Cell variance divided by sigma8^2, + # which is the common normalization + transfers = totalvar_mod / sigma8square_model + + # model ratio: formula for posterior distribution on unknown + # multiplicative factor in multivariate Gaussian likelihood + FOT = np.sum(transfers[:, 0] * totalvar_obs[:, 0] / totalvar_var[:, 0]) + FTT = np.sum(transfers[:, 0] * transfers[:, 0] / totalvar_var[:, 0]) + # mean and stddev of multiplicative factor + sigma8square_fit = FOT / FTT + sigma8square_error = FTT ** -0.5 + + return sigma8square_fit, sigma8square_error, sigma8square_model + + # solve for the gaussian posterior distribution on sigma8^2 + sigma8square_fit, sigma8square_error, sigma8square_model = solve_for_multiplicative_factor( + spuriousdensitypowers, self.density_tomography_model, fskys, self.lmin, self.power_multiplier + ) + + results_sigma8_square_bias = (sigma8square_fit - sigma8square_model) / sigma8square_error + if not self.convert_to_sigma8: + return results_sigma8_square_bias + + else: + + # turn result into bias on sigma8, + # via change of variable and simple propagation of uncertainty. + sigma8_fit = sigma8square_fit ** 0.5 + sigma8_model = sigma8square_model ** 0.5 + sigma8_error = 0.5 * sigma8square_error * sigma8_fit / sigma8square_fit + results_sigma8_bias = (sigma8_fit - sigma8_model) / sigma8_error + print(sigma8square_model, sigma8square_fit, sigma8square_error, results_sigma8_square_bias) + print(sigma8_model, sigma8_fit, sigma8_error, results_sigma8_bias) + return results_sigma8_bias + + + +class MultibandMeanzBiasMetric(BaseMetric): + """ + Run as summary metric on MultibandExgalM5. + + Parameters + ---------- + year : `int`, optional + The year of the survey to calculate the bias. + This is used to derive the dm/dz derivative used to translate m5 rms into dz rms. + + meanz_tomograph_model : `dict` + dictionary containing models calculated for fiducial N(z): + + meanz: numpy.float + the meanz within a tomographic bin at a given band + dz_dm5: numpy.float + the absolute value of the derivative of z in a bin as a function + of the depth, m5 magnitude. This is later used to translate fluctuations + in depth to fluctuations in tomographic redshift + + Returns + ------- + result : `float` array + The ratio of this bias to the desired DESC y1 upper bound on the bias, and the ratio + between the clbias and the y10 DESC SRD requirement. + Desired values are less than 1 by Y10. + + Notes + ----- + This is a summary metric to be run on the results + of the MultibandExgalM5. + + MultibandExgalM5 provides the m5 depth in all LSST bands given a specific slice. + + This summary metric takes those depths and reads the derivatives from the tomogrpahic model, + computes the bias in shear signal Cl and then computes the ratio of that bias to the Y1 goal and + Y10 science requirement. + """ + + def __init__(self, + meanz_tomography_model, + filter_list=["u", "g", "r", "i", "z", "y"], + year=10, n_filters=6, **kwargs): + + super().__init__(col="metricdata", **kwargs) + # Set mask_val, so that we receive metric_values.filled(mask_val) + self.mask_val = hp.UNSEEN + self.badval = hp.UNSEEN + self.rmsMetric = RmsMetric() + self.year = year + self.filter_list = filter_list + + def run(self, data_slice, slice_point=None): + + def compute_dzfromdm(zbins, model_z, band_ind, year): + + """This computes the dm/dz relationship calibrated from simulations + by Jeff Newmann and Qianjun Hang, which forms the meanz_tomographic_model. + + Parameters + ---------- + zbins : `int` + The number of tomographic bins considered. For now this is zbins < 5 + filter : `str` + The assumed filter band + model_z: dict + The meanz_tomographic_model assumed in this work + + Returns + ------- + dzdminterp : `float` array + The interpolated value of the derivative dz/dm + meanzinterp : `float` array + The meanz in each tomographic bin. + + """ + import pandas as pd + + filter_list = ["u", "g", "r", "i", "z", "y"] + meanzinterp=np.zeros(zbins) + dzdminterp=np.zeros(zbins) + + for z in range(zbins): + meanzinterp[z]=model_z["year%i"%(year+1)]["meanz"][z][filter] + dzdminterp[z] =model_z["year%i"%(year+1)]["dz_dm5"][z][filter] + + return dzdminterp, meanzinterp + + ## Not entirely sure if we should include these here or in a separate module/file + def use_zbins(meanz_vals, figure_9_mean_z=np.array([0.2, 0.4, 0.7, 1.0]), figure_9_width=0.2): + """This computes which redshift bands are within the range + specified in https://arxiv.org/pdf/2305.15406.pdf and can safely be used + to compute what Cl bias result from z fluctuations caused by rms variations in the m5. + + + Parameters + ---------- + meanz_vals : `float` array + Array of meanz values to be used. + + Returns + ------- + use_bins : `boolean` array + An array of boolean values of length meanz_vals + + """ + max_z_use = np.max(figure_9_mean_z) + 2 * figure_9_width + use_bins = meanz_vals < max_z_use + + return use_bins + + def compute_Clbias(meanz_vals, scatter_mean_z_values): + """This computes the Cl bias + that results z fluctuations caused by rms variations in the m5. + + + + Parameters + ---------- + meanz_vals : `float` array + Array of meanz values to be used. + + scatter_mean_z_values : `float` array + Array of rms values of the z fluctuations + + + Returns + ------- + clbiasvals : `float` array + An array of values of the clbias + + mean_z_values_use : `float` array + An array of the meanz values that are within the interpolation range of 2305.15406 + + Notes + ------ + This interpolates from the Figure 9 in https://arxiv.org/pdf/2305.15406.pdf + + """ + import numpy as np + + figure_9_mean_z = np.array([0.2, 0.4, 0.7, 1.0]) + figure_9_Clbias = np.array([1e-3, 2e-3, 5e-3, 1.1e-2]) + figure_9_width = 0.2 + figure_9_mean_z_scatter = 0.02 + + mzvals = np.array([float(mz) for mz in meanz_vals]) + sctz = np.array([float(sz) for sz in scatter_mean_z_values]) + + fit_res = np.polyfit(figure_9_mean_z, figure_9_Clbias, 2) + poly_fit = np.poly1d(fit_res) + use_bins = use_zbins(meanz_vals, figure_9_mean_z, figure_9_width) + + mean_z_values_use = mzvals[use_bins] + sctz_use = sctz[use_bins] + + Clbias = poly_fit(mean_z_values_use) + rescale_fac = sctz_use / figure_9_mean_z_scatter + Clbias *= rescale_fac + fit_res_bias = np.polyfit(mean_z_values_use, Clbias, 1) + poly_fit_bias = np.poly1d(fit_res_bias) + + clbiasvals = poly_fit_bias(mean_z_values_use) + return clbiasvals, mean_z_values_use + + result = np.empty(1, dtype=[("name", np.str_, 20), ("y1ratio", float), ("y10ratio", float)]) + result["name"][0] = "MultibandMeanzBiasMetric" + + # Technically don't need this for now (isn't used in previous one) + # need to define an array of bad values for the masked pixels + badval_arr = np.repeat(self.badval, len(self.filter_list)) + # converts the input recarray to an array + data_slice_list = [ + badval_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + lengths = np.array([len(x) for x in data_slice_list]) + # should be (nbins, npix) + data_slice_arr = np.asarray(data_slice_list, dtype=float).T + data_slice_arr[~np.isfinite(data_slice_arr)] = ( + hp.UNSEEN + ) # need to work with TotalPowerMetric and healpix + + # measure rms in each bin. + # The original metric returns an array at each slice_point (of the + # original slicer) -- so there is a bit of "rearrangement" that + # has to happen to be able to pass a np.array with right dtype + # (i.e. dtype = [("metricdata", float)]) to each call to + # the rmsMetric `run` methods. + rmsarray = np.array( + [ + self.rmsMetric.run(np.core.records.fromrecords(x, dtype=[("metricdata", float)])) + for x in data_slice_arr + ] + ) # rms values + + totdz = 0 + avmeanz = 0 + totclbias = 0 + + for i, filt in enumerate(self.filter_list): + dzdminterp, meanzinterp = compute_dzfromdm(self.zbins, filt, self.year) + stdz = [float(np.abs(dz)) * float(rmsarray[i]) for dz in dzdminterp] + + clbias, meanz_use = compute_Clbias(meanzinterp, stdz) + + totdz += [float(st ** 2) for st in stdz] + totclbias += clbias + avmeanz += meanzinterp + + # These should be included in self rather than hard coded + y10_req = 0.003 + y1_goal = 0.013 + + # clbiastot = np.max(clbias) # if adding doesn't work over z range -- CHECK + y10ratio = totclbias / y10_req + y1ratio = totclbias / y1_goal + + result["y1ratio"] = y1ratio + result["y10ratio"] = y10ratio + + +# Let's make code that pulls out the northern/southern galactic regions, and gets statistics of the footprint by region. +def _is_ngp(ra, dec): + c = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame="icrs") + lat = c.galactic.b.deg + return lat >= 0 + + +def get_stats_by_region(use_map, nside, maskval=0, region="all"): + if region not in ["all", "north", "south"]: + raise ValueError("Invalid region %s" % region) + to_use = use_map > maskval + + if region != "all": + # Find the north/south part of the map as requested + npix = hp.nside2npix(nside) + ra, dec = hp.pix2ang(hp.npix2nside(npix), range(npix), lonlat=True) + ngp_mask = _is_ngp(ra, dec) + if region == "north": + reg_mask = ngp_mask + else: + reg_mask = ~ngp_mask + to_use = to_use & reg_mask + + # Calculate the desired stats + reg_mad = median_abs_deviation(use_map[to_use]) + reg_median = np.median(use_map[to_use]) + reg_std = np.std(use_map[to_use]) + + # Return the values + return (reg_mad, reg_median, reg_std) + + +def stripiness_test_statistic(data_slice, nside): + """ + A utility to find whether a particular routine has stripey features in the exposure time map. + """ + # Analyze the exposure time map to get MAD, median, std in north/south. + mad = {} + med = {} + frac_scatter = {} + regions = ["north", "south"] + for region in regions: + mad[region], med[region], _ = get_stats_by_region(data_slice, nside, region=region) + frac_scatter[region] = mad[region] / med[region] + test_statistic = frac_scatter["north"] / frac_scatter["south"] - 1 + return test_statistic + + +class StripinessMetric(BaseMetric): + """ + Run as summary metric on NestedRIZExptimeExgalM5Metric. + + This metric uses maps of the combined RIZ exposure time to identify stripes or residual rolling features, + for the UniformAreaFoMFractionMetric. + The number returned quantifies the stripiness of the field. + In UniformAreaFoMFractionMetric it is a test statistic: if it is outside of +-0.7/np.sqrt(year) then + the RIZ exposure time map is considered to have stripes. + + Points of contact / contributors: Rachel Mandelbaum, Boris Leistedt + + Parameters + ---------- + year: `int` + year of observation, in order to adjust the cut on the depth fluctuations for the stipe detection + nside: `int` + must be the nside at which the base metric NestedRIZExptimeExgalM5Metric is calculated at + verbose: bool, optional + if true, will display the segmentation maps and the areas and FOMs of the two regions found + Returns + ------- + result: `float` + A number quantifying the stripiness of the RIZ exposure time. + """ + + def __init__( + self, + year, + nside, + verbose=True, + **kwargs, + ): + self.year = year + self.mask_val = hp.UNSEEN + self.verbose = verbose + self.nside = nside + names = ["exgal_m5", "riz_exptime"] + types = [float] * 2 + self.mask_val_arr = np.zeros(1, dtype=list(zip(names, types))) + self.mask_val_arr["exgal_m5"] = self.mask_val + self.mask_val_arr["riz_exptime"] = self.mask_val + + super().__init__(col="metricdata", **kwargs) + + def run(self, data_slice, slice_point=None): + data_slice_list = [ + self.mask_val_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + data_slice_arr = np.asarray(data_slice_list, dtype=self.mask_val_arr.dtype) + # a bit of gymnastics to make sure all bad values (nan, -666) are recast as hp.UNSEEN + ind = data_slice_arr["riz_exptime"] == -666 + ind |= ~np.isfinite(data_slice_arr["riz_exptime"]) + data_slice_arr["riz_exptime"][ind.ravel()] = hp.UNSEEN + # sanity check + nside = hp.npix2nside(data_slice["metricdata"].size) + assert nside == self.nside + + test_statistic = stripiness_test_statistic(data_slice_arr["riz_exptime"].ravel(), nside) + + return test_statistic + + +class UniformAreaFoMFractionMetric(BaseMetric): + """ + Run as summary metric on NestedRIZExptimeExgalM5Metric. + + This metric uses maps of the combined RIZ exposure time and i-band depth maps (with a consistent set of area cuts) + to identify potential reductions in cosmological constraining power due to substantial large-scale power + in non-uniform coadds at a particular data release. + The RIZ exposure time map is used to identify whether there are residual rolling features. + Under the hood, this runs the StripinessMetric. + If not, the metric returns 1. If there are such features, then the region is segmented into similar-depth regions + and the one with the largest cosmological constraining power is presumed to be used for science. + In that case, the metric returns the 3x2pt FoM (StaticProbesFoMEmulatorMetric, + quantifying weak lensing and large-scale structure constraining power) for the largest of those regions, + divided by that for the full region if it had been usable. + + + + Points of contact / contributors: Rachel Mandelbaum, Boris Leistedt + + Parameters + ---------- + year: `int` + year of observation, in order to adjust the cut on the depth fluctuations for the stipe detection + nside: `int` + must be the nside at which the base metric NestedRIZExptimeExgalM5Metric is calculated at + verbose: bool, optional + if true, will display the segmentation maps and the areas and FOMs of the two regions found + Returns + ------- + result: `float` + The ratio of the FOMs of the two areas found + """ + + def __init__( + self, + year, + nside, + verbose=True, + **kwargs, + ): + self.year = year + self.mask_val = hp.UNSEEN + self.verbose = verbose + self.nside = nside + names = ["exgal_m5", "riz_exptime"] + types = [float] * 2 + self.mask_val_arr = np.zeros(1, dtype=list(zip(names, types))) + self.mask_val_arr["exgal_m5"] = self.mask_val + self.mask_val_arr["riz_exptime"] = self.mask_val + + self.threebyTwoSummary = StaticProbesFoMEmulatorMetric( + nside=nside, metric_name="3x2ptFoM", col="exgal_m5" + ) + super().__init__(col="metricdata", **kwargs) + + def run(self, data_slice, slice_point=None): + data_slice_list = [ + self.mask_val_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + data_slice_arr = np.asarray(data_slice_list, dtype=self.mask_val_arr.dtype) + # a bit of gymnastics to make sure all bad values (nan, -666) are recast as hp.UNSEEN + ind = data_slice_arr["riz_exptime"] == -666 + ind |= ~np.isfinite(data_slice_arr["riz_exptime"]) + ind = data_slice_arr["exgal_m5"] == -666 + ind |= ~np.isfinite(data_slice_arr["exgal_m5"]) + data_slice_arr["exgal_m5"][ind.ravel()] = hp.UNSEEN + data_slice_arr["riz_exptime"][ind.ravel()] = hp.UNSEEN + # sanity check + nside = hp.npix2nside(data_slice["metricdata"].size) + assert nside == self.nside + + # Check for stripiness + threshold = 0.7 / np.sqrt(self.year) + test_statistic = stripiness_test_statistic(data_slice_arr["riz_exptime"].ravel(), nside) + + if np.abs(test_statistic) < np.abs(threshold): + stripes = False + else: + stripes = True + + def apply_clustering(clustering_data): + # A thin wrapper around sklearn routines (can swap out which one we are using systematically). + # We fix parameters like `n_clusters` since realistically we know for rolling that we should expect 2 clusters. + # from sklearn.cluster import SpectralClustering + # clustering = SpectralClustering(n_clusters=2, assign_labels='discretize', random_state=0).fit(clustering_data) + + clustering = KMeans(n_clusters=2, random_state=0, n_init="auto").fit(clustering_data) + labels = clustering.labels_ + 1 + return labels + + def expand_labels(depth_map, labels, maskval=0): + # A utility to apply the labels from a masked version of a depth map back to the entire depth map. + expanded_labels = np.zeros(hp.nside2npix(nside)) + cutval = maskval + 0.1 + expanded_labels[depth_map > cutval] = labels + return expanded_labels + + def get_area_stats(depth_map, labels, maskval=0, n_clusters=2): + # A routine to get some statistics of the clustering: area fractions, median map values + expanded_labels = expand_labels(depth_map, labels, maskval=maskval) + cutval = maskval + 0.1 + area_frac = [] + med_val = [] + for i in range(n_clusters): + new_map = depth_map.copy() + new_map[expanded_labels != i + 1] = maskval + this_area_frac = len(new_map[new_map > cutval]) / len(depth_map[depth_map > cutval]) + this_med_val = np.median(new_map[new_map > cutval]) + area_frac.append(this_area_frac) + med_val.append(this_med_val) + return area_frac, med_val + + def show_clusters(depth_map, labels, maskval=0, n_clusters=2, min=500, max=3000): + # A routine to show the clusters found by the unsupervised clustering algorithm (start with original map then 2 clusters). + expanded_labels = expand_labels(depth_map, labels, maskval=maskval) + hp.visufunc.mollview(depth_map, min=min, max=max) + for i in range(n_clusters): + new_map = depth_map.copy() + new_map[expanded_labels != i + 1] = maskval + hp.visufunc.mollview(new_map, min=min, max=max) + return get_area_stats(depth_map, labels, maskval=maskval, n_clusters=n_clusters) + + def make_clustering_dataset(depth_map, maskval=0, priority_fac=0.9, nside=64): + # A utility routine to get a dataset for unsupervised clustering. Note: + # - We want the unmasked regions of the depth map only. + # - We assume masked regions are set to `maskval`, and cut 0.1 magnitudes above that. + # - We really want it to look at depth fluctuations. So, we have to rescale the + # RA/dec dimensions to avoid them being prioritized because their values are larger and + # have more variation than depth. Currently we rescale RA/dec such that their + # standard deviations are 1-priority_fac times the standard deviation of the depth map. + # That's why priority_fac is a tunable parameter; it should be between 0 and 1 + if priority_fac < 0 or priority_fac >= 1: + raise ValueError("priority_fac must lie between 0 and 1") + theta, phi = hp.pixelfunc.pix2ang(nside, ipix=np.arange(hp.nside2npix(nside))) + # theta is 0 at the north pole, pi/2 at equator, pi at south pole; phi maps to RA + ra = np.rad2deg(phi) + dec = np.rad2deg(0.5 * np.pi - theta) + + # Make a 3D numpy array containing the unmasked regions, including a rescaling factor to prioritize the depth + n_unmasked = len(depth_map[depth_map > 0.1]) + my_data = np.zeros((n_unmasked, 3)) + cutval = 0.1 + maskval + my_data[:, 0] = ( + ra[depth_map > cutval] + * (1 - priority_fac) + * np.std(depth_map[depth_map > cutval]) + / np.std(ra[depth_map > cutval]) + ) + my_data[:, 1] = ( + dec[depth_map > cutval] + * (1 - priority_fac) + * np.std(depth_map[depth_map > cutval]) + / np.std(dec[depth_map > cutval]) + ) + my_data[:, 2] = depth_map[depth_map > cutval] + return my_data + + if not stripes: + return 1 + else: + # Do the clustering if we got to this point + if self.verbose: + print("Verbose mode - Carrying out the clustering exercise for this map") + clustering_data = make_clustering_dataset(data_slice_arr["riz_exptime"].ravel(), nside=nside) + labels = apply_clustering(clustering_data) + area_frac, med_val = get_area_stats(data_slice_arr["riz_exptime"].ravel(), labels) + if self.verbose: + print("Verbose mode - showing original map and clusters identified for this map") + show_clusters(data_slice_arr["riz_exptime"].ravel(), labels) + print("Area fractions", area_frac) + print("Median exposure time values", med_val) + print("Median exposure time ratio", np.max(med_val) / np.min(med_val)) + print("Verbose mode - proceeding with area cuts") + # Get the FoM without/with cuts. We want to check the FoM for each area, if we're doing cuts, and + # return the higher one. This will typically be for the larger area, but not necessarily, if the smaller area + # is deeper. + expanded_labels = expand_labels(data_slice_arr["riz_exptime"].ravel(), labels) + my_hpid_1 = expanded_labels == 1 # np.where(expanded_labels == 1)[0] + my_hpid_2 = expanded_labels == 2 # np.where(expanded_labels == 2)[0] + # should this be labels or expanded_labels?! + + # copies need to be recarrays + data_slice_subset_2 = deepcopy(data_slice_arr) + data_slice_subset_1 = deepcopy(data_slice_arr) + data_slice_subset_1[my_hpid_2] = self.mask_val_arr + data_slice_subset_2[my_hpid_1] = self.mask_val_arr + fom1 = self.threebyTwoSummary.run(data_slice_subset_1) + fom2 = self.threebyTwoSummary.run(data_slice_subset_2) + fom = np.max((fom1, fom2)) + fom_total = self.threebyTwoSummary.run(data_slice_arr) + if self.verbose: + print("FOMs:", fom1, fom2, fom, fom_total) + return fom / fom_total + + +class uDropoutsNumbersShallowestDepth(BaseMetric): + """ + Run as summary metric on MultibandExgalM5. + + + + Parameters + ---------- + year : `int`, optional + The year of the survey to calculate the bias. + This is used to derive the dm/dz derivative used to translate m5 rms into dz rms. + + Returns + ------- + result : `float` array + The ratio of this bias to the desired DESC y1 upper bound on the bias, and the ratio + between the clbias and the y10 DESC SRD requirement. + Desired values are less than 1 by Y10. + + + Notes + ----- + This is a summary metric to be run on the results + of the MultibandExgalM5. + + MultibandExgalM5 provides the m5 depth in all LSST bands given a specific slice. + + """ + + def __init__(self, + filter_list=["u", "g", "r", "i", "z", "y"], + **kwargs): + super().__init__(col="metricdata", percentile=5, **kwargs) + # Set mask_val, so that we receive metric_values.filled(mask_val) + self.mask_val = hp.UNSEEN + self.badval = hp.UNSEEN + self.percentileMetric = PercentileMetric(percentile=5, name='percentile') + self.filter_list = filter_list + + def run(self, data_slice, slice_point=None): + # Technically don't need this for now (isn't used in previous one) + # need to define an array of bad values for the masked pixels + badval_arr = np.repeat(self.badval, len(self.filter_list)) + # converts the input recarray to an array + data_slice_list = [ + badval_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + lengths = np.array([len(x) for x in data_slice_list]) + # should be (nbins, npix) + data_slice_arr = np.asarray(data_slice_list, dtype=float).T + data_slice_arr[~np.isfinite(data_slice_arr)] = ( + hp.UNSEEN + ) # need to work with TotalPowerMetric and healpix + + # measure rms in each bin. + # The original metric returns an array at each slice_point (of the + # original slicer) -- so there is a bit of "rearrangement" that + # has to happen to be able to pass a np.array with right dtype + # (i.e. dtype = [("metricdata", float)]) to each call to + # the rmsMetric `run` methods. + percentiles = np.array( + [ + self.percentileMetric.run(np.core.records.fromrecords(x, dtype=[("metricdata", float)])) + for x in data_slice_arr + ] + ) # 6 numbers + + from astropy.modeling.models import Schechter1D + from astropy.cosmology import Planck18 as cosmo + import astropy.units as u + + def schecter_lf( + m_grid: np.ndarray, + log_phi_star: float = -2.84, + M_star: float = -20.91, + alpha: float = -1.68, + redshift: float = 3, + ) -> np.ndarray: + """Schecter Luminosity Function on grid of apparent magnitudes. + + Defaults are for z~3 u-dropout Luminosity Function from Table 6 + of Harikane et al. 2022. + + Parameters + ---------- + m_grid: np.ndarray + Array of apparent AB magnitudes on which to calculate the + luminosity function. + log_phi_star: float, default=-2.84 + Natural log of phi_star, the normalization of the luminosity + function in units of mag^-1 Mpc^-3 + M_star: float, default=-20.91 + The characteristic absolute magnitude where the power-law form + of the luminosity function cuts off. + alpha: float, default=-1.68 + The power law index, also known as the faint-end slope. + redshift: float, default=3 + Redshift used for converting apparent magnitudes into absolute + magnitudes. + + Returns + ------- + np.ndarray + Number density in units mag^-1 Mpc^-3 + """ + # Convert observed magnitudes to absolute + DL = cosmo.luminosity_distance(redshift).to(u.pc).value # Lum. Dist. in pc + M_grid = m_grid - 5 * np.log10(DL / 10) + 2.5 * np.log10(1 + redshift) + + # Calculate luminosity function in absolute magnitudes + schechter = Schechter1D(10 ** log_phi_star, M_star, alpha) + + return schechter(M_grid) + + def number_density( + m_grid: np.ndarray, + LF: np.ndarray, + redshift: float = 3, + dz: float = 1, + ) -> float: + """Calculate number density per deg^2. + + Parameters + ---------- + m_grid: np.ndarray + Array of apparent AB magnitudes corresponding to the + luminosity function values. + LF: np.ndarray + The luminosity function evaluated on m_grid, in units + mag^-1 Mpc^-3. + redshift: float, default=3 + The central redshift used for evaluating comoving quantities. + dz: float, default=1 + The full width of the redshift bin + + Returns + ------- + float + The total number density of galaxies in units deg^-2. + """ + # Calculate comoving depth of redshift bin (Mpc) + chi_far = cosmo.comoving_distance(redshift + dz / 2) + chi_near = cosmo.comoving_distance(redshift - dz / 2) + dchi = chi_far - chi_near + + # Calculate number density (mag^-1 Mpc^-2) + n_dm = (LF / u.Mpc ** 3) * dchi + + # Convert to mag^-1 deg^-2 + deg_per_Mpc = cosmo.arcsec_per_kpc_comoving(redshift).to(u.deg / u.Mpc) + n_dm /= deg_per_Mpc ** 2 + + # Integrate the luminosity function + n = np.trapz(n_dm, m_grid) + + return n.value + + u5 = percentiles[self.filter_list == 'u'] # uband m5 + u3 = u5 + 2.5 * np.log10(5 / 3) # convert to uband 3-sigma + g_cut = u3 - 1.5 # cut on g band + # Calculate LF up to g_cut + m_grid = np.linspace(20, g_cut) + LF = schecter_lf(m_grid) + n_deg2 = number_density(m_grid, LF) # number of objects per sq deg + + return n_deg2 + + +class uDropoutsArea(BaseMetric): + """ + Run as summary metric on MultibandExgalM5. + + Parameters + ---------- + year : `int`, optional + The year of the survey to calculate the bias. + This is used to derive the dm/dz derivative used to translate m5 rms into dz rms. + + Returns + ------- + result : `float` array + The ratio of this bias to the desired DESC y1 upper bound on the bias, and the ratio + between the clbias and the y10 DESC SRD requirement. + Desired values are less than 1 by Y10. + + + Notes + ----- + This is a summary metric to be run on the results + of the MultibandExgalM5. + + MultibandExgalM5 provides the m5 depth in all LSST bands given a specific slice. + + """ + + def __init__(self, + filter_list=["u", "g", "r", "i", "z", "y"], + **kwargs): + super().__init__(col="metricdata", percentile=5, **kwargs) + # Set mask_val, so that we receive metric_values.filled(mask_val) + self.mask_val = hp.UNSEEN + self.badval = hp.UNSEEN + self.percentileMetric = PercentileMetric(percentile=5, name='percentile') + self.filter_list = filter_list + + def run(self, data_slice, slice_point=None): + # Technically don't need this for now (isn't used in previous one) + # need to define an array of bad values for the masked pixels + badval_arr = np.repeat(self.badval, len(self.filter_list)) + # converts the input recarray to an array + data_slice_list = [ + badval_arr if isinstance(x, float) else x for x in data_slice["metricdata"].tolist() + ] + lengths = np.array([len(x) for x in data_slice_list]) + # should be (nbins, npix) + data_slice_arr = np.asarray(data_slice_list, dtype=float).T + data_slice_arr[~np.isfinite(data_slice_arr)] = ( + hp.UNSEEN + ) # need to work with TotalPowerMetric and healpix + + # measure rms in each bin. + # The original metric returns an array at each slice_point (of the + # original slicer) -- so there is a bit of "rearrangement" that + # has to happen to be able to pass a np.array with right dtype + # (i.e. dtype = [("metricdata", float)]) to each call to + # the rmsMetric `run` methods. + percentiles = np.array( + [ + self.percentileMetric.run(np.core.records.fromrecords(x, dtype=[("metricdata", float)])) + for x in data_slice_arr + ] + ) # 6 numbers + + diff --git a/rubin_sim/maf/metrics/summary_metrics.py b/rubin_sim/maf/metrics/summary_metrics.py index eac41369..a8c205ea 100644 --- a/rubin_sim/maf/metrics/summary_metrics.py +++ b/rubin_sim/maf/metrics/summary_metrics.py @@ -5,15 +5,11 @@ "IdentityMetric", "NormalizeMetric", "ZeropointMetric", - "TotalPowerMetric", - "StaticProbesFoMEmulatorMetricSimple", ) -import warnings import healpy as hp import numpy as np -from scipy import interpolate from .base_metric import BaseMetric @@ -203,125 +199,3 @@ def run(self, data_slice, slice_point=None): return result[0] else: return result - - -class TotalPowerMetric(BaseMetric): - """ - Calculate the total power in the angular power spectrum between lmin/lmax. - """ - - def __init__( - self, col="metricdata", lmin=100.0, lmax=300.0, remove_dipole=True, mask_val=np.nan, **kwargs - ): - self.lmin = lmin - self.lmax = lmax - self.remove_dipole = remove_dipole - super(TotalPowerMetric, self).__init__(col=col, mask_val=mask_val, **kwargs) - - def run(self, data_slice, slice_point=None): - # Calculate the power spectrum. - if self.remove_dipole: - cl = hp.anafast(hp.remove_dipole(data_slice[self.colname], verbose=False)) - else: - cl = hp.anafast(data_slice[self.colname]) - ell = np.arange(np.size(cl)) - condition = np.where((ell <= self.lmax) & (ell >= self.lmin))[0] - totalpower = np.sum(cl[condition] * (2 * ell[condition] + 1)) - return totalpower - - -class StaticProbesFoMEmulatorMetricSimple(BaseMetric): - """This calculates the Figure of Merit for the combined - static probes (3x2pt, i.e., Weak Lensing, LSS, Clustering). - This FoM is purely statistical and does not factor in systematics. - - This version of the emulator was used to generate the results in - https://ui.adsabs.harvard.edu/abs/2018arXiv181200515L/abstract - - A newer version is being created. This version has been renamed - Simple in anticipation of the newer, more sophisticated metric - replacing it. - - Note that this is truly a summary metric and should be run on the output of - Exgalm5_with_cuts. - """ - - def __init__(self, nside=128, year=10, col=None, **kwargs): - """ - Args: - nside (int): healpix resolution - year (int): year of the FoM emulated values, - can be one of [1, 3, 6, 10] - col (str): column name of metric data. - """ - self.nside = nside - super().__init__(col=col, **kwargs) - if col is None: - self.col = "metricdata" - self.year = year - - def run(self, data_slice, slice_point=None): - """ - Args: - data_slice (ndarray): Values passed to metric by the slicer, - which the metric will use to calculate metric values - at each slice_point. - slice_point (Dict): Dictionary of slice_point metadata passed - to each metric. - Returns: - float: Interpolated static-probe statistical Figure-of-Merit. - Raises: - ValueError: If year is not one of the 4 for which a FoM is calculated - """ - # Chop off any outliers - good_pix = np.where(data_slice[self.col] > 0)[0] - - # Calculate area and med depth from - area = hp.nside2pixarea(self.nside, degrees=True) * np.size(good_pix) - median_depth = np.median(data_slice[self.col][good_pix]) - - # FoM is calculated at the following values - if self.year == 1: - areas = [7500, 13000, 16000] - depths = [24.9, 25.2, 25.5] - fom_arr = [ - [1.212257e02, 1.462689e02, 1.744913e02], - [1.930906e02, 2.365094e02, 2.849131e02], - [2.316956e02, 2.851547e02, 3.445717e02], - ] - elif self.year == 3: - areas = [10000, 15000, 20000] - depths = [25.5, 25.8, 26.1] - fom_arr = [ - [1.710645e02, 2.246047e02, 2.431472e02], - [2.445209e02, 3.250737e02, 3.516395e02], - [3.173144e02, 4.249317e02, 4.595133e02], - ] - - elif self.year == 6: - areas = [10000, 15000, 2000] - depths = [25.9, 26.1, 26.3] - fom_arr = [ - [2.346060e02, 2.414678e02, 2.852043e02], - [3.402318e02, 3.493120e02, 4.148814e02], - [4.452766e02, 4.565497e02, 5.436992e02], - ] - - elif self.year == 10: - areas = [10000, 15000, 20000] - depths = [26.3, 26.5, 26.7] - fom_arr = [ - [2.887266e02, 2.953230e02, 3.361616e02], - [4.200093e02, 4.292111e02, 4.905306e02], - [5.504419e02, 5.624697e02, 6.441837e02], - ] - else: - warnings.warn("FoMEmulator is not defined for this year") - return self.badval - - # Interpolate FoM to the actual values for this sim - areas = [[i] * 3 for i in areas] - depths = [depths] * 3 - f = interpolate.interp2d(areas, depths, fom_arr, bounds_error=False) - fom = f(area, median_depth)[0] - return fom diff --git a/rubin_sim/maf/metrics/tomography_models.py b/rubin_sim/maf/metrics/tomography_models.py new file mode 100644 index 00000000..0064164a --- /dev/null +++ b/rubin_sim/maf/metrics/tomography_models.py @@ -0,0 +1,1235 @@ +__all__ = ("DENSITY_TOMOGRAPHY_MODEL","MEANZ_TOMOGRAPHY_MODEL") + +import numpy as np + +"""A dictionary of TOMOGRAPHY models for use with the +cosmology summary metrics (TomographicClusteringSigma8bias, MultibandMeanzBiasMetric). + +This dictionary is derived from work shown in +https://github.com/ixkael/ObsStrat/blob/ +meanz_uniformity_maf/code/meanz_uniformity/ +romanrubinmock_for_sigma8tomoghy.ipynb +""" + +# The dictionary DENSITY_TOMOGRAPHY_MODEL contains the current model. +# It is intended for use with various DESC cosmology metrics, in +# particular the TomographicClusteringSigma8bias metric. +# +# The first set of keys are the years (year1, ..., year10), +# since this would change typical depth and galaxy catalog cuts. +# In what follows we have 5 tomographic bins. +# +# The second nested dictionary has the following: +# `sigma8square_model` : +# The fiducial sigma8^2 value used in CCL for the theory predictions. +# `poly1d_coefs_loglog` : +# a polynomial (5th degree) describing the angular power spectra +# (in log log space) in the 5 tomographic bins considered, +# thus has shape (5, 6) +# `lmax` : +# the lmax limits to sum the Cells over when calculating sigma8 +# for each tomographic bin. thus is it of shape (5, ) +# `dlogN_dm5` : +# the derivatives of logN wrt m5 calculated in Qianjun & Jeff's simulations. +# It is an array of 5 dictionaries (5 = the tomographic bins) +# +# Each dictionary must have keys that are the lsst bands. +# If some are missing they are ignored in the linear model. +# These badnpasses are the ones which will be fed to +# NestedLinearMultibandModelMetric. +# Everything else above is going into the modeling. +# +# The notebook I used to make this dictionary is +# https://github.com/ixkael/ObsStrat/blob/meanz_uniformity_maf/code/ +# meanz_uniformity/romanrubinmock_for_sigma8tomography.ipynb + +# The MEANZ_TOMOGRAPHY_MODEL contains a nested dictionary of: + +# 'meanz' whichh is the redshift per band +# 'dz_dm5' the absolute value of the derivative of z in a bin as a function +# of the depth, m5 magnitude +# this is later used to translate fluctuations in depth to fluctuations in tomographic +# redshift + +# The first set of keys are the years (year1, ..., year10), +# since this would change typical depth and galaxy catalog cuts. +# In what follows we have 5 tomographic bins. + +# the mean z values for a given year over 5 z bins were taken from Jeff Newman and Qianjun Hang's +# code which originally produced results for 7 bands (ugrizY and combined) +# +# We are assuming a fixed delta density index (we pulled from index 5 +# in the is for zero m_5 shift) in Qianjun's simulation + +# Each dictionary must have keys that are the lsst bands. +# If some are missing they are ignored in the linear model. +# These badnpasses are the ones which will be fed to +# NestedLinearMultibandModelMetric. + +DENSITY_TOMOGRAPHY_MODEL = { + "year1": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.94546823e-04, + 2.62494452e-02, + -2.27597655e-01, + 4.67266497e-01, + 3.38850627e-01, + -1.14120412e01, + ], + [ + 3.98978657e-04, + 6.69012313e-03, + -1.35005872e-01, + 3.72800783e-01, + 3.97123211e-01, + -1.29908869e01, + ], + [ + 1.03940293e-03, + -5.18759967e-03, + -6.85178825e-02, + 2.69468206e-01, + 4.60031701e-01, + -1.41217016e01, + ], + [ + 1.42986845e-03, + -1.27862974e-02, + -2.26263084e-02, + 1.85536763e-01, + 5.12886349e-01, + -1.48861030e01, + ], + [ + 1.73667255e-03, + -1.89601890e-02, + 1.80214246e-02, + 9.50637447e-02, + 5.58980246e-01, + -1.55053201e01, + ], + ] + ), + "lmax": np.array([64, 94, 124, 146, 164]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.043962974958662145, + "g": 0.02066177878696687, + "r": 0.0018693130274623476, + "i": 0.04276675705247136, + "z": 0.08638732823394031, + "y": 0.0657505127473751, + "ugrizy": 0.2472461161462007, + }, + { + "cst": 0.0, + "u": 0.056364745687503465, + "g": 0.034752650841693544, + "r": -0.0166573314175779, + "i": 0.03094286532175563, + "z": 0.06502345706021054, + "y": 0.05768411667472194, + "ugrizy": 0.22095351609058403, + }, + { + "cst": 0.0, + "u": 0.038530883579453466, + "g": 0.04538858565704199, + "r": 0.033787326976188484, + "i": -0.040733526061764155, + "z": 0.025969868799574303, + "y": 0.03702233496870881, + "ugrizy": 0.1250655075680508, + }, + { + "cst": 0.0, + "u": -0.02161852795176933, + "g": 0.0222336809395544, + "r": 0.05908765697086321, + "i": -0.058749555322422424, + "z": 0.0188797224519434, + "y": 0.035785840616014225, + "ugrizy": 0.046115493661878317, + }, + { + "cst": 0.0, + "u": -0.07624090622672856, + "g": -0.060659802088035716, + "r": -0.01676696950407968, + "i": -0.04141021293176227, + "z": -0.08739685057046752, + "y": -0.033876500857633024, + "ugrizy": -0.2762112449739721, + }, + ], + }, + "year2": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.72116403e-04, + 2.59654127e-02, + -2.26957718e-01, + 4.69432713e-01, + 3.39213385e-01, + -1.13989048e01, + ], + [ + 4.22801028e-04, + 6.36607807e-03, + -1.34108644e-01, + 3.74841312e-01, + 3.96986022e-01, + -1.29600665e01, + ], + [ + 1.07354523e-03, + -5.73388904e-03, + -6.60135335e-02, + 2.67377905e-01, + 4.62220942e-01, + -1.40738990e01, + ], + [ + 1.46948747e-03, + -1.34766061e-02, + -1.89609170e-02, + 1.80340598e-01, + 5.17052663e-01, + -1.48679098e01, + ], + [ + 1.79954923e-03, + -1.99827467e-02, + 2.28427111e-02, + 9.06204330e-02, + 5.69080814e-01, + -1.55353679e01, + ], + ] + ), + "lmax": np.array([64, 95, 123, 147, 166]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.055492947394854004, + "g": 0.012699380270242842, + "r": 0.007530749024132768, + "i": 0.03513585435458864, + "z": 0.07561295665384882, + "y": 0.0736592982969794, + "ugrizy": 0.24374481489942412, + }, + { + "cst": 0.0, + "u": 0.046332184492532485, + "g": 0.029608214159172044, + "r": -0.006075560506575478, + "i": 0.03333244846260029, + "z": 0.06027754177286897, + "y": 0.06506093356954053, + "ugrizy": 0.2084760154143301, + }, + { + "cst": 0.0, + "u": 0.03334437333732201, + "g": 0.032573705405151664, + "r": 0.01704700362928206, + "i": -0.02614338304111813, + "z": 0.03924255373249242, + "y": 0.041759198393945374, + "ugrizy": 0.13123919055353145, + }, + { + "cst": 0.0, + "u": 0.0030051697578186262, + "g": 0.03783532078713637, + "r": 0.06049801557289073, + "i": -0.061409368487154156, + "z": 0.03808865014429153, + "y": 0.02731197586031816, + "ugrizy": 0.10545655695602359, + }, + { + "cst": 0.0, + "u": -0.02882987071992067, + "g": -0.013641440195713459, + "r": 0.03418483342474093, + "i": -0.018501441200469277, + "z": -0.050114709867968066, + "y": 0.021268476475702833, + "ugrizy": -0.049794732397725236, + }, + ], + }, + "year3": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.47335451e-04, + 2.55922442e-02, + -2.25366760e-01, + 4.68243951e-01, + 3.39807736e-01, + -1.14004586e01, + ], + [ + 4.37952430e-04, + 6.13359097e-03, + -1.33135431e-01, + 3.74413364e-01, + 3.97397465e-01, + -1.29291594e01, + ], + [ + 1.08534795e-03, + -5.92068724e-03, + -6.51754697e-02, + 2.66734891e-01, + 4.63243671e-01, + -1.40474513e01, + ], + [ + 1.47812824e-03, + -1.35845536e-02, + -1.88487898e-02, + 1.82317694e-01, + 5.17486035e-01, + -1.48384185e01, + ], + [ + 1.82831353e-03, + -2.03855872e-02, + 2.40741011e-02, + 9.27502779e-02, + 5.73850795e-01, + -1.55657849e01, + ], + ] + ), + "lmax": np.array([64, 94, 123, 147, 168]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.04571336657551873, + "g": 0.01002695291717099, + "r": 0.006381748317232131, + "i": 0.03346386406073478, + "z": 0.07000026154199343, + "y": 0.07122103158922388, + "ugrizy": 0.2276375817331448, + }, + { + "cst": 0.0, + "u": 0.05405844464453291, + "g": 0.03578384485292078, + "r": -0.006971742709660895, + "i": 0.03055814323956491, + "z": 0.0608644383751681, + "y": 0.06285483616682519, + "ugrizy": 0.2191177469764914, + }, + { + "cst": 0.0, + "u": 0.03701180196157078, + "g": 0.02283319745010296, + "r": 0.023135003462698395, + "i": -0.020108849022374618, + "z": 0.03313299331060819, + "y": 0.042722067239431505, + "ugrizy": 0.12769514475136706, + }, + { + "cst": 0.0, + "u": 0.005411661204687025, + "g": 0.03102947466687921, + "r": 0.04758492397282435, + "i": -0.05756589965508176, + "z": 0.038139167579015414, + "y": 0.028170175608809505, + "ugrizy": 0.09819722293690358, + }, + { + "cst": 0.0, + "u": -0.007052479985258033, + "g": 0.0072989635471763055, + "r": 0.04502123045630644, + "i": -0.0073176602744123255, + "z": -0.025917249340435374, + "y": 0.035858309542520166, + "ugrizy": 0.042230607766813, + }, + ], + }, + "year4": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.41568743e-04, + 2.55297105e-02, + -2.25339653e-01, + 4.69342516e-01, + 3.39385437e-01, + -1.13829013e01, + ], + [ + 4.47724999e-04, + 5.99613130e-03, + -1.32666973e-01, + 3.74708889e-01, + 3.97389350e-01, + -1.29064781e01, + ], + [ + 1.09566923e-03, + -6.07450334e-03, + -6.45513085e-02, + 2.66522400e-01, + 4.63653700e-01, + -1.40182124e01, + ], + [ + 1.48898237e-03, + -1.37812174e-02, + -1.77616026e-02, + 1.80633787e-01, + 5.19149830e-01, + -1.48277616e01, + ], + [ + 1.84812024e-03, + -2.07364151e-02, + 2.60418879e-02, + 8.93719930e-02, + 5.77645934e-01, + -1.55728084e01, + ], + ] + ), + "lmax": np.array([63, 94, 123, 148, 169]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.054293436011708406, + "g": 0.009952046153107954, + "r": 0.0054797065106343325, + "i": 0.02576260292778629, + "z": 0.06525553880320857, + "y": 0.06423711415381689, + "ugrizy": 0.20811399377168, + }, + { + "cst": 0.0, + "u": 0.049744205203719596, + "g": 0.034776420860264834, + "r": -0.008954957603949931, + "i": 0.03315665664175325, + "z": 0.06433009373174128, + "y": 0.06440730378199838, + "ugrizy": 0.22417139228019037, + }, + { + "cst": 0.0, + "u": 0.04304128101754591, + "g": 0.017432851239669426, + "r": 0.023084985340137268, + "i": -0.01232102883069664, + "z": 0.03998793202306738, + "y": 0.044236619617342786, + "ugrizy": 0.1392270140374441, + }, + { + "cst": 0.0, + "u": 0.009135822126354363, + "g": 0.036515151515151424, + "r": 0.04077785494592716, + "i": -0.056015513609789944, + "z": 0.04361845936433841, + "y": 0.029657965796579696, + "ugrizy": 0.10045949490109377, + }, + { + "cst": 0.0, + "u": -0.001320660763385228, + "g": 0.009954114661184083, + "r": 0.05948333012553189, + "i": 8.127849262724018e-05, + "z": -0.02151012665531872, + "y": 0.0434040047114252, + "ugrizy": 0.08833836494032939, + }, + ], + }, + "year5": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.35347360e-04, + 2.54490508e-02, + -2.25123057e-01, + 4.69736089e-01, + 3.39586115e-01, + -1.13755955e01, + ], + [ + 4.53471097e-04, + 5.90580865e-03, + -1.32288337e-01, + 3.74594679e-01, + 3.97577407e-01, + -1.29025614e01, + ], + [ + 1.10564516e-03, + -6.23727532e-03, + -6.37987371e-02, + 2.65929839e-01, + 4.64485800e-01, + -1.40110253e01, + ], + [ + 1.49887589e-03, + -1.39372292e-02, + -1.70839226e-02, + 1.80323328e-01, + 5.19727794e-01, + -1.47989678e01, + ], + [ + 1.87614685e-03, + -2.11800406e-02, + 2.81877828e-02, + 8.68286284e-02, + 5.79267039e-01, + -1.55754670e01, + ], + ] + ), + "lmax": np.array([63, 95, 123, 148, 169]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.049239020937134025, + "g": 0.01383328060395141, + "r": -0.0007656653768851196, + "i": 0.028065747821641045, + "z": 0.06050907785720702, + "y": 0.06377977731569703, + "ugrizy": 0.20957610639279134, + }, + { + "cst": 0.0, + "u": 0.055603049706020896, + "g": 0.03635061259366789, + "r": -0.0031966661282920375, + "i": 0.03213295163976452, + "z": 0.06832536553023698, + "y": 0.06457919740205376, + "ugrizy": 0.23517375380435723, + }, + { + "cst": 0.0, + "u": 0.03475293260889745, + "g": 0.017981923705511698, + "r": 0.018654370753291738, + "i": -0.016107465625364203, + "z": 0.03550081300116625, + "y": 0.048124231680847056, + "ugrizy": 0.13850667560597368, + }, + { + "cst": 0.0, + "u": 0.01511037258240387, + "g": 0.033244669574385696, + "r": 0.04459876391783935, + "i": -0.057050518892029986, + "z": 0.04070290636155964, + "y": 0.022294059577929574, + "ugrizy": 0.09371570861821565, + }, + { + "cst": 0.0, + "u": 0.009871915026554134, + "g": 0.0110669656320889, + "r": 0.05829720129030573, + "i": 0.0026049915531623793, + "z": -0.022746170354716745, + "y": 0.04106604585185824, + "ugrizy": 0.09054948216340614, + }, + ], + }, + "year6": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.29402173e-04, + 2.53591187e-02, + -2.24744974e-01, + 4.69511416e-01, + 3.39616928e-01, + -1.13734220e01, + ], + [ + 4.57098537e-04, + 5.86340062e-03, + -1.32257952e-01, + 3.75362916e-01, + 3.97431135e-01, + -1.28861252e01, + ], + [ + 1.10975980e-03, + -6.30581393e-03, + -6.34236372e-02, + 2.65250360e-01, + 4.64594967e-01, + -1.39905973e01, + ], + [ + 1.50470306e-03, + -1.40535883e-02, + -1.63140359e-02, + 1.78531058e-01, + 5.20616617e-01, + -1.47817593e01, + ], + [ + 1.87684022e-03, + -2.11967381e-02, + 2.82475769e-02, + 8.70224981e-02, + 5.80844971e-01, + -1.55660655e01, + ], + ] + ), + "lmax": np.array([63, 94, 123, 147, 169]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.04731277533039653, + "g": 0.007373242386006846, + "r": 0.008974537732206396, + "i": 0.021952456469413933, + "z": 0.0631987046576748, + "y": 0.06698033256378032, + "ugrizy": 0.20345569174304204, + }, + { + "cst": 0.0, + "u": 0.05541230971237573, + "g": 0.03307274442076383, + "r": -0.011785139084349847, + "i": 0.03544678226556244, + "z": 0.06044332758921493, + "y": 0.06316287559760261, + "ugrizy": 0.2248884735953701, + }, + { + "cst": 0.0, + "u": 0.03224724233127598, + "g": 0.018519321739837606, + "r": 0.02629040528053873, + "i": -0.01879121908695766, + "z": 0.043381666839462575, + "y": 0.0465917194286579, + "ugrizy": 0.1360573901416739, + }, + { + "cst": 0.0, + "u": 0.022608143953455426, + "g": 0.029168555210221787, + "r": 0.04195762939125768, + "i": -0.051460395769436555, + "z": 0.03947548537932443, + "y": 0.0191135636135187, + "ugrizy": 0.10394405843289192, + }, + { + "cst": 0.0, + "u": 0.008432569998479788, + "g": 0.015244199381057874, + "r": 0.0564324179119396, + "i": 0.0034085056574326596, + "z": -0.0231010021253421, + "y": 0.049696845834184995, + "ugrizy": 0.10479317991975469, + }, + ], + }, + "year7": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.18947651e-04, + 2.52048402e-02, + -2.24128689e-01, + 4.69247373e-01, + 3.39996105e-01, + -1.13691749e01, + ], + [ + 4.66194938e-04, + 5.71142011e-03, + -1.31505217e-01, + 3.74462215e-01, + 3.97920687e-01, + -1.28754286e01, + ], + [ + 1.11096134e-03, + -6.32112700e-03, + -6.34065660e-02, + 2.65513757e-01, + 4.64794474e-01, + -1.39828608e01, + ], + [ + 1.50696925e-03, + -1.40997435e-02, + -1.60608483e-02, + 1.78259845e-01, + 5.21324331e-01, + -1.47778620e01, + ], + [ + 1.89997066e-03, + -2.16155858e-02, + 3.08453940e-02, + 8.11771421e-02, + 5.82883411e-01, + -1.55574865e01, + ], + ] + ), + "lmax": np.array([63, 94, 123, 148, 170]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.04691474328780031, + "g": 0.010732781728918056, + "r": 0.010261971860717898, + "i": 0.022190117396450985, + "z": 0.05969082742663293, + "y": 0.06903541851329315, + "ugrizy": 0.2089634228500948, + }, + { + "cst": 0.0, + "u": 0.0586609387003185, + "g": 0.02926172329832247, + "r": -0.00867144813333158, + "i": 0.03238973478100006, + "z": 0.06403490756697708, + "y": 0.061258549036290445, + "ugrizy": 0.21531623412108158, + }, + { + "cst": 0.0, + "u": 0.037673080027243955, + "g": 0.024329841045213924, + "r": 0.024028335574766118, + "i": -0.004173864313724429, + "z": 0.04098643494507197, + "y": 0.04706774658791262, + "ugrizy": 0.16192259431907577, + }, + { + "cst": 0.0, + "u": 0.014161985933060993, + "g": 0.02818398259140731, + "r": 0.03724357335900577, + "i": -0.05102903641687596, + "z": 0.04037709323317632, + "y": 0.019008759510214236, + "ugrizy": 0.08457753101369372, + }, + { + "cst": 0.0, + "u": 0.021186578579099866, + "g": 0.01743166214581306, + "r": 0.06093757433043147, + "i": -0.004343243303855121, + "z": -0.026688237135998293, + "y": 0.03848277935815396, + "ugrizy": 0.1090269688631813, + }, + ], + }, + "year8": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.20978717e-04, + 2.52386810e-02, + -2.24297998e-01, + 4.69472361e-01, + 3.39895949e-01, + -1.13578618e01, + ], + [ + 4.67843599e-04, + 5.69127210e-03, + -1.31454331e-01, + 3.74573469e-01, + 3.97784625e-01, + -1.28655285e01, + ], + [ + 1.11711699e-03, + -6.42150038e-03, + -6.29326395e-02, + 2.65065150e-01, + 4.65281866e-01, + -1.39739921e01, + ], + [ + 1.51463261e-03, + -1.42349222e-02, + -1.53268735e-02, + 1.77158327e-01, + 5.22084943e-01, + -1.47629374e01, + ], + [ + 1.89029193e-03, + -2.14244957e-02, + 2.94913874e-02, + 8.48496470e-02, + 5.82823818e-01, + -1.55582137e01, + ], + ] + ), + "lmax": np.array([63, 94, 123, 148, 170]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.04786137493054727, + "g": 0.009655534393294174, + "r": -0.0004343582618342657, + "i": 0.022915274064171095, + "z": 0.059461135879046355, + "y": 0.0620954964843791, + "ugrizy": 0.19278674763853387, + }, + { + "cst": 0.0, + "u": 0.061847229214042274, + "g": 0.03327254021698463, + "r": -0.003037483153393612, + "i": 0.032612371741557, + "z": 0.06330901593731424, + "y": 0.06700098294204117, + "ugrizy": 0.22947625091188906, + }, + { + "cst": 0.0, + "u": 0.033831383811744026, + "g": 0.019055419055419107, + "r": 0.02176497146894746, + "i": -0.005173853421466327, + "z": 0.04318093054632704, + "y": 0.05310544905099789, + "ugrizy": 0.14598595925075922, + }, + { + "cst": 0.0, + "u": 0.014811863932028239, + "g": 0.029658558690905158, + "r": 0.03239827338991352, + "i": -0.0542981307448116, + "z": 0.03624916138976098, + "y": 0.017430252211788642, + "ugrizy": 0.0816796281398052, + }, + { + "cst": 0.0, + "u": 0.008756393562914043, + "g": 0.010614732066291235, + "r": 0.05656537211028906, + "i": 0.001649522457501062, + "z": -0.02285214586340533, + "y": 0.04088282647966966, + "ugrizy": 0.10694124411573667, + }, + ], + }, + "year9": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.23432734e-04, + 2.52943590e-02, + -2.24707794e-01, + 4.70555869e-01, + 3.39390877e-01, + -1.13448965e01, + ], + [ + 4.68298841e-04, + 5.69359050e-03, + -1.31559784e-01, + 3.75124408e-01, + 3.97498952e-01, + -1.28555902e01, + ], + [ + 1.11687532e-03, + -6.40919987e-03, + -6.30785495e-02, + 2.65647380e-01, + 4.64932561e-01, + -1.39661254e01, + ], + [ + 1.51771172e-03, + -1.42863837e-02, + -1.50790664e-02, + 1.76934887e-01, + 5.22416665e-01, + -1.47607945e01, + ], + [ + 1.90353259e-03, + -2.16360592e-02, + 3.05135359e-02, + 8.36649036e-02, + 5.84231771e-01, + -1.55657045e01, + ], + ] + ), + "lmax": np.array([63, 94, 123, 148, 170]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.044390298549263456, + "g": 0.011808056276693719, + "r": 0.006910724770876357, + "i": 0.018412740395619607, + "z": 0.05811594349550721, + "y": 0.06252677636632357, + "ugrizy": 0.18925712972055841, + }, + { + "cst": 0.0, + "u": 0.06119330455885738, + "g": 0.024830399951491312, + "r": -0.006768194405460846, + "i": 0.030065619863371994, + "z": 0.06255158775650581, + "y": 0.07147566736996605, + "ugrizy": 0.22306568518819123, + }, + { + "cst": 0.0, + "u": 0.03789658719578499, + "g": 0.027326499247561357, + "r": 0.022781496880732575, + "i": -0.007626568451506756, + "z": 0.04268771512656232, + "y": 0.04633870731223419, + "ugrizy": 0.1571925190810045, + }, + { + "cst": 0.0, + "u": 0.01457970225016212, + "g": 0.027698293355913416, + "r": 0.03530620058462831, + "i": -0.04816414607472254, + "z": 0.033660844914206074, + "y": 0.022692062998703265, + "ugrizy": 0.08856228929221631, + }, + { + "cst": 0.0, + "u": 0.01976694240068723, + "g": 0.01451018240142195, + "r": 0.04892114102580124, + "i": 0.004035038256796502, + "z": -0.024440454014018773, + "y": 0.03989380894491724, + "ugrizy": 0.10691935463909878, + }, + ], + }, + "year10": { + "sigma8square_model": 0.8**2.0, + "poly1d_coefs_loglog": np.array( + [ + [ + -7.16061620e-04, + 2.51744502e-02, + -2.24107456e-01, + 4.69700181e-01, + 3.39695775e-01, + -1.13492481e01, + ], + [ + 4.69488083e-04, + 5.67705574e-03, + -1.31512276e-01, + 3.75205008e-01, + 3.97691903e-01, + -1.28530855e01, + ], + [ + 1.11838359e-03, + -6.43543666e-03, + -6.29190558e-02, + 2.65259815e-01, + 4.65238973e-01, + -1.39566261e01, + ], + [ + 1.52245292e-03, + -1.43703455e-02, + -1.45911836e-02, + 1.75990818e-01, + 5.22698913e-01, + -1.47457598e01, + ], + [ + 1.92138091e-03, + -2.19184749e-02, + 3.20303750e-02, + 8.07319417e-02, + 5.85683971e-01, + -1.55499003e01, + ], + ] + ), + "lmax": np.array([63, 95, 123, 148, 170]), + "dlogN_dm5": [ + { + "cst": 0.0, + "u": 0.05535070282190964, + "g": 0.011193545843687215, + "r": 0.006737582346111497, + "i": 0.025833960371206527, + "z": 0.06136534095644072, + "y": 0.06324681820708049, + "ugrizy": 0.2074608064455781, + }, + { + "cst": 0.0, + "u": 0.05615851073325299, + "g": 0.029973080589983655, + "r": -0.008700108117848502, + "i": 0.030847889011579494, + "z": 0.06273045722713866, + "y": 0.0650054579514396, + "ugrizy": 0.2228857063055601, + }, + { + "cst": 0.0, + "u": 0.03889599549031938, + "g": 0.02232895038144942, + "r": 0.02262813226603151, + "i": -0.0020788337259347776, + "z": 0.04185161853072303, + "y": 0.046998102923906174, + "ugrizy": 0.1611649091746466, + }, + { + "cst": 0.0, + "u": 0.01329424791323612, + "g": 0.025204737395854025, + "r": 0.031181310573597302, + "i": -0.046678413257774214, + "z": 0.04040377859096919, + "y": 0.020907507991371956, + "ugrizy": 0.08643457382953187, + }, + { + "cst": 0.0, + "u": 0.010769393019973022, + "g": 0.018559094751972344, + "r": 0.04731613285883742, + "i": -0.002550702324093634, + "z": -0.025182273951223154, + "y": 0.04311795497390813, + "ugrizy": 0.08119502007816377, + }, + ], + }, +} + + + +MEANZ_TOMOGRAPHY_MODEL = { + 'year1': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.0061936031125660335, 'g': 0.0017477457637154574, 'r': 0.00857817319840909, 'i': 0.002261268563912217, 'z': 0.003615973972222958, 'y': 0.001392612510731922}, + {'cst': 0.0, 'u': 0.004345685997944008, 'g': 0.0001066449832057932, 'r': 0.02111794222593694, 'i': 0.006836212378911876, 'z': 0.006273169567217962, 'y': 0.005162933422399783}, + {'cst': 0.0, 'u': 0.014700677476920066, 'g': 0.0003565533900656737, 'r': 0.0017381491514960698, 'i': 0.009309643264830638, 'z': 0.016059289172704643, 'y': 0.0047741817275022145}, + {'cst': 0.0, 'u': 0.024205912205433812, 'g': 0.012361098286786311, 'r': 0.01880919231214671, 'i': 0.03358810565134138, 'z': 0.023483929856863206, 'y': 0.014918766325311002}, + {'cst': 0.0, 'u': 0.0321597986461887, 'g': 0.03136635025995545, 'r': 0.021574815075134906, 'i': 0.0014333459879824232, 'z': 0.03291250497198226, 'y': 0.06394244184953382}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year2': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.004338970181070121, 'g': 0.0016768047683343727, 'r': 0.0060592089529286535, 'i': 0.0011217592017519536, 'z': 0.004208995658444005, 'y': 0.0029192155014597277}, + {'cst': 0.0, 'u': 0.006308529577304101, 'g': 0.004489671830412379, 'r': 0.0149282556690614, 'i': 0.005663599202125596, 'z': 0.006618542577755489, 'y': 0.005328390112670172}, + {'cst': 0.0, 'u': 0.014010536894606778, 'g': 0.006531813328583142, 'r': 0.001064844420727941, 'i': 0.007693182966817698, 'z': 0.014105928956172072, 'y': 0.008872500627931887}, + {'cst': 0.0, 'u': 0.01841603610100139, 'g': 0.011810269454770199, 'r': 0.01614766678988043, 'i': 0.03186190002388773, 'z': 0.02311972156473258, 'y': 0.014271489579513739}, + {'cst': 0.0, 'u': 0.021739743015393127, 'g': 0.025149521657059408, 'r': 0.014078456177089594, 'i': 0.0018102347442217066, 'z': 0.028058043868573437, 'y': 0.04699521850017141}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year3': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.004859417694659237, 'g': 0.002180837502931259, 'r': 0.005832952837817742, 'i': 0.0009531182753823333, 'z': 0.0023033760076126217, 'y': 0.0015007675788997344}, + {'cst': 0.0, 'u': 0.004609599247070641, 'g': 0.005031031522752048, 'r': 0.015624451866627847, 'i': 0.005869943031945628, 'z': 0.006884672593261995, 'y': 0.0056058690161543}, + {'cst': 0.0, 'u': 0.0144032500579391, 'g': 0.005722985926304894, 'r': 0.0017250290713755655, 'i': 0.007451163225555915, 'z': 0.013392732698667628, 'y': 0.006364500709516715}, + {'cst': 0.0, 'u': 0.017533605831362774, 'g': 0.008251147477987851, 'r': 0.012096195758780895, 'i': 0.030878807791409498, 'z': 0.022084047620954738, 'y': 0.011990493375150112}, + {'cst': 0.0, 'u': 0.022143260252450846, 'g': 0.022411705429521, 'r': 0.014556293159154483, 'i': 0.0010309571383136837, 'z': 0.020234465337130497, 'y': 0.03985137154196112}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year4': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.003985240196509952, 'g': 0.0027532606644155378, 'r': 0.005483338187710845, 'i': 0.0007528805958253858, 'z': 0.003461810214170973, 'y': 0.0021065880721598085}, + {'cst': 0.0, 'u': 0.006118747084253836, 'g': 0.004545432277725995, 'r': 0.014095107550755952, 'i': 0.005018878296676052, 'z': 0.00706919993983708, 'y': 0.004919183921287182}, + {'cst': 0.0, 'u': 0.013473093290497944, 'g': 0.006327245488813308, 'r': 0.000752143146261697, 'i': 0.0075992674636556215, 'z': 0.015226643705576057, 'y': 0.006045347823248635}, + {'cst': 0.0, 'u': 0.01637953497116186, 'g': 0.009238657231159889, 'r': 0.010311599441057185, 'i': 0.02979323222276616, 'z': 0.02034221470767409, 'y': 0.011741402879354814}, + {'cst': 0.0, 'u': 0.023010046933870196, 'g': 0.019579037549717713, 'r': 0.009522904153614719, 'i': 0.0009216742185484149, 'z': 0.019082329197965237, 'y': 0.03606026413099511}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year5': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.004322174236571098, 'g': 0.000925503896899105, 'r': 0.006001372775289819, 'i': 0.0003851728264808069, 'z': 0.0037290713855516423, 'y': 0.0018691538664168846}, + {'cst': 0.0, 'u': 0.004104191848902968, 'g': 0.0036709692434323334, 'r': 0.015247997305358245, 'i': 0.005108499501041228, 'z': 0.007003013302855577, 'y': 0.0031784313282287586}, + {'cst': 0.0, 'u': 0.014034469176035529, 'g': 0.005501576418784268, 'r': 0.00039697969301331296, 'i': 0.006822303576260055, 'z': 0.01354573898685822, 'y': 0.004348395091403874}, + {'cst': 0.0, 'u': 0.015760374944149894, 'g': 0.009708817471841392, 'r': 0.009439249899976899, 'i': 0.027437747734854297, 'z': 0.01955710723508006, 'y': 0.011708679565838107}, + {'cst': 0.0, 'u': 0.01916921601729655, 'g': 0.017258039968740434, 'r': 0.009435982710271604, 'i': 0.001148713967482839, 'z': 0.018870876477435586, 'y': 0.035129232319458734}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year6': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.00401463803084737, 'g': 0.00227896032952828, 'r': 0.00414517011439339, 'i': 0.00012034334150151153, 'z': 0.0031041475428703466, 'y': 0.0022068346243767634}, + {'cst': 0.0, 'u': 0.004279660869635816, 'g': 0.004974508298049341, 'r': 0.01441969958455425, 'i': 0.005300713532368238, 'z': 0.006705268515480176, 'y': 0.0026288264795089867}, + {'cst': 0.0, 'u': 0.011695908885343312, 'g': 0.006836776262757935, 'r': 0.0008162813403802943, 'i': 0.00941240196835918, 'z': 0.014164973202726755, 'y': 0.005960504500210602}, + {'cst': 0.0, 'u': 0.01287769687115447, 'g': 0.008060701157048416, 'r': 0.00794477580985339, 'i': 0.02660463236164379, 'z': 0.020185492682803233, 'y': 0.011994471806209214}, + {'cst': 0.0, 'u': 0.018022239870737054, 'g': 0.013750767748987586, 'r': 0.0069591222407607446, 'i': 0.0001375124388373366, 'z': 0.01780189567775608, 'y': 0.035301070492626214}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year7': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.0046113242356673875, 'g': 0.001030768662063273, 'r': 0.005149200542821557, 'i': 0.00043617919435715296, 'z': 0.0021465207198366754, 'y': 0.003167482139907165}, + {'cst': 0.0, 'u': 0.003679450398532821, 'g': 0.006178780562933841, 'r': 0.014510390666100199, 'i': 0.005447842562172617, 'z': 0.0066705116666000265, 'y': 0.004619373515126778}, + {'cst': 0.0, 'u': 0.013195868783163119, 'g': 0.006731851559839656, 'r': 0.0013538026851885776, 'i': 0.008233517186360746, 'z': 0.013937828257940694, 'y': 0.007153756573591662}, + {'cst': 0.0, 'u': 0.014532421647278295, 'g': 0.009230187136861906, 'r': 0.006799179893831536, 'i': 0.02661610610698625, 'z': 0.019198295578364243, 'y': 0.012168017632711078}, + {'cst': 0.0, 'u': 0.020163521342320887, 'g': 0.01627643852466436, 'r': 0.0052828216701703975, 'i': 0.002156641147824091, 'z': 0.018451714697112643, 'y': 0.03176559180596637}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year8': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.0034867436464466823, 'g': 0.002391141146475631, 'r': 0.006141460633276389, 'i': 0.0003930247234630615, 'z': 0.001538742083076441, 'y': 0.0020875799389788955}, + {'cst': 0.0, 'u': 0.004651513338892838, 'g': 0.004121379848735558, 'r': 0.01399200393728178, 'i': 0.0047569905088256405, 'z': 0.0058610452249306735, 'y': 0.005329404898274454}, + {'cst': 0.0, 'u': 0.012434514701489941, 'g': 0.006678104495364717, 'r': 0.0005887192070473445, 'i': 0.006891274800122103, 'z': 0.014133567077072695, 'y': 0.005598358228319747}, + {'cst': 0.0, 'u': 0.014334526831459799, 'g': 0.006369695121937281, 'r': 0.0052034988737709505, 'i': 0.026018772671167534, 'z': 0.0200727488322856, 'y': 0.01237458403141792}, + {'cst': 0.0, 'u': 0.015596769593596307, 'g': 0.015692069235335986, 'r': 0.004465109773747347, 'i': 0.0006260066313149872, 'z': 0.017607954960679975, 'y': 0.030327041562193614}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year9': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.004186776052779476, 'g': 0.0016596736105282271, 'r': 0.005646086736528863, 'i': 0.00013094597170195023, 'z': 0.0025711413209097613, 'y': 0.001942007590166414}, + {'cst': 0.0, 'u': 0.004140140794723791, 'g': 0.006036640456360094, 'r': 0.013765288589381483, 'i': 0.00373324246197823, 'z': 0.0063827481880761195, 'y': 0.0040198472711352315}, + {'cst': 0.0, 'u': 0.011217853811063638, 'g': 0.006298778337846079, 'r': 0.00010721557965193624, 'i': 0.009306425563727428, 'z': 0.01384228149777034, 'y': 0.005957723690415474}, + {'cst': 0.0, 'u': 0.013233271264130663, 'g': 0.007701581322318664, 'r': 0.005654670214168566, 'i': 0.025923646930363253, 'z': 0.020148845207008867, 'y': 0.012030159744817224}, + {'cst': 0.0, 'u': 0.01858456186277238, 'g': 0.015107369783887435, 'r': 0.002706504896526665, 'i': 0.0007447286074517631, 'z': 0.015925483270651648, 'y': 0.028646916816506118}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, + 'year10': { + 'dz_dm5': [ + {'cst': 0.0, 'u': 0.003354880770259003, 'g': 0.001458765181158561, 'r': 0.0055856686517475745, 'i': 0.0009030647353411612, 'z': 0.0032068067215789953, 'y': 0.0019167284161829108}, + {'cst': 0.0, 'u': 0.0049884576772368005, 'g': 0.005428041699247426, 'r': 0.013727547832691051, 'i': 0.005432212663811623, 'z': 0.006042559409003354, 'y': 0.00434503390924071}, + {'cst': 0.0, 'u': 0.01032419109971061, 'g': 0.00525952527032317, 'r': 0.0005514979865284969, 'i': 0.008857060937095687, 'z': 0.013099492539662139, 'y': 0.005312102829780564}, + {'cst': 0.0, 'u': 0.013229857206590936, 'g': 0.007145690671123007, 'r': 0.004670048309991343, 'i': 0.024555656028367286, 'z': 0.019766085030518556, 'y': 0.010660490885596226}, + {'cst': 0.0, 'u': 0.014555683932784726, 'g': 0.01291152464999209, 'r': 0.0009053946982488664, 'i': 0.002285602665316371, 'z': 0.016921396981975284, 'y': 0.028678557472965396}, + ], + 'meanz': [ + {'cst': 0.0, 'u': 0.1412097514129552, 'g': 0.1406497428388376, 'r': 0.1408687434840175, 'i': 0.14179709129841037, 'z': 0.1425811131555499, 'y': 0.14130581564879285}, + {'cst': 0.0, 'u': 0.31887802365309226, 'g': 0.3190376986628495, 'r': 0.3209289780716591, 'i': 0.32072880008997834, 'z': 0.32187019107155196, 'y': 0.3201021125646}, + {'cst': 0.0, 'u': 0.4978288897927948, 'g': 0.4979842029016838, 'r': 0.4978477649299082, 'i': 0.4976698193750615, 'z': 0.4959446673587056, 'y': 0.4966033026751485}, + {'cst': 0.0, 'u': 0.6950850709405102, 'g': 0.6928492900731955, 'r': 0.6927398597939921, 'i': 0.6930562085764114, 'z': 0.6944189402802977, 'y': 0.6923577239768269}, + {'cst': 0.0, 'u': 0.8682973716127108, 'g': 0.8704231456789853, 'r': 0.8727016782075796, 'i': 0.870275847683205, 'z': 0.8705515709345578, 'y': 0.8734045000404805}, + ], + }, +} \ No newline at end of file diff --git a/rubin_sim/maf/metrics/uniformity_metrics.py b/rubin_sim/maf/metrics/uniformity_metrics.py new file mode 100644 index 00000000..92667993 --- /dev/null +++ b/rubin_sim/maf/metrics/uniformity_metrics.py @@ -0,0 +1,446 @@ +__all__ = ("SingleLinearMultibandModelMetric", "NestedLinearMultibandModelMetric") + +import numpy as np +import healpy as hp + +from .base_metric import BaseMetric +from .exgal_m5 import ExgalM5 +from .weak_lensing_systematics_metric import RIZDetectionCoaddExposureTime, ExgalM5WithCuts + + +class SingleLinearMultibandModelMetric(BaseMetric): + """Calculate a single linear combination of depths. + + This is useful for calculating density or redshift fluctuations + from depth fluctuations on the sky, + for a single redshift bins (thus it is a single metric). + For multiple bins, see NestedLinearMultibandModelMetric. + + Parameters + ---------- + arr_of_model_dicts : `dict` + Linear coefficients to be applied to the M5 depth values. + Keys should be filter names + 'cst' if there is a constant offset. + post_processing_fn : lambda + Post processing function to apply to the linear combination of inputs. + Default `lambda x: x` simply returns the output. + extinction_cut : `float`, optional + E(B-V) cut on extinction (0.2 by default). + n_filters : `int`, optional + Cut on the number of filters required (6 by default). + min_depth_cut : `dict`, optional + Cut the areas of low depths, based on cuts in this dict. + Keys should be filter names. Default is no cuts. + max_depth_cut : `dict`, optional + Cut the areas of high depths, based on cuts in this dict. + Keys should be filter names. Default is no cuts. + mean_depth : `dict`, optional + Mean depths, which will be subtracted from the inputs + before applying model. + Keys should be filter names. Default is zero in each band. + In a lot of cases, instead of feeding mean_depths, + one can remove the monopole in the constructed healpix maps. + m5_col : `str`, optional + Column name for the m5 depth. + filter_col : `str`, optional + Column name for the filter. + units : `str`, optional + Label for "units" in the output, for use in plots. + badval : `float`, optional + Value to return for metric failure. + Specified here so that exgalm5 uses the same value. + + Returns + ------- + result : `float` + Linear combination of the M5 depths (minus mean depth, if provided). + result = np.sum([ + model_dict[band] * M5_depth[band] + for band in ["u", "g", "r", "i", "z", "y"] + ]) + if post_processing_fn is provided: + result => post_processing_fn(result) + """ + + def __init__( + self, + model_dict, + post_processing_fn=lambda x: x, + extinction_cut=0.2, + n_filters=6, + min_depth_cut=None, + max_depth_cut=None, + mean_depth=None, + m5_col="fiveSigmaDepth", + filter_col="filter", + units="mag", + badval=np.nan, + **kwargs, + ): + if min_depth_cut is None: + min_depth_cut = {} + self.min_depth_cut = min_depth_cut + if max_depth_cut is None: + max_depth_cut = {} + self.max_depth_cut = max_depth_cut + if mean_depth is None: + # Set mean_depth to 0 if not specified + mean_depth = dict([(f, 0) for f in "ugrizy"]) + self.mean_depth = mean_depth + + self.n_filters = n_filters + self.extinction_cut = extinction_cut + if "cst" in model_dict: + self.cst = model_dict["cst"] + else: + self.cst = 0.0 + lsst_filters = ["u", "g", "r", "i", "z", "y"] + self.bands = [x for x in model_dict.keys() if x in lsst_filters] + self.model_arr = np.array([model_dict[x] for x in self.bands])[:, None] + + self.m5_col = m5_col + self.filter_col = filter_col + self.post_processing_fn = post_processing_fn + + self.exgal_m5 = ExgalM5(m5_col=m5_col, filter_col=filter_col, badval=badval) + super().__init__( + col=[m5_col, filter_col], units=units, maps=self.exgal_m5.maps, badval=badval, **kwargs + ) + + def run(self, data_slice, slice_point=None): + """ """ + + if slice_point["ebv"] > self.extinction_cut: + return self.badval + + n_filters = len(set(data_slice[self.filter_col])) + if n_filters < self.n_filters: + return self.badval + + # add extinction_cut and depth cuts? and n_filters cuts + depths = np.vstack( + [ + self.exgal_m5.run(data_slice[data_slice[self.filter_col] == lsst_filter], slice_point) + for lsst_filter in self.bands + ] + ).T + for i, lsst_filter in enumerate(self.bands): + if lsst_filter in self.mean_depth: + depths[i] -= self.mean_depth[lsst_filter] + # Don't raise Exceptions in the middle of metrics - + # assert depths.shape[0] >= 1 + if depths.shape[0] < 1: + return self.badval + + if np.any(depths == self.badval): + return self.badval + + for i, lsst_filter in enumerate(self.bands): + if ( + lsst_filter in self.min_depth_cut + and depths[0, i] < self.min_depth_cut[lsst_filter] - self.mean_depth[lsst_filter] + ): + return self.badval + if ( + lsst_filter in self.max_depth_cut + and depths[0, i] > self.max_depth_cut[lsst_filter] - self.mean_depth[lsst_filter] + ): + return self.badval + + val = self.post_processing_fn(self.cst + np.dot(depths, self.model_arr)) + return val + + +class NestedLinearMultibandModelMetric(BaseMetric): + """Calculate multiple linear combinations of depths. + + This is useful for calculating density or redshift fluctuations + from depth fluctuations on the sky, + for multiple redshift bins (thus it is a nested metric). + For a single bin, see LinearMultibandModelMetric. + + Points of contact / contributors: Boris Leistedt + + Parameters + ---------- + arr_of_model_dicts : `list` [ `dicts` ] + Array of linear coefficients to be applied to the M5 depth values. + Keys should be filter names + 'cst' if there is a constant offset. + post_processing_fn : lambda + Post processing function to apply to the linear combination of inputs. + Default `lambda x: x` simply returns the output. + extinction_cut : `float`, optional + E(B-V) cut on extinction (0.2 by default). + n_filters : `int`, optional + Cut on the number of filters required (6 by default). + min_depth_cut : `dict`, optional + Cut the areas of low depths, based on cuts in this dict. + Keys should be filter names. Default is no cuts. + max_depth_cut : `dict`, optional + Cut the areas of high depths, based on cuts in this dict. + Keys should be filter names. Default is no cuts. + mean_depth : `dict`, optional + Mean depths, which will be subtracted from the inputs + before applying model. + Keys should be filter names. Default is zero in each band. + In a lot of cases, instead of feeding mean_depths, + one can remove the monopole in the constructed healpix maps. + m5_col : `str`, optional + Column name for the m5 depth. + filter_col : `str`, optional + Column name for the filter. + units : `str`, optional + Label for "units" in the output, for use in plots. + badval : `float`, optional + Value to return for metric failure. + Specified here so that exgalm5 uses the same value. + + Returns + ------- + result : `list` [ `float` ] + List is the same size as arr_of_model_dicts. + For each element, return a linear combination of the M5 depths + (minus mean depth, if provided). + for i, model_dict in enumerate(arr_of_model_dicts): + result[i] = np.sum([ + model_dict[band] * M5_depth[band] + for band in ["u", "g", "r", "i", "z", "y"] + ]) + if post_processing_fn is provided: + result[i] => post_processing_fn(result[i]) + """ + + def __init__( + self, + arr_of_model_dicts, + post_processing_fn=lambda x: x, + extinction_cut=0.2, + n_filters=6, + min_depth_cut=None, + max_depth_cut=None, + mean_depth=None, + m5_col="fiveSigmaDepth", + filter_col="filter", + units="mag", + badval=np.nan, + **kwargs, + ): + self.arr_of_model_dicts = arr_of_model_dicts + self.n_filters = n_filters + self.extinction_cut = extinction_cut + if min_depth_cut is None: + min_depth_cut = {} + self.min_depth_cut = min_depth_cut + if max_depth_cut is None: + max_depth_cut = {} + self.max_depth_cut = max_depth_cut + if mean_depth is None: + mean_depth = {} + self.mean_depth = mean_depth + + self.m5_col = m5_col + self.badval = badval + self.mask_val = hp.UNSEEN + self.filter_col = filter_col + self.post_processing_fn = post_processing_fn + + self.n_bins = len(self.arr_of_model_dicts) + self.bad_val_arr = np.repeat(self.badval, self.n_bins) + + self.exgal_m5 = ExgalM5( + m5_col=m5_col, + filter_col=filter_col, + badval=self.badval, + ) + super().__init__( + col=[m5_col, filter_col], + badval=self.badval, + units=units, + maps=self.exgal_m5.maps, + metric_dtype="object", + **kwargs, + ) + + def run(self, data_slice, slice_point=None): + # apply extinction and n_filters cut + if slice_point["ebv"] > self.extinction_cut: + return self.bad_val_arr + + n_filters = len(set(data_slice[self.filter_col])) + if n_filters < self.n_filters: + return self.bad_val_arr + + lsst_filters = ["u", "g", "r", "i", "z", "y"] + + # initialize dictionary of outputs + # types = [float]*n_bins + result_arr = np.zeros((self.n_bins,), dtype=float) + + for binnumber, model_dict in enumerate(self.arr_of_model_dicts): + + # put together a vector of depths + depths = np.vstack( + [ + self.exgal_m5.run(data_slice[data_slice[self.filter_col] == lsst_filter], slice_point) + for lsst_filter in lsst_filters + ] + ).T + + # if there are depth cuts, apply them + for i, lsst_filter in enumerate(lsst_filters): + if lsst_filter in self.min_depth_cut and depths[0, i] < self.min_depth_cut[lsst_filter]: + depths[0, i] = self.badval + if lsst_filter in self.max_depth_cut and depths[0, i] > self.max_depth_cut[lsst_filter]: + depths[0, i] = self.badval + + # if provided, subtract the mean + for i, lsst_filter in enumerate(lsst_filters): + if lsst_filter in self.mean_depth: + depths[i] -= self.mean_depth[lsst_filter] + + # now assemble the results + model_arr = np.array([model_dict[x] if x in model_dict else 0.0 for x in lsst_filters]) + cst = model_dict["cst"] if "cst" in model_dict else 0.0 + + result_arr[binnumber] = self.post_processing_fn(cst + np.dot(depths, model_arr)) + + # returns array where each element of the original + # input arr_of_model_dicts contains a scalar value resulting + # from a linear combination of the six depths + return result_arr + + +class MultibandExgalM5(BaseMetric): + """Calculate multiple linear combinations of depths.""" + + def __init__( + self, + extinction_cut=0.2, + n_filters=6, + m5_col="fiveSigmaDepth", + filter_col="filter", + units="mag", + badval=np.nan, + **kwargs, + ): + self.n_filters = n_filters + self.extinction_cut = extinction_cut + + self.m5_col = m5_col + self.badval = badval + self.filter_col = filter_col + + self.lsst_filters = ["u", "g", "r", "i", "z", "y"] + self.bad_val_arr = np.repeat(self.badval, len(self.lsst_filters)) + + self.exgal_m5 = ExgalM5( + m5_col=m5_col, + filter_col=filter_col, + badval=self.badval, + ) + super().__init__( + col=[m5_col, filter_col], + badval=self.badval, + units=units, + maps=self.exgal_m5.maps, + metric_dtype="object", + **kwargs, + ) + + def run(self, data_slice, slice_point=None): + # apply extinction and n_filters cut + if slice_point["ebv"] > self.extinction_cut: + return self.bad_val_arr + + n_filters = len(set(data_slice[self.filter_col])) + if n_filters < self.n_filters: + return self.bad_val_arr + + depths = np.vstack( + [ + self.exgal_m5.run(data_slice[data_slice[self.filter_col] == lsst_filter], slice_point) + for lsst_filter in self.lsst_filters + ] + ).T + + return depths.ravel() + + +class NestedRIZExptimeExgalM5Metric(BaseMetric): + """ + This is a simple wrapper metric which returns values of both RIZDetectionCoaddExposureTime and ExgalM5WithCuts + in a recarray with cvolumns ["exgal_m5", "riz_exptime"] + + Points of contact / contributors: Rachel Mandelbaum, Boris Leistedt + + Parameters + ---------- + extinction_cut : `float`, optional + E(B-V) cut on extinction (0.2 by default) (for ExgalM5WithCuts). + n_filters : `int`, optional + Cut on the number of filters required (6 by default) (for ExgalM5WithCuts). + lsst_filter : `str`, optional + The filter choice for ExgalM5WithCuts + exptime_col: `str`, optional + Column name for the exposure time (for RIZDetectionCoaddExposureTime). + m5_col : `str`, optional + Column name for the m5 depth (for ExgalM5WithCuts). + filter_col : `str`, optional + Column name for the filter. + units : `str`, optional + Label for "units" in the output, for use in plots. + badval : `float`, optional + Value to return for metric failure. + + Returns + ------- + result : `recarray`, np.array(dtype=[("exgal_m5", float), ("riz_exptime", float)]) + The values of RIZDetectionCoaddExposureTime and ExgalM5WithCuts + + """ + + def __init__( + self, + m5_col="fiveSigmaDepth", + filter_col="filter", + exptime_col="visitExposureTime", + extinction_cut=0.2, + n_filters=6, + depth_cut=24, + metric_name="new_nested_FoM", + lsst_filter="i", + badval=np.nan, + **kwargs, + ): + maps = ["DustMap"] + self.m5_col = m5_col + self.filter_col = filter_col + self.exptime_col = exptime_col + self.lsst_filter = lsst_filter + self.n_filters = n_filters + + cols = [self.m5_col, self.filter_col, self.exptime_col] + super().__init__(cols, metric_name=metric_name, maps=maps, badval=badval, **kwargs) + self.riz_exptime_metric = RIZDetectionCoaddExposureTime(ebvlim=extinction_cut) + self.exgalm5_metric = ExgalM5WithCuts( + m5_col=m5_col, + filter_col=filter_col, + extinction_cut=extinction_cut, + depth_cut=depth_cut, + lsst_filter=lsst_filter, + badval=badval, + n_filters=n_filters, + ) + + self.metric_dtype = "object" + + def run(self, data_slice, slice_point): + + # set up array to hold the results + names = ["exgal_m5", "riz_exptime"] + types = [float] * 2 + result = np.zeros(1, dtype=list(zip(names, types))) + result["exgal_m5"] = self.exgalm5_metric.run(data_slice, slice_point) # if using ExgalM5WithCuts + result["riz_exptime"] = self.riz_exptime_metric.run(data_slice, slice_point) + + return result diff --git a/tests/maf/test_cosmology_summaries.py b/tests/maf/test_cosmology_summaries.py new file mode 100644 index 00000000..3d53629f --- /dev/null +++ b/tests/maf/test_cosmology_summaries.py @@ -0,0 +1,20 @@ +import unittest + +import healpy as hp +import numpy as np + +import rubin_sim.maf.metrics as metrics + + +class TestCosmologySummaryMetrics(unittest.TestCase): + + def test_total_power_metric(self): + nside = 64 + data = np.ones(hp.nside2npix(nside), dtype=list(zip(["testcol"], ["float"]))) + metric = metrics.TotalPowerMetric(col="testcol") + result = metric.run(data) + np.testing.assert_equal(result, 0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/maf/test_summarymetrics.py b/tests/maf/test_summarymetrics.py index 3ae2750c..8ed84be9 100644 --- a/tests/maf/test_summarymetrics.py +++ b/tests/maf/test_summarymetrics.py @@ -77,13 +77,6 @@ def test_zeropoint_metric(self): result = metric.run(data) np.testing.assert_equal(result, np.ones(10, float) + 5.5) - def test_total_power_metric(self): - nside = 128 - data = np.ones(12 * nside**2, dtype=list(zip(["testcol"], ["float"]))) - metric = metrics.TotalPowerMetric(col="testcol") - result = metric.run(data) - np.testing.assert_equal(result, 0.0) - if __name__ == "__main__": unittest.main()