diff --git a/examples/ssi/config_des_mag.yml b/examples/ssi/config_des_mag.yml new file mode 100755 index 00000000..efb099df --- /dev/null +++ b/examples/ssi/config_des_mag.yml @@ -0,0 +1,49 @@ +# Values in this section are accessible to all the different stages. +# They can be overridden by individual stages though. +global: + # This is read by many stages that read complete + # catalog data, and tells them how many rows to read + # at once + chunk_rows: 10000 + # These mapping options are also read by a range of stages + pixelization: healpix + nside: 256 + sparse: true # Generate sparse maps - faster if using small areas + + +TXIngestSSIMatchedDESBalrogMag: + name: TXIngestSSIMatchedDESBalrogMag + +TXMeanLensSelectorSSI: + name: TXMeanLensSelectorSSI + # Mag cuts + lens_zbin_edges: [0.20, 0.40, 0.55, 0.70, 0.85, 0.95, 1.05] + selection_type: DESmaglim + maglim_band: i + bright_limit: 17.5 + a: 4 + b: 18 + extra_cols: [EXTENDED_CLASS_SOF, FLAGS_GOLD] + +TXMeanLensSelectorSSIMag: + name: TXMeanLensSelectorSSIMag + # Mag cuts + lens_zbin_edges: [0.20, 0.40, 0.55, 0.70, 0.85, 0.95, 1.05] + selection_type: DESmaglim + maglim_band: i + bright_limit: 17.5 + a: 4 + b: 18 + extra_cols: [EXTENDED_CLASS_SOF, FLAGS_GOLD] + +TXLensCatalogSplitterSSI: + name: TXLensCatalogSplitterSSI + +TXLensCatalogSplitterSSIMag: + name: TXLensCatalogSplitterSSIMag + +TXSSIMagnification: + name: TXSSIMagnification + applied_magnification: 1.02 + n_patches: 500 + bootstrap_error: true diff --git a/examples/ssi/pipeline_des_mag.yml b/examples/ssi/pipeline_des_mag.yml new file mode 100755 index 00000000..78cfac2b --- /dev/null +++ b/examples/ssi/pipeline_des_mag.yml @@ -0,0 +1,74 @@ +# Stages to run +stages: + - name: TXIngestSSIMatchedDESBalrog # process the unmagnified SSI inputs + - name: TXIngestSSIMatchedDESBalrogMag # process the magnified SSI inputs + classname: TXIngestSSIMatchedDESBalrog + aliases: + balrog_matched_catalog: balrog_matched_catalog_mag + matched_ssi_photometry_catalog: matched_ssi_photometry_catalog_mag + + - name: TXMeanLensSelectorSSI #make unmagnified SSI lens samples + classname: TXMeanLensSelector + aliases: + photometry_catalog: matched_ssi_photometry_catalog + lens_photoz_pdfs: balrog_lens_photoz_pdfs + - name: TXLensCatalogSplitterSSI + classname: TXLensCatalogSplitter + aliases: + photometry_catalog: matched_ssi_photometry_catalog + lens_photoz_pdfs: balrog_lens_photoz_pdfs + + - name: TXMeanLensSelectorSSIMag #make magnified SSI lens samples + classname: TXMeanLensSelector + aliases: + photometry_catalog: matched_ssi_photometry_catalog_mag + lens_tomography_catalog_unweighted: lens_tomography_catalog_unweighted_mag + lens_photoz_pdfs: balrog_lens_photoz_pdfs_mag + - name: TXLensCatalogSplitterSSIMag + classname: TXLensCatalogSplitter + aliases: + photometry_catalog: matched_ssi_photometry_catalog_mag + lens_tomography_catalog_unweighted: lens_tomography_catalog_unweighted_mag + binned_lens_catalog_unweighted: binned_lens_catalog_unweighted_mag + lens_photoz_pdfs: balrog_lens_photoz_pdfs_mag + + - name: TXSSIMagnification # compute magnification coefficients for lens sample + aliases: + binned_lens_catalog_nomag: binned_lens_catalog_unweighted + binned_lens_catalog_mag: binned_lens_catalog_unweighted_mag + + +modules: txpipe + +# Where to put outputs +output_dir: data/example/outputs_ssi_des_mag/ + +# How to run the pipeline: mini, parsl, or cwl +launcher: + name: mini + interval: 1.0 + +# Where to run the pipeline: cori-interactive, cori-batch, or local +site: + name: local + max_threads: 2 + +# configuration settings +config: examples/ssi/config_des_mag.yml + +inputs: + # See README for paths to download these files + balrog_matched_catalog: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_matched_catalog_sof_run2a_v1.4.fits + balrog_matched_catalog_mag: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_matched_catalog_sof_run2a-mag_v1.4.fits + balrog_lens_photoz_pdfs: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_matched_catalog_sof_run2a_v1.4_dnf.h5 + balrog_lens_photoz_pdfs_mag: /global/cfs/cdirs/lsst/groups/WL/users/jelvinpo/balrog_txpipe_example/balrog_matched_catalog_sof_run2a-mag_v1.4_dnf.h5 + + fiducial_cosmology: data/fiducial_cosmology.yml + +# if supported by the launcher, restart the pipeline where it left off +# if interrupted +resume: false +# where to put output logs for individual stages +log_dir: data/example/outputs_ssi_des_mag/logs/ +# where to put an overall parsl pipeline log +pipeline_log: data/example/outputs_ssi_des_mag/log.txt diff --git a/txpipe/ingest/__init__.py b/txpipe/ingest/__init__.py index 9ea2a2bb..64b85659 100644 --- a/txpipe/ingest/__init__.py +++ b/txpipe/ingest/__init__.py @@ -2,4 +2,4 @@ from .mocks import TXCosmoDC2Mock, TXGaussianSimsMock from .gcr import TXMetacalGCRInput, TXIngestStars from .redmagic import TXIngestRedmagic -from .ssi import TXIngestSSIGCR, TXMatchSSI \ No newline at end of file +from .ssi import TXIngestSSIGCR, TXMatchSSI, TXIngestSSIMatched, TXIngestSSIMatchedDESBalrog \ No newline at end of file diff --git a/txpipe/ingest/ssi.py b/txpipe/ingest/ssi.py index 708cf202..546496b3 100644 --- a/txpipe/ingest/ssi.py +++ b/txpipe/ingest/ssi.py @@ -1,6 +1,7 @@ from ..base_stage import PipelineStage from ..data_types import ( HDFFile, + FitsFile, ) from ..utils import ( nanojansky_to_mag_ab, @@ -274,3 +275,110 @@ def finalize_output(self, outfile, group, ntot): col.resize((ntot,)) outfile.close() return + + + +class TXIngestSSIMatched(PipelineStage): + """ + Base-stage for ingesting a matched SSI catalog + + This stage will just read in a file in a given format and output to a + HDF file that TXPIPE can use + + """ + + name = "TXIngestSSIMatched" + + outputs = [ + ("matched_ssi_photometry_catalog", HDFFile), + ] + + config_options = { + "chunk_rows": 100_000, + "magnification":0, # magnification label + } + + +class TXIngestSSIMatchedDESBalrog(TXIngestSSIMatched): + """ + Class for ingesting a matched "SSI" catalog from DES (AKA Balrog) + """ + + name = "TXIngestSSIMatchedDESBalrog" + + inputs = [ + ("balrog_matched_catalog", FitsFile), + ] + + def run(self): + """ + Run the analysis for this stage. + """ + print('Ingesting DES Balrog matched catalog') + + #get some basic onfo about the input file + f = self.open_input("balrog_matched_catalog") + n = f[1].get_nrows() + dtypes = f[1].get_rec_dtype()[0] + f.close() + + print(f'{n} objects in matched catalog') + + chunk_rows = self.config["chunk_rows"] + + #we will only load a subset of columns to save space + column_names = { + "bal_id": "bal_id", # Unique identifier for object (created during balrog process) + "true_bdf_mag_deredden": "inj_mag", # Magnitude of the original deep field object, dereddened + "true_id": "inj_id", # Original coadd_obj_id of deep field object + "meas_id": "id", # Coadd_object_id of injection + "meas_ra": "ra", # measured RA of the injection + "meas_dec": "dec", # measured DEC of the injection + "meas_cm_mag_deredden": "mag", # measured magnitude of the injection + "meas_cm_T": "cm_T", # measured size parameter T (x^2+y^2) + "meas_EXTENDED_CLASS_SOF": "EXTENDED_CLASS_SOF", # Star galaxy classifier (0,1=star, 2,3=Galaxy) + "meas_FLAGS_GOLD_SOF_ONLY": "FLAGS_GOLD", # Measured flags (short version) + } + cols = list(column_names.keys()) + + #set up the output file columns + output = self.open_output("matched_ssi_photometry_catalog") + g = output.create_group("photometry") + for col in cols: + dtype = dtypes[col] + + if "_mag" in col: + #per band + dtype = dtype.subdtype[0] + for b in "griz": + g.create_dataset(column_names[col]+f"_{b}", (n,), dtype=dtype) + + #also create an empty array for the mag errors + #TO DO: compute mag errors from flux errors + g.create_dataset(column_names[col]+f"_err_{b}", (n,), dtype=dtype) + else: + g.create_dataset(column_names[col], (n,), dtype=dtype) + + #iterate over the input file and save to the output columns + for (s, e, data) in self.iterate_fits("balrog_matched_catalog", 1, cols, chunk_rows): + print(s,e,n) + for col in cols: + if "_mag" in col: + for iband,b in enumerate("griz"): + g[column_names[col]+f"_{b}"][s:e] = data[col][:,iband] + else: + g[column_names[col]][s:e] = data[col] + + # Set up any dummy columns with sentinal values + # that were not in the original files + dummy_columns = { + "redshift_true":10.0, + } + for col_name in dummy_columns.keys(): + g.create_dataset(col_name, data=np.full(n,dummy_columns[col_name])) + + output.close() + + + + diff --git a/txpipe/lens_selector.py b/txpipe/lens_selector.py index 8d1ff871..d8020a8f 100755 --- a/txpipe/lens_selector.py +++ b/txpipe/lens_selector.py @@ -49,6 +49,7 @@ class TXBaseLensSelector(PipelineStage): "selection_type": "boss", "maglim_band": "i", "maglim_limit": 24.1, + "extra_cols": [""], } @@ -219,6 +220,8 @@ def select_lens(self, phot_data): s = self.select_lens_boss(phot_data) elif t == "maglim": s= self.select_lens_maglim(phot_data) + elif t == "DESmaglim": + s= self.select_lens_DESmaglim(phot_data) else: raise ValueError(f"Unknown lens selection type {t} - expected boss or maglim") ntot = s.size @@ -281,6 +284,30 @@ def select_lens_maglim(self, phot_data): s = (mag_i < limit).astype(np.int8) return s + def select_lens_DESmaglim(self, phot_data): + band = self.config["maglim_band"] + bright_limit = self.config["bright_limit"] + # mag < a*zphot + b + a = self.config["a"] + b = self.config["b"] + + z = phot_data["z"] + mag_i = phot_data[f"mag_{band}"] + sg = phot_data["EXTENDED_CLASS_SOF"] + flags = phot_data["FLAGS_GOLD"] + + cut1 = (mag_i > bright_limit) #cut very bright gaalxies + cut2 = (mag_i < a*z + b) #the z-dependant mag cut + cut3 = (sg == 3) #star/galaxy separator + + # flags, selects objects flagged by bits 2^1, 2^2, 2^3, 2^4, 2^5, 2^6 but not 2^7 + # from Rodriguez-Monroy et al + cut4 = (np.bitwise_and(flags,126) == 0) + + + s = (cut1 & cut2 & cut3 & cut4).astype(np.int8) + return s + def calculate_tomography(self, pz_data, phot_data, lens_gals): nbin = len(self.config["lens_zbin_edges"]) - 1 @@ -318,6 +345,8 @@ def data_iterator(self): print(f"We are cheating and using the true redshift.") chunk_rows = self.config["chunk_rows"] phot_cols = ["mag_i", "mag_r", "mag_g", "redshift_true"] + extra_cols = [c for c in self.config["extra_cols"] if c] + phot_cols += extra_cols # Input data. These are iterators - they lazily load chunks # of the data one by one later when we do the for loop. # This code can be run in parallel, and different processes will @@ -347,6 +376,8 @@ def data_iterator(self): phot_cols = ["mag_i", "mag_r", "mag_g"] z_cols = ["zmean"] rename = {"zmean": "z"} + extra_cols = [c for c in self.config["extra_cols"] if c] + phot_cols += extra_cols it = self.combined_iterators( chunk_rows, @@ -360,7 +391,6 @@ def data_iterator(self): return rename_iterated(it, rename) - class TXModeLensSelector(TXBaseLensSelector): """ Select lens objects based on best-fit redshifts and BOSS criteria @@ -379,6 +409,8 @@ def data_iterator(self): phot_cols = ["mag_i", "mag_r", "mag_g"] z_cols = ["zmode"] rename = {"zmode": "z"} + extra_cols = [c for c in self.config["extra_cols"] if c] + phot_cols += extra_cols it = self.combined_iterators( chunk_rows, @@ -411,6 +443,8 @@ class TXRandomForestLensSelector(TXBaseLensSelector): def data_iterator(self): chunk_rows = self.config["chunk_rows"] phot_cols = ["mag_u", "mag_g", "mag_r", "mag_i", "mag_z", "mag_y"] + extra_cols = [c for c in self.config["extra_cols"] if c] + phot_cols += extra_cols for s, e, data in self.iterate_hdf( "photometry_catalog", "photometry", phot_cols, chunk_rows