From cac4c6c3a931cd2a8b36615f588242a5220bdcff Mon Sep 17 00:00:00 2001 From: empEvil Date: Tue, 4 Nov 2025 13:40:47 +0000 Subject: [PATCH 01/78] inital commit of anacal ingester --- txpipe/ingest/anacal.py | 124 ++++++++++++++++++++++++++++++ txpipe/utils/__init__.py | 2 +- txpipe/utils/calibration_tools.py | 2 + 3 files changed, 127 insertions(+), 1 deletion(-) create mode 100644 txpipe/ingest/anacal.py diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py new file mode 100644 index 000000000..84a9668ab --- /dev/null +++ b/txpipe/ingest/anacal.py @@ -0,0 +1,124 @@ +from base import TXIngestCatalogFits +from ..data_types import ShearCatalog, PhotometryCatalog, HDFFile, FileCollection, FitsFile +from .lsst import process_photometry_data, process_shear_data +from ceci.config import StageParameter +import numpy as np +from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, +from ..utils.hdf_tools import h5py_shorten, repack + +class TXIngestAnacal(TXIngestCatalogFits): + """ + Ingest an anacal catalog!, + """ + + name = "TXIngestAnacal" + input = [ + ("Aanacal_catalog", FitsFile) + ] + outputs = [ + ("photometry_catalog", PhotometryCatalog), + ("shear_catalog", ShearCatalog), + ("exposures", HDFFile), + ("survey_propety_maps", FileCollection), + ] + config_options = { + "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(str, "grizy", msg="string of flux bands"), + "scale": StageParameter(str, "gauss2", msg="scale radius for the convolution with Gaussian PSF") + } + + def run(self): + tracts = self.config["tracts"] + file_path = self.config["file_path"] + + n, dtypes = self.get_meta(f"{file_path}/anacal_anacal_table.fits") + cols = self.setup_input("shear_catalog") + prefix = self.config["prefix"] + + 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" + #self.write_output(shear_outfile, "shear", shear_data) + + 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"] + cols = ( + [ + "ra", + "dec", + "wsel", + "wdet", + f"{prefix}_e1", + f"{prefix}_e2", + ]) + cols += ["dwsel"+ suffix for suffix in ["_dg1", "_dg2"]] + cols += [ + prefix +delta + suffix + for delta in ["_de1", "_de2"] + for suffix in ["_dg1", "_dg2"] + ] + bands = self.config["bands"] + for i in range(4): + cols += [band + "_flux_gauss" + i for band in bands] + cols += [band + "_flux_gauss" + i + "_err" for band in bands] + cols += [ + band + "_dflux_gauss" + i + suffix + for band in bands + for suffix in ["_dg1", "_dg2"] + ] + return cols + + def process_anacal_shear_data(data): + bands = self.config["bands"] + s = self.config["scale"] + output = { + "ra": data["ra"], + "dec": data["dec"], + "weight":data["wsel"], + "weight_detection": data["wdet"], + "weight_dg1": data["dwsel_dg1"], + "weight_dg2": data["dwsel_dg2"], + } + for band in bands: + f = data[f"{band}_flux_{s}"] + f_err = data[f"{band}_flux_{s}_err"] + output[f"mag_{band}"] = nanojansky_to_mag_ab(f) + output[f"mag_err_{band}"] = nanojansky_err_to_mag_ab(f, f_err) + + if band == "i": + output["s2n"] = f / f_err + + for d in ["_dg1", "_dg2"]: + dd = data[f"{band}_dflux_{s}"+d] + output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) + + return output + + def setup_output(self, tag, group, first_chunk, n): + f = self.open_output(tag) + g = f.create_group(group) + + for name, col in first_chunk.items(): + g.create_dataset(name, shape=(n,), dtype=col.dtype) + return f + + def write_output(self, outfile, group, data, start, end): + g = outfile[group] + for name, col in data.items(): + # replace masked values with nans + if np.ma.isMaskedArray(col): + col = col.filled(np.nan) + g[name][start:end] = col + diff --git a/txpipe/utils/__init__.py b/txpipe/utils/__init__.py index 6f02561ce..371489bfc 100755 --- a/txpipe/utils/__init__.py +++ b/txpipe/utils/__init__.py @@ -11,7 +11,7 @@ HSCCalibrator, ) from .splitters import Splitter, DynamicSplitter -from .calibration_tools import read_shear_catalog_type, band_variants, metacal_variants, metadetect_variants +from .calibration_tools import read_shear_catalog_type, band_variants, metacal_variants, metadetect_variants, anacal_variants from .calibration_tools import MetacalCalculator, LensfitCalculator, MeanShearInBins from .conversion import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, moments_to_shear, mag_ab_to_nanojansky from .timer import Timer diff --git a/txpipe/utils/calibration_tools.py b/txpipe/utils/calibration_tools.py index 333f7049d..eec482723 100755 --- a/txpipe/utils/calibration_tools.py +++ b/txpipe/utils/calibration_tools.py @@ -22,6 +22,8 @@ def metacal_variants(*names): def metadetect_variants(*names): return [f"{group}/{name}" for group in ["00", "1p", "1m", "2p", "2m"] for name in names] +def anacal_variants(*names): + return [name + suffix for suffix in ["_dg1", "_dg2"] for name in names] def band_variants(bands, *names, shear_catalog_type="metacal"): if shear_catalog_type == "metacal": From 52093b0cbbc7de6c7aacd0163584e9b71ba4722a Mon Sep 17 00:00:00 2001 From: empEvil Date: Tue, 4 Nov 2025 16:56:58 +0000 Subject: [PATCH 02/78] fixed, first attempt --- examples/Anacal_test/config.yml | 5 ++++ examples/Anacal_test/pipeline.yml | 20 ++++++++++++++ txpipe/ingest/__init__.py | 1 + txpipe/ingest/anacal.py | 46 ++++++++++++++----------------- 4 files changed, 47 insertions(+), 25 deletions(-) create mode 100644 examples/Anacal_test/config.yml create mode 100644 examples/Anacal_test/pipeline.yml diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml new file mode 100644 index 000000000..eb16e9f98 --- /dev/null +++ b/examples/Anacal_test/config.yml @@ -0,0 +1,5 @@ +global: + chunk_rows: 100000 + sparse: true # Generate sparse maps - faster if using small areas + pixelization: healpix + nside: 512 \ No newline at end of file diff --git a/examples/Anacal_test/pipeline.yml b/examples/Anacal_test/pipeline.yml new file mode 100644 index 000000000..1434cdb63 --- /dev/null +++ b/examples/Anacal_test/pipeline.yml @@ -0,0 +1,20 @@ +stages: + - name: TXIngestAnacal + +output_dir: data/example/output_anacal + +launcher: + name: mini + interval: 1.0 + +site: + name: local + max_threads: 2 + +modules: > + txpipe + +config: examples/Anacal_test/config.yml + +inputs: + anacal_catalog: /Users/eskemollerpedersen/Research/TXPipe/Catalog_examples/Anacal/anacal_table.fits diff --git a/txpipe/ingest/__init__.py b/txpipe/ingest/__init__.py index 0b381c8d7..2dafdfec9 100644 --- a/txpipe/ingest/__init__.py +++ b/txpipe/ingest/__init__.py @@ -11,3 +11,4 @@ ) from .dp1 import TXIngestDataPreview1 from .legacy import TXIngestDESY3Gold +from .anacal import TXIngestAnacal \ No newline at end of file diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 84a9668ab..3c2bf936d 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -1,9 +1,9 @@ -from base import TXIngestCatalogFits +from .base import TXIngestCatalogFits from ..data_types import ShearCatalog, PhotometryCatalog, HDFFile, FileCollection, FitsFile from .lsst import process_photometry_data, process_shear_data from ceci.config import StageParameter import numpy as np -from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, +from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab from ..utils.hdf_tools import h5py_shorten, repack class TXIngestAnacal(TXIngestCatalogFits): @@ -12,14 +12,11 @@ class TXIngestAnacal(TXIngestCatalogFits): """ name = "TXIngestAnacal" - input = [ - ("Aanacal_catalog", FitsFile) + inputs = [ + ("anacal_catalog", FitsFile) ] outputs = [ - ("photometry_catalog", PhotometryCatalog), ("shear_catalog", ShearCatalog), - ("exposures", HDFFile), - ("survey_propety_maps", FileCollection), ] config_options = { "tracts": StageParameter(str, "", msg="Comma-separated list of tracts to use (empty for all)."), @@ -30,13 +27,12 @@ class TXIngestAnacal(TXIngestCatalogFits): def run(self): tracts = self.config["tracts"] - file_path = self.config["file_path"] - n, dtypes = self.get_meta(f"{file_path}/anacal_anacal_table.fits") - cols = self.setup_input("shear_catalog") + n, dtypes = self.get_meta("anacal_catalog") + cols = self.setup_input() prefix = self.config["prefix"] - file = self.open_input("Anacal_catalog") + file = self.open_input("anacal_catalog") data = file[1][cols] shear_data = self.process_anacal_shear_data(data) @@ -70,30 +66,30 @@ def setup_input(self): for suffix in ["_dg1", "_dg2"] ] bands = self.config["bands"] - for i in range(4): - cols += [band + "_flux_gauss" + i for band in bands] - cols += [band + "_flux_gauss" + i + "_err" for band in bands] + for i in [0, 2, 4]: + cols += [band + "_flux_gauss" + f"{i}" for band in bands] + cols += [band + "_flux_gauss" + f"{i}" + "_err" for band in bands] cols += [ - band + "_dflux_gauss" + i + suffix + band + "_dflux_gauss" + f"{i}" + suffix for band in bands for suffix in ["_dg1", "_dg2"] ] return cols - def process_anacal_shear_data(data): + def process_anacal_shear_data(self, data): bands = self.config["bands"] s = self.config["scale"] output = { - "ra": data["ra"], - "dec": data["dec"], - "weight":data["wsel"], - "weight_detection": data["wdet"], - "weight_dg1": data["dwsel_dg1"], - "weight_dg2": data["dwsel_dg2"], + "ra": data["ra"][:], + "dec": data["dec"][:], + "weight":data["wsel"][:], + "weight_detection": data["wdet"][:], + "weight_dg1": data["dwsel_dg1"][:], + "weight_dg2": data["dwsel_dg2"][:], } for band in bands: - f = data[f"{band}_flux_{s}"] - f_err = data[f"{band}_flux_{s}_err"] + f = data[f"{band}_flux_{s}"][:] + f_err = data[f"{band}_flux_{s}_err"][:] output[f"mag_{band}"] = nanojansky_to_mag_ab(f) output[f"mag_err_{band}"] = nanojansky_err_to_mag_ab(f, f_err) @@ -101,7 +97,7 @@ def process_anacal_shear_data(data): output["s2n"] = f / f_err for d in ["_dg1", "_dg2"]: - dd = data[f"{band}_dflux_{s}"+d] + dd = data[f"{band}_dflux_{s}"+d][:] output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) return output From abba4699afd1e0260d0ce8b7a3f3851c359fe3bb Mon Sep 17 00:00:00 2001 From: empEvil Date: Fri, 17 Apr 2026 14:27:50 +0100 Subject: [PATCH 03/78] fixing the ingestion part --- txpipe/ingest/anacal.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 3c2bf936d..4f6b90b70 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -53,16 +53,18 @@ def setup_input(self): cols = ( [ "ra", - "dec", + "dec", "wsel", "wdet", f"{prefix}_e1", f"{prefix}_e2", + f"{prefix}_m0", + f"{prefix}_m2" ]) cols += ["dwsel"+ suffix for suffix in ["_dg1", "_dg2"]] cols += [ prefix +delta + suffix - for delta in ["_de1", "_de2"] + for delta in ["_de1", "_de2", "_dm0", "_dm2"] for suffix in ["_dg1", "_dg2"] ] bands = self.config["bands"] @@ -79,6 +81,7 @@ def setup_input(self): def process_anacal_shear_data(self, data): bands = self.config["bands"] s = self.config["scale"] + prefix = self.config["prefix"] output = { "ra": data["ra"][:], "dec": data["dec"][:], @@ -86,7 +89,14 @@ def process_anacal_shear_data(self, data): "weight_detection": data["wdet"][:], "weight_dg1": data["dwsel_dg1"][:], "weight_dg2": data["dwsel_dg2"][:], + "g1": data[f"{prefix}_e1"][:], + "g2": data[f"{prefix}_e2"][:], + "m0": data[f"{prefix}_m0"][:], + "m2": data[f"{prefix}_m2"][:], } + for delta in ["de1", "de2", "dm0", "dm2"]: + output[f"{delta}_dg1"] = data[f"{prefix}_{delta}_dg1"] + output[f"{delta}_dg2"] = data[f"{prefix}_{delta}_dg2"] for band in bands: f = data[f"{band}_flux_{s}"][:] f_err = data[f"{band}_flux_{s}_err"][:] From fd42ff6574e74594fdb85233b7d8245037ae5363 Mon Sep 17 00:00:00 2001 From: empEvil Date: Fri, 17 Apr 2026 14:30:46 +0100 Subject: [PATCH 04/78] updating test configs --- examples/Anacal_test/config.yml | 6 +++++- examples/Anacal_test/pipeline.yml | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index eb16e9f98..434e91749 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -2,4 +2,8 @@ global: chunk_rows: 100000 sparse: true # Generate sparse maps - faster if using small areas pixelization: healpix - nside: 512 \ No newline at end of file + nside: 512 + +TXIngestAnacal: + bands: griz + \ No newline at end of file diff --git a/examples/Anacal_test/pipeline.yml b/examples/Anacal_test/pipeline.yml index 1434cdb63..a197ceb91 100644 --- a/examples/Anacal_test/pipeline.yml +++ b/examples/Anacal_test/pipeline.yml @@ -17,4 +17,4 @@ modules: > config: examples/Anacal_test/config.yml inputs: - anacal_catalog: /Users/eskemollerpedersen/Research/TXPipe/Catalog_examples/Anacal/anacal_table.fits + anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits \ No newline at end of file From 1c99217c299bc5ec0a33e459964ed8a5ea951e6f Mon Sep 17 00:00:00 2001 From: empEvil Date: Fri, 17 Apr 2026 14:33:11 +0100 Subject: [PATCH 05/78] more input edits --- examples/Anacal_test/pipeline.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/Anacal_test/pipeline.yml b/examples/Anacal_test/pipeline.yml index a197ceb91..b9ba1ed8f 100644 --- a/examples/Anacal_test/pipeline.yml +++ b/examples/Anacal_test/pipeline.yml @@ -17,4 +17,6 @@ modules: > config: examples/Anacal_test/config.yml inputs: - anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits \ No newline at end of file + anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits + +resume: False \ No newline at end of file From 8126a27139317e44a694ec553f1533ed2ac7e74b Mon Sep 17 00:00:00 2001 From: empEvil Date: Fri, 17 Apr 2026 14:38:17 +0100 Subject: [PATCH 06/78] typo fix --- txpipe/ingest/anacal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 4f6b90b70..b70e285ec 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -95,8 +95,8 @@ def process_anacal_shear_data(self, data): "m2": data[f"{prefix}_m2"][:], } for delta in ["de1", "de2", "dm0", "dm2"]: - output[f"{delta}_dg1"] = data[f"{prefix}_{delta}_dg1"] - output[f"{delta}_dg2"] = data[f"{prefix}_{delta}_dg2"] + output[f"{delta}_dg1"] = data[f"{prefix}_{delta}_dg1"][:] + output[f"{delta}_dg2"] = data[f"{prefix}_{delta}_dg2"][:] for band in bands: f = data[f"{band}_flux_{s}"][:] f_err = data[f"{band}_flux_{s}_err"][:] From 96914db05a9c184bad67f7bedf4fd49451305395 Mon Sep 17 00:00:00 2001 From: empEvil Date: Fri, 17 Apr 2026 15:27:34 +0100 Subject: [PATCH 07/78] change naming --- txpipe/ingest/anacal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index b70e285ec..77144abbc 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -89,8 +89,8 @@ def process_anacal_shear_data(self, data): "weight_detection": data["wdet"][:], "weight_dg1": data["dwsel_dg1"][:], "weight_dg2": data["dwsel_dg2"][:], - "g1": data[f"{prefix}_e1"][:], - "g2": data[f"{prefix}_e2"][:], + "e1": data[f"{prefix}_e1"][:], + "e2": data[f"{prefix}_e2"][:], "m0": data[f"{prefix}_m0"][:], "m2": data[f"{prefix}_m2"][:], } From 09bf4e30671c05fc3d58b1954a3ca498e910c66a Mon Sep 17 00:00:00 2001 From: empEvil Date: Mon, 27 Apr 2026 15:09:36 +0100 Subject: [PATCH 08/78] WIP on anacal calibrator --- .../calibration_calculators.py | 141 ++++++++++++++++++ txpipe/shear_calibration/calibrators.py | 28 ++++ 2 files changed, 169 insertions(+) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index f786e8215..f812c9f7f 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -849,4 +849,145 @@ 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, selectora): + """ + 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 selector can take further *args and **kwargs, passed in when adding + data. + + Parameters + ---------- + selector: function + Function that selects objects + delta_gamma: float + The difference in applied g between 1p and 1m variants + """ + from parallel_statistics import ParallelMean + super().__init__(selector) + + self.cal_bias_means = ParallelMean(size=5) + #self.sel_bias_means = ParallelMean(size=8) + + 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 + """ + + select = self.selector(data, *args, **kwargs) + + + e1 = data["e1"] + e2 = data["e2"] + weight = data["weight"] + weight_dg1 = data["weight_dg1"] + weight_dg2 = data["weight_dg2"] + de1_dg1 = data["de1_dg1"] + de1_dg2 = data["de1_dg2"] + de2_dg1 = data["de2_dg1"] + de2_dg2 = data["de2_dg2"] + + n = e1[select].size + + self.count += n + self.sum_weights += np.sum(weight[select]) + self.sum_sq_weights += np.sum(weight[select]**2) + + w00 = weight[select] + R00 = weight_dg1[select] * e1[select] + weight[select] * de1_dg1[select] + R01 = weight_dg2[select] * e1[select] + weight[select] * de1_dg2[select] + R10 = weight_dg1[select] * e2[select] + weight[select] * de2_dg1[select] + R11 = weight_dg2[select] * e2[select] + weight[select] * de2_dg2[select] + + self.cal_bias_means.add_data(0, R00, w00) + self.cal_bias_means.add_data(1, R01, w00) + self.cal_bias_means.add_data(2, R10, w00) + self.cal_bias_means.add_data(3, R11, w00) + self.cal_bias_means.add_data(4, w00, w00) + + self.shear_stats.add_data(0, e1[select], w00) + self.shear_stats.add_data(1, e2[select], w00) + + 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. + """ + # collect all the things we need + if comm is not None: + if allgather: + count = comm.allreduce(self.count) + sum_weights = comm.allreduce(self.sum_weights) + sum_sq_weights = comm.allreduce(self.sum_sq_weights) + else: + count = comm.reduce(self.count) + sum_weights = comm.reduce(self.sum_weights) + sum_sq_weights = comm.reduce(self.sum_sq_weights) + else: + count = self.count + sum_weights = self.sum_weights + sum_sq_weights = self.sum_sq_weights + + # Collect the mean values we need + mode = "allgather" if allgather else "gather" + _, R = self.cal_bias_means.collect(comm, mode) + _, mean_e, var_e = self.shear_stats.collect(comm, mode) + + # Unpack the flat mean R: + R_mean = np.zeros((2, 2)) + R_mean[0, 0] = R[0] + R_mean[0, 1] = R[1] + R_mean[1, 0] = R[2] + R_mean[1, 1] = R[3] + w00 = R[4] + + if sum_weights is None: + Neff = None + else: + Neff = sum_weights**2 / sum_sq_weights + + calibrator = AnaCalibrator(R_mean, mean_e, w00, mu_is_weighted=False) + sigma_e = calibrator.calibrate_variance_to_sigma_e(var_e) + sigma = calibrator.calibrate_sigma(np.sqrt(var_e)) + bin_stats = BinStats(count, Neff, calibrator.mu, sigma_e, sigma, calibrator) + return bin_stats + + diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index ec8f1fe16..594bedd1a 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -630,3 +630,31 @@ def calibrate_sigma(self, sigma): """ return np.array(sigma) / (2 * self.R) / (1 + self.K) + +def AnaCalibrator(Calibrator): + """Stores information needed to calibrate a Anacal shear method""" + def __init__(self, R, mu, weights, mu_is_weigthed=True): + self.R = R + self.Rinv = np.linalg.inv(R) + self.w = weights + if mu_is_weigthed: + self.mu = np.array(mu) + else: + self.mu = self.Rinv @ (weights * mu ) + + def apply(self, g1, g2 substract_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) + """ + """ From 21103dd81cc4a64dd3a6a6a8dd564b60052e3e9e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Mon, 27 Apr 2026 20:06:20 +0100 Subject: [PATCH 09/78] WIP simplifying the calibration --- .../calibration_calculators.py | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index f812c9f7f..70c99c012 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -855,7 +855,7 @@ class AnaCalCalculator(CalibrationCalculator): See the CalibrationCalculator class for the use and contents of this class """ - def __init__(self, selectora): + def __init__(self, selector): """ Initialize the Calibrator using the funtion you will use to select objects. That function should take at least one argument, @@ -875,7 +875,7 @@ def __init__(self, selectora): from parallel_statistics import ParallelMean super().__init__(selector) - self.cal_bias_means = ParallelMean(size=5) + self.response_means = ParallelMean(size=4) #self.sel_bias_means = ParallelMean(size=8) def add_data(self, data, *args, **kwargs): @@ -905,8 +905,6 @@ def add_data(self, data, *args, **kwargs): weight_dg1 = data["weight_dg1"] weight_dg2 = data["weight_dg2"] de1_dg1 = data["de1_dg1"] - de1_dg2 = data["de1_dg2"] - de2_dg1 = data["de2_dg1"] de2_dg2 = data["de2_dg2"] n = e1[select].size @@ -915,17 +913,17 @@ def add_data(self, data, *args, **kwargs): self.sum_weights += np.sum(weight[select]) self.sum_sq_weights += np.sum(weight[select]**2) - w00 = weight[select] - R00 = weight_dg1[select] * e1[select] + weight[select] * de1_dg1[select] - R01 = weight_dg2[select] * e1[select] + weight[select] * de1_dg2[select] - R10 = weight_dg1[select] * e2[select] + weight[select] * de2_dg1[select] - R11 = weight_dg2[select] * e2[select] + weight[select] * de2_dg2[select] + wsel = weight[select] + de1_dg1_sub = de1_dg1[select] + de2_dg2_sub = de2_dg2[select] + dwsel_dg1 = weight_dg1[select] + dwsel_dg2 = weight_dg2[select] + + self.response_means.add_data(0, de1_dg1_sub, wsel) + self.response_means.add_data(1, de2_dg2_sub, wsel) + self.response_means.add_data(2, dwsel_dg1, wsel) + self.response_means.add_data(3, dwsel_dg2, wsel) - self.cal_bias_means.add_data(0, R00, w00) - self.cal_bias_means.add_data(1, R01, w00) - self.cal_bias_means.add_data(2, R10, w00) - self.cal_bias_means.add_data(3, R11, w00) - self.cal_bias_means.add_data(4, w00, w00) self.shear_stats.add_data(0, e1[select], w00) self.shear_stats.add_data(1, e2[select], w00) From 5881da1823e9630b67b87d9c9880b2b77d3d9356 Mon Sep 17 00:00:00 2001 From: empEvil Date: Tue, 28 Apr 2026 12:32:59 +0100 Subject: [PATCH 10/78] added calibration calculator --- txpipe/ingest/anacal.py | 10 ++- .../calibration_calculators.py | 84 +++++++++++++++---- txpipe/shear_calibration/calibrators.py | 15 +++- 3 files changed, 86 insertions(+), 23 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 77144abbc..9deab1263 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -56,6 +56,7 @@ def setup_input(self): "dec", "wsel", "wdet", + "mask_value", f"{prefix}_e1", f"{prefix}_e2", f"{prefix}_m0", @@ -85,16 +86,17 @@ def process_anacal_shear_data(self, data): output = { "ra": data["ra"][:], "dec": data["dec"][:], - "weight":data["wsel"][:], + "weight": data["wsel"][:], + "mask_value": data["mask_value"][:], "weight_detection": data["wdet"][:], "weight_dg1": data["dwsel_dg1"][:], "weight_dg2": data["dwsel_dg2"][:], "e1": data[f"{prefix}_e1"][:], "e2": data[f"{prefix}_e2"][:], - "m0": data[f"{prefix}_m0"][:], - "m2": data[f"{prefix}_m2"][:], + "m00": data[f"{prefix}_m00"][:], + "m20": data[f"{prefix}_m20"][:], } - for delta in ["de1", "de2", "dm0", "dm2"]: + 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"][:] for band in bands: diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 70c99c012..02cf3a493 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -855,7 +855,7 @@ class AnaCalCalculator(CalibrationCalculator): See the CalibrationCalculator class for the use and contents of this class """ - def __init__(self, selector): + 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, @@ -874,8 +874,13 @@ def __init__(self, selector): """ from parallel_statistics import ParallelMean super().__init__(selector) - + + self.delta_gamma = delta_gamma self.response_means = ParallelMean(size=4) + self.sel_response_e1p1 = 0 + self.sel_response_e1m1 = 0 + self.sel_response_e2p2 = 0 + self.sel_response_e2m2 = 0 #self.sel_bias_means = ParallelMean(size=8) def add_data(self, data, *args, **kwargs): @@ -898,7 +903,6 @@ def add_data(self, data, *args, **kwargs): select = self.selector(data, *args, **kwargs) - e1 = data["e1"] e2 = data["e2"] weight = data["weight"] @@ -907,16 +911,41 @@ def add_data(self, data, *args, **kwargs): de1_dg1 = data["de1_dg1"] de2_dg2 = data["de2_dg2"] + m00 = data["m00"] + m20 = data["m20"] + + dm00_dg1 = data["dm00_dg1"] + dm00_dg2 = data["dm00_dg2"] + dm20_dg1 = data["dm20_dg1"] + dm20_dg2 = data["dm20_dg2"] + mask_value = data["mask_value"] + n = e1[select].size + # Record the count for this chunk, for summation later self.count += n self.sum_weights += np.sum(weight[select]) self.sum_sq_weights += np.sum(weight[select]**2) + # Next we calculate the cuts per variation masks wsel = weight[select] + mask = mask_value[select] < 40 + + mask_p1 = mask & self.get_submask(m00, m20, dm00_dg1, dm20_dg1, +1) + mask_m1 = mask & self.get_submask(m00, m20, dm00_dg1, dm20_dg1, -1) + mask_p2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, +1) + mask_m2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, -1) + + self.sel_response_e1p1 += np.sum(weight[mask_p1]*e1[mask_p1]) + self.sel_response_e1m1 += np.sum(weight[mask_m1]*e1[mask_m1]) + self.sel_response_e2p2 += np.sum(weight[mask_p2]*e2[mask_p2]) + self.sel_response_e2m2 += np.sum(weight[mask_m2]*e2[mask_m2]) + + # Next we find the means part needed for the Shape response, + # and the Weight-bias response de1_dg1_sub = de1_dg1[select] de2_dg2_sub = de2_dg2[select] - dwsel_dg1 = weight_dg1[select] + dwsel_dg1 = weight_dg1[select] dwsel_dg2 = weight_dg2[select] self.response_means.add_data(0, de1_dg1_sub, wsel) @@ -924,12 +953,17 @@ def add_data(self, data, *args, **kwargs): self.response_means.add_data(2, dwsel_dg1, wsel) self.response_means.add_data(3, dwsel_dg2, wsel) - - self.shear_stats.add_data(0, e1[select], w00) - self.shear_stats.add_data(1, e2[select], w00) + self.shear_stats.add_data(0, e1[select], wsel) + self.shear_stats.add_data(1, e2[select], wsel) return select + def get_submask(self, m00, m20, dm00_dg, dm20_dg, sign): + m0 = m00 + sign * self.delta_gamma * dm00_dg + m2 = m20 + sign * self.delta_gamma * dm20_dg + return (m0 + m2)/ m0 > 0.1 + + def collect(self, comm=None, allgather=False) -> BinStats: """ Finalize and sum up all the response values, and return a BinStats @@ -954,34 +988,54 @@ def collect(self, comm=None, allgather=False) -> BinStats: count = comm.allreduce(self.count) sum_weights = comm.allreduce(self.sum_weights) sum_sq_weights = comm.allreduce(self.sum_sq_weights) + sum_sel_response_e1p1 = comm.allreduce(self.sel_response_e1p1) + sum_sel_response_e1m1 = comm.allreduce(self.sel_response_e1m1) + sum_sel_response_e2p2 = comm.allreduce(self.sel_response_e2p2) + sum_sel_response_e2m2 = comm.allreduce(self.sel_response_e2m2) else: count = comm.reduce(self.count) sum_weights = comm.reduce(self.sum_weights) sum_sq_weights = comm.reduce(self.sum_sq_weights) + sum_sel_response_e1p1 = comm.reduce(self.sel_response_e1p1) + sum_sel_response_e1m1 = comm.reduce(self.sel_response_e1m1) + sum_sel_response_e2p2 = comm.reduce(self.sel_response_e2p2) + sum_sel_response_e2m2 = comm.reduce(self.sel_response_e2m2) else: count = self.count sum_weights = self.sum_weights sum_sq_weights = self.sum_sq_weights + sum_sel_response_e1p1 = self.sel_response_e1p1 + sum_sel_response_e1m1 = self.sel_response_e1m1 + sum_sel_response_e2p2 = self.sel_response_e2p2 + sum_sel_response_e2m2 = self.sel_response_e2m2 # Collect the mean values we need mode = "allgather" if allgather else "gather" - _, R = self.cal_bias_means.collect(comm, mode) + _, R = self.response_means.collect(comm, mode) _, mean_e, var_e = self.shear_stats.collect(comm, mode) # Unpack the flat mean R: - R_mean = np.zeros((2, 2)) - R_mean[0, 0] = R[0] - R_mean[0, 1] = R[1] - R_mean[1, 0] = R[2] - R_mean[1, 1] = R[3] - w00 = R[4] + mean_de1_dg1 = R[0] + mean_de2_dg2 = R[1] + mean_dwsel_dg1_e1 = R[2] + mean_dwsel_dg2_e2 = R[3] + + # Reducing down the responses + R_shape = 0.5 * (mean_de1_dg1 + mean_de2_dg2) + R_weight = 0.5 * (mean_dwsel_dg1_e1 + mean_dwsel_dg2_e2) + + R_sel_1 = (sum_sel_response_e1p1 - sum_sel_response_e1m1) / (2.0 * self.delta_gamma) / count + R_sel_2 = (sum_sel_response_e2p2 - sum_sel_response_e2m2) / (2.0 * self.delta_gamma) / count + R_sel = 0.5 * (R_sel_1 + R_sel_2) + + R_total = R_shape + R_weight + R_sel if sum_weights is None: Neff = None else: Neff = sum_weights**2 / sum_sq_weights - calibrator = AnaCalibrator(R_mean, mean_e, w00, mu_is_weighted=False) + calibrator = AnaCalibrator(R_total, mean_e, w00, mu_is_weighted=False) sigma_e = calibrator.calibrate_variance_to_sigma_e(var_e) sigma = calibrator.calibrate_sigma(np.sqrt(var_e)) bin_stats = BinStats(count, Neff, calibrator.mu, sigma_e, sigma, calibrator) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index 594bedd1a..ae562c974 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -631,7 +631,7 @@ def calibrate_sigma(self, sigma): return np.array(sigma) / (2 * self.R) / (1 + self.K) -def AnaCalibrator(Calibrator): +class AnaCalibrator(Calibrator): """Stores information needed to calibrate a Anacal shear method""" def __init__(self, R, mu, weights, mu_is_weigthed=True): self.R = R @@ -640,9 +640,9 @@ def __init__(self, R, mu, weights, mu_is_weigthed=True): if mu_is_weigthed: self.mu = np.array(mu) else: - self.mu = self.Rinv @ (weights * mu ) + self.mu = self.Rinv @ (weights * mu) - def apply(self, g1, g2 substract_mean=True): + def apply(self, g1, g2, substract_mean=True): """ Calibrate a set of shears using the response matrix and mean shear substraction @@ -657,4 +657,11 @@ def apply(self, g1, g2 substract_mean=True): subtract_mean: bool whether to subtract mean shear (default True) """ - """ + if not substract_mean: + g1, g2 = self.Rinv @ self.w * [g1, g2] + elif np.isscalar(g1): + g1, g2 = self.Rinv @ self.w * [g1, g2] - self.mu + else: + g1, g2 = self.Rinv @ self.w * [g1, g2] - self.mu[:, np.newaxis] + return g1, g2 + From 7966f16320ac313fd49f2910a70b547f856c09e0 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 5 May 2026 16:09:44 +0100 Subject: [PATCH 11/78] WIP --- txpipe/shear_calibration/__init__.py | 11 +++- txpipe/shear_calibration/calibrators.py | 51 ++++++++++++++---- txpipe/source_selection/Anacal.py | 70 +++++++++++++++++++++++++ 3 files changed, 122 insertions(+), 10 deletions(-) create mode 100644 txpipe/source_selection/Anacal.py 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/calibrators.py b/txpipe/shear_calibration/calibrators.py index ae562c974..1b280f9bf 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -631,18 +631,16 @@ def calibrate_sigma(self, sigma): return np.array(sigma) / (2 * self.R) / (1 + self.K) -class AnaCalibrator(Calibrator): +class AnaCalibrator(MetaCalibrator): """Stores information needed to calibrate a Anacal shear method""" - def __init__(self, R, mu, weights, mu_is_weigthed=True): + def __init__(self, R, mu, mu_is_weigthed=True): self.R = R - self.Rinv = np.linalg.inv(R) - self.w = weights if mu_is_weigthed: self.mu = np.array(mu) else: - self.mu = self.Rinv @ (weights * mu) + assert("Anacal needs an already calibrated mu.") - def apply(self, g1, g2, substract_mean=True): + def apply(self, g1, g2, weights, substract_mean=True): """ Calibrate a set of shears using the response matrix and mean shear substraction @@ -658,10 +656,45 @@ def apply(self, g1, g2, substract_mean=True): whether to subtract mean shear (default True) """ if not substract_mean: - g1, g2 = self.Rinv @ self.w * [g1, g2] + g1, g2 = weights * [g1, g2] / self.R elif np.isscalar(g1): - g1, g2 = self.Rinv @ self.w * [g1, g2] - self.mu + g1, g2 = weights * [g1, g2] / self.R - self.mu else: - g1, g2 = self.Rinv @ self.w * [g1, g2] - self.mu[:, np.newaxis] + g1, g2 = weights * [g1, g2] / self.R - self.mu[:, np.newaxis] return g1, g2 + + @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 AnacaltCalibrators, one per bin + """ + import h5py + + with h5py.File(tomo_file, "r") as f: + R = f["response/R"][:] + n = len(R) + + mu1 = f["counts/mean_e1"][:] + mu2 = f["counts/mean_e2"][:] + + calibrators = [cls(R[i], [mu1[i], mu2[i]]) for i in range(n)] + return calibrators + + def save(self, outfile, i): + outfile["respponse/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/Anacal.py b/txpipe/source_selection/Anacal.py new file mode 100644 index 000000000..cedbb8a52 --- /dev/null +++ b/txpipe/source_selection/Anacal.py @@ -0,0 +1,70 @@ +from .base import TXSourceSelectorBase +from .base import select_weak_lensing_sample, select_tomographic_weak_lensing_sample +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 subcallss is designed for anacal-type catalogs like those + that DESC plans to produce for Ruin data. + """ + + name = "TXSourceSelectorAnaCal" + + config_options = { + **TXSourceSelectorBase.config_options, + "delta_gamma": StageParameter( + float, + required=True, + msg= "Delta gmamma value for hte AnaCal response calculation" + ), + } + + 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", + "mask_value", + "e1", + "e2", + "m00", + "m20", + "de1_dg1", + "de2_dg2", + "dm00_dg1", + "dm00_dg2", + "dm20_dg1", + "dm20_dg2" + ] + shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="Anacal") + + if self.config["input_pz"]: + shear_cols += ["mean_z"] + elif self.config["true_z"]: + shear_cols += ["redshift_true"] + + chunk_rows = self.config["chunk_rows"] + return self.iterate_hdf("shear_catalog", "shear", shear_cols, chunk_rows) + + def setup_output(self): + """ + Prepare the output columns for the response values generated bby Anacal + """ + outfile = super().setup_output() + n = outfile["count/counts"].size + group = outfile.create_group("response") + group.create_dataset("R", (n, 1, 1), dtype="f") + return outfile + From 6fb51b9407336871c2316bc2cbfccfe8fad8f1b4 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 5 May 2026 16:27:54 +0100 Subject: [PATCH 12/78] initial version of source selector --- txpipe/source_selection/Anacal.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index cedbb8a52..d02f0a631 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -68,3 +68,14 @@ def setup_output(self): group.create_dataset("R", (n, 1, 1), dtype="f") return outfile + def setup_response_calculators(self, nbin_source): + delta_gamma = self.config["delta_gamma"] + calculators = [AnaCalCalculator(select_tomographic_weak_lensing_sample, delta_gamma) for i in range(nbin_source)] + calculators.append(AnaCalCalculator(select_weak_lensing_sample, delta_gamma)) + return calculators + + def write_tomography(self, outfile, start, end, source_bin, R): + super().write_tomogrpahy(outfile, start, end, source_bin, R) + group = outfile["response"] + group["R"][start:end] = R + \ No newline at end of file From 5c444a180f2191ca907d0b20a9a14c76203d3ddf Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Wed, 27 May 2026 17:25:14 +0100 Subject: [PATCH 13/78] adding missing stage to init --- txpipe/__init__.py | 1 + 1 file changed, 1 insertion(+) 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 From 3cc70174532c0b0aa5aea5222588d6c40f411332 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Wed, 27 May 2026 17:27:08 +0100 Subject: [PATCH 14/78] more import --- txpipe/source_selection/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/txpipe/source_selection/__init__.py b/txpipe/source_selection/__init__.py index 4fac2c893..29ad90cd4 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 From a7b56a6196c431af9fa56dba6445cf3139c3ea9b Mon Sep 17 00:00:00 2001 From: empEvil Date: Thu, 28 May 2026 13:49:35 +0100 Subject: [PATCH 15/78] added butler version to Anacal WIP --- txpipe/ingest/anacal.py | 114 +++++++++++++++++++++++++++++++++-- txpipe/ingest/dp1_details.py | 57 ++++++++++++++++++ 2 files changed, 167 insertions(+), 4 deletions(-) create mode 100644 txpipe/ingest/dp1_details.py diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 9deab1263..899593834 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -1,6 +1,7 @@ from .base import TXIngestCatalogFits from ..data_types import ShearCatalog, PhotometryCatalog, HDFFile, FileCollection, FitsFile from .lsst import process_photometry_data, process_shear_data +from .dp1_details import DP1_COSMOLOGY_FIELDS, DP1_TRACTS, DP1_COSMOLOGY_TRACTS, DP1_FIELD_CENTERS, DP1_SURVEY_PROPERTIES from ceci.config import StageParameter import numpy as np from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab @@ -8,7 +9,7 @@ class TXIngestAnacal(TXIngestCatalogFits): """ - Ingest an anacal catalog!, + Ingest an anacal catalog!, """ name = "TXIngestAnacal" @@ -19,6 +20,14 @@ class TXIngestAnacal(TXIngestCatalogFits): ("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."), + "cosmology_tracts_only": StageParameter(bool, True, msg="Use only cosmology tracts."), + "select_field": StageParameter(str, "", msg="Field to select (overrides cosmology_tracts_only)."), + "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(str, "grizy", msg="string of flux bands"), @@ -26,6 +35,93 @@ class TXIngestAnacal(TXIngestCatalogFits): } 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: + raise ImportError(error_msg) + + # 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"] + try: + butler = Butler(butler_config_file, collections=collections) + except: + raise RuntimeError(error_msg) + + if self.config["select_field"]: + tracts = DP1_TRACTS[self.config["select_field"]] + elif self.config["cosmology_tracts_only"]: + tracts = DP1_COSMOLOGY_TRACTS + else: + tracts = ALL_TRACTS + + n = self.get_catalog_size(butler, "ShearObject") + #shear_outfile = self.open_output("shear_catalog") + #group = shear_outfile.create_group("shear") + #shear_outfile["shear"].attrs["catalog_type"] = "Anacal" + + created_files = False + data_set_refs = butler.query_datasets('object_shear_all') + n_chunks = len(data_set_refs) + input_columns = self.get_input_columns() + + shear_start = 0 + for i, ref in enumerate(data_set_refs): + tract = ref.dataId["tract"] + if tract not in tracts: + print(f"Skipping chunk {i + 1} / {n_chunks} since tract {tract} is not selected") + continue + + d = butler.get('object_shear_all', + 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_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 {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): tracts = self.config["tracts"] n, dtypes = self.get_meta("anacal_catalog") @@ -47,7 +143,6 @@ def run(self): shear_outfile.close() - def setup_input(self): prefix = self.config["prefix"] cols = ( @@ -107,11 +202,11 @@ def process_anacal_shear_data(self, data): if band == "i": output["s2n"] = f / f_err - + for d in ["_dg1", "_dg2"]: dd = data[f"{band}_dflux_{s}"+d][:] output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) - + return output def setup_output(self, tag, group, first_chunk, n): @@ -130,3 +225,14 @@ def write_output(self, outfile, group, data, start, end): col = col.filled(np.nan) g[name][start:end] = col + 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 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", +} From fcc9c07a0ed4f188accceef50313eeb1192ad7b4 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:12:39 +0100 Subject: [PATCH 16/78] fixing butler calls --- txpipe/ingest/anacal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 899593834..ef40073a4 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -72,7 +72,7 @@ def butler_run(self): else: tracts = ALL_TRACTS - n = self.get_catalog_size(butler, "ShearObject") + n = self.get_catalog_size(butler, "deep_coadd_cell_anacal_merged") #shear_outfile = self.open_output("shear_catalog") #group = shear_outfile.create_group("shear") #shear_outfile["shear"].attrs["catalog_type"] = "Anacal" @@ -89,7 +89,7 @@ def butler_run(self): print(f"Skipping chunk {i + 1} / {n_chunks} since tract {tract} is not selected") continue - d = butler.get('object_shear_all', + d = butler.get("deep_coadd_cell_anacal_merged", #This name might change, we should double check with Xiangchong / PO dataId=ref.dataId, parameters={"columns": input_columns} ) From 75f3e4a7b1203ea6fff65aabe50b27c14018b833 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:24:14 +0100 Subject: [PATCH 17/78] fixin input columns --- txpipe/ingest/anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index ef40073a4..19fed9575 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -80,7 +80,7 @@ def butler_run(self): created_files = False data_set_refs = butler.query_datasets('object_shear_all') n_chunks = len(data_set_refs) - input_columns = self.get_input_columns() + input_columns = self.setup_input() shear_start = 0 for i, ref in enumerate(data_set_refs): From 71f004fcb225a0f349d4f15433160980632a3caf Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:30:17 +0100 Subject: [PATCH 18/78] fixing ingestion --- examples/Anacal_test/config.yml | 3 ++- txpipe/ingest/anacal.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index 434e91749..1a2560176 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -6,4 +6,5 @@ global: TXIngestAnacal: bands: griz - \ No newline at end of file + collections: "u/xiangchl/dp1/a360_anacal" + cosmology_tracts_only: False \ No newline at end of file diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 19fed9575..3aee2648d 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -1,7 +1,7 @@ from .base import TXIngestCatalogFits from ..data_types import ShearCatalog, PhotometryCatalog, HDFFile, FileCollection, FitsFile from .lsst import process_photometry_data, process_shear_data -from .dp1_details import DP1_COSMOLOGY_FIELDS, DP1_TRACTS, DP1_COSMOLOGY_TRACTS, DP1_FIELD_CENTERS, DP1_SURVEY_PROPERTIES +from .dp1_details import DP1_COSMOLOGY_FIELDS, DP1_TRACTS, DP1_COSMOLOGY_TRACTS, DP1_FIELD_CENTERS, DP1_SURVEY_PROPERTIES, ALL_TRACTS from ceci.config import StageParameter import numpy as np from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab From dcb80cc1d36b383c636de8fc21ddfb49f4473153 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:36:12 +0100 Subject: [PATCH 19/78] removing missing detail --- txpipe/ingest/anacal.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 3aee2648d..da817230e 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -150,7 +150,6 @@ def setup_input(self): "ra", "dec", "wsel", - "wdet", "mask_value", f"{prefix}_e1", f"{prefix}_e2", @@ -183,7 +182,6 @@ def process_anacal_shear_data(self, data): "dec": data["dec"][:], "weight": data["wsel"][:], "mask_value": data["mask_value"][:], - "weight_detection": data["wdet"][:], "weight_dg1": data["dwsel_dg1"][:], "weight_dg2": data["dwsel_dg2"][:], "e1": data[f"{prefix}_e1"][:], From 901b0c20ce9d293fbdf6e01af56d233877712863 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:38:13 +0100 Subject: [PATCH 20/78] fixing prefix --- examples/Anacal_test/config.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index 1a2560176..e171c3497 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -7,4 +7,5 @@ global: TXIngestAnacal: bands: griz collections: "u/xiangchl/dp1/a360_anacal" - cosmology_tracts_only: False \ No newline at end of file + cosmology_tracts_only: False + prefix: "fpsf1" \ No newline at end of file From b39c0d265be174d683bbb76d003bedb84b0ab9d2 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:39:38 +0100 Subject: [PATCH 21/78] typo --- examples/Anacal_test/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index e171c3497..c2f61fe3e 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -8,4 +8,4 @@ TXIngestAnacal: bands: griz collections: "u/xiangchl/dp1/a360_anacal" cosmology_tracts_only: False - prefix: "fpsf1" \ No newline at end of file + prefix: "fpfs1" \ No newline at end of file From cb0d00282df05950bccde6d38e110ddad600480d Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:41:09 +0100 Subject: [PATCH 22/78] fixing inputs --- txpipe/ingest/anacal.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index da817230e..c4908715e 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -153,13 +153,13 @@ def setup_input(self): "mask_value", f"{prefix}_e1", f"{prefix}_e2", - f"{prefix}_m0", - f"{prefix}_m2" + f"{prefix}_m00", + f"{prefix}_m20" ]) cols += ["dwsel"+ suffix for suffix in ["_dg1", "_dg2"]] cols += [ prefix +delta + suffix - for delta in ["_de1", "_de2", "_dm0", "_dm2"] + for delta in ["_de1", "_de2", "_dm00", "_dm20"] for suffix in ["_dg1", "_dg2"] ] bands = self.config["bands"] From c6f2a3bf3f9deb8f3ec5699aa4b9c17e2bb562c7 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:50:21 +0100 Subject: [PATCH 23/78] fixing flux reading --- examples/Anacal_test/config.yml | 3 ++- txpipe/ingest/anacal.py | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index c2f61fe3e..87104b493 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -8,4 +8,5 @@ TXIngestAnacal: bands: griz collections: "u/xiangchl/dp1/a360_anacal" cosmology_tracts_only: False - prefix: "fpfs1" \ No newline at end of file + prefix: "fpfs1" + scale: "fpfs1" \ No newline at end of file diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index c4908715e..ac9b090cd 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -145,6 +145,7 @@ def file_run(self): def setup_input(self): prefix = self.config["prefix"] + scale = self.config["scale"] cols = ( [ "ra", @@ -163,14 +164,18 @@ def setup_input(self): for suffix in ["_dg1", "_dg2"] ] bands = self.config["bands"] - for i in [0, 2, 4]: - cols += [band + "_flux_gauss" + f"{i}" for band in bands] - cols += [band + "_flux_gauss" + f"{i}" + "_err" for band in bands] - cols += [ - band + "_dflux_gauss" + f"{i}" + suffix - for band in bands - for suffix in ["_dg1", "_dg2"] - ] + cols += [band + "_flux_" + scale for band in bands] + cols += [band + "_flux_" + scale + "_err" for band in bands] + cols += [band + "_dflux" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] + + # for i in [0, 2, 4]: + # cols += [band + "_flux_" + scale + f"{i}" for band in bands] + # cols += [band + "_flux_" + scale + f"{i}" + "_err" for band in bands] + # cols += [ + # band + "_dflux_" + prefix + f"{i}" + suffix + # for band in bands + # for suffix in ["_dg1", "_dg2"] + # ] return cols def process_anacal_shear_data(self, data): From ba9fe191cb140c4f12d2fa32fb0eb1183b217a5f Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:51:38 +0100 Subject: [PATCH 24/78] typo --- txpipe/ingest/anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index ac9b090cd..53073fc0c 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -166,7 +166,7 @@ def setup_input(self): bands = self.config["bands"] cols += [band + "_flux_" + scale for band in bands] cols += [band + "_flux_" + scale + "_err" for band in bands] - cols += [band + "_dflux" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] + cols += [band + "_dflux_" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] # for i in [0, 2, 4]: # cols += [band + "_flux_" + scale + f"{i}" for band in bands] From cbce432608d883121c907b0bda87dddaefe5f95e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 16:56:20 +0100 Subject: [PATCH 25/78] another typo --- txpipe/ingest/anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 53073fc0c..fc8b04e81 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -99,7 +99,7 @@ def butler_run(self): print(f"Skipping chunk {i + 1} / {n_chunks} since it is empty") continue - shear_data = self.process_anacal_data(d) + shear_data = self.process_anacal_shear_data(d) if not created_files: created_files = True shear_outfile = self.setup_output("shear_catalog", "shear", From 276d8054221adc19e9328a6ab67b327a0fa9d791 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Thu, 28 May 2026 17:09:35 +0100 Subject: [PATCH 26/78] fixing sourceselection --- txpipe/source_selection/Anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index d02f0a631..52097c959 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -63,7 +63,7 @@ def setup_output(self): Prepare the output columns for the response values generated bby Anacal """ outfile = super().setup_output() - n = outfile["count/counts"].size + n = outfile["tomography/bin"].size group = outfile.create_group("response") group.create_dataset("R", (n, 1, 1), dtype="f") return outfile From 99eec0b9933ba5d03a7541223863e34bd403b7b4 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 14:30:53 +0100 Subject: [PATCH 27/78] fixing the calculator --- txpipe/shear_calibration/calibration_calculators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 025c7aa56..ddf5db4cb 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -900,7 +900,7 @@ def add_data(self, data, *args, **kwargs): sel: array The indicies of the objects selected from this chunk of data """ - + data = _DataWrapper(data, "") select = self.selector(data, *args, **kwargs) e1 = data["e1"] From c5744a75781b1302a3c6f3e714edafbac95da785 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:41:36 +0100 Subject: [PATCH 28/78] adding specific selectors --- txpipe/shear_calibration/calibrators.py | 2 ++ txpipe/source_selection/Anacal.py | 38 ++++++++++++++++++++++--- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index 1b280f9bf..d3819dacd 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") diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index 52097c959..1c16b107e 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -46,7 +46,8 @@ def data_iterator(self): "dm00_dg1", "dm00_dg2", "dm20_dg1", - "dm20_dg2" + "dm20_dg2", + "s2n" ] shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="Anacal") @@ -70,12 +71,41 @@ def setup_output(self): def setup_response_calculators(self, nbin_source): delta_gamma = self.config["delta_gamma"] - calculators = [AnaCalCalculator(select_tomographic_weak_lensing_sample, delta_gamma) for i in range(nbin_source)] - calculators.append(AnaCalCalculator(select_weak_lensing_sample, 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_tomogrpahy(outfile, start, end, source_bin, R) group = outfile["response"] group["R"][start:end] = R - \ No newline at end of file + +def select_anacal_weak_lensing_sample(data, config, calling_from_select=False): + s2n_cut = config["s2n_cut"] + verbose = config["verbose"] + + flag = data["mask_value"] + s2n = data["s2n"] + + n0 = len(flag) + sel = flag == 0 + f1 = sel.sum() / n0 + + sel &= s2n > s2n_cut + f2 = sel.sum() / n0 + + sel &= data["zbin"] >= 0 + f3 = sel.sum() / n0 + + if verbose and calling_from_select: + print(f"Tomo selection {f1:.2%} flag, {f2:.2%} SNR, ", end ="") + elif verbose: + print(f"2D selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.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 From a795deacdadfcf463855f2a6866e452a86dd30df Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:44:59 +0100 Subject: [PATCH 29/78] adding missing columns --- txpipe/source_selection/Anacal.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index 1c16b107e..e21fdf535 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -47,7 +47,9 @@ def data_iterator(self): "dm00_dg2", "dm20_dg1", "dm20_dg2", - "s2n" + "s2n", + "weight_dg1", + "weight_dg2" ] shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="Anacal") From 30b4e45b3761ecf0067282eb6d8a26f848db7d76 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:50:48 +0100 Subject: [PATCH 30/78] fixing the calibration --- txpipe/shear_calibration/calibration_calculators.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index ddf5db4cb..b96c3f0f2 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -929,6 +929,8 @@ def add_data(self, data, *args, **kwargs): # Next we calculate the cuts per variation masks wsel = weight[select] + e1_sel = e1[select] + e2_sel = e2[select] mask = mask_value[select] < 40 mask_p1 = mask & self.get_submask(m00, m20, dm00_dg1, dm20_dg1, +1) @@ -936,10 +938,10 @@ def add_data(self, data, *args, **kwargs): mask_p2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, +1) mask_m2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, -1) - self.sel_response_e1p1 += np.sum(weight[mask_p1]*e1[mask_p1]) - self.sel_response_e1m1 += np.sum(weight[mask_m1]*e1[mask_m1]) - self.sel_response_e2p2 += np.sum(weight[mask_p2]*e2[mask_p2]) - self.sel_response_e2m2 += np.sum(weight[mask_m2]*e2[mask_m2]) + self.sel_response_e1p1 += np.sum(weight[mask_p1] * e1_sel[mask_p1]) + self.sel_response_e1m1 += np.sum(weight[mask_m1] * e1_sel[mask_m1]) + self.sel_response_e2p2 += np.sum(weight[mask_p2] * e2_sel[mask_p2]) + self.sel_response_e2m2 += np.sum(weight[mask_m2] * e2_sel[mask_m2]) # Next we find the means part needed for the Shape response, # and the Weight-bias response From 1f29f6aaeae7ae88e10bec35f911197170e4b232 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:52:19 +0100 Subject: [PATCH 31/78] more fixes --- txpipe/shear_calibration/calibration_calculators.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index b96c3f0f2..1cf770e21 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -933,10 +933,10 @@ def add_data(self, data, *args, **kwargs): e2_sel = e2[select] mask = mask_value[select] < 40 - mask_p1 = mask & self.get_submask(m00, m20, dm00_dg1, dm20_dg1, +1) - mask_m1 = mask & self.get_submask(m00, m20, dm00_dg1, dm20_dg1, -1) - mask_p2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, +1) - mask_m2 = mask & self.get_submask(m00, m20, dm00_dg2, dm20_dg2, -1) + mask_p1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], +1) + mask_m1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], -1) + mask_p2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], +1) + mask_m2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], -1) self.sel_response_e1p1 += np.sum(weight[mask_p1] * e1_sel[mask_p1]) self.sel_response_e1m1 += np.sum(weight[mask_m1] * e1_sel[mask_m1]) From 9ad3c5e288cf58dc8504b4a13315982b29df1c6a Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:55:07 +0100 Subject: [PATCH 32/78] more fixes --- txpipe/shear_calibration/calibration_calculators.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 1cf770e21..1aa858f52 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -938,10 +938,10 @@ def add_data(self, data, *args, **kwargs): mask_p2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], +1) mask_m2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], -1) - self.sel_response_e1p1 += np.sum(weight[mask_p1] * e1_sel[mask_p1]) - self.sel_response_e1m1 += np.sum(weight[mask_m1] * e1_sel[mask_m1]) - self.sel_response_e2p2 += np.sum(weight[mask_p2] * e2_sel[mask_p2]) - self.sel_response_e2m2 += np.sum(weight[mask_m2] * e2_sel[mask_m2]) + self.sel_response_e1p1 += np.sum(wsel[mask_p1] * e1_sel[mask_p1]) + self.sel_response_e1m1 += np.sum(wsel[mask_m1] * e1_sel[mask_m1]) + self.sel_response_e2p2 += np.sum(wsel[mask_p2] * e2_sel[mask_p2]) + self.sel_response_e2m2 += np.sum(wsel[mask_m2] * e2_sel[mask_m2]) # Next we find the means part needed for the Shape response, # and the Weight-bias response From 2459f909fb42d2968994ced1b8d490159a6709a9 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 15:57:32 +0100 Subject: [PATCH 33/78] typo --- txpipe/source_selection/Anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index e21fdf535..78ccba308 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -78,7 +78,7 @@ def setup_response_calculators(self, nbin_source): return calculators def write_tomography(self, outfile, start, end, source_bin, R): - super().write_tomogrpahy(outfile, start, end, source_bin, R) + super().write_tomography(outfile, start, end, source_bin, R) group = outfile["response"] group["R"][start:end] = R From 7bbc45092d3d92ee2b085cfa2a9d8978719ed2d4 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:00:13 +0100 Subject: [PATCH 34/78] fixing a missing import --- txpipe/shear_calibration/calibration_calculators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 1aa858f52..8aad96d4c 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: From f727a960dc0cca6667d4e02fe51cae613d2a324e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:02:16 +0100 Subject: [PATCH 35/78] fixed a typo --- txpipe/shear_calibration/calibration_calculators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 8aad96d4c..a3c38fbe9 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -1037,7 +1037,7 @@ def collect(self, comm=None, allgather=False) -> BinStats: else: Neff = sum_weights**2 / sum_sq_weights - calibrator = AnaCalibrator(R_total, mean_e, w00, mu_is_weighted=False) + 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)) bin_stats = BinStats(count, Neff, calibrator.mu, sigma_e, sigma, calibrator) From 06e1820d2425e0f01a5f80c4ed336021df3d989f Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:03:43 +0100 Subject: [PATCH 36/78] fixed a typo --- txpipe/shear_calibration/calibrators.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index d3819dacd..f8bbf35dd 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -635,9 +635,9 @@ def calibrate_sigma(self, sigma): class AnaCalibrator(MetaCalibrator): """Stores information needed to calibrate a Anacal shear method""" - def __init__(self, R, mu, mu_is_weigthed=True): + def __init__(self, R, mu, mu_is_weighted=True): self.R = R - if mu_is_weigthed: + if mu_is_weighted: self.mu = np.array(mu) else: assert("Anacal needs an already calibrated mu.") From 05ac192680d069fa9a0155e707ab471c7c8d56af Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:07:49 +0100 Subject: [PATCH 37/78] fixing the calibrator --- txpipe/shear_calibration/calibrators.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index f8bbf35dd..a37a50adf 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -640,6 +640,7 @@ def __init__(self, R, mu, mu_is_weighted=True): 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, substract_mean=True): @@ -665,6 +666,12 @@ def apply(self, g1, g2, weights, substract_mean=True): 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): """ From 089c8bb4b77b3f3f98ebb76979e75d98113ceca5 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:10:36 +0100 Subject: [PATCH 38/78] typo --- txpipe/shear_calibration/calibrators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index a37a50adf..3c62bd3d0 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -703,7 +703,7 @@ def load(cls, tomo_file): return calibrators def save(self, outfile, i): - outfile["respponse/R"][i] = self.R + outfile["response/R"][i] = self.R outfile["counts/mean_e1"][i] = self.mu[0] outfile["counts/mean_e2"][i] = self.mu[1] From 64e95ec2c0e758f254eb64270fe3878a0afe00f2 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:16:36 +0100 Subject: [PATCH 39/78] changed save --- txpipe/shear_calibration/calibrators.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index 3c62bd3d0..d269cffd2 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -703,7 +703,12 @@ def load(cls, tomo_file): return calibrators def save(self, outfile, i): - outfile["response/R"][i] = self.R - outfile["counts/mean_e1"][i] = self.mu[0] - outfile["counts/mean_e2"][i] = self.mu[1] + 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] From 803cee298513b7f6458070ec5f93ede6b088ffec Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 29 May 2026 16:27:06 +0100 Subject: [PATCH 40/78] adding 2d --- txpipe/source_selection/Anacal.py | 1 + 1 file changed, 1 insertion(+) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index 78ccba308..bc8917474 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -69,6 +69,7 @@ def setup_output(self): 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): From e761385abe9d143e52d5c6b6891b57cea183b4d9 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Mon, 1 Jun 2026 16:30:26 +0100 Subject: [PATCH 41/78] trying to fix calibration --- txpipe/shear_calibration/calibrators.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index d269cffd2..b623e942f 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -693,14 +693,17 @@ def load(cls, tomo_file): import h5py with h5py.File(tomo_file, "r") as f: - R = f["response/R"][:] - n = len(R) - 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)] - return calibrators + calibrator2d = cls(R_2d, [mu1_2d, mu2_2d]) + return calibrators, calibrator2d def save(self, outfile, i): if i == "2d": From 62ad6a606ab18085da217352029a8fe9de873ccb Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Mon, 1 Jun 2026 16:41:32 +0100 Subject: [PATCH 42/78] adding anacal definitions --- txpipe/data_types.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/txpipe/data_types.py b/txpipe/data_types.py index 3326e8d2d..8ed93ca9b 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 = {} From 82413b6fdff383fb5eb5d5415bcc685534845874 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Mon, 1 Jun 2026 16:47:16 +0100 Subject: [PATCH 43/78] fixing typos --- txpipe/shear_calibration/calibrators.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index b623e942f..fcbf92228 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -643,7 +643,7 @@ def __init__(self, R, mu, mu_is_weighted=True): self.mu = np.array(mu) / R assert("Anacal needs an already calibrated mu.") - def apply(self, g1, g2, weights, substract_mean=True): + def apply(self, g1, g2, weights, subtract_mean=True): """ Calibrate a set of shears using the response matrix and mean shear substraction @@ -658,7 +658,7 @@ def apply(self, g1, g2, weights, substract_mean=True): subtract_mean: bool whether to subtract mean shear (default True) """ - if not substract_mean: + if not subtract_mean: g1, g2 = weights * [g1, g2] / self.R elif np.isscalar(g1): g1, g2 = weights * [g1, g2] / self.R - self.mu From 20363b3393ff0c75fd128943683f8d8e8ed3e25e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Mon, 1 Jun 2026 16:51:07 +0100 Subject: [PATCH 44/78] adding anacal specifics --- txpipe/calibrate.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/txpipe/calibrate.py b/txpipe/calibrate.py index 37342ce09..f7836ed3e 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) From 8a900fe4d6bad5a498302a94818a7bb5bc4282e5 Mon Sep 17 00:00:00 2001 From: empEvil Date: Wed, 1 Jul 2026 14:46:05 +0100 Subject: [PATCH 45/78] adding tests for anacal's calibrator --- txpipe/test/test_cal.py | 61 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index c721067f8..ca46cb613 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,60 @@ def core_metadet(comm): assert stats.source_count == N * nproc +def core_anacal(comm): + delta_gamma = .02 + nproc = 1 if comm is None else comm.size + N = 10 + + + R_shape_true = 0.4 # known test value + R_weight_true = 0.3 + + base_data = { + "e1": np.random.normal(0, 0.1, size=N), + "e2": np.random.normal(0, 0.1, size=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), + "mask_value": 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) + assert np.allclose(stats.calibrator.R, R_weight_true) + + # 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) + assert np.allclose(stats.calibrator.R, R_shape_true + R_weight_true) + + def test_metacalibrator_serial(): core_metacal(None) @@ -168,6 +225,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 +436,4 @@ def test_null_mean(): test_metadetect_parallel() test_mean_shear_no_weights() test_mean_shear_weights() + test_anacal_serial() From 1ec31296ba8064621bac89f37cbb506ffa00dcd9 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Fri, 10 Jul 2026 13:53:55 +0100 Subject: [PATCH 46/78] small upates --- examples/Anacal_test/config.yml | 2 +- txpipe/ingest/anacal.py | 20 +++++++++++--------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/Anacal_test/config.yml b/examples/Anacal_test/config.yml index 87104b493..764be5626 100644 --- a/examples/Anacal_test/config.yml +++ b/examples/Anacal_test/config.yml @@ -6,7 +6,7 @@ global: TXIngestAnacal: bands: griz - collections: "u/xiangchl/dp1/a360_anacal" + collections: "u/xiangchl/dp1/a360_anacal2" cosmology_tracts_only: False prefix: "fpfs1" scale: "fpfs1" \ No newline at end of file diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index fc8b04e81..1982bfbbf 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -25,6 +25,7 @@ class TXIngestAnacal(TXIngestCatalogFits): "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)."), "collections": StageParameter(str, "LSSTComCam/DP1", msg="Butler collections to use."), @@ -52,8 +53,8 @@ def butler_run(self): ) try: from lsst.daf.butler import Butler - except: - raise ImportError(error_msg) + 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. @@ -62,8 +63,8 @@ def butler_run(self): collections = self.config["collections"] try: butler = Butler(butler_config_file, collections=collections) - except: - raise RuntimeError(error_msg) + except Exception as e: + raise RuntimeError(error_msg) from e if self.config["select_field"]: tracts = DP1_TRACTS[self.config["select_field"]] @@ -72,13 +73,14 @@ def butler_run(self): else: tracts = ALL_TRACTS - n = self.get_catalog_size(butler, "deep_coadd_cell_anacal_merged") + object_name = self.config["butler_object_name"] + n = self.get_catalog_size(butler, object_name) #shear_outfile = self.open_output("shear_catalog") #group = shear_outfile.create_group("shear") #shear_outfile["shear"].attrs["catalog_type"] = "Anacal" created_files = False - data_set_refs = butler.query_datasets('object_shear_all') + data_set_refs = butler.query_datasets(object_name) n_chunks = len(data_set_refs) input_columns = self.setup_input() @@ -89,7 +91,7 @@ def butler_run(self): print(f"Skipping chunk {i + 1} / {n_chunks} since tract {tract} is not selected") continue - d = butler.get("deep_coadd_cell_anacal_merged", #This name might change, we should double check with Xiangchong / PO + d = butler.get(object_name, dataId=ref.dataId, parameters={"columns": input_columns} ) @@ -206,8 +208,8 @@ def process_anacal_shear_data(self, data): if band == "i": output["s2n"] = f / f_err - for d in ["_dg1", "_dg2"]: - dd = data[f"{band}_dflux_{s}"+d][:] + for d in ["dg1", "dg2"]: + dd = data[f"{band}_dflux_{s}_"+d][:] output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) return output From 20f554ed2683ebd27a3b88646912b46a2f5ff984 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 15:25:38 +0100 Subject: [PATCH 47/78] adding s2n derivatives --- txpipe/ingest/anacal.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 1982bfbbf..e7cba764a 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -7,6 +7,7 @@ from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab from ..utils.hdf_tools import h5py_shorten, repack + class TXIngestAnacal(TXIngestCatalogFits): """ Ingest an anacal catalog!, @@ -211,6 +212,9 @@ def process_anacal_shear_data(self, data): for d in ["dg1", "dg2"]: dd = data[f"{band}_dflux_{s}_"+d][:] output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) + if band == "i": + output[f"ds2n_{d}"] = dd/f_err + return output From baaadffa9e401eb983b7b5f888b10689672c9e2e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 15:30:14 +0100 Subject: [PATCH 48/78] fixing response flux to magnitude --- txpipe/ingest/anacal.py | 6 +++--- txpipe/utils/__init__.py | 2 +- txpipe/utils/conversion.py | 9 +++++++++ 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index e7cba764a..40922f0f1 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -4,7 +4,7 @@ from .dp1_details import DP1_COSMOLOGY_FIELDS, DP1_TRACTS, DP1_COSMOLOGY_TRACTS, DP1_FIELD_CENTERS, DP1_SURVEY_PROPERTIES, ALL_TRACTS from ceci.config import StageParameter import numpy as np -from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab +from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, anacal_mag_response from ..utils.hdf_tools import h5py_shorten, repack @@ -211,10 +211,10 @@ def process_anacal_shear_data(self, data): for d in ["dg1", "dg2"]: dd = data[f"{band}_dflux_{s}_"+d][:] - output[f"mag_{band}_{d}"] = nanojansky_to_mag_ab(dd) + output[f"mag_{band}_{d}"] = anacal_mag_response(f, dd) if band == "i": output[f"ds2n_{d}"] = dd/f_err - + return output diff --git a/txpipe/utils/__init__.py b/txpipe/utils/__init__.py index c6c3e3925..56cc25bbd 100755 --- a/txpipe/utils/__init__.py +++ b/txpipe/utils/__init__.py @@ -3,7 +3,7 @@ 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, anacal_mag_response 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..53cbfc49a 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -21,3 +21,12 @@ def moments_to_shear(Ixx, Iyy, Ixy): e1 = (Ixx - Iyy) / b e2 = 2 * Ixy / b return e1, e2 + + +def anacal_mag_response(flux, response): + mag_resp = np.divide( -2.5 / np.log(10) * response, + flux, + out=np.full_like(f, np.nan), + where=f > 0 + ) + return mag_resp \ No newline at end of file From cd1a7fe7645c5a5557a0917934e5c384512dc48a Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 16:43:45 +0100 Subject: [PATCH 49/78] fixing mask_threshold --- txpipe/shear_calibration/calibration_calculators.py | 7 +++++-- txpipe/source_selection/Anacal.py | 10 +++++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index a3c38fbe9..d4e0564d7 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -855,7 +855,7 @@ class AnaCalCalculator(CalibrationCalculator): See the CalibrationCalculator class for the use and contents of this class """ - def __init__(self, selector, delta_gamma): + def __init__(self, selector, delta_gamma, mask_threshold): """ Initialize the Calibrator using the funtion you will use to select objects. That function should take at least one argument, @@ -871,11 +871,14 @@ def __init__(self, selector, delta_gamma): Function that selects objects delta_gamma: float The difference in applied g between 1p and 1m variants + mask_threshold: float + The value of the mask at which to cut """ from parallel_statistics import ParallelMean super().__init__(selector) self.delta_gamma = delta_gamma + self.mask_threshold = mask_threshold self.response_means = ParallelMean(size=4) self.sel_response_e1p1 = 0 self.sel_response_e1m1 = 0 @@ -931,7 +934,7 @@ def add_data(self, data, *args, **kwargs): wsel = weight[select] e1_sel = e1[select] e2_sel = e2[select] - mask = mask_value[select] < 40 + mask = mask_value[select] < self.mask_threshold mask_p1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], +1) mask_m1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], -1) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index bc8917474..cc748c31f 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -22,6 +22,7 @@ class TXSourceSelectorAnacal(TXSourceSelectorBase): required=True, msg= "Delta gmamma value for hte AnaCal response calculation" ), + "mask_threshold": StageParameter(float, 40, msg = "mask threshold, for when to mask objects." ), } def data_iterator(self): @@ -48,6 +49,8 @@ def data_iterator(self): "dm20_dg1", "dm20_dg2", "s2n", + "ds2n_dg1", + "ds2n_dg2", "weight_dg1", "weight_dg2" ] @@ -74,8 +77,9 @@ def setup_output(self): 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)) + mask_threshold = self.config["mask_threshold"] + calculators = [AnaCalCalculator(select_anacal_tomographic_weak_lensing_sample, delta_gamma, mask_threshold) for i in range(nbin_source)] + calculators.append(AnaCalCalculator(select_anacal_weak_lensing_sample, delta_gamma, mask_threshold)) return calculators def write_tomography(self, outfile, start, end, source_bin, R): @@ -87,7 +91,7 @@ def select_anacal_weak_lensing_sample(data, config, calling_from_select=False): s2n_cut = config["s2n_cut"] verbose = config["verbose"] - flag = data["mask_value"] + flag = data["mask_value"] < config["mask_threshold"] s2n = data["s2n"] n0 = len(flag) From d80dd630dd00b537a26c89406bc98db4909894a7 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 16:44:41 +0100 Subject: [PATCH 50/78] fixing folder --- examples/{Anacal_test => anacal}/config.yml | 0 examples/{Anacal_test => anacal}/pipeline.yml | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename examples/{Anacal_test => anacal}/config.yml (100%) rename examples/{Anacal_test => anacal}/pipeline.yml (100%) diff --git a/examples/Anacal_test/config.yml b/examples/anacal/config.yml similarity index 100% rename from examples/Anacal_test/config.yml rename to examples/anacal/config.yml diff --git a/examples/Anacal_test/pipeline.yml b/examples/anacal/pipeline.yml similarity index 100% rename from examples/Anacal_test/pipeline.yml rename to examples/anacal/pipeline.yml From 5a66b750b4b7799ffa813234d53e4ae8322f4f82 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 16:51:36 +0100 Subject: [PATCH 51/78] fixing catalog name --- txpipe/calibrate.py | 2 +- txpipe/ingest/anacal.py | 17 ++--------------- txpipe/shear_calibration/calibrators.py | 2 +- 3 files changed, 4 insertions(+), 17 deletions(-) diff --git a/txpipe/calibrate.py b/txpipe/calibrate.py index f7836ed3e..a5227826c 100644 --- a/txpipe/calibrate.py +++ b/txpipe/calibrate.py @@ -181,7 +181,7 @@ 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": + 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/ingest/anacal.py b/txpipe/ingest/anacal.py index 40922f0f1..3bb76d75f 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -76,9 +76,6 @@ def butler_run(self): object_name = self.config["butler_object_name"] n = self.get_catalog_size(butler, object_name) - #shear_outfile = self.open_output("shear_catalog") - #group = shear_outfile.create_group("shear") - #shear_outfile["shear"].attrs["catalog_type"] = "Anacal" created_files = False data_set_refs = butler.query_datasets(object_name) @@ -108,7 +105,7 @@ def butler_run(self): shear_outfile = self.setup_output("shear_catalog", "shear", shear_data, n) - shear_outfile["shear"].attrs["catalog_type"] = "Anacal" + 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, @@ -136,8 +133,7 @@ def file_run(self): 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" - #self.write_output(shear_outfile, "shear", shear_data) + shear_outfile["shear"].attrs["catalog_type"] = "anacal" print("Trimming shear columns:") for col in shear_data.keys(): @@ -171,14 +167,6 @@ def setup_input(self): cols += [band + "_flux_" + scale + "_err" for band in bands] cols += [band + "_dflux_" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] - # for i in [0, 2, 4]: - # cols += [band + "_flux_" + scale + f"{i}" for band in bands] - # cols += [band + "_flux_" + scale + f"{i}" + "_err" for band in bands] - # cols += [ - # band + "_dflux_" + prefix + f"{i}" + suffix - # for band in bands - # for suffix in ["_dg1", "_dg2"] - # ] return cols def process_anacal_shear_data(self, data): @@ -215,7 +203,6 @@ def process_anacal_shear_data(self, data): if band == "i": output[f"ds2n_{d}"] = dd/f_err - return output def setup_output(self, tag, group, first_chunk, n): diff --git a/txpipe/shear_calibration/calibrators.py b/txpipe/shear_calibration/calibrators.py index fcbf92228..eab644082 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -58,7 +58,7 @@ def load(cls, tomo_file, null=False): subcls = LensfitCalibrator elif cat_type == "hsc": subcls = HSCCalibrator - elif cat_type == "Anacal": + elif cat_type == "anacal": subcls = AnaCalibrator else: raise ValueError(f"Unknown catalog type {cat_type} in tomo file") From 69f08ef40d4138ff4b81a6953dfbeece92189908 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 16:55:04 +0100 Subject: [PATCH 52/78] minor fixes --- txpipe/ingest/anacal.py | 2 +- txpipe/source_selection/Anacal.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 3bb76d75f..ee4c4b9ca 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -32,7 +32,7 @@ class TXIngestAnacal(TXIngestCatalogFits): "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(str, "grizy", msg="string of flux bands"), + "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") } diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index cc748c31f..ccd035ce6 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -20,7 +20,7 @@ class TXSourceSelectorAnacal(TXSourceSelectorBase): "delta_gamma": StageParameter( float, required=True, - msg= "Delta gmamma value for hte AnaCal response calculation" + msg= "Delta gamma value for the AnaCal response calculation" ), "mask_threshold": StageParameter(float, 40, msg = "mask threshold, for when to mask objects." ), } @@ -66,7 +66,7 @@ def data_iterator(self): def setup_output(self): """ - Prepare the output columns for the response values generated bby Anacal + Prepare the output columns for the response values generated by Anacal """ outfile = super().setup_output() n = outfile["tomography/bin"].size From a0275a97ca8415fb530e3dc87601f8e1f7c69dc4 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 16:56:43 +0100 Subject: [PATCH 53/78] fixing doctstring --- txpipe/ingest/anacal.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index ee4c4b9ca..ca821938f 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -10,7 +10,9 @@ class TXIngestAnacal(TXIngestCatalogFits): """ - Ingest an anacal catalog!, + 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" From 89b6543b2a4b3e9165d1f31250c651777f12c03e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:00:23 +0100 Subject: [PATCH 54/78] fixing error message --- txpipe/ingest/anacal.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index ca821938f..5f92a40c4 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -64,10 +64,13 @@ def butler_run(self): 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_msg) from e + raise RuntimeError(error_msg2) from e if self.config["select_field"]: tracts = DP1_TRACTS[self.config["select_field"]] From 3a72a90670ed47794c964e45ecb68e9ef9fffee2 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:03:47 +0100 Subject: [PATCH 55/78] adding select_tracts option --- txpipe/ingest/anacal.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 5f92a40c4..949b46f89 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -31,6 +31,7 @@ class TXIngestAnacal(TXIngestCatalogFits): "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 "), @@ -74,6 +75,8 @@ def butler_run(self): 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: From 7346ae7e43ddcec51944bbaf93d7065f8c822e39 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:12:30 +0100 Subject: [PATCH 56/78] trying to fix the calibration of the detection response --- .../calibration_calculators.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index d4e0564d7..2b0cd293e 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -879,11 +879,13 @@ def __init__(self, selector, delta_gamma, mask_threshold): self.delta_gamma = delta_gamma self.mask_threshold = mask_threshold - self.response_means = ParallelMean(size=4) + self.response_means = ParallelMean(size=2) self.sel_response_e1p1 = 0 self.sel_response_e1m1 = 0 self.sel_response_e2p2 = 0 self.sel_response_e2m2 = 0 + self.det_response_e1 = 0 + self.det_response_e2 = 0 #self.sel_bias_means = ParallelMean(size=8) def add_data(self, data, *args, **kwargs): @@ -955,8 +957,9 @@ def add_data(self, data, *args, **kwargs): self.response_means.add_data(0, de1_dg1_sub, wsel) self.response_means.add_data(1, de2_dg2_sub, wsel) - self.response_means.add_data(2, dwsel_dg1, wsel) - self.response_means.add_data(3, dwsel_dg2, wsel) + + self.det_response_e1 += np.sum(dwsel_dg1 * e1_sel) + self.det_response_e2 += np.sum(dwsel_dg2 * e2_sel) self.shear_stats.add_data(0, e1[select], wsel) self.shear_stats.add_data(1, e2[select], wsel) @@ -997,6 +1000,8 @@ def collect(self, comm=None, allgather=False) -> BinStats: sum_sel_response_e1m1 = comm.allreduce(self.sel_response_e1m1) sum_sel_response_e2p2 = comm.allreduce(self.sel_response_e2p2) sum_sel_response_e2m2 = comm.allreduce(self.sel_response_e2m2) + sum_det_response_e1 = comm.allreduce(self.det_response_e1) + sum_det_response_e2 = comm.allreduce(self.det_response_e2) else: count = comm.reduce(self.count) sum_weights = comm.reduce(self.sum_weights) @@ -1005,6 +1010,8 @@ def collect(self, comm=None, allgather=False) -> BinStats: sum_sel_response_e1m1 = comm.reduce(self.sel_response_e1m1) sum_sel_response_e2p2 = comm.reduce(self.sel_response_e2p2) sum_sel_response_e2m2 = comm.reduce(self.sel_response_e2m2) + sum_det_response_e1 = comm.reduce(self.det_respponse_e1) + sum_det_response_e2 = comm.reduce(self.det_response_e2) else: count = self.count sum_weights = self.sum_weights @@ -1013,6 +1020,8 @@ def collect(self, comm=None, allgather=False) -> BinStats: sum_sel_response_e1m1 = self.sel_response_e1m1 sum_sel_response_e2p2 = self.sel_response_e2p2 sum_sel_response_e2m2 = self.sel_response_e2m2 + sum_det_response_e1 = self.det_response_e1 + sum_det_response_e2 = self.det_response_e2 # Collect the mean values we need mode = "allgather" if allgather else "gather" @@ -1027,7 +1036,7 @@ def collect(self, comm=None, allgather=False) -> BinStats: # Reducing down the responses R_shape = 0.5 * (mean_de1_dg1 + mean_de2_dg2) - R_weight = 0.5 * (mean_dwsel_dg1_e1 + mean_dwsel_dg2_e2) + R_weight = 0.5 * (sum_det_response_e1 + sum_det_response_e2) / sum_weights R_sel_1 = (sum_sel_response_e1p1 - sum_sel_response_e1m1) / (2.0 * self.delta_gamma) / count R_sel_2 = (sum_sel_response_e2p2 - sum_sel_response_e2m2) / (2.0 * self.delta_gamma) / count From 51317983b48b533bf30cf296bdeb9ddd53387120 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:13:16 +0100 Subject: [PATCH 57/78] removing unused things --- txpipe/shear_calibration/calibration_calculators.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 2b0cd293e..234a95cdf 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -1031,8 +1031,6 @@ def collect(self, comm=None, allgather=False) -> BinStats: # Unpack the flat mean R: mean_de1_dg1 = R[0] mean_de2_dg2 = R[1] - mean_dwsel_dg1_e1 = R[2] - mean_dwsel_dg2_e2 = R[3] # Reducing down the responses R_shape = 0.5 * (mean_de1_dg1 + mean_de2_dg2) From 0e5890a1388c89da4f27f69e392da270f25aef7a Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:20:02 +0100 Subject: [PATCH 58/78] fixing threshold flagging --- txpipe/source_selection/Anacal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index ccd035ce6..89aeeb53f 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -95,7 +95,7 @@ def select_anacal_weak_lensing_sample(data, config, calling_from_select=False): s2n = data["s2n"] n0 = len(flag) - sel = flag == 0 + sel = flag f1 = sel.sum() / n0 sel &= s2n > s2n_cut From 32d49796f3c7cbb6ad1ebe305dc8d56b7860950d Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:26:21 +0100 Subject: [PATCH 59/78] fixing missing input in test --- txpipe/test/test_cal.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index ca46cb613..c4d905636 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -166,6 +166,7 @@ def core_metadet(comm): def core_anacal(comm): delta_gamma = .02 + mask_threshold = 40 nproc = 1 if comm is None else comm.size N = 10 @@ -192,7 +193,7 @@ def core_anacal(comm): # 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 = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) cal.add_data(data) stats = cal.collect(comm, allgather=True) assert np.allclose(stats.calibrator.R, R_shape_true) @@ -200,7 +201,7 @@ def core_anacal(comm): # 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 = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) cal.add_data(data) stats = cal.collect(comm, allgather=True) assert np.allclose(stats.calibrator.R, R_weight_true) @@ -212,7 +213,7 @@ def core_anacal(comm): "weight_dg1": np.full(N, R_weight_true), "weight_dg2": np.full(N, R_weight_true), } - cal = AnaCalCalculator(select_all_anacal, delta_gamma) + cal = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) cal.add_data(data) stats = cal.collect(comm, allgather=True) assert np.allclose(stats.calibrator.R, R_shape_true + R_weight_true) From 8d3a0be5aac31d7894c37efefec94bbaf05ec6a0 Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Tue, 14 Jul 2026 17:44:50 +0100 Subject: [PATCH 60/78] correcting test --- txpipe/test/test_cal.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index c4d905636..1fbe0ba0d 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -204,7 +204,9 @@ def core_anacal(comm): cal = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) cal.add_data(data) stats = cal.collect(comm, allgather=True) - assert np.allclose(stats.calibrator.R, R_weight_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, @@ -216,7 +218,9 @@ def core_anacal(comm): cal = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) cal.add_data(data) stats = cal.collect(comm, allgather=True) - assert np.allclose(stats.calibrator.R, R_shape_true + R_weight_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(): From e98a09a1a5a8232fae269d41ac0051fa42aa2a34 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Tue, 14 Jul 2026 22:48:32 -0700 Subject: [PATCH 61/78] 1. unify shear_catalog_type to "anacal". There was 'Anacal' and 'anacal' 2. correct for the typo in calibration_calibrators.py, det_respponse_e1 3. in conversion.py f was not defined, changed to flux --- txpipe/data_types.py | 2 +- txpipe/shear_calibration/calibration_calculators.py | 2 +- txpipe/source_selection/Anacal.py | 2 +- txpipe/utils/conversion.py | 8 ++++---- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/txpipe/data_types.py b/txpipe/data_types.py index 8ed93ca9b..b5a65e1b0 100755 --- a/txpipe/data_types.py +++ b/txpipe/data_types.py @@ -108,7 +108,7 @@ 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": + elif self.catalog_type == "anacal": shear_cols = ["e1", "e2", "ra", "dec", "weight"] rename = {"e1":"g1", "e2":"g2"} else: diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 234a95cdf..4eb5a6c30 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -1010,7 +1010,7 @@ def collect(self, comm=None, allgather=False) -> BinStats: sum_sel_response_e1m1 = comm.reduce(self.sel_response_e1m1) sum_sel_response_e2p2 = comm.reduce(self.sel_response_e2p2) sum_sel_response_e2m2 = comm.reduce(self.sel_response_e2m2) - sum_det_response_e1 = comm.reduce(self.det_respponse_e1) + sum_det_response_e1 = comm.reduce(self.det_response_e1) sum_det_response_e2 = comm.reduce(self.det_response_e2) else: count = self.count diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index 89aeeb53f..bdcc914a6 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -54,7 +54,7 @@ def data_iterator(self): "weight_dg1", "weight_dg2" ] - shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="Anacal") + shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="anacal") if self.config["input_pz"]: shear_cols += ["mean_z"] diff --git a/txpipe/utils/conversion.py b/txpipe/utils/conversion.py index 53cbfc49a..2c9d43a72 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -24,9 +24,9 @@ def moments_to_shear(Ixx, Iyy, Ixy): def anacal_mag_response(flux, response): - mag_resp = np.divide( -2.5 / np.log(10) * response, - flux, - out=np.full_like(f, np.nan), - where=f > 0 + mag_resp = np.divide( -2.5 / np.log(10) * response, + flux, + out=np.full_like(flux, np.nan, dtype=np.float64), + where=flux > 0, ) return mag_resp \ No newline at end of file From 2312eb1241cccc8107a56dc30cd638667f4faa05 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Tue, 14 Jul 2026 23:27:27 -0700 Subject: [PATCH 62/78] fix R_sel normalization and guard nanojansky_err_to_mag_ab MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two equivalent bookkeeping conventions give an unbiased shear: (A) = Σ(wsel·e) / N R = Σ(de/dg·wsel + dw/dg·e) / N → R_sel normalized by count. (B) = Σ(wsel·e) / N / R = Σ(de/dg·wsel + dw/dg·e) / N / → R_sel normalized by sum(wsel). The rest of the code — R_shape, R_weight, and every mean-e computation — uses convention (B), so R_sel must match with a sum_weights denominator, not count. Also: nanojansky_err_to_mag_ab returns NaN for flux <= 0 via np.divide(where=flux>0) instead of producing negative sigma_mag and RuntimeWarnings. --- txpipe/shear_calibration/calibration_calculators.py | 4 ++-- txpipe/utils/conversion.py | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 4eb5a6c30..cae035c3e 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -1036,8 +1036,8 @@ def collect(self, comm=None, allgather=False) -> BinStats: R_shape = 0.5 * (mean_de1_dg1 + mean_de2_dg2) R_weight = 0.5 * (sum_det_response_e1 + sum_det_response_e2) / sum_weights - R_sel_1 = (sum_sel_response_e1p1 - sum_sel_response_e1m1) / (2.0 * self.delta_gamma) / count - R_sel_2 = (sum_sel_response_e2p2 - sum_sel_response_e2m2) / (2.0 * self.delta_gamma) / count + R_sel_1 = (sum_sel_response_e1p1 - sum_sel_response_e1m1) / (2.0 * self.delta_gamma) / sum_weights + R_sel_2 = (sum_sel_response_e2p2 - sum_sel_response_e2m2) / (2.0 * self.delta_gamma) / sum_weights R_sel = 0.5 * (R_sel_1 + R_sel_2) R_total = R_shape + R_weight + R_sel diff --git a/txpipe/utils/conversion.py b/txpipe/utils/conversion.py index 2c9d43a72..e06ad2a53 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -13,7 +13,11 @@ 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): From 04dc546c12987ba26dc8a743d0704273e70636cc Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Tue, 14 Jul 2026 23:49:22 -0700 Subject: [PATCH 63/78] update the examples/anacal so that they are consistent --- examples/anacal/config.yml | 16 +++++++++++++--- examples/anacal/pipeline.yml | 2 +- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/examples/anacal/config.yml b/examples/anacal/config.yml index 764be5626..ced5deb26 100644 --- a/examples/anacal/config.yml +++ b/examples/anacal/config.yml @@ -5,8 +5,18 @@ global: nside: 512 TXIngestAnacal: - bands: griz - collections: "u/xiangchl/dp1/a360_anacal2" + # DP1-v2 AnaCal outputs for the a360 cluster field (mergePatches per-tract). + butler_config_file: /global/cfs/cdirs/lsst/production/gen3/rubin/DP1/repo/butler.yaml + collections: "u/xiangchl/dp1-v2/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". a360 was + # measured on griz only (no u/y coadds). prefix: "fpfs1" - scale: "fpfs1" \ No newline at end of file + scale: "gauss2" + bands: "griz" diff --git a/examples/anacal/pipeline.yml b/examples/anacal/pipeline.yml index b9ba1ed8f..d35a3abbc 100644 --- a/examples/anacal/pipeline.yml +++ b/examples/anacal/pipeline.yml @@ -14,7 +14,7 @@ site: modules: > txpipe -config: examples/Anacal_test/config.yml +config: examples/anacal/config.yml inputs: anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits From e112deaad622d9db006972324e503140dd3c65be Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Wed, 15 Jul 2026 20:15:26 -0700 Subject: [PATCH 64/78] =?UTF-8?q?1.=20Extended=20the=20examples/anacal/=20?= =?UTF-8?q?pipeline.=20Added=20TXSourceSelectorAnaCal=20=20=20=20as=20a=20?= =?UTF-8?q?second=20stage=20after=20ingest,=20with=20a=20matching=20config?= =?UTF-8?q?.yml=20block:=20=20=20=20delta=5Fgamma:=200.01,=20mask=5Fthresh?= =?UTF-8?q?old:=2040,=20s2n=5Fcut:=205.0,=20T=5Fcut:=200.1,=20=20=20=205-b?= =?UTF-8?q?in=20source=5Fzbin=5Fedges,=20input=5Fpz:=20true.=202.=20Ingest?= =?UTF-8?q?ed=20photo-z=20point=20estimates.=20TXIngestAnacal=20now=20read?= =?UTF-8?q?s=20=20=20=20zmode=5F{0,1p,1m,2p,2m}=20from=20the=20merged=20ca?= =?UTF-8?q?talog=20and=20exposes=20them=20on=20the=20=20=20=20shear=20cata?= =?UTF-8?q?log=20as=20mean=5Fz,=20mean=5Fz=5F{1p,1m,2p,2m}=20=E2=80=94=20b?= =?UTF-8?q?aseline=20for=20=20=20=20tomographic=20binning,=20the=20four=20?= =?UTF-8?q?variants=20for=20the=20R=5Fsel=20bin-migration=20=20=20=20term.?= =?UTF-8?q?=203.=20Built=20=C2=B1=CE=B3=20variants=20of=20every=20quantity?= =?UTF-8?q?=20the=20selector=20cuts=20on.=20In=20=20=20=20TXSourceSelector?= =?UTF-8?q?AnaCal.apply=5Fsimple=5Fredshift=5Fcut,=20injected=20=20=20=20s?= =?UTF-8?q?2n=5F{1p,1m,2p,2m}=20=3D=20s2n=20=C2=B1=20=CE=B4=CE=B3=20=C2=B7?= =?UTF-8?q?=20ds2n=5Fdg{1,2},=20and=20analogous=20variants=20=20=20=20of?= =?UTF-8?q?=20m00,=20m20=20(from=20dm00=5Fdg,=20dm20=5Fdg)=20and=20zbin=20?= =?UTF-8?q?(from=20mean=5Fz=5F*).=20All=20=20=20=20four=20cut=20quantities?= =?UTF-8?q?=20now=20respond=20to=20shear=20via=20a=20single=20uniform=20?= =?UTF-8?q?=20=20=20injection=20step.=204.=20Unified=20the=20calculator's?= =?UTF-8?q?=20cut=20logic=20under=20=5FDataWrapper.=20Deleted=20the=20=20?= =?UTF-8?q?=20=20ad-hoc=20get=5Fsubmask=20size-cut=20machinery=20in=20AnaC?= =?UTF-8?q?alCalculator.add=5Fdata.=20=20=20=20The=20selector=20now=20appl?= =?UTF-8?q?ies=20mask=20+=20s2n=20+=20size=20+=20zbin=20in=20one=20place;?= =?UTF-8?q?=20the=20=20=20=20calculator=20just=20runs=20the=20same=20selec?= =?UTF-8?q?tor=20on=20five=20wrappers=20(baseline=20+=20=20=20=20=5F1p/=5F?= =?UTF-8?q?1m/=5F2p/=5F2m),=20and=20the=20wrapper's=20suffix-lookup=20rout?= =?UTF-8?q?es=20the=20shifted=20=20=20=20columns=20automatically.=20AnaCal?= =?UTF-8?q?Calculator.=5F=5Finit=5F=5F=20dropped=20its=20=20=20=20mask=5Ft?= =?UTF-8?q?hreshold=20argument=20as=20a=20consequence.=205.=20Consistent?= =?UTF-8?q?=20response=20formulae.=20Rewrote=20R=5Fshape,=20R=5Fweight,=20?= =?UTF-8?q?R=5Fsel=5F{1,2}=20=20=20=20in=20the=20calculator=20so=20all=20t?= =?UTF-8?q?hree=20share=20the=20same=20code=20shape:=20=CE=A3(w=20=C2=B7?= =?UTF-8?q?=20X)=20/=20=20=20=20sum=5Fweights=20for=20R=5Fshape=20and=20R?= =?UTF-8?q?=5Fweight,=20and=20[=CE=A3(w=C2=B7e)|sel=5F+=20=E2=88=92=20=20?= =?UTF-8?q?=20=20=CE=A3(w=C2=B7e)|sel=5F=E2=88=92]=20/=20(2=20=CE=B4=CE=B3?= =?UTF-8?q?=20=C2=B7=20sum=5Fweights)=20for=20R=5Fsel.=20R=5Fshape=20switc?= =?UTF-8?q?hed=20from=20=20=20=20a=20ParallelMean=20to=20a=20plain=20scala?= =?UTF-8?q?r=20sum=20for=20parity;=20det=5Fresponse=5F*=20=20=20=20renamed?= =?UTF-8?q?=20to=20weight=5Fresponse=5F*=20to=20match=20the=20"R=5Fweight"?= =?UTF-8?q?=20label;=20a=20small=20=20=20=20=5Freduce=20helper=20replaced?= =?UTF-8?q?=2030=20lines=20of=20MPI=20reduce=20boilerplate.=206.=20Promote?= =?UTF-8?q?d=20the=20size=20threshold=20to=20config.=20The=20AnaCal=20size?= =?UTF-8?q?=20cut,=20=20=20=20previously=20hard-coded=20as=20(m00=20+=20m2?= =?UTF-8?q?0)=20/=20m00=20>=200.1,=20is=20now=20>=20T=5Fcut=20=E2=80=94=20?= =?UTF-8?q?=20=20=20driven=20by=20the=20base=20selector's=20T=5Fcut=20conf?= =?UTF-8?q?ig=20option=20that=20was=20previously=20=20=20=20unused=20for?= =?UTF-8?q?=20AnaCal.=20TXSourceSelectorAnaCal.config=5Foptions=20gives=20?= =?UTF-8?q?it=20a=20=20=20=20default=20of=200.1=20alongside=20a=20new=20s2?= =?UTF-8?q?n=5Fcut=20default=20of=205.0,=20so=20a=20minimal=20=20=20=20yam?= =?UTF-8?q?l=20block=20still=20runs;=20both=20can=20be=20overridden=20in?= =?UTF-8?q?=20the=20yaml=20for=20=20=20=20tighter/looser=20cuts.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/anacal/config.yml | 36 +++ examples/anacal/pipeline.yml | 15 +- txpipe/ingest/anacal.py | 40 ++- .../calibration_calculators.py | 264 ++++++++---------- txpipe/shear_calibration/calibrators.py | 4 +- txpipe/source_selection/Anacal.py | 117 ++++++-- txpipe/test/test_cal.py | 10 +- 7 files changed, 308 insertions(+), 178 deletions(-) diff --git a/examples/anacal/config.yml b/examples/anacal/config.yml index ced5deb26..2b07a400e 100644 --- a/examples/anacal/config.yml +++ b/examples/anacal/config.yml @@ -20,3 +20,39 @@ TXIngestAnacal: 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 (i-band flux / flux_err). Default in code is 5.0. + s2n_cut: 5.0 + + # AnaCal size cut: (m00 + m20) / m00 > T_cut. Default in code is 0.1. + T_cut: 0.1 + + # 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 diff --git a/examples/anacal/pipeline.yml b/examples/anacal/pipeline.yml index d35a3abbc..d346c33f7 100644 --- a/examples/anacal/pipeline.yml +++ b/examples/anacal/pipeline.yml @@ -1,5 +1,6 @@ stages: - name: TXIngestAnacal + - name: TXSourceSelectorAnaCal output_dir: data/example/output_anacal @@ -17,6 +18,18 @@ modules: > 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-v2 published + # per-field catalog to satisfy the graph builder. anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits -resume: False \ No newline at end of file + # 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/ingest/anacal.py b/txpipe/ingest/anacal.py index 949b46f89..1c8a67b12 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -10,9 +10,9 @@ 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). + 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" @@ -44,7 +44,7 @@ def run(self): self.butler_run() else: self.file_run() - + print("repacking files") repack(self.get_output("shear_catalog")) @@ -59,7 +59,7 @@ def butler_run(self): 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. @@ -147,9 +147,16 @@ def file_run(self): for col in shear_data.keys(): print(" ", col) h5py_shorten(shear_outfile["shear"], col, len(shear_data["ra"])) - + shear_outfile.close() + # 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") + def setup_input(self): prefix = self.config["prefix"] scale = self.config["scale"] @@ -166,7 +173,7 @@ def setup_input(self): ]) cols += ["dwsel"+ suffix for suffix in ["_dg1", "_dg2"]] cols += [ - prefix +delta + suffix + prefix +delta + suffix for delta in ["_de1", "_de2", "_dm00", "_dm20"] for suffix in ["_dg1", "_dg2"] ] @@ -175,8 +182,13 @@ def setup_input(self): cols += [band + "_flux_" + scale + "_err" for band in bands] cols += [band + "_dflux_" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] + # 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 self._PZ_SUFFIXES] + return cols - + def process_anacal_shear_data(self, data): bands = self.config["bands"] s = self.config["scale"] @@ -209,7 +221,17 @@ def process_anacal_shear_data(self, data): dd = data[f"{band}_dflux_{s}_"+d][:] output[f"mag_{band}_{d}"] = anacal_mag_response(f, dd) if band == "i": - output[f"ds2n_{d}"] = dd/f_err + output[f"ds2n_{d}"] = dd/f_err + + # zmode_0 → mean_z (baseline photo-z used by TXSourceSelectorAnaCal + # in input_pz mode for tomographic binning). + # zmode_{1p,1m,2p,2m} → mean_z_{1p,1m,2p,2m} (shifted variants + # used by the AnaCal calculator's ±γ selection response — the + # _DataWrapper suffix lookup routes them into the selector when it + # runs on the shifted samples). + output["mean_z"] = data["zmode_0"][:] + for suf in ("1p", "1m", "2p", "2m"): + output[f"mean_z_{suf}"] = data[f"zmode_{suf}"][:] return output diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index cae035c3e..cafb29d03 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -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. """ @@ -849,48 +849,57 @@ 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, mask_threshold): + 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, + 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 selector can take further *args and **kwargs, passed in when adding - data. - + 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 + 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 - The difference in applied g between 1p and 1m variants - mask_threshold: float - The value of the mask at which to cut + 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 - self.mask_threshold = mask_threshold - self.response_means = ParallelMean(size=2) - self.sel_response_e1p1 = 0 - self.sel_response_e1m1 = 0 - self.sel_response_e2p2 = 0 - self.sel_response_e2m2 = 0 - self.det_response_e1 = 0 - self.det_response_e2 = 0 - #self.sel_bias_means = ParallelMean(size=8) + # Every accumulator below is Σ(per-source quantity). In collect() + # each of R_shape, R_weight, R_sel is computed as + # 0.5 * (accum_e1 + accum_e2) / sum_weights + # (with a 2·delta_gamma factor for R_sel because it's a finite + # difference). Same code shape for all three. + self.shape_response_e1 = 0.0 # Σ w · de1/dg1 + self.shape_response_e2 = 0.0 # Σ w · de2/dg2 + self.weight_response_e1 = 0.0 # Σ (dw/dg1) · e1 + self.weight_response_e2 = 0.0 # Σ (dw/dg2) · e2 + self.sel_response_e1p1 = 0.0 # Σ w · e1 · 𝟙[sel_1p] + self.sel_response_e1m1 = 0.0 + self.sel_response_e2p2 = 0.0 + self.sel_response_e2m2 = 0.0 def add_data(self, data, *args, **kwargs): """Select objects from a new chunk of data and tally their responses - + Parameters ---------- data: dict @@ -899,79 +908,68 @@ def add_data(self, data, *args, **kwargs): 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 """ - data = _DataWrapper(data, "") - select = self.selector(data, *args, **kwargs) - - e1 = data["e1"] - e2 = data["e2"] - weight = data["weight"] - weight_dg1 = data["weight_dg1"] - weight_dg2 = data["weight_dg2"] - de1_dg1 = data["de1_dg1"] - de2_dg2 = data["de2_dg2"] - - m00 = data["m00"] - m20 = data["m20"] - - dm00_dg1 = data["dm00_dg1"] - dm00_dg2 = data["dm00_dg2"] - dm20_dg1 = data["dm20_dg1"] - dm20_dg2 = data["dm20_dg2"] - mask_value = data["mask_value"] + # 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 (mask_value, 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") - n = e1[select].size + # 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) - # Record the count for this chunk, for summation later - self.count += n - self.sum_weights += np.sum(weight[select]) - self.sum_sq_weights += np.sum(weight[select]**2) + e1 = data_00["e1"] + e2 = data_00["e2"] + weight = data_00["weight"] + weight_dg1 = data_00["weight_dg1"] + weight_dg2 = data_00["weight_dg2"] + de1_dg1 = data_00["de1_dg1"] + de2_dg2 = data_00["de2_dg2"] - # Next we calculate the cuts per variation masks - wsel = weight[select] + # Baseline slices shared by every accumulator below. + w_sel = weight[select] e1_sel = e1[select] e2_sel = e2[select] - mask = mask_value[select] < self.mask_threshold - - mask_p1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], +1) - mask_m1 = mask & self.get_submask(m00[select], m20[select], dm00_dg1[select], dm20_dg1[select], -1) - mask_p2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], +1) - mask_m2 = mask & self.get_submask(m00[select], m20[select], dm00_dg2[select], dm20_dg2[select], -1) - - self.sel_response_e1p1 += np.sum(wsel[mask_p1] * e1_sel[mask_p1]) - self.sel_response_e1m1 += np.sum(wsel[mask_m1] * e1_sel[mask_m1]) - self.sel_response_e2p2 += np.sum(wsel[mask_p2] * e2_sel[mask_p2]) - self.sel_response_e2m2 += np.sum(wsel[mask_m2] * e2_sel[mask_m2]) - # Next we find the means part needed for the Shape response, - # and the Weight-bias response - de1_dg1_sub = de1_dg1[select] - de2_dg2_sub = de2_dg2[select] - dwsel_dg1 = weight_dg1[select] - dwsel_dg2 = weight_dg2[select] - - self.response_means.add_data(0, de1_dg1_sub, wsel) - self.response_means.add_data(1, de2_dg2_sub, wsel) - - self.det_response_e1 += np.sum(dwsel_dg1 * e1_sel) - self.det_response_e2 += np.sum(dwsel_dg2 * e2_sel) - - self.shear_stats.add_data(0, e1[select], wsel) - self.shear_stats.add_data(1, e2[select], wsel) + # Bookkeeping: raw + weighted counts. + self.count += e1_sel.size + self.sum_weights += np.sum(w_sel) + self.sum_sq_weights += np.sum(w_sel ** 2) + + # All three response-term numerators are accumulated in exactly the + # same shape (Σ over the chunk of a per-source product). collect() + # divides each by sum_weights so they match convention (B): + # weighted-mean normalization consistent with mean_e = Σ(w·e)/Σw. + self.shape_response_e1 += np.sum(w_sel * de1_dg1[select]) + self.shape_response_e2 += np.sum(w_sel * de2_dg2[select]) + self.weight_response_e1 += np.sum(weight_dg1[select] * e1_sel) + self.weight_response_e2 += np.sum(weight_dg2[select] * e2_sel) + self.sel_response_e1p1 += np.sum(weight[sel_1p] * e1[sel_1p]) + self.sel_response_e1m1 += np.sum(weight[sel_1m] * e1[sel_1m]) + self.sel_response_e2p2 += np.sum(weight[sel_2p] * e2[sel_2p]) + self.sel_response_e2m2 += np.sum(weight[sel_2m] * e2[sel_2m]) + + self.shear_stats.add_data(0, e1_sel, w_sel) + self.shear_stats.add_data(1, e2_sel, w_sel) return select - def get_submask(self, m00, m20, dm00_dg, dm20_dg, sign): - m0 = m00 + sign * self.delta_gamma * dm00_dg - m2 = m20 + sign * self.delta_gamma * dm20_dg - return (m0 + m2)/ m0 > 0.1 - - def collect(self, comm=None, allgather=False) -> BinStats: """ Finalize and sum up all the response values, and return a BinStats @@ -990,68 +988,46 @@ def collect(self, comm=None, allgather=False) -> BinStats: bin_stats: BinStats An object containing the final calibration and statistics for this bin. """ - # collect all the things we need - if comm is not None: - if allgather: - count = comm.allreduce(self.count) - sum_weights = comm.allreduce(self.sum_weights) - sum_sq_weights = comm.allreduce(self.sum_sq_weights) - sum_sel_response_e1p1 = comm.allreduce(self.sel_response_e1p1) - sum_sel_response_e1m1 = comm.allreduce(self.sel_response_e1m1) - sum_sel_response_e2p2 = comm.allreduce(self.sel_response_e2p2) - sum_sel_response_e2m2 = comm.allreduce(self.sel_response_e2m2) - sum_det_response_e1 = comm.allreduce(self.det_response_e1) - sum_det_response_e2 = comm.allreduce(self.det_response_e2) - else: - count = comm.reduce(self.count) - sum_weights = comm.reduce(self.sum_weights) - sum_sq_weights = comm.reduce(self.sum_sq_weights) - sum_sel_response_e1p1 = comm.reduce(self.sel_response_e1p1) - sum_sel_response_e1m1 = comm.reduce(self.sel_response_e1m1) - sum_sel_response_e2p2 = comm.reduce(self.sel_response_e2p2) - sum_sel_response_e2m2 = comm.reduce(self.sel_response_e2m2) - sum_det_response_e1 = comm.reduce(self.det_response_e1) - sum_det_response_e2 = comm.reduce(self.det_response_e2) - else: - count = self.count - sum_weights = self.sum_weights - sum_sq_weights = self.sum_sq_weights - sum_sel_response_e1p1 = self.sel_response_e1p1 - sum_sel_response_e1m1 = self.sel_response_e1m1 - sum_sel_response_e2p2 = self.sel_response_e2p2 - sum_sel_response_e2m2 = self.sel_response_e2m2 - sum_det_response_e1 = self.det_response_e1 - sum_det_response_e2 = self.det_response_e2 - - # Collect the mean values we need + # 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) + S_e1 = _reduce(self.shape_response_e1) + S_e2 = _reduce(self.shape_response_e2) + W_e1 = _reduce(self.weight_response_e1) + W_e2 = _reduce(self.weight_response_e2) + sel_e1p1 = _reduce(self.sel_response_e1p1) + sel_e1m1 = _reduce(self.sel_response_e1m1) + sel_e2p2 = _reduce(self.sel_response_e2p2) + sel_e2m2 = _reduce(self.sel_response_e2m2) + + # shear_stats has its own MPI reduction inside .collect() mode = "allgather" if allgather else "gather" - _, R = self.response_means.collect(comm, mode) _, mean_e, var_e = self.shear_stats.collect(comm, mode) - # Unpack the flat mean R: - mean_de1_dg1 = R[0] - mean_de2_dg2 = R[1] - - # Reducing down the responses - R_shape = 0.5 * (mean_de1_dg1 + mean_de2_dg2) - R_weight = 0.5 * (sum_det_response_e1 + sum_det_response_e2) / sum_weights - - R_sel_1 = (sum_sel_response_e1p1 - sum_sel_response_e1m1) / (2.0 * self.delta_gamma) / sum_weights - R_sel_2 = (sum_sel_response_e2p2 - sum_sel_response_e2m2) / (2.0 * self.delta_gamma) / sum_weights - R_sel = 0.5 * (R_sel_1 + R_sel_2) + # All three responses in one consistent shape: + # term = 0.5 * (accum_e1 + accum_e2) / sum_weights + # (R_sel picks up the extra 2·delta_gamma from the finite diff). + R_shape = 0.5 * (S_e1 + S_e2) / sum_weights + R_weight = 0.5 * (W_e1 + W_e2) / sum_weights + R_sel_1 = ( + (sel_e1p1 - sel_e1m1) / (2.0 * self.delta_gamma) / sum_weights + ) + R_sel_2 = ( + (sel_e2p2 - sel_e2m2) / (2.0 * self.delta_gamma) / sum_weights + ) + R_sel = 0.5 * (R_sel_1 + R_sel_2) + R_total = R_shape + R_weight + R_sel + + Neff = None if sum_weights is None else sum_weights ** 2 / sum_sq_weights - R_total = R_shape + R_weight + R_sel - - if sum_weights is None: - Neff = None - else: - Neff = 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)) - bin_stats = BinStats(count, Neff, calibrator.mu, sigma_e, sigma, calibrator) - return bin_stats - - - + 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 eab644082..14d96a8ee 100644 --- a/txpipe/shear_calibration/calibrators.py +++ b/txpipe/shear_calibration/calibrators.py @@ -634,7 +634,7 @@ def calibrate_sigma(self, sigma): class AnaCalibrator(MetaCalibrator): - """Stores information needed to calibrate a Anacal shear method""" + """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: @@ -688,7 +688,7 @@ def load(cls, tomo_file): Returns ------- cals: list - A set of AnacaltCalibrators, one per bin + A set of AnaCalibrators, one per bin """ import h5py diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index bdcc914a6..7ce9b8ac9 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -7,10 +7,10 @@ class TXSourceSelectorAnacal(TXSourceSelectorBase): """ - Source selection and tomography for AnaCal catalogs + Source selection and tomography for AnaCal catalogs. - This selector subcallss is designed for anacal-type catalogs like those - that DESC plans to produce for Ruin data. + This selector subclass is designed for anacal-type catalogs like those + that DESC plans to produce for Rubin data. """ name = "TXSourceSelectorAnaCal" @@ -18,11 +18,26 @@ class TXSourceSelectorAnacal(TXSourceSelectorBase): config_options = { **TXSourceSelectorBase.config_options, "delta_gamma": StageParameter( - float, + float, required=True, - msg= "Delta gamma value for the AnaCal response calculation" + 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="AnaCal size cut: (m00 + m20) / m00 > T_cut.", ), - "mask_threshold": StageParameter(float, 40, msg = "mask threshold, for when to mask objects." ), } def data_iterator(self): @@ -57,12 +72,66 @@ def data_iterator(self): shear_cols += band_variants(bands, "mag", "mag_err", shear_catalog_type="anacal") if self.config["input_pz"]: - shear_cols += ["mean_z"] + # Baseline photo-z + the four shifted variants used to derive + # zbin_{1p,1m,2p,2m} for the selection-response calculation. + shear_cols += ["mean_z", "mean_z_1p", "mean_z_1m", "mean_z_2p", "mean_z_2m"] elif self.config["true_z"]: shear_cols += ["redshift_true"] chunk_rows = self.config["chunk_rows"] return self.iterate_hdf("shear_catalog", "shear", shear_cols, chunk_rows) + + def apply_simple_redshift_cut(self, shear_data): + """Extend the base tomographic binning to emit ±γ shifted variants + of every quantity the AnaCal selector cuts on: zbin (from mean_z), + s2n (from the analytic derivative ds2n_dg), and the AnaCal size + moments m00, m20 (from dm00_dg, dm20_dg). + + Once these are in ``pz_data`` (which the base + ``calculate_tomography`` merges into ``data`` before dispatch), + ``AnaCalCalculator`` needs no special-case size-cut logic — the + _DataWrapper suffix lookup routes every cut through the same code + path uniformly. + + zbin variants are only produced under ``input_pz`` (truth-z has no + shear response); s2n and moment variants are always emitted since + the selector applies those cuts in every mode. + """ + pz_data = super().apply_simple_redshift_cut(shear_data) + + dg = self.config["delta_gamma"] + + # zbin_{1p,1m,2p,2m}: only under input_pz (from shifted photo-z). + if self.config["input_pz"]: + edges = self.config["source_zbin_edges"] + for suf in ("1p", "1m", "2p", "2m"): + zz = shear_data[f"mean_z_{suf}"] + b = np.full(len(zz), -1, dtype=int) + for zi in range(len(edges) - 1): + m = (zz >= edges[zi]) & (zz < edges[zi + 1]) + b[m] = zi + pz_data[f"zbin_{suf}"] = b + + # s2n_{1p,1m,2p,2m}: shifted S/N via ds2n_dg{1,2}. + s2n = shear_data["s2n"] + pz_data["s2n_1p"] = s2n + dg * shear_data["ds2n_dg1"] + pz_data["s2n_1m"] = s2n - dg * shear_data["ds2n_dg1"] + pz_data["s2n_2p"] = s2n + dg * shear_data["ds2n_dg2"] + pz_data["s2n_2m"] = s2n - dg * shear_data["ds2n_dg2"] + + # m00 and m20 variants: needed for the size cut (m00+m20)/m00 > 0.1 + # under ±γ. Kept identical in form to the metacal variant scheme. + m00 = shear_data["m00"] + m20 = shear_data["m20"] + pz_data["m00_1p"] = m00 + dg * shear_data["dm00_dg1"] + pz_data["m00_1m"] = m00 - dg * shear_data["dm00_dg1"] + pz_data["m00_2p"] = m00 + dg * shear_data["dm00_dg2"] + pz_data["m00_2m"] = m00 - dg * shear_data["dm00_dg2"] + pz_data["m20_1p"] = m20 + dg * shear_data["dm20_dg1"] + pz_data["m20_1m"] = m20 - dg * shear_data["dm20_dg1"] + pz_data["m20_2p"] = m20 + dg * shear_data["dm20_dg2"] + pz_data["m20_2m"] = m20 - dg * shear_data["dm20_dg2"] + return pz_data def setup_output(self): """ @@ -77,9 +146,8 @@ def setup_output(self): def setup_response_calculators(self, nbin_source): delta_gamma = self.config["delta_gamma"] - mask_threshold = self.config["mask_threshold"] - calculators = [AnaCalCalculator(select_anacal_tomographic_weak_lensing_sample, delta_gamma, mask_threshold) for i in range(nbin_source)] - calculators.append(AnaCalCalculator(select_anacal_weak_lensing_sample, delta_gamma, mask_threshold)) + 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): @@ -88,26 +156,43 @@ def write_tomography(self, outfile, start, end, source_bin, R): group["R"][start:end] = R def select_anacal_weak_lensing_sample(data, config, calling_from_select=False): + """Baseline AnaCal weak-lensing selection. + + Applies four cuts, each looked up through the ``_DataWrapper`` so the + ±γ shifted variants (m00_1p, s2n_1p, zbin_1p, ...) are used + automatically when the caller wraps ``data`` with a suffix: + + 1. mask_value < mask_threshold (shear-independent, no variants). + 2. s2n > s2n_cut (variants from ds2n_dg{1,2}). + 3. AnaCal size cut (m00 + m20) / m00 > T_cut + (variants from dm00_dg, dm20_dg). + 4. zbin >= 0 (variants from shifted mean_z). + """ s2n_cut = config["s2n_cut"] + T_cut = config["T_cut"] verbose = config["verbose"] flag = data["mask_value"] < config["mask_threshold"] - s2n = data["s2n"] - n0 = len(flag) sel = flag f1 = sel.sum() / n0 - sel &= s2n > s2n_cut + sel &= data["s2n"] > s2n_cut f2 = sel.sum() / n0 - sel &= data["zbin"] >= 0 + # AnaCal size cut, was previously duplicated in AnaCalCalculator.get_submask. + m00 = data["m00"] + m20 = data["m20"] + sel &= (m00 + m20) / m00 > T_cut f3 = sel.sum() / n0 + sel &= data["zbin"] >= 0 + f4 = sel.sum() / n0 + if verbose and calling_from_select: - print(f"Tomo selection {f1:.2%} flag, {f2:.2%} SNR, ", end ="") + print(f"Tomo selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.2%} size, ", end="") elif verbose: - print(f"2D selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.2%} any z bin") + print(f"2D selection {f1:.2%} flag, {f2:.2%} SNR, {f3:.2%} size, {f4:.2%} any z bin") print("total 2D", sel.sum()) return sel diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index 1fbe0ba0d..46be1dff3 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -165,12 +165,10 @@ def core_metadet(comm): def core_anacal(comm): - delta_gamma = .02 - mask_threshold = 40 + delta_gamma = 0.02 nproc = 1 if comm is None else comm.size N = 10 - R_shape_true = 0.4 # known test value R_weight_true = 0.3 @@ -193,7 +191,7 @@ def core_anacal(comm): # 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, mask_threshold) + 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) @@ -201,7 +199,7 @@ def core_anacal(comm): # 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, mask_threshold) + 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"]) + @@ -215,7 +213,7 @@ def core_anacal(comm): "weight_dg1": np.full(N, R_weight_true), "weight_dg2": np.full(N, R_weight_true), } - cal = AnaCalCalculator(select_all_anacal, delta_gamma, mask_threshold) + 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"]) + From 6cf485e5f5bd3d4d421d0a4e232d34d2c0ae4bd6 Mon Sep 17 00:00:00 2001 From: empEvil Date: Thu, 16 Jul 2026 14:38:37 +0100 Subject: [PATCH 65/78] Moving the PZ_SUFFIXES for clarity --- txpipe/ingest/anacal.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 1c8a67b12..5ef127f7c 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -7,6 +7,12 @@ from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, anacal_mag_response 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): """ @@ -150,13 +156,6 @@ def file_run(self): shear_outfile.close() - # 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") - def setup_input(self): prefix = self.config["prefix"] scale = self.config["scale"] @@ -185,7 +184,7 @@ def setup_input(self): # 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 self._PZ_SUFFIXES] + cols += [f"zmode_{s}" for s in PZ_SUFFIXES] return cols From 77654ac6d3c115fa042689df5e693390fc6db4ee Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Wed, 15 Jul 2026 20:38:35 -0700 Subject: [PATCH 66/78] update test to use consistent delta_gamma for anacal --- txpipe/test/test_cal.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index 46be1dff3..4ed4057a5 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -165,7 +165,11 @@ def core_metadet(comm): def core_anacal(comm): - delta_gamma = 0.02 + # 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 From 0f2d95fac90c118d7e0c32056fd104493425b641 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Thu, 16 Jul 2026 20:46:04 -0700 Subject: [PATCH 67/78] 1. use ParallelMean; 2. use consistent definition with paper e_mea = e_raw * wsel and g= e_mea / R_total wher = sum (R_shape + R_detect) / Ngal + R_sel; and weight is set to 1. This ensures that treecorr get the correct correlation / ^2 --- txpipe/ingest/anacal.py | 39 ++++-- .../calibration_calculators.py | 122 ++++++++++-------- txpipe/source_selection/Anacal.py | 55 +++++--- txpipe/utils/conversion.py | 12 +- 4 files changed, 134 insertions(+), 94 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 5ef127f7c..e3121f0c8 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -192,18 +192,35 @@ def process_anacal_shear_data(self, data): bands = self.config["bands"] s = self.config["scale"] prefix = self.config["prefix"] + # 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"][:] output = { - "ra": data["ra"][:], - "dec": data["dec"][:], - "weight": data["wsel"][:], - "mask_value": data["mask_value"][:], - "weight_dg1": data["dwsel_dg1"][:], - "weight_dg2": data["dwsel_dg2"][:], - "e1": data[f"{prefix}_e1"][:], - "e2": data[f"{prefix}_e2"][:], - "m00": data[f"{prefix}_m00"][:], - "m20": data[f"{prefix}_m20"][:], - } + "ra": data["ra"][:], + "dec": data["dec"][:], + "weight": np.ones_like(wsel), # uniform 1 for treecorr + "wsel": wsel, # raw wsel (for R_shape) + "mask_value": data["mask_value"][:], + "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": data[f"{prefix}_m00"][:], + "m20": data[f"{prefix}_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"][:] diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index cafb29d03..03a88c1e5 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -880,22 +880,23 @@ def __init__(self, selector, delta_gamma): 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 - # Every accumulator below is Σ(per-source quantity). In collect() - # each of R_shape, R_weight, R_sel is computed as - # 0.5 * (accum_e1 + accum_e2) / sum_weights - # (with a 2·delta_gamma factor for R_sel because it's a finite - # difference). Same code shape for all three. - self.shape_response_e1 = 0.0 # Σ w · de1/dg1 - self.shape_response_e2 = 0.0 # Σ w · de2/dg2 - self.weight_response_e1 = 0.0 # Σ (dw/dg1) · e1 - self.weight_response_e2 = 0.0 # Σ (dw/dg2) · e2 - self.sel_response_e1p1 = 0.0 # Σ w · e1 · 𝟙[sel_1p] - self.sel_response_e1m1 = 0.0 - self.sel_response_e2p2 = 0.0 - self.sel_response_e2m2 = 0.0 + # 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 @@ -934,39 +935,54 @@ def add_data(self, data, *args, **kwargs): sel_2p = self.selector(data_2p, *args, **kwargs) sel_2m = self.selector(data_2m, *args, **kwargs) - e1 = data_00["e1"] + # 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"] - weight = data_00["weight"] + 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 shared by every accumulator below. - w_sel = weight[select] - e1_sel = e1[select] + # Baseline slices. + w_sel = wsel[select] + e1_sel = e1[select] # already wsel · e_raw e2_sel = e2[select] - # Bookkeeping: raw + weighted counts. + # 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 three response-term numerators are accumulated in exactly the - # same shape (Σ over the chunk of a per-source product). collect() - # divides each by sum_weights so they match convention (B): - # weighted-mean normalization consistent with mean_e = Σ(w·e)/Σw. - self.shape_response_e1 += np.sum(w_sel * de1_dg1[select]) - self.shape_response_e2 += np.sum(w_sel * de2_dg2[select]) - self.weight_response_e1 += np.sum(weight_dg1[select] * e1_sel) - self.weight_response_e2 += np.sum(weight_dg2[select] * e2_sel) - self.sel_response_e1p1 += np.sum(weight[sel_1p] * e1[sel_1p]) - self.sel_response_e1m1 += np.sum(weight[sel_1m] * e1[sel_1m]) - self.sel_response_e2p2 += np.sum(weight[sel_2p] * e2[sel_2p]) - self.sel_response_e2m2 += np.sum(weight[sel_2m] * e2[sel_2m]) - - self.shear_stats.add_data(0, e1_sel, w_sel) - self.shear_stats.add_data(1, e2_sel, w_sel) + # 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 @@ -998,32 +1014,24 @@ def _reduce(x): count = _reduce(self.count) sum_weights = _reduce(self.sum_weights) sum_sq_weights = _reduce(self.sum_sq_weights) - S_e1 = _reduce(self.shape_response_e1) - S_e2 = _reduce(self.shape_response_e2) - W_e1 = _reduce(self.weight_response_e1) - W_e2 = _reduce(self.weight_response_e2) - sel_e1p1 = _reduce(self.sel_response_e1p1) - sel_e1m1 = _reduce(self.sel_response_e1m1) - sel_e2p2 = _reduce(self.sel_response_e2p2) - sel_e2m2 = _reduce(self.sel_response_e2m2) - - # shear_stats has its own MPI reduction inside .collect() + + # 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) - # All three responses in one consistent shape: - # term = 0.5 * (accum_e1 + accum_e2) / sum_weights - # (R_sel picks up the extra 2·delta_gamma from the finite diff). - R_shape = 0.5 * (S_e1 + S_e2) / sum_weights - R_weight = 0.5 * (W_e1 + W_e2) / sum_weights - R_sel_1 = ( - (sel_e1p1 - sel_e1m1) / (2.0 * self.delta_gamma) / sum_weights - ) - R_sel_2 = ( - (sel_e2p2 - sel_e2m2) / (2.0 * self.delta_gamma) / sum_weights - ) - R_sel = 0.5 * (R_sel_1 + R_sel_2) - R_total = R_shape + R_weight + R_sel + # 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 diff --git a/txpipe/source_selection/Anacal.py b/txpipe/source_selection/Anacal.py index 7ce9b8ac9..14699855a 100644 --- a/txpipe/source_selection/Anacal.py +++ b/txpipe/source_selection/Anacal.py @@ -47,14 +47,17 @@ def data_iterator(self): 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 "mask_value", - "e1", + "e1", # e_meas ≡ wsel · e_raw (pre-multiplied) "e2", + "e1_raw", # raw e — needed by R_detect numerator + "e2_raw", "m00", "m20", "de1_dg1", @@ -82,23 +85,31 @@ def data_iterator(self): return self.iterate_hdf("shear_catalog", "shear", shear_cols, chunk_rows) def apply_simple_redshift_cut(self, shear_data): - """Extend the base tomographic binning to emit ±γ shifted variants - of every quantity the AnaCal selector cuts on: zbin (from mean_z), - s2n (from the analytic derivative ds2n_dg), and the AnaCal size - moments m00, m20 (from dm00_dg, dm20_dg). + """Override of the base hook that TXSourceSelectorBase.run() calls + under ``input_pz`` / ``true_z``. We keep the base name so + polymorphism fires, but the real work — injecting ±γ shifted + variants of every quantity the AnaCal selector cuts on — lives in + ``add_sheared_variant_columns`` for a more descriptive name. + """ + pz_data = super().apply_simple_redshift_cut(shear_data) + self.add_sheared_variant_columns(pz_data, shear_data) + return pz_data + + def add_sheared_variant_columns(self, pz_data, shear_data): + """Inject ±γ shifted variants of every quantity the AnaCal selector + cuts on: zbin (from mean_z), s2n (from ds2n_dg), and the AnaCal + size moments m00, m20 (from dm00_dg, dm20_dg). - Once these are in ``pz_data`` (which the base - ``calculate_tomography`` merges into ``data`` before dispatch), - ``AnaCalCalculator`` needs no special-case size-cut logic — the - _DataWrapper suffix lookup routes every cut through the same code - path uniformly. + The additions land in ``pz_data`` — the base + ``calculate_tomography`` will merge that dict into ``data`` before + the calculator sees it, so ``_DataWrapper(data, "_1p")["s2n"]`` + transparently routes to the injected ``s2n_1p`` column. This is + why the calculator needs no per-cut special casing. zbin variants are only produced under ``input_pz`` (truth-z has no shear response); s2n and moment variants are always emitted since the selector applies those cuts in every mode. """ - pz_data = super().apply_simple_redshift_cut(shear_data) - dg = self.config["delta_gamma"] # zbin_{1p,1m,2p,2m}: only under input_pz (from shifted photo-z). @@ -119,8 +130,8 @@ def apply_simple_redshift_cut(self, shear_data): pz_data["s2n_2p"] = s2n + dg * shear_data["ds2n_dg2"] pz_data["s2n_2m"] = s2n - dg * shear_data["ds2n_dg2"] - # m00 and m20 variants: needed for the size cut (m00+m20)/m00 > 0.1 - # under ±γ. Kept identical in form to the metacal variant scheme. + # m00 and m20 variants: needed for the size cut (m00+m20)/m00 > T_cut + # under ±γ. Kept identical in form to the metacal variant scheme. m00 = shear_data["m00"] m20 = shear_data["m20"] pz_data["m00_1p"] = m00 + dg * shear_data["dm00_dg1"] @@ -131,8 +142,7 @@ def apply_simple_redshift_cut(self, shear_data): pz_data["m20_1m"] = m20 - dg * shear_data["dm20_dg1"] pz_data["m20_2p"] = m20 + dg * shear_data["dm20_dg2"] pz_data["m20_2m"] = m20 - dg * shear_data["dm20_dg2"] - return pz_data - + def setup_output(self): """ Prepare the output columns for the response values generated by Anacal @@ -146,10 +156,15 @@ def setup_output(self): 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)) + 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"] diff --git a/txpipe/utils/conversion.py b/txpipe/utils/conversion.py index e06ad2a53..7510b45c1 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -28,9 +28,9 @@ def moments_to_shear(Ixx, Iyy, Ixy): def anacal_mag_response(flux, response): - mag_resp = np.divide( -2.5 / np.log(10) * response, - flux, - out=np.full_like(flux, np.nan, dtype=np.float64), - where=flux > 0, - ) - return mag_resp \ No newline at end of file + return np.divide( + -2.5 / np.log(10) * response, + flux, + out=np.full_like(flux, np.nan, dtype=np.float64), + where=flux > 0, + ) \ No newline at end of file From 1812a901fb89e94530d035e9d8a808b9a43c91b3 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Thu, 16 Jul 2026 21:13:19 -0700 Subject: [PATCH 68/78] add necessary columns in test_cal.py --- txpipe/test/test_cal.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index 4ed4057a5..01198e4fc 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -176,9 +176,19 @@ def core_anacal(comm): 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": np.random.normal(0, 0.1, size=N), - "e2": np.random.normal(0, 0.1, size=N), + "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), From a21049eac87b5a3015c3b6990e26f7b697ff81a1 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Thu, 16 Jul 2026 21:44:15 -0700 Subject: [PATCH 69/78] separate the ingestor from selector and calibrator in examples/anacal/ --- examples/anacal/ingest.yml | 28 ++++++++++++++++++++++++++++ examples/anacal/pipeline.yml | 9 ++++----- 2 files changed, 32 insertions(+), 5 deletions(-) create mode 100644 examples/anacal/ingest.yml diff --git a/examples/anacal/ingest.yml b/examples/anacal/ingest.yml new file mode 100644 index 000000000..eb715f2ac --- /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-v2 published + # per-field catalog to satisfy the graph builder. + anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits + +resume: False diff --git a/examples/anacal/pipeline.yml b/examples/anacal/pipeline.yml index d346c33f7..a62a14b90 100644 --- a/examples/anacal/pipeline.yml +++ b/examples/anacal/pipeline.yml @@ -1,5 +1,7 @@ +# AnaCal selector + calibrator pipeline — reads the shear catalog produced +# by ``examples/anacal/ingest.yml``. + stages: - - name: TXIngestAnacal - name: TXSourceSelectorAnaCal output_dir: data/example/output_anacal @@ -18,10 +20,7 @@ modules: > 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-v2 published - # per-field catalog to satisfy the graph builder. - anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits + 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 From a4682351594a53bc170d8f6f117f947963d6b55c Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Wed, 22 Jul 2026 00:32:53 -0700 Subject: [PATCH 70/78] 1. add selection on magnitude in different bands (with selection response); 2. add selection on |e| (with selection response); 3. the magnitude is calculated from anacal code with smooth truncation at 40; --- examples/anacal/config.yml | 26 ++- examples/anacal/ingest.yml | 4 +- txpipe/ingest/anacal.py | 219 ++++++++++++------ txpipe/source_selection/__init__.py | 2 +- .../source_selection/{Anacal.py => anacal.py} | 194 +++++++++++----- 5 files changed, 317 insertions(+), 128 deletions(-) rename txpipe/source_selection/{Anacal.py => anacal.py} (50%) diff --git a/examples/anacal/config.yml b/examples/anacal/config.yml index 2b07a400e..d8707a390 100644 --- a/examples/anacal/config.yml +++ b/examples/anacal/config.yml @@ -5,9 +5,9 @@ global: nside: 512 TXIngestAnacal: - # DP1-v2 AnaCal outputs for the a360 cluster field (mergePatches per-tract). + # 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-v2/a360_anacal2" + 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). @@ -15,7 +15,10 @@ TXIngestAnacal: cosmology_tracts_only: False # Merged catalog carries band-combined shape moments under the "fpfs1_*" - # prefix and Gaussian-aperture fluxes at scale "gauss2". a360 was + # 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" @@ -32,12 +35,25 @@ TXSourceSelectorAnaCal: # AnaCal reports a bit-packed flag; 1 means "any bit set" (strictest). mask_threshold: 40 - # AnaCal S/N cut (i-band flux / flux_err). Default in code is 5.0. + # AnaCal S/N cut (flux / flux_err). Default in code is 5.0. s2n_cut: 5.0 - # AnaCal size cut: (m00 + m20) / m00 > T_cut. Default in code is 0.1. + # 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] diff --git a/examples/anacal/ingest.yml b/examples/anacal/ingest.yml index eb715f2ac..87701e68e 100644 --- a/examples/anacal/ingest.yml +++ b/examples/anacal/ingest.yml @@ -21,8 +21,8 @@ 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-v2 published + # declared input to have a path — pointed at the DP1-v3 published # per-field catalog to satisfy the graph builder. - anacal_catalog: /pscratch/sd/x/xiangchl/data/DP1/catalogs/anacal_catalog_a360.fits + anacal_catalog: /global/cfs/cdirs/desc-wl/projects/anacal/DP1-v3/catalogs/anacal_catalog_a360.fits resume: False diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index e3121f0c8..853c84170 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -1,12 +1,15 @@ from .base import TXIngestCatalogFits -from ..data_types import ShearCatalog, PhotometryCatalog, HDFFile, FileCollection, FitsFile -from .lsst import process_photometry_data, process_shear_data -from .dp1_details import DP1_COSMOLOGY_FIELDS, DP1_TRACTS, DP1_COSMOLOGY_TRACTS, DP1_FIELD_CENTERS, DP1_SURVEY_PROPERTIES, ALL_TRACTS +from ..data_types import ShearCatalog, FitsFile +from .dp1_details import ( + DP1_TRACTS, + DP1_COSMOLOGY_TRACTS, + ALL_TRACTS, +) from ceci.config import StageParameter import numpy as np -from ..utils import nanojansky_err_to_mag_ab, nanojansky_to_mag_ab, anacal_mag_response 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 @@ -14,6 +17,7 @@ # 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 @@ -29,20 +33,49 @@ class TXIngestAnacal(TXIngestCatalogFits): ("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") + "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", + ), } def run(self): @@ -100,13 +133,17 @@ def butler_run(self): for i, ref in enumerate(data_set_refs): tract = ref.dataId["tract"] if tract not in tracts: - print(f"Skipping chunk {i + 1} / {n_chunks} since tract {tract} is not selected") + 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} - ) + d = butler.get( + object_name, + dataId=ref.dataId, + parameters={"columns": input_columns}, + ) chunk_size = len(d) if chunk_size == 0: @@ -116,16 +153,20 @@ def butler_run(self): 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 = 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 {shear_start:,} - {shear_end:,}") + 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:") @@ -136,17 +177,16 @@ def butler_run(self): shear_outfile.close() def file_run(self): - tracts = self.config["tracts"] - n, dtypes = self.get_meta("anacal_catalog") cols = self.setup_input() - prefix = self.config["prefix"] 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 = self.setup_output( + "shear_catalog", "shear", shear_data, n, + ) shear_outfile["shear"].attrs["catalog_type"] = "anacal" print("Trimming shear columns:") @@ -159,27 +199,51 @@ def file_run(self): def setup_input(self): prefix = self.config["prefix"] scale = self.config["scale"] - cols = ( - [ - "ra", - "dec", - "wsel", - "mask_value", - f"{prefix}_e1", - f"{prefix}_e2", - f"{prefix}_m00", - f"{prefix}_m20" - ]) - cols += ["dwsel"+ suffix for suffix in ["_dg1", "_dg2"]] + cols = [ + "ra", + "dec", + "wsel", + "mask_value", + 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"] - ] + prefix + delta + suffix + for delta in ["_de1", "_de2", "_dm00", "_dm20"] + for suffix in ["_dg1", "_dg2"] + ] bands = self.config["bands"] - cols += [band + "_flux_" + scale for band in bands] - cols += [band + "_flux_" + scale + "_err" for band in bands] - cols += [band + "_dflux_" + scale + suffix for band in bands for suffix in ["_dg1", "_dg2"]] + # 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| s2n_cut (variants from ds2n_dg{1,2}). - 3. AnaCal size cut (m00 + m20) / m00 > T_cut - (variants from dm00_dg, dm20_dg). - 4. zbin >= 0 (variants from shifted mean_z). + 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"] + t_cut = config["T_cut"] + esq_max = config["emax"] ** 2 verbose = config["verbose"] flag = data["mask_value"] < config["mask_threshold"] @@ -195,24 +266,41 @@ def select_anacal_weak_lensing_sample(data, config, calling_from_select=False): sel &= data["s2n"] > s2n_cut f2 = sel.sum() / n0 - # AnaCal size cut, was previously duplicated in AnaCalCalculator.get_submask. + # Band-combined size cut on fpfs1 moments. m00 = data["m00"] m20 = data["m20"] - sel &= (m00 + m20) / m00 > T_cut + sel &= (m00 + m20) / m00 > t_cut f3 = sel.sum() / n0 - sel &= data["zbin"] >= 0 + 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, ", end="") + 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, {f4:.2%} any z bin") + 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 = select_anacal_weak_lensing_sample( + data, config, calling_from_select=True, + ) sel &= zbin == bin_index return sel From 9085096c7c3bc825f76f3562f85177aaf3389fa1 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Wed, 22 Jul 2026 00:39:53 -0700 Subject: [PATCH 71/78] removed unused anacal code in utils/conversion.py --- txpipe/utils/__init__.py | 7 ++++++- txpipe/utils/conversion.py | 11 +---------- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/txpipe/utils/__init__.py b/txpipe/utils/__init__.py index 56cc25bbd..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, anacal_mag_response +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 7510b45c1..66ced1b5c 100644 --- a/txpipe/utils/conversion.py +++ b/txpipe/utils/conversion.py @@ -24,13 +24,4 @@ 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 - - -def anacal_mag_response(flux, response): - return np.divide( - -2.5 / np.log(10) * response, - flux, - out=np.full_like(flux, np.nan, dtype=np.float64), - where=flux > 0, - ) \ No newline at end of file + return e1, e2 \ No newline at end of file From ae7e7c29c62eca4f1a4f6795e3e137dfb342b89e Mon Sep 17 00:00:00 2001 From: "Eske M. Pedersen" Date: Wed, 29 Jul 2026 19:00:36 +0100 Subject: [PATCH 72/78] moving variant columns from selector to ingestor --- txpipe/ingest/anacal.py | 66 +++++++++++++++++++++++++++---- txpipe/source_selection/anacal.py | 53 ++++--------------------- 2 files changed, 66 insertions(+), 53 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 853c84170..3794fee2b 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -76,6 +76,10 @@ class TXIngestAnacal(TXIngestCatalogFits): 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." + ) } def run(self): @@ -256,6 +260,7 @@ def process_anacal_shear_data(self, data): bands = self.config["bands"] s = self.config["scale"] prefix = self.config["prefix"] + dg = self.config["delta_gamma"] # 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" @@ -270,6 +275,9 @@ def process_anacal_shear_data(self, data): 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"][:], @@ -282,8 +290,8 @@ def process_anacal_shear_data(self, data): "e2": wsel * e2_raw, "e1_raw": e1_raw, # raw shape (for R_detect) "e2_raw": e2_raw, - "m00": data[f"{prefix}_m00"][:], - "m20": data[f"{prefix}_m20"][:], + "m00": m00, + "m20": m20, } for delta in ["de1", "de2", "dm00", "dm20"]: output[f"{delta}_dg1"] = data[f"{prefix}_{delta}_dg1"][:] @@ -292,9 +300,12 @@ def process_anacal_shear_data(self, data): # 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. - output["s2n"] = data["lsst_i_s2n_fpfs1"][:] - output["ds2n_dg1"] = data["lsst_i_ds2n_fpfs1_dg1"][:] - output["ds2n_dg2"] = data["lsst_i_ds2n_fpfs1_dg2"][:] + 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 @@ -317,9 +328,12 @@ def process_anacal_shear_data(self, data): # Band-combined shape magnitude + shear derivatives — feeds the # |e| T_cut - # under ±γ. Kept identical in form to the metacal variant scheme. - m00 = shear_data["m00"] - m20 = shear_data["m20"] - pz_data["m00_1p"] = m00 + dg * shear_data["dm00_dg1"] - pz_data["m00_1m"] = m00 - dg * shear_data["dm00_dg1"] - pz_data["m00_2p"] = m00 + dg * shear_data["dm00_dg2"] - pz_data["m00_2m"] = m00 - dg * shear_data["dm00_dg2"] - pz_data["m20_1p"] = m20 + dg * shear_data["dm20_dg1"] - pz_data["m20_1m"] = m20 - dg * shear_data["dm20_dg1"] - pz_data["m20_2p"] = m20 + dg * shear_data["dm20_dg2"] - pz_data["m20_2m"] = m20 - dg * shear_data["dm20_dg2"] - - # esq variants: shape-magnitude cut |e| Date: Wed, 29 Jul 2026 19:17:53 +0100 Subject: [PATCH 73/78] fixing g1/g2 vs e1/e2 --- txpipe/diagnostics.py | 2 ++ txpipe/psf_diagnostics.py | 5 +++++ 2 files changed, 7 insertions(+) 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/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] From a84dae6d503cfb4bb8ca6ef759922531c376c7c3 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Wed, 22 Jul 2026 12:52:42 -0400 Subject: [PATCH 74/78] add calibrator to example --- examples/anacal/config.yml | 5 +++++ examples/anacal/pipeline.yml | 1 + 2 files changed, 6 insertions(+) diff --git a/examples/anacal/config.yml b/examples/anacal/config.yml index d8707a390..0212dc9dc 100644 --- a/examples/anacal/config.yml +++ b/examples/anacal/config.yml @@ -72,3 +72,8 @@ TXSourceSelectorAnaCal: 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/pipeline.yml b/examples/anacal/pipeline.yml index a62a14b90..7b0bd341c 100644 --- a/examples/anacal/pipeline.yml +++ b/examples/anacal/pipeline.yml @@ -3,6 +3,7 @@ stages: - name: TXSourceSelectorAnaCal + - name: TXShearCalibration output_dir: data/example/output_anacal From 956f9cc229da229529807e412720084dd3f9579f Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Sun, 30 Aug 2026 21:35:02 -0700 Subject: [PATCH 75/78] Ingest anacal mask column as n_mask_base The DP2 anacal merged catalog carries the Gaussian-weighted masked fraction as n_mask_base, not mask_value. Read and expose it under that name in the anacal ingest, the anacal source selector, and the calibration test fixture. --- txpipe/ingest/anacal.py | 6 +++--- txpipe/shear_calibration/calibration_calculators.py | 2 +- txpipe/source_selection/anacal.py | 6 +++--- txpipe/test/test_cal.py | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 3794fee2b..2b5ac5165 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -207,7 +207,7 @@ def setup_input(self): "ra", "dec", "wsel", - "mask_value", + "n_mask_base", f"{prefix}_e1", f"{prefix}_e2", f"{prefix}_m00", @@ -277,13 +277,13 @@ def process_anacal_shear_data(self, data): 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) - "mask_value": data["mask_value"][:], + "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 diff --git a/txpipe/shear_calibration/calibration_calculators.py b/txpipe/shear_calibration/calibration_calculators.py index 03a88c1e5..688de9655 100755 --- a/txpipe/shear_calibration/calibration_calculators.py +++ b/txpipe/shear_calibration/calibration_calculators.py @@ -919,7 +919,7 @@ def add_data(self, data, *args, **kwargs): # 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 (mask_value, e1, weight, ...) + # Columns without a matching variant (n_mask_base, e1, weight, ...) # fall back to baseline. data_00 = _DataWrapper(data, "") data_1p = _DataWrapper(data, "_1p") diff --git a/txpipe/source_selection/anacal.py b/txpipe/source_selection/anacal.py index 036732575..4bfec82b3 100644 --- a/txpipe/source_selection/anacal.py +++ b/txpipe/source_selection/anacal.py @@ -72,7 +72,7 @@ def data_iterator(self): "dec", "weight", "wsel", # raw wsel — needed by R_shape numerator - "mask_value", + "n_mask_base", "e1", # e_meas ≡ wsel · e_raw (pre-multiplied) "e2", "e1_raw", # raw e — needed by R_detect numerator @@ -204,7 +204,7 @@ def select_anacal_weak_lensing_sample( ±γ 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. mask_value < mask_threshold (shear-independent, no variants). + 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,...} = @@ -221,7 +221,7 @@ def select_anacal_weak_lensing_sample( esq_max = config["emax"] ** 2 verbose = config["verbose"] - flag = data["mask_value"] < config["mask_threshold"] + flag = data["n_mask_base"] < config["mask_threshold"] n0 = len(flag) sel = flag f1 = sel.sum() / n0 diff --git a/txpipe/test/test_cal.py b/txpipe/test/test_cal.py index 01198e4fc..0d691f343 100644 --- a/txpipe/test/test_cal.py +++ b/txpipe/test/test_cal.py @@ -200,7 +200,7 @@ def core_anacal(comm): "dm00_dg2": np.zeros(N), "dm20_dg1": np.zeros(N), "dm20_dg2": np.zeros(N), - "mask_value": np.zeros(N) + "n_mask_base": np.zeros(N) } # case 1: pure shape response From bdb46cc700450fd665d84d7d59d5721189d7f859 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Sun, 30 Aug 2026 21:36:59 -0700 Subject: [PATCH 76/78] Add per-band PSF properties to anacal ingest Emit psf_g1/psf_g2/psf_T_mean for each band in the config bands from that band's PSF HSM second moments (lsst__ext_shapeHSM_HsmPsfMoments_*), NaN-ing out failed HSM fits via the flag column. Controlled by the new add_psf_properties option (default True). --- txpipe/ingest/anacal.py | 56 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 2b5ac5165..4668b66ed 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -79,7 +79,14 @@ class TXIngestAnacal(TXIngestCatalogFits): "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): @@ -249,6 +256,19 @@ def setup_input(self): # TXSourceSelectorAnaCal. cols += ["esq", "desq_dg1", "desq_dg2"] + # 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). @@ -261,6 +281,9 @@ def process_anacal_shear_data(self, data): 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" @@ -341,10 +364,39 @@ def process_anacal_shear_data(self, data): # used by the AnaCal calculator's ±γ selection response — the # _DataWrapper suffix lookup routes them into the selector when it # runs on the shifted samples). + # Photo-z point estimates from the merged catalog itself: + # zmode_0 → mean_z (baseline for tomographic binning), and the four + # shifted variants zmode_{1p,1m,2p,2m} → mean_z_{...} (photoZPipe + # built them at dg=0.01, so the selector must use delta_gamma=0.01). output["mean_z"] = data["zmode_0"][:] for suf in ("1p", "1m", "2p", "2m"): output[f"mean_z_{suf}"] = data[f"zmode_{suf}"][:] - + + # PSF ellipticity + size per band, from each band's PSF HSM second + # moments: psf_g1_{b}=(xx-yy)/(xx+yy), psf_g2_{b}=2xy/(xx+yy), + # psf_T_mean_{b}=xx+yy. Failed HSM fits (flag set) are NaN'd. + if self.config["add_psf_properties"]: + for b in bands: + base = f"lsst_{b}_ext_shapeHSM_HsmPsfMoments" + pxx = np.asarray(data[f"{base}_xx"][:], dtype=np.float64) + pyy = np.asarray(data[f"{base}_yy"][:], dtype=np.float64) + pxy = np.asarray(data[f"{base}_xy"][:], dtype=np.float64) + ptr = pxx + pyy + with np.errstate(divide="ignore", invalid="ignore"): + psf_g1 = (pxx - pyy) / ptr + psf_g2 = 2.0 * pxy / ptr + bad = ~np.isfinite(ptr) | (ptr <= 0) + flag_col = f"{base}_flag" + if flag_col in colnames: + bad |= np.asarray(data[flag_col][:], dtype=bool) + psf_g1[bad] = np.nan + psf_g2[bad] = np.nan + ptr = ptr.copy() + ptr[bad] = np.nan + output[f"psf_g1_{b}"] = psf_g1 + output[f"psf_g2_{b}"] = psf_g2 + output[f"psf_T_mean_{b}"] = ptr + # shifted values that are needed for selection output["s2n_1p"] = s2n + dg * ds2n_dg1 output["s2n_1m"] = s2n - dg * ds2n_dg1 From 0009f706ecccb1d4a355dc1e9a3ae62761ce437b Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Sun, 30 Aug 2026 21:37:53 -0700 Subject: [PATCH 77/78] Ingest all tracts in the collection when no selection is given The anacal ingest fell back to the DP1 ALL_TRACTS list when no field, tract list, or cosmology-only flag was set, which selects nothing on a non-DP1 (e.g. DP2) collection. Use tracts=None to disable the filter and ingest every tract present in the collection instead. --- txpipe/ingest/anacal.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 4668b66ed..913a8a834 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -3,7 +3,6 @@ from .dp1_details import ( DP1_TRACTS, DP1_COSMOLOGY_TRACTS, - ALL_TRACTS, ) from ceci.config import StageParameter import numpy as np @@ -130,7 +129,11 @@ def butler_run(self): elif self.config["cosmology_tracts_only"]: tracts = DP1_COSMOLOGY_TRACTS else: - tracts = ALL_TRACTS + # 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) @@ -143,7 +146,7 @@ def butler_run(self): shear_start = 0 for i, ref in enumerate(data_set_refs): tract = ref.dataId["tract"] - if tract not in tracts: + 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" From 1f4ad0f3e36e7d9bf00acd4d8c3fb74fe15b78d7 Mon Sep 17 00:00:00 2001 From: Xiangchong Li Date: Sun, 30 Aug 2026 22:31:43 -0700 Subject: [PATCH 78/78] add extinction column from dm pipeline --- txpipe/ingest/anacal.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/txpipe/ingest/anacal.py b/txpipe/ingest/anacal.py index 913a8a834..e94657f9b 100644 --- a/txpipe/ingest/anacal.py +++ b/txpipe/ingest/anacal.py @@ -259,6 +259,9 @@ def setup_input(self): # TXSourceSelectorAnaCal. cols += ["esq", "desq_dg1", "desq_dg2"] + # Per-band extinction lsst_a_ 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. @@ -352,6 +355,12 @@ def process_anacal_shear_data(self, data): 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