diff --git a/examples/anacal/config.yml b/examples/anacal/config.yml new file mode 100644 index 000000000..0212dc9dc --- /dev/null +++ b/examples/anacal/config.yml @@ -0,0 +1,79 @@ +global: + chunk_rows: 100000 + sparse: true # Generate sparse maps - faster if using small areas + pixelization: healpix + nside: 512 + +TXIngestAnacal: + # DP1-v3 AnaCal outputs for the a360 cluster field (per-tract). + butler_config_file: /global/cfs/cdirs/lsst/production/gen3/rubin/DP1/repo/butler.yaml + collections: "u/xiangchl/dp1-v3/a360_anacal2" + + # a360 sits in the "Low Ecliptic Latitude Field" (LELF). DP1_TRACTS["LELF"] + # contains the two live tracts 10463, 10464 (plus four never-processed). + select_field: "LELF" + cosmology_tracts_only: False + + # Merged catalog carries band-combined shape moments under the "fpfs1_*" + # prefix and Gaussian-aperture fluxes at scale "gauss2". Per-band + # columns are survey-prefixed as ``lsst_{band}_...`` (v3 schema) and + # carry pre-computed magnitudes + dmag/dg responses, so no + # nanojansky→mag conversion happens inside the ingest. a360 was + # measured on griz only (no u/y coadds). + prefix: "fpfs1" + scale: "gauss2" + bands: "griz" + +TXSourceSelectorAnaCal: + # AnaCal shear-response step size for selection bias correction — must + # match the δγ used to build the on-disk photo-z distortions + # (zmode_{1p,1m,2p,2m}) in xlens' photoZPipe (currently ±0.01 in + # DISTORTIONS). + delta_gamma: 0.01 + + # mask_value threshold — sources with mask_value >= this get dropped. + # AnaCal reports a bit-packed flag; 1 means "any bit set" (strictest). + mask_threshold: 40 + + # AnaCal S/N cut (flux / flux_err). Default in code is 5.0. + s2n_cut: 5.0 + + # Band-combined size cut (m00 + m20) / m00 > T_cut on fpfs1 moments. + # Default 0.1. + T_cut: 0.1 + + # Band-combined shape magnitude cut |e| < emax (selector applies + # esq < emax**2 on the ``esq`` column). + emax: 0.5 + + # Per-band AB mag upper bounds (uses pre-computed mag_{b} from v3 + # ingest). Default 50 is above the xlens smooth-truncation cap + # (mag=40) so no cut is applied; set per band to gate on brightness. + g_hi_cut: 50.0 + r_hi_cut: 50.0 + i_hi_cut: 50.0 + z_hi_cut: 50.0 + + # Bands available in the a360 ingest. + bands: [g, r, i, z] + + # 5 tomographic bins — will need edited to fit your science n(z). + source_zbin_edges: [0.3, 0.6, 0.9, 1.2, 1.5, 2.0] + + # ``input_pz: true`` reads ``mean_z`` from the shear catalog and skips the + # random-forest classifier. TXIngestAnacal exposes ``mean_z`` (from + # ``zmode_0`` in the merged catalog) plus the four shifted variants + # ``mean_z_{1p,1m,2p,2m}`` (from ``zmode_{1p,1m,2p,2m}``) — the selector + # consumes the baseline column for tomographic binning and the four + # variants for the selection response computation. + # If ``input_pz: false``, a shear_tomography_classifier pickle + + # spectroscopic_catalog are required (see pipeline.yml inputs). + input_pz: true + true_z: false + verbose: True + chunk_rows: 100000 + +TXShearCalibration: + shear_catalog_type: anacal + subtract_mean_shear: true + chunk_rows: 100000 diff --git a/examples/anacal/ingest.yml b/examples/anacal/ingest.yml new file mode 100644 index 000000000..87701e68e --- /dev/null +++ b/examples/anacal/ingest.yml @@ -0,0 +1,28 @@ +# AnaCal ingest — runs TXIngestAnacal only. Writes the shear catalog +# into ``data/example/anacal_inputs`` which is then consumed by +# ``examples/anacal/pipeline.yml``. + +stages: + - name: TXIngestAnacal + +output_dir: data/example/anacal_inputs + +launcher: + name: mini + interval: 1.0 + +site: + name: local + +modules: > + txpipe + +config: examples/anacal/config.yml + +inputs: + # Butler-mode ingest doesn't read this FITS, but ceci requires every + # declared input to have a path — pointed at the DP1-v3 published + # per-field catalog to satisfy the graph builder. + anacal_catalog: /global/cfs/cdirs/desc-wl/projects/anacal/DP1-v3/catalogs/anacal_catalog_a360.fits + +resume: False diff --git a/examples/anacal/pipeline.yml b/examples/anacal/pipeline.yml new file mode 100644 index 000000000..7b0bd341c --- /dev/null +++ b/examples/anacal/pipeline.yml @@ -0,0 +1,35 @@ +# AnaCal selector + calibrator pipeline — reads the shear catalog produced +# by ``examples/anacal/ingest.yml``. + +stages: + - name: TXSourceSelectorAnaCal + - name: TXShearCalibration + +output_dir: data/example/output_anacal + +launcher: + name: mini + interval: 1.0 + +site: + name: local + max_threads: 2 + +modules: > + txpipe + +config: examples/anacal/config.yml + +inputs: + shear_catalog: data/example/anacal_inputs/shear_catalog.hdf5 + + # Selector-side inputs. With ``input_pz: true`` in config.yml the + # selector reads ``mean_z`` from the shear catalog (populated by + # TXIngestAnacal from ``zmode_0``), so these two are only exercised + # under ``input_pz: false`` (random-forest tomography). Kept here to + # satisfy the ceci graph builder — ceci checks every declared input + # has a path, even when the stage doesn't open it. + spectroscopic_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/spec_catalog_a360.hdf5 + shear_tomography_classifier: /pscratch/sd/x/xiangchl/data/DP1/catalogs/shear_tomography_classifier_a360.pkl + +resume: False diff --git a/txpipe/__init__.py b/txpipe/__init__.py index 8164ef4d5..5f174171e 100644 --- a/txpipe/__init__.py +++ b/txpipe/__init__.py @@ -10,6 +10,7 @@ TXSourceSelectorMetacal, TXSourceSelectorLensfit, TXSourceSelectorMetadetect, + TXSourceSelectorAnacal, ) from .lens_selector import TXMeanLensSelector from .photoz_stack import TXPhotozStack, TXPhotozPlot, TXTruePhotozStack diff --git a/txpipe/calibrate.py b/txpipe/calibrate.py index 37342ce09..a5227826c 100644 --- a/txpipe/calibrate.py +++ b/txpipe/calibrate.py @@ -181,6 +181,8 @@ def run(self): # therefore, we add dec to split data into these fields. # You can choose not to by setting dec_cut = 90 in the config, for example. d["g1"], d["g2"] = cal.apply(d["g1"], d["g2"], d["dec"], subtract_mean=subtract_mean_shear) + elif cat_type == "anacal": + d["g1"], d["g2"] = cal.apply(d["g1"], d["g2"], d["weight"], subtract_mean=subtract_mean_shear) else: d["g1"], d["g2"] = cal.apply(d["g1"], d["g2"], subtract_mean=subtract_mean_shear) diff --git a/txpipe/data_types.py b/txpipe/data_types.py index 3326e8d2d..b5a65e1b0 100755 --- a/txpipe/data_types.py +++ b/txpipe/data_types.py @@ -108,6 +108,9 @@ def get_primary_catalog_names(self, true_shear=False): elif self.catalog_type == "metadetect": shear_cols = ["00/g1", "00/g2", "00/ra", "00/dec", "00/weight"] rename = {c: c[3:] for c in shear_cols} + elif self.catalog_type == "anacal": + shear_cols = ["e1", "e2", "ra", "dec", "weight"] + rename = {"e1":"g1", "e2":"g2"} else: shear_cols = ["g1", "g2", "ra", "dec", "weight"] rename = {} diff --git a/txpipe/diagnostics.py b/txpipe/diagnostics.py index adbb5b1fa..56fab8932 100644 --- a/txpipe/diagnostics.py +++ b/txpipe/diagnostics.py @@ -223,6 +223,8 @@ def run(self): "weight", ) shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="metadetect") + elif cat_type == "anacal": + raise NameError("anacal does not work with this atm.") else: shear_cols = [ "dec", diff --git a/txpipe/ingest/__init__.py b/txpipe/ingest/__init__.py index 38573dda4..583476bed 100644 --- a/txpipe/ingest/__init__.py +++ b/txpipe/ingest/__init__.py @@ -10,4 +10,5 @@ TXIngestSSIDetectionDESBalrog, ) from .dp1 import TXIngestDataPreview1 +from .anacal import TXIngestAnacal from .stage3 import TXIngestDESY3Gold diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py new file mode 100644 index 000000000..e94657f9b --- /dev/null +++ b/txpipe/ingest/anacal.py @@ -0,0 +1,484 @@ +from .base import TXIngestCatalogFits +from ..data_types import ShearCatalog, FitsFile +from .dp1_details import ( + DP1_TRACTS, + DP1_COSMOLOGY_TRACTS, +) +from ceci.config import StageParameter +import numpy as np +from ..utils.hdf_tools import h5py_shorten, repack + + +# Suffixes on the merged catalog's photo-z point-estimate columns +# (zmode_0, zmode_1p, zmode_1m, zmode_2p, zmode_2m). These become +# ``mean_z`` / ``mean_z_{1p,1m,2p,2m}`` in the ingested shear catalog +# and are consumed by TXSourceSelectorAnaCal for the tomographic +# bin-migration term of R_sel via _DataWrapper suffix lookup. +PZ_SUFFIXES = ("0", "1p", "1m", "2p", "2m") + + +class TXIngestAnacal(TXIngestCatalogFits): + """ + Ingestion of an anacal catalog, generated from actual Rubin data. This + stage, will take an anacal catalog, from either the butler, or a file + (parquet), and ingest it into TXPipe format (HDF5). + """ + + name = "TXIngestAnacal" + inputs = [ + ("anacal_catalog", FitsFile) + ] + outputs = [ + ("shear_catalog", ShearCatalog), + ] + config_options = { + "use_butler": StageParameter( + bool, True, + msg="Should be left on, unless you got an external file, " + "in that case knock yourself out!", + ), + "butler_config_file": StageParameter( + str, + "/global/cfs/cdirs/lsst/production/gen3/rubin/DP1/repo/butler.yaml", + msg="Path to the LSST butler config file.", + ), + "butler_object_name": StageParameter( + str, "deep_coadd_cell_anacal_merged", + ), + "cosmology_tracts_only": StageParameter( + bool, True, msg="Use only cosmology tracts.", + ), + "select_field": StageParameter( + str, "", + msg="Field to select (overrides cosmology_tracts_only).", + ), + "select_tracts": StageParameter( + list, [], + msg="list of tracts (overrides cosmology_tracts_only, but " + "not select_field).", + ), + "collections": StageParameter( + str, "LSSTComCam/DP1", msg="Butler collections to use.", + ), + "tracts": StageParameter( + str, "", + msg="Comma-separated list of tracts to use (empty for all).", + ), + "prefix": StageParameter( + str, "fpfs", + msg="prefix indicating the method used to calculate the ", + ), + "bands": StageParameter( + list, ["g", "r", "i", "z", "y"], msg="string of flux bands", + ), + "scale": StageParameter( + str, "gauss2", + msg="scale radius for the convolution with Gaussian PSF", + ), + "delta_gamma": StageParameter( + float, 0.2, + msg="delta gamma value used for the analytical shearing." + ), + "add_psf_properties": StageParameter( + bool, True, + msg="Emit per-band psf_g1/psf_g2/psf_T_mean (one set per band in " + "'bands') from the PSF HSM second moments " + "(lsst__ext_shapeHSM_HsmPsfMoments_*). Set False for " + "catalogs that predate doPsfHsmMoments.", + ), + } + + def run(self): + if self.config["use_butler"]: + self.butler_run() + else: + self.file_run() + + print("repacking files") + repack(self.get_output("shear_catalog")) + + def butler_run(self): + error_msg = ( + "The LSST Science Pipelines are not installed in this environment, " + "or are not configured correctly to access the data. " + "See the note in the file example/dp1/ingest.yml for how to set " + "this up on NERSC." + ) + try: + from lsst.daf.butler import Butler + except Exception as e: + raise ImportError(error_msg) from e + + # Configure and create the butler. There are several ways to do this, + # Here we use a central collective butler yaml file from NERSC. + + butler_config_file = self.config["butler_config_file"] + collections = self.config["collections"] + error_msg2 = error_msg + ( + ' Or there is a typo in the collection you have looked for.' + ) + try: + butler = Butler(butler_config_file, collections=collections) + except Exception as e: + raise RuntimeError(error_msg2) from e + + if self.config["select_field"]: + tracts = DP1_TRACTS[self.config["select_field"]] + elif self.config["select_tracts"]: + tracts = self.config["select_tracts"] + elif self.config["cosmology_tracts_only"]: + tracts = DP1_COSMOLOGY_TRACTS + else: + # No explicit selection: ingest every tract present in the + # collection (works for any survey, incl. DP2). Previously this + # fell back to the DP1 ALL_TRACTS list, which selects nothing on + # a DP2 collection. + tracts = None + + object_name = self.config["butler_object_name"] + n = self.get_catalog_size(butler, object_name) + + created_files = False + data_set_refs = butler.query_datasets(object_name) + n_chunks = len(data_set_refs) + input_columns = self.setup_input() + + shear_start = 0 + for i, ref in enumerate(data_set_refs): + tract = ref.dataId["tract"] + if tracts is not None and tract not in tracts: + print( + f"Skipping chunk {i + 1} / {n_chunks} since tract " + f"{tract} is not selected" + ) + continue + + d = butler.get( + object_name, + dataId=ref.dataId, + parameters={"columns": input_columns}, + ) + chunk_size = len(d) + + if chunk_size == 0: + print(f"Skipping chunk {i + 1} / {n_chunks} since it is empty") + continue + + shear_data = self.process_anacal_shear_data(d) + if not created_files: + created_files = True + shear_outfile = self.setup_output( + "shear_catalog", "shear", shear_data, n, + ) + shear_outfile["shear"].attrs["catalog_type"] = "anacal" + + shear_end = shear_start + len(shear_data["ra"]) + self.write_output( + shear_outfile, "shear", shear_data, shear_start, shear_end, + ) + + print( + f"Processing chunk {i + 1} / {n_chunks} into rows " + f"{shear_start:,} - {shear_end:,}" + ) + shear_start = shear_end + + print("Trimming shear columns:") + for col in shear_data.keys(): + print(" ", col) + h5py_shorten(shear_outfile["shear"], col, shear_end) + + shear_outfile.close() + + def file_run(self): + n, dtypes = self.get_meta("anacal_catalog") + cols = self.setup_input() + + file = self.open_input("anacal_catalog") + data = file[1][cols] + shear_data = self.process_anacal_shear_data(data) + + shear_outfile = self.setup_output( + "shear_catalog", "shear", shear_data, n, + ) + shear_outfile["shear"].attrs["catalog_type"] = "anacal" + + print("Trimming shear columns:") + for col in shear_data.keys(): + print(" ", col) + h5py_shorten(shear_outfile["shear"], col, len(shear_data["ra"])) + + shear_outfile.close() + + def setup_input(self): + prefix = self.config["prefix"] + scale = self.config["scale"] + cols = [ + "ra", + "dec", + "wsel", + "n_mask_base", + f"{prefix}_e1", + f"{prefix}_e2", + f"{prefix}_m00", + f"{prefix}_m20", + ] + cols += ["dwsel" + suffix for suffix in ["_dg1", "_dg2"]] + cols += [ + prefix + delta + suffix + for delta in ["_de1", "_de2", "_dm00", "_dm20"] + for suffix in ["_dg1", "_dg2"] + ] + bands = self.config["bands"] + # i-band S/N + shear response, always from the fpfs1 family. + # ``scale`` below only picks the flux used for magnitudes; the + # brightness cut consumes the pre-computed fpfs1 S/N so it is + # consistent regardless of the mag scale. + cols += [ + "lsst_i_s2n_fpfs1", + "lsst_i_ds2n_fpfs1_dg1", + "lsst_i_ds2n_fpfs1_dg2", + ] + # Pre-computed AB magnitudes + shear responses on the DP1-v3 + # merged catalog (xlens.add_magnitude_columns at the fixed + # MAG_ZERO_AB zeropoint) — passed through so downstream stages + # can consume mags without redoing the nanojansky→mag math. + # Both ``mag`` and its error ``mag_err`` carry shear responses. + for b in bands: + cols += [ + f"lsst_{b}_mag_{scale}", + f"lsst_{b}_dmag_{scale}_dg1", + f"lsst_{b}_dmag_{scale}_dg2", + f"lsst_{b}_mag_{scale}_err", + f"lsst_{b}_dmag_{scale}_err_dg1", + f"lsst_{b}_dmag_{scale}_err_dg2", + ] + # Band-combined shape magnitude ``esq = e1^2 + e2^2`` and its + # shear derivatives, emitted by xlens.MergePipe on the + # WCS-corrected fpfs1 shape. Consumed by the |e| from the merged catalog. + cols += self._extinction_source_columns() + + # PSF second moments -> per-band psf_g1/psf_g2/psf_T_mean (see + # process_anacal_shear_data), one set per band in ``bands``. The + # flag column lets us NaN out failed HSM measurements. + if self.config["add_psf_properties"]: + for b in bands: + base = f"lsst_{b}_ext_shapeHSM_HsmPsfMoments" + cols += [ + f"{base}_xx", + f"{base}_yy", + f"{base}_xy", + f"{base}_flag", + ] + + # zmode_0 → mean_z; zmode_1p, zmode_1m, zmode_2p, zmode_2m → the + # metacal-style shifted variants (built with dg=0.01 in xlens' + # photoZPipe, so TXSourceSelectorAnaCal must use delta_gamma=0.01). + cols += [f"zmode_{s}" for s in PZ_SUFFIXES] + + return cols + + def process_anacal_shear_data(self, data): + bands = self.config["bands"] + s = self.config["scale"] + prefix = self.config["prefix"] + dg = self.config["delta_gamma"] + # Column names work for both an astropy Table (butler mode) and a + # numpy/FITS structured array (file mode). + colnames = list(getattr(data, "colnames", None) or data.dtype.names) + # The dm computed e1/e2 columns store the pre-multiplied observable + # e_meas = wsel · e_raw, and "weight" is uniformly set to 1. + # This way downstream GGCorrelation with weight_column="weight" + # computes xi_e = Σ (wsel_i e_i)(wsel_j e_j) / N_pairs instead of a + # ⟨wsel wsel⟩-weighted mean of raw shapes. + # xi_g = xi_e / ^2 + + # The raw shapes and wsel are + # still exposed as separate columns (wsel, e1_raw, e2_raw) so + # TXSourceSelectorAnaCal can compute R_shape (⟨wsel · de/dg⟩) and + # R_detect (⟨(dwsel/dg) · e_raw⟩). + wsel = data["wsel"][:] + e1_raw = data[f"{prefix}_e1"][:] + e2_raw = data[f"{prefix}_e2"][:] + m00 = data[f"{prefix}_m00"][:] + m20 = data[f"{prefix}_m20"][:] + + output = { + "ra": data["ra"][:], + "dec": data["dec"][:], + "weight": np.ones_like(wsel), # uniform 1 for treecorr + "wsel": wsel, # raw wsel (for R_shape) + "n_mask_base": data["n_mask_base"][:], + "weight_dg1": data["dwsel_dg1"][:], + "weight_dg2": data["dwsel_dg2"][:], + "e1": wsel * e1_raw, # e_meas ≡ wsel · e_raw + "e2": wsel * e2_raw, + "e1_raw": e1_raw, # raw shape (for R_detect) + "e2_raw": e2_raw, + "m00": m00, + "m20": m20, + } + for delta in ["de1", "de2", "dm00", "dm20"]: + output[f"{delta}_dg1"] = data[f"{prefix}_{delta}_dg1"][:] + output[f"{delta}_dg2"] = data[f"{prefix}_{delta}_dg2"][:] + + # i-band S/N + shear response — passed through from the + # pre-computed fpfs1 columns. ``scale`` only picks the flux + # family for magnitudes, not S/N. + s2n = data["lsst_i_s2n_fpfs1"][:] + ds2n_dg1 = data["lsst_i_ds2n_fpfs1_dg1"][:] + ds2n_dg2 = data["lsst_i_ds2n_fpfs1_dg2"][:] + output["s2n"] = s2n + output["ds2n_dg1"] = ds2n_dg1 + output["ds2n_dg2"] = ds2n_dg2 + + # Per-band AB magnitudes come pre-computed on the v3 merged + # catalog (xlens.add_magnitude_columns writes them at the fixed + # MAG_ZERO_AB zeropoint with smooth truncation and the analytic + # dmag/dg shear responses), so TXPipe simply forwards them + # rather than redoing nanojansky→mag inside the ingest. ``mag`` + # and ``mag_err`` both carry ``d*/dg{1,2}`` shear derivatives — + # exposed as ``dmag_{band}_dg{c}`` / ``dmag_err_{band}_dg{c}`` on + # the shear catalog so downstream stages can build ±γ variants + # of every quantity a mag-based cut consumes. + for band in bands: + b = f"lsst_{band}" + output[f"mag_{band}"] = data[f"{b}_mag_{s}"][:] + output[f"mag_err_{band}"] = data[f"{b}_mag_{s}_err"][:] + for d in ("dg1", "dg2"): + output[f"dmag_{band}_{d}"] = data[f"{b}_dmag_{s}_{d}"][:] + output[f"dmag_err_{band}_{d}"] = ( + data[f"{b}_dmag_{s}_err_{d}"][:] + ) + + # Per-band extinction a_, from the merged catalog's + # lsst_a_ column (the mags above are NOT dereddened, so these + # give downstream a place to apply extinction). + for band in bands: + output[f"a_{band}"] = data[f"lsst_a_{band}"][:] + + # Band-combined shape magnitude + shear derivatives — feeds the + # |e|`` for each band in the config ``bands``.""" + return [f"lsst_a_{band}" for band in self.config["bands"]] + + def get_catalog_size(self, butler, dataset_type): + import pyarrow.parquet + + n = 0 + for ref in butler.query_datasets(dataset_type): + uri = butler.getURI(ref) + if not uri.path.endswith(".parq"): + raise ValueError( + f"Some data in dataset {dataset_type} was not in " + f"parquet format: {uri.path}" + ) + with pyarrow.parquet.ParquetFile(uri.path) as f: + n += f.metadata.num_rows + return n diff --git a/txpipe/ingest/dp1_details.py b/txpipe/ingest/dp1_details.py new file mode 100644 index 000000000..ad217e36e --- /dev/null +++ b/txpipe/ingest/dp1_details.py @@ -0,0 +1,57 @@ +# All TRACT INFORMATION SHOULD BE MOVED ELSEWHERE +DP1_COSMOLOGY_FIELDS = [ + "EDFS", + "ECDFS", + "LGLF", +] + + +DP1_TRACTS = { + # Euclid Deep Field South + "EDFS": [2393, 2234, 2235, 2394], + # Extended Chandra Deep Field South + "ECDFS": [5062, 5063, 5064, 4848, 4849], + # Low Galactic Latitude Field / Rubin_SV_095_-25 + "LGLF": [5305, 5306, 5525, 5526], + # Fornax Dwarf Spheroidal Galaxy + "FDSG": [4016, 4217, 4218, 4017], + # Low Ecliptic Latitude Field / Rubin_SV_38_7 + "LELF": [10464, 10221, 10222, 10704, 10705, 10463], + # Seagull Nebula + "Seagull": [7850, 7849, 7610, 7611], + # 47 Tuc Globular Cluster + "47Tuc": [531, 532, 453, 454], +} + +DP1_COSMOLOGY_TRACTS = sum([DP1_TRACTS[_field] for _field in DP1_COSMOLOGY_FIELDS], []) +ALL_TRACTS = sum(DP1_TRACTS.values(), []) + + +# In case useful later: +DP1_FIELD_CENTERS = { + "47 Tuc Globular Cluster": (6.02, -72.08), + "Low Ecliptic Latitude Field": (37.86, 6.98), + "Fornax Dwarf Spheroidal Galaxy": (40.00, -34.45), + "Extended Chandra Deep Field South": (53.13, -28.10), + "Euclid Deep Field South": (59.10, -48.73), + "Low Galactic Latitude Field": (95.00, -25.00), + "Seagull Nebula": (106.23, -10.51), +} + + +DP1_SURVEY_PROPERTIES = { + "deepCoadd_exposure_time_consolidated_map_sum": "Total exposure time accumulated per sky position (second)", + "deepCoadd_epoch_consolidated_map_min": "Earliest observation epoch (MJD)", + "deepCoadd_epoch_consolidated_map_max": "Latest observation epoch (MJD)", + "deepCoadd_epoch_consolidated_map_mean": "Mean observation epoch (MJD)", + "deepCoadd_psf_size_consolidated_map_weighted_mean": "Weighted mean of PSF characteristic width as computed from the determinant radius (pixel)", + "deepCoadd_psf_e1_consolidated_map_weighted_mean": "Weighted mean of PSF ellipticity component e1", + "deepCoadd_psf_e2_consolidated_map_weighted_mean": "Weighted mean of PSF ellipticity component e2", + "deepCoadd_psf_maglim_consolidated_map_weighted_mean": "Weighted mean of PSF flux 5σ magnitude limit (magAB)", + "deepCoadd_sky_background_consolidated_map_weighted_mean": "Weighted mean of background light level from the sky (nJy)", + "deepCoadd_sky_noise_consolidated_map_weighted_mean": "Weighted mean of standard deviation of the sky level (nJy)", + "deepCoadd_dcr_dra_consolidated_map_weighted_mean": "Weighted mean of DCR-induced astrometric shift in right ascension direction, expressed as a proportionality factor", + "deepCoadd_dcr_ddec_consolidated_map_weighted_mean": "Weighted mean of DCR-induced astrometric shift in declination direction, expressed as a proportionality factor", + "deepCoadd_dcr_e1_consolidated_map_weighted_mean": "Weighted mean of DCR-induced change in PSF ellipticity (e1), expressed as a proportionality factor", + "deepCoadd_dcr_e2_consolidated_map_weighted_mean": "Weighted mean of DCR-induced change in PSF ellipticity (e2), expressed as a proportionality factor", +} diff --git a/txpipe/psf_diagnostics.py b/txpipe/psf_diagnostics.py index 6de43adf8..87c4619ab 100644 --- a/txpipe/psf_diagnostics.py +++ b/txpipe/psf_diagnostics.py @@ -1134,6 +1134,11 @@ def load_galaxies(self): g2 = g["g2"][:][mask] weight = g["weight"][:][mask] + elif cat_type == "anacal": + g1 = g["e1"][:][mask] + g2 = g["e2"][:][mask] + weight = g["weight"][:][mask] + else: g1 = g["g1"][:][mask] g2 = g["g2"][:][mask] diff --git a/txpipe/shear_calibration/__init__.py b/txpipe/shear_calibration/__init__.py index 89f32f186..fcbb8f529 100644 --- a/txpipe/shear_calibration/__init__.py +++ b/txpipe/shear_calibration/__init__.py @@ -5,9 +5,18 @@ LensfitCalibrator, HSCCalibrator, MetaDetectCalibrator, + AnaCalibrator ) -from .calibration_calculators import MetacalCalculator, LensfitCalculator, HSCCalculator, MetaDetectCalculator, MockCalculator +from .calibration_calculators import ( + MetacalCalculator, + LensfitCalculator, + HSCCalculator, + MetaDetectCalculator, + MockCalculator, + AnaCalCalculator + ) + from .mean_shear_in_bins import MeanShearInBins from .names import band_variants, metacal_variants, metadetect_variants, META_VARIANTS from .utils import BinStats diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 9430f6315..688de9655 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -1,6 +1,6 @@ import numpy as np from .names import META_VARIANTS -from .calibrators import MetaCalibrator, LensfitCalibrator, HSCCalibrator, MetaDetectCalibrator, NullCalibrator +from .calibrators import MetaCalibrator, LensfitCalibrator, HSCCalibrator, MetaDetectCalibrator, NullCalibrator, AnaCalibrator from .utils import BinStats class _DataWrapper: @@ -58,7 +58,7 @@ class CalibrationCalculator: The selection function does not need to know about all these variants. The calculator will wrap the data dictionary passed in in a special class that chooses variant columns when they are looked up. So your selection function can just ask for, e.g., - "T", "s2n", or "mag_r", and the calculator will make sure it gets the right variant of + "T", "s2n", or "mag_r", and the calculator will make sure it gets the right variant of that column for each selection. The final results are in the form of a BinStats object, which contains: @@ -314,8 +314,8 @@ class MetaDetectCalculator(CalibrationCalculator): """A calibration and stats calculator for metadetect catalogs. See the CalibrationCalculator class for the use and contents of this class, - but note that the attributes of the class are different because we have to - keep track of 5 variants separately, and the shear_stats keeps track of the + but note that the attributes of the class are different because we have to + keep track of 5 variants separately, and the shear_stats keeps track of the mean and std dev of all of the variants. """ @@ -850,3 +850,192 @@ def collect(self, comm=None, allgather=False) -> BinStats: return bin_stats +class AnaCalCalculator(CalibrationCalculator): + """Calibration and stats calculator for AnaCal catalogs + + See the CalibrationCalculator class for the use and contents of this class + """ + + def __init__(self, selector, delta_gamma): + """ + Initialize the Calibrator using the funtion you will use to select + objects. That function should take at least one argument, + the chunk of data to select on. It should look up the original + names of the columns to select on; the calculator wraps the data + with _DataWrapper so the same call also picks up ±γ variants when + needed (mask + s2n + size + zbin cuts all go through this single + code path). + + The selector can take further *args and **kwargs, passed in when + adding data. + + Parameters + ---------- + selector: function + Function that selects objects. Must apply every cut whose + shear response contributes to the selection bias (mask, s2n, + size, zbin) — the calculator no longer applies any cut of its + own. + delta_gamma: float + Half-difference in applied g between the ±1 variants + (each variant is at ±delta_gamma from baseline). + """ + from parallel_statistics import ParallelMean + super().__init__(selector) + + self.delta_gamma = delta_gamma + # e_meas ≡ wsel · e is the observable. All accumulators + # use uniform weight=1 and per-source values are pre-multiplied by wsel + # where the chain-rule term needs it, so each collect() returns a + # plain sample mean ⟨wsel · X⟩ (never divided by ⟨wsel⟩). + # Note that wsel is smooth detection weight and 96% of them are 1. + # Ensemble shear estimator: + # γ = ⟨e_meas⟩ / R_total, R_total = R_shape + R_detect + R_sel + # 0: ⟨wsel · de1/dg1⟩, 1: ⟨wsel · de2/dg2⟩ + self.shape_response = ParallelMean(size=2) + # 0: ⟨(dwsel/dg1)·e1⟩, 1: ⟨(dwsel/dg2)·e2⟩ + self.detect_response = ParallelMean(size=2) + # 0..3: ⟨wsel · e · 𝟙[sel_±]⟩ + self.sel_response = ParallelMean(size=4) + + def add_data(self, data, *args, **kwargs): + """Select objects from a new chunk of data and tally their responses + + Parameters + ---------- + data: dict + Dictionary of data columns to select on and add + *args + Positional arguments to be passed to the selection function + **kwargs + Keyword arguments to be passed to the selection function. + + Returns + ------- + sel: array + The indicies of the objects selected from this chunk of data + """ + # Baseline + four shifted views on the same chunk. _DataWrapper + # routes ``data_1p["s2n"] -> data["s2n_1p"]``, + # ``data_1p["m00"] -> data["m00_1p"]``, + # ``data_1p["zbin"] -> data["zbin_1p"]``, etc. + # Columns without a matching variant (n_mask_base, e1, weight, ...) + # fall back to baseline. + data_00 = _DataWrapper(data, "") + data_1p = _DataWrapper(data, "_1p") + data_1m = _DataWrapper(data, "_1m") + data_2p = _DataWrapper(data, "_2p") + data_2m = _DataWrapper(data, "_2m") + + # A single selector call per variant applies mask + s2n + size + + # zbin cuts uniformly. No per-cut special casing lives here now. + select = self.selector(data_00, *args, **kwargs) + sel_1p = self.selector(data_1p, *args, **kwargs) + sel_1m = self.selector(data_1m, *args, **kwargs) + sel_2p = self.selector(data_2p, *args, **kwargs) + sel_2m = self.selector(data_2m, *args, **kwargs) + + # Convention-A column plan (see TXIngestAnacal): + # e1, e2 – e_meas ≡ wsel · e_raw (pre-multiplied, feed treecorr) + # e1_raw, e2_raw – raw shape (needed for R_detect numerator) + # wsel – raw wsel (needed for R_shape numerator, Neff, μ) + # weight – uniform 1 (feed treecorr's weight_column) + e1 = data_00["e1"] # = wsel · e_raw + e2 = data_00["e2"] + e1_raw = data_00["e1_raw"] + e2_raw = data_00["e2_raw"] + wsel = data_00["wsel"] + weight_dg1 = data_00["weight_dg1"] + weight_dg2 = data_00["weight_dg2"] + de1_dg1 = data_00["de1_dg1"] + de2_dg2 = data_00["de2_dg2"] + + # Baseline slices. + w_sel = wsel[select] + e1_sel = e1[select] # already wsel · e_raw + e2_sel = e2[select] + + # Bookkeeping: raw count + wsel sums (used for Kish's Neff). + self.count += e1_sel.size + self.sum_weights += np.sum(w_sel) + self.sum_sq_weights += np.sum(w_sel ** 2) + + # All accumulators use weight=1 (uniform). Per-source values are + # pre-multiplied by wsel where the chain-rule term needs it, so each + # ParallelMean returns ⟨wsel · X⟩ over the fed sample. + + # R_shape numerator: ⟨wsel · de_raw/dg⟩ over the baseline sample. + self.shape_response.add_data(0, w_sel * de1_dg1[select]) + self.shape_response.add_data(1, w_sel * de2_dg2[select]) + # R_detect numerator: ⟨(dwsel/dg) · e_raw⟩ over the baseline sample. + self.detect_response.add_data(0, weight_dg1[select] * e1_raw[select]) + self.detect_response.add_data(1, weight_dg2[select] * e2_raw[select]) + # R_sel numerator: ⟨wsel · e_raw⟩ over each ±γ variant sample — + # equivalently ⟨e_meas⟩ since e1, e2 are already pre-multiplied. + # (Spin-2 argument: N_1p, N_1m, N_2p, N_2m converge in the large-N + # ensemble, so per-subset means telescope cleanly under the finite + # difference; residual bias ~ (δN/N)·⟨e⟩, zero in ensemble.) + self.sel_response.add_data(0, e1[sel_1p]) + self.sel_response.add_data(1, e1[sel_1m]) + self.sel_response.add_data(2, e2[sel_2p]) + self.sel_response.add_data(3, e2[sel_2m]) + + # μ numerator: mean of e_meas (already pre-multiplied) over baseline. + self.shear_stats.add_data(0, e1_sel) + self.shear_stats.add_data(1, e2_sel) + + return select + + def collect(self, comm=None, allgather=False) -> BinStats: + """ + Finalize and sum up all the response values, and return a BinStats + obejct that collections calibration and statistics. + + Parameters + ---------- + comm: MPI Communicator + If supplied, all processors response values will be combined together. + All processes will return the same final value + allgather: bool + If True, the response values will be returned for all the processors. + + Returns + ------- + bin_stats: BinStats + An object containing the final calibration and statistics for this bin. + """ + # One helper for every scalar accumulator; comm.reduce returns + # None on non-root ranks, which we let propagate through Neff. + def _reduce(x): + if comm is None: + return x + return comm.allreduce(x) if allgather else comm.reduce(x) + + count = _reduce(self.count) + sum_weights = _reduce(self.sum_weights) + sum_sq_weights = _reduce(self.sum_sq_weights) + + # Each ParallelXxx does its own MPI reduction inside .collect(). + mode = "allgather" if allgather else "gather" + _, shape_means = self.shape_response.collect(comm, mode) # (2,) means + _, detect_means = self.detect_response.collect(comm, mode) # (2,) means + _, sel_means = self.sel_response.collect(comm, mode) # (4,) means + _, mean_e, var_e = self.shear_stats.collect(comm, mode) + + # Convention A: every mean is a plain sample mean ⟨wsel · X⟩ over + # the baseline (shape, detect, μ) or ±γ variant (sel) sample. + # No ⟨wsel⟩ denominator anywhere — R_total already contains wsel + # via the pre-multiplication in add_data. + R_shape = 0.5 * (shape_means[0] + shape_means[1]) + R_detect = 0.5 * (detect_means[0] + detect_means[1]) + R_sel_1 = (sel_means[0] - sel_means[1]) / (2.0 * self.delta_gamma) + R_sel_2 = (sel_means[2] - sel_means[3]) / (2.0 * self.delta_gamma) + R_sel = 0.5 * (R_sel_1 + R_sel_2) + R_total = R_shape + R_detect + R_sel + + Neff = None if sum_weights is None else sum_weights ** 2 / sum_sq_weights + + calibrator = AnaCalibrator(R_total, mean_e, mu_is_weighted=False) + sigma_e = calibrator.calibrate_variance_to_sigma_e(var_e) + sigma = calibrator.calibrate_sigma(np.sqrt(var_e)) + return BinStats(count, Neff, calibrator.mu, sigma_e, sigma, calibrator) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index ec8f1fe16..14d96a8ee 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -58,6 +58,8 @@ def load(cls, tomo_file, null=False): subcls = LensfitCalibrator elif cat_type == "hsc": subcls = HSCCalibrator + elif cat_type == "anacal": + subcls = AnaCalibrator else: raise ValueError(f"Unknown catalog type {cat_type} in tomo file") @@ -630,3 +632,86 @@ def calibrate_sigma(self, sigma): """ return np.array(sigma) / (2 * self.R) / (1 + self.K) + +class AnaCalibrator(MetaCalibrator): + """Stores information needed to calibrate an AnaCal shear method.""" + def __init__(self, R, mu, mu_is_weighted=True): + self.R = R + if mu_is_weighted: + self.mu = np.array(mu) + else: + self.mu = np.array(mu) / R + assert("Anacal needs an already calibrated mu.") + + def apply(self, g1, g2, weights, subtract_mean=True): + """ + Calibrate a set of shears using the response matrix and + mean shear substraction + Parameters + ---------- + g1: array or float + Shear 1 component + + g2: array or float + Shear 2 component + + subtract_mean: bool + whether to subtract mean shear (default True) + """ + if not subtract_mean: + g1, g2 = weights * [g1, g2] / self.R + elif np.isscalar(g1): + g1, g2 = weights * [g1, g2] / self.R - self.mu + else: + g1, g2 = weights * [g1, g2] / self.R - self.mu[:, np.newaxis] + return g1, g2 + + def calibrate_variance_to_sigma_e(self, var_e): + return np.sqrt(0.5 *np.sum(var_e)) / self.R + + def calibrate_sigma(self, sigma): + return np.array(sigma) / self.R + + @classmethod + def load(cls, tomo_file): + """ + Make a set of AnaCal calibrators using the info in a tomography file. + + You can use the parent Calibrator.load to automatically + load the correct subclass. + + Parameters + ---------- + tomo_file: str + A tomography file name the cal factors are read from + + Returns + ------- + cals: list + A set of AnaCalibrators, one per bin + """ + import h5py + + with h5py.File(tomo_file, "r") as f: + mu1 = f["counts/mean_e1"][:] + mu2 = f["counts/mean_e2"][:] + n = len(mu1) + R = f["response/R"][:n] + R_2d = f["response/R_2d"][0] + mu1_2d = f["counts/mean_e1_2d"][0] + mu2_2d = f["counts/mean_e2_2d"][0] + + calibrators = [cls(R[i], [mu1[i], mu2[i]]) for i in range(n)] + calibrator2d = cls(R_2d, [mu1_2d, mu2_2d]) + return calibrators, calibrator2d + + def save(self, outfile, i): + if i == "2d": + outfile["response/R_2d"][0] = self.R + outfile["counts/mean_e1_2d"][0] = self.mu[0] + outfile["counts/mean_e2_2d"][0] = self.mu[1] + else: + outfile["response/R"][i] = self.R + outfile["counts/mean_e1"][i] = self.mu[0] + outfile["counts/mean_e2"][i] = self.mu[1] + diff --git a/txpipe/source_selection/__init__.py b/txpipe/source_selection/__init__.py index 4fac2c893..f5c6461d8 100644 --- a/txpipe/source_selection/__init__.py +++ b/txpipe/source_selection/__init__.py @@ -4,3 +4,4 @@ from .metacal import TXSourceSelectorMetacal from .metadetect import TXSourceSelectorMetadetect from .tomography import TXSourceTomography +from .anacal import TXSourceSelectorAnacal diff --git a/txpipe/source_selection/anacal.py b/txpipe/source_selection/anacal.py new file mode 100644 index 000000000..4bfec82b3 --- /dev/null +++ b/txpipe/source_selection/anacal.py @@ -0,0 +1,269 @@ +from .base import TXSourceSelectorBase +from ..shear_calibration import AnaCalCalculator, band_variants +import numpy as np +from ceci.config import StageParameter + + +class TXSourceSelectorAnacal(TXSourceSelectorBase): + """ + Source selection and tomography for AnaCal catalogs. + + This selector subclass is designed for anacal-type catalogs like those + that DESC plans to produce for Rubin data. + """ + + name = "TXSourceSelectorAnaCal" + + config_options = { + **TXSourceSelectorBase.config_options, + "delta_gamma": StageParameter( + float, + required=True, + msg="Delta gamma value for the AnaCal response calculation", + ), + "mask_threshold": StageParameter( + float, 40, + msg="mask threshold, for when to mask objects.", + ), + # AnaCal-specific defaults for the two base cuts. These override + # the base class's required-with-no-default declarations so a + # minimal yaml block still runs; users can still set them + # explicitly in the yaml if they want tighter/looser cuts. + "s2n_cut": StageParameter( + float, 5.0, + msg="AnaCal S/N cut (i-band flux / flux_err).", + ), + "T_cut": StageParameter( + float, 0.1, + msg="Band-combined size cut (m00 + m20) / m00 > T_cut on the " + "fpfs1 moments emitted by xlens.MergePipe.", + ), + # Band-combined shape magnitude cut |e| < emax on the ``esq`` + # column emitted by xlens' MergePipe (esq = e1**2 + e2**2 on the + # WCS-corrected fpfs1 shape). The selector cuts on esq < emax**2 + # internally so this stays a familiar |e| threshold. + "emax": StageParameter( + float, 0.5, + msg="Band-combined shape magnitude cut |e| < emax " + "(applied as esq < emax**2 on the merge-stage esq column).", + ), + # Per-band AB mag upper bounds — one StageParameter per band, + # matching the ``r_hi_cut`` / ``i_hi_cut`` convention in + # lens_selector.py. Default 50 is above the xlens + # smooth-truncation cap (MAG_CAP=40), so the default effectively + # applies no cut; tighten in the yaml to gate on brightness. + "g_hi_cut": StageParameter(float, 50.0, msg="Upper g-mag cut."), + "r_hi_cut": StageParameter(float, 50.0, msg="Upper r-mag cut."), + "i_hi_cut": StageParameter(float, 50.0, msg="Upper i-mag cut."), + "z_hi_cut": StageParameter(float, 50.0, msg="Upper z-mag cut."), + } + + def data_iterator(self): + """ + This iterator returns chunks of data in dictionaries one by one. + + We call to a parent class method to do the main iteration; the work + here is just choosing which columns to read. + """ + + bands = self.config["bands"] + shear_cols = [ + "ra", + "dec", + "weight", + "wsel", # raw wsel — needed by R_shape numerator + "n_mask_base", + "e1", # e_meas ≡ wsel · e_raw (pre-multiplied) + "e2", + "e1_raw", # raw e — needed by R_detect numerator + "e2_raw", + "m00", + "m20", + "de1_dg1", + "de2_dg2", + "dm00_dg1", + "dm00_dg2", + "dm20_dg1", + "dm20_dg2", + "s2n", + "ds2n_dg1", + "ds2n_dg2", + "weight_dg1", + "weight_dg2", + ] + # ±γ shifted variants of every quantity the selector cuts on, + # precomputed by TXIngestAnacal — consumed via _DataWrapper by + # AnaCalCalculator (data["s2n_1p"], data["m00_1p"], ...). + for suf in ("1p", "1m", "2p", "2m"): + shear_cols += [f"s2n_{suf}", f"m00_{suf}", f"m20_{suf}", f"esq_{suf}"] + shear_cols += [f"mag_{b}_{suf}" for b in bands] + shear_cols += band_variants( + bands, "mag", "mag_err", shear_catalog_type="anacal", + ) + # Per-band shear-response derivatives — ``dmag_{b}_dg{c}`` is + # what ``add_sheared_variant_columns`` uses to build the ±γ + # variants of the per-band mag cut. ``dmag_err_{b}_dg{c}`` is + # carried alongside for downstream code that wants a shear + # response on the magnitude error. + for b in bands: + shear_cols += [ + f"dmag_{b}_dg1", f"dmag_{b}_dg2", + f"dmag_err_{b}_dg1", f"dmag_err_{b}_dg2", + ] + # Band-combined shape magnitude and its shear derivatives, + # emitted by xlens.MergePipe — feeds the |e|= edges[zi]) & (zz < edges[zi + 1]) + b[m] = zi + pz_data[f"zbin_{suf}"] = b + + def setup_output(self): + """ + Prepare the output columns for the response values generated by Anacal + """ + outfile = super().setup_output() + n = outfile["tomography/bin"].size + group = outfile.create_group("response") + group.create_dataset("R", (n, 1, 1), dtype="f") + group.create_dataset("R_2d", (1,), dtype="f") + return outfile + + def setup_response_calculators(self, nbin_source): + delta_gamma = self.config["delta_gamma"] + calculators = [ + AnaCalCalculator( + select_anacal_tomographic_weak_lensing_sample, delta_gamma, + ) + for i in range(nbin_source) + ] + calculators.append( + AnaCalCalculator( + select_anacal_weak_lensing_sample, delta_gamma, + ) + ) + return calculators + + def write_tomography(self, outfile, start, end, source_bin, r): + super().write_tomography(outfile, start, end, source_bin, r) + group = outfile["response"] + group["R"][start:end] = r + + +def select_anacal_weak_lensing_sample( + data, config, calling_from_select=False): + """Baseline AnaCal weak-lensing selection. + + Applies six cuts, each looked up through the ``_DataWrapper`` so the + ±γ shifted variants (esq_1p, mag_g_1p, m00_1p, s2n_1p, zbin_1p, ...) + are used automatically when the caller wraps ``data`` with a suffix: + + 1. n_mask_base < mask_threshold (shear-independent, no variants). + 2. s2n > s2n_cut (variants from ds2n_dg{1,2}). + 3. Band-combined size cut (m00 + m20) / m00 > T_cut on the fpfs1 + moments (variants automatically use m00_{1p,...}/m20_{1p,...} = + m + dg * dm_dg{c}). + 4. Band-combined shape magnitude esq < emax**2 + (variants from desq_dg{1,2}, both derived from the WCS-corrected + fpfs1 shape by xlens.MergePipe). + 5. Per-band mag_{b} < {b}_hi_cut for every band in config.bands + (variants from mag_{b}_dg{1,2}). + 6. zbin >= 0 (variants from shifted mean_z). + """ + s2n_cut = config["s2n_cut"] + t_cut = config["T_cut"] + esq_max = config["emax"] ** 2 + verbose = config["verbose"] + + flag = data["n_mask_base"] < config["mask_threshold"] + n0 = len(flag) + sel = flag + f1 = sel.sum() / n0 + + sel &= data["s2n"] > s2n_cut + f2 = sel.sum() / n0 + + # Band-combined size cut on fpfs1 moments. + m00 = data["m00"] + m20 = data["m20"] + sel &= (m00 + m20) / m00 > t_cut + f3 = sel.sum() / n0 + + sel &= data["esq"] < esq_max + f4 = sel.sum() / n0 + + for b in config["bands"]: + sel &= data[f"mag_{b}"] < config[f"{b}_hi_cut"] + f5 = sel.sum() / n0 + + sel &= data["zbin"] >= 0 + f6 = sel.sum() / n0 + + if verbose and calling_from_select: + print( + f"Tomo selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.2%} size, " + f"{f4:.2%} |e|, {f5:.2%} mag, ", + end="", + ) + elif verbose: + print( + f"2D selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.2%} size, " + f"{f4:.2%} |e|, {f5:.2%} mag, {f6:.2%} any z bin" + ) + print("total 2D", sel.sum()) + return sel + + +def select_anacal_tomographic_weak_lensing_sample(data, config, bin_index): + zbin = data["zbin"] + sel = select_anacal_weak_lensing_sample( + data, config, calling_from_select=True, + ) + sel &= zbin == bin_index + return sel diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index c721067f8..0d691f343 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -4,6 +4,7 @@ MetaDetectCalculator, MetaCalibrator, NullCalibrator, + AnaCalCalculator ) import numpy as np @@ -22,6 +23,8 @@ def select_all_where(data): # we just want to select everything here too return np.where(data["g2"] * 0 == 0) +def select_all_anacal(data): + return np.repeat(True, data["e2"].size) @@ -161,6 +164,77 @@ def core_metadet(comm): assert stats.source_count == N * nproc +def core_anacal(comm): + # Half-step convention: each ±1 variant is at ±delta_gamma from baseline, + # matching xlens' photoZPipe DISTORTIONS (±0.01) and the on-disk zmode + # variants in the merged catalog. Not to be confused with metacal, where + # delta_gamma above is the full 1p−1m separation (0.02 there). + delta_gamma = 0.01 + nproc = 1 if comm is None else comm.size + N = 10 + + R_shape_true = 0.4 # known test value + R_weight_true = 0.3 + + # Convention A: TXIngestAnacal exposes the pre-multiplied observable + # "e1"/"e2" (= wsel · e_raw), plus separate "e1_raw"/"e2_raw" and "wsel" + # columns. With wsel=1 uniformly (as in this test) the pre-multiplied + # and raw shapes coincide, and Σweight_dg·e / Σweight equals the plain + # ⟨weight_dg · e_raw⟩ that the new R_detect accumulator computes. + e1 = np.random.normal(0, 0.1, size=N) + e2 = np.random.normal(0, 0.1, size=N) + base_data = { + "e1": e1, + "e2": e2, + "e1_raw": e1, + "e2_raw": e2, + "wsel": np.ones(N), + "weight": np.ones(N), + "weight_dg1": np.zeros(N), + "weight_dg2": np.zeros(N), + "de1_dg1": np.zeros(N), + "de2_dg2": np.zeros(N), + "m00": np.ones(N), + "m20": np.ones(N), + "dm00_dg1": np.zeros(N), + "dm00_dg2": np.zeros(N), + "dm20_dg1": np.zeros(N), + "dm20_dg2": np.zeros(N), + "n_mask_base": np.zeros(N) + } + + # case 1: pure shape response + data = {**base_data, "de1_dg1": np.full(N, R_shape_true), "de2_dg2": np.full(N, R_shape_true)} + cal = AnaCalCalculator(select_all_anacal, delta_gamma) + cal.add_data(data) + stats = cal.collect(comm, allgather=True) + assert np.allclose(stats.calibrator.R, R_shape_true) + assert stats.source_count == N * nproc + + # case 2: pure weight-bias response + data = {**base_data, "weight_dg1": np.full(N, R_weight_true), "weight_dg2": np.full(N, R_weight_true)} + cal = AnaCalCalculator(select_all_anacal, delta_gamma) + cal.add_data(data) + stats = cal.collect(comm, allgather=True) + expected_R_weight = 0.5 * (np.sum(data["weight_dg1"] * data["e1"]) + + np.sum(data["weight_dg2"] * data["e2"])) / np.sum(data["weight"]) + assert np.allclose(stats.calibrator.R, expected_R_weight) + + # Case 3: both contributions add correctly + data = {**base_data, + "de1_dg1": np.full(N, R_shape_true), + "de2_dg2": np.full(N, R_shape_true), + "weight_dg1": np.full(N, R_weight_true), + "weight_dg2": np.full(N, R_weight_true), + } + cal = AnaCalCalculator(select_all_anacal, delta_gamma) + cal.add_data(data) + stats = cal.collect(comm, allgather=True) + expected_R_weight = 0.5 * (np.sum(data["weight_dg1"] * data["e1"]) + + np.sum(data["weight_dg2"] * data["e2"])) / np.sum(data["weight"]) + assert np.allclose(stats.calibrator.R, R_shape_true + expected_R_weight) + + def test_metacalibrator_serial(): core_metacal(None) @@ -168,6 +242,9 @@ def test_metacalibrator_serial(): def test_metadetect_serial(): core_metadet(None) +def test_anacal_serial(): + core_anacal(None) + def test_metadetect_parallel(): mockmpi.mock_mpiexec(2, core_metadet) @@ -376,3 +453,4 @@ def test_null_mean(): test_metadetect_parallel() test_mean_shear_no_weights() test_mean_shear_weights() + test_anacal_serial() diff --git a/txpipe/utils/__init__.py b/txpipe/utils/__init__.py index c6c3e3925..fbbe3a9b9 100755 --- a/txpipe/utils/__init__.py +++ b/txpipe/utils/__init__.py @@ -3,7 +3,12 @@ from .misc import array_hash, unique_list, hex_escape, rename_iterated, read_shear_catalog_type from .healpix import dilated_healpix_map from .splitters import Splitter, DynamicSplitter -from .conversion import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, moments_to_shear, mag_ab_to_nanojansky +from .conversion import ( + nanojansky_err_to_mag_ab, + nanojansky_to_mag_ab, + moments_to_shear, + mag_ab_to_nanojansky, +) from .timer import Timer from .debuggable_dask import import_dask from .mpi_utils import in_place_reduce, mpi_reduce_large \ No newline at end of file diff --git a/txpipe/utils/conversion.py b/txpipe/utils/conversion.py index bb601e810..66ced1b5c 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -13,11 +13,15 @@ def mag_ab_to_nanojansky(mag): def nanojansky_err_to_mag_ab(flux, flux_err): - return 2.5 / np.log(10) * (flux_err / flux) + return np.divide(2.5 / np.log(10) * flux_err, + flux, + out=np.full_like(flux, np.nan, dtype=np.float64), + where=flux > 0, + ) def moments_to_shear(Ixx, Iyy, Ixy): b = Ixx + Iyy + 2 * np.sqrt(Ixx * Iyy - Ixy**2) e1 = (Ixx - Iyy) / b e2 = 2 * Ixy / b - return e1, e2 + return e1, e2 \ No newline at end of file