Merge pull request #358 from LSSTDESC/ssi_des_inputs

SSI DES inputs

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

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

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
+from .ssi import TXIngestSSIGCR, TXMatchSSI, TXIngestSSIMatched, TXIngestSSIMatchedDESBalrog

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()
+
+
+
+

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