Merge 98a8396 into 9326f1f

bin/picca_export.py

@@ -258,15 +258,22 @@ def main(cmdargs):
 
     # Check if we need blinding and apply it
     if 'BLIND' in data_name or blinding != 'none':
-        if blinding == 'corr_yshift':
-            userprint("Blinding using strategy corr_yshift.")
+        blinding_dir = '/global/cfs/projectdirs/desi/science/lya/y1-kp6/blinding/'
+        blinding_templates = {'desi_m2': {'standard': 'm2_blinding_v1.2_standard_29_03_2022.h5',
+                                          'grid': 'm2_blinding_v1.2_regular_grid_29_03_2022.h5'},
+                              'desi_y1': {'standard': 'y1_blinding_v2_standard_17_12_2022.h5',
+                                          'grid': 'y1_blinding_v2_regular_grid_17_12_2022.h5'},
+                              'desi_y3': {'standard': 'y3_blinding_v3_standard_18_12_2022.h5',
+                                          'grid': 'y3_blinding_v3_regular_grid_18_12_2022.h5'}}
+
+        if blinding in blinding_templates:
+            userprint(f"Blinding using seed for {blinding}")
         else:
-            raise ValueError("Expected blinding to be 'corr_yshift' or 'minimal'."
-                             " Found {}.".format(blinding))
+            raise ValueError(f"Expected blinding to be one of {blinding_templates.keys()}."
+                             f" Found {blinding}.")
 
         if args.blind_corr_type is None:
-            raise ValueError("Blinding strategy 'corr_yshift' requires"
-                             " argument --blind_corr_type.")
+            raise ValueError("Blinding requires argument --blind_corr_type.")
 
         # Check type of correlation and get size and regular binning
         if args.blind_corr_type in ['lyaxlya', 'lyaxlyb']:
@@ -282,12 +289,10 @@ def main(cmdargs):
 
         if corr_size == len(xi):
             # Read the blinding file and get the right template
-            blinding_filename = ('/global/cfs/projectdirs/desi/science/lya/y1-kp6/'
-                                 'blinding/y1_blinding_v1.2_standard_29_03_2022.h5')
+            blinding_filename = blinding_dir + blinding_templates[blinding]['standard']
         else:
             # Read the regular grid blinding file and get the right template
-            blinding_filename = ('/global/cfs/projectdirs/desi/science/lya/y1-kp6/'
-                                 'blinding/y1_blinding_v1.2_regular_grid_29_03_2022.h5')
+            blinding_filename = blinding_dir + blinding_templates[blinding]['grid']
 
         if not os.path.isfile(blinding_filename):
             raise RuntimeError("Missing blinding file. Make sure you are running at"

py/picca/delta_extraction/data_catalogues/desi_data.py

@@ -14,20 +14,21 @@
     accepted_options as accepted_options_quasar_catalogue)
 from picca.delta_extraction.quasar_catalogues.desi_quasar_catalogue import (
     defaults as defaults_quasar_catalogue)
-from picca.delta_extraction.utils import ACCEPTED_BLINDING_STRATEGIES
+from picca.delta_extraction.utils import (
+    ACCEPTED_BLINDING_STRATEGIES, UNBLINDABLE_STRATEGIES)
 from picca.delta_extraction.utils_pk1d import spectral_resolution_desi, exp_diff_desi
 from picca.delta_extraction.utils import update_accepted_options, update_default_options
 
 
 accepted_options = update_accepted_options(accepted_options, accepted_options_quasar_catalogue)
 accepted_options = update_accepted_options(
     accepted_options,
-    ["blinding", "use non-coadded spectra", "wave solution"])
+    ["unblind", "use non-coadded spectra", "wave solution"])
 
 defaults = update_default_options(defaults, {
     "delta lambda": 0.8,
     "delta log lambda": 3e-4,
-    "blinding": "corr_yshift",
+    "unblind": False,
     "use non-coadded spectra": False,
     "wave solution": "lin",
 })
@@ -104,7 +105,7 @@ def __init__(self, config):
         super().__init__(config)
 
         # load variables from config
-        self.blinding = None
+        self.unblind = None
         self.use_non_coadded_spectra = None
         self.__parse_config(config)
 
@@ -118,12 +119,13 @@ def __init__(self, config):
         # read data
         t0 = time.time()
         self.logger.progress("Reading data")
-        is_mock, is_sv = self.read_data()
+        is_mock = self.read_data()
         t1 = time.time()
         self.logger.progress(f"Time spent reading data: {t1-t0}")
 
         # set blinding
-        self.set_blinding(is_mock, is_sv)
+        self.blinding = None
+        self.set_blinding(is_mock)
 
     def __parse_config(self, config):
         """Parse the configuration options
@@ -138,14 +140,9 @@ def __parse_config(self, config):
         DataError upon missing required variables
         """
         # instance variables
-        self.blinding = config.get("blinding")
-        if self.blinding is None:
-            raise DataError("Missing argument 'blinding' required by DesiData")
-        if self.blinding not in ACCEPTED_BLINDING_STRATEGIES:
-            raise DataError(
-                "Unrecognized blinding strategy. Accepted strategies "
-                f"are {ACCEPTED_BLINDING_STRATEGIES}. "
-                f"Found '{self.blinding}'")
+        self.unblind = config.getboolean("unblind")
+        if self.unblind is None:
+            raise DataError("Missing argument 'unblind' required by DesiData")
 
         self.use_non_coadded_spectra = config.getboolean(
             "use non-coadded spectra")
@@ -166,17 +163,14 @@ def read_data(self):
         is_mock: bool
         True if mocks are read, False otherwise
 
-        is_sv: bool
-        True if all the read data belong to SV. False otherwise
-
         Raise
         -----
         DataError if no quasars were found
         """
         raise DataError(
             "Function 'read_data' was not overloaded by child class")
 
-    def set_blinding(self, is_mock, is_sv):
+    def set_blinding(self, is_mock):
         """Set the blinding in Forest.
 
         Update the stored value if necessary.
@@ -185,9 +179,6 @@ def set_blinding(self, is_mock, is_sv):
         ----------
         is_mock: boolean
         True if reading mocks, False otherwise
-
-        is_sv: boolean
-        True if reading SV data only, False otherwise
         """
         # blinding checks
         if is_mock:
@@ -196,19 +187,24 @@ def set_blinding(self, is_mock, is_sv):
                                     "being ignored as mocks should not be "
                                     "blinded. 'none' blinding engaged")
                 self.blinding = "none"
-        elif is_sv:
-            if self.blinding != "none":
-                self.logger.warning(f"Selected blinding, {self.blinding} is "
-                                    "being ignored as SV data should not be "
-                                    "blinded. 'none' blinding engaged")
-                self.blinding = "none"
-        # TODO: remove this when we are ready to unblind
         else:
-            if self.blinding != "corr_yshift":
-                self.logger.warning(f"Selected blinding, {self.blinding} is "
-                                    "being ignored as data should be blinded. "
-                                    "'corr_yshift' blinding engaged")
-                self.blinding = "corr_yshift"
+            if all(self.catalogue["LASTNIGHT"] < 20210520):
+                # sv data, no blinding
+                self.blinding = "none"
+            elif all(self.catalogue["LASTNIGHT"] < 20210801):
+                self.blinding = "desi_m2"
+            elif all(self.catalogue["LASTNIGHT"] < 20220801):
+                self.blinding = "desi_y1"
+            else:
+                self.blinding = "desi_y3"
+
+            if self.unblind:
+                if self.blinding not in UNBLINDABLE_STRATEGIES:
+                    raise DataError(
+                        "In DesiData: Requested unblinding but data requires blinding strategy "
+                        f"{blinding_strategy} and this strategy do not support "
+                        "unblinding. If you believe this is an error, contact "
+                        "picca developers")
 
         # set blinding strategy
         Forest.blinding = self.blinding

py/picca/delta_extraction/data_catalogues/desi_healpix.py

@@ -92,24 +92,18 @@ def read_data(self):
         is_mock: bool
         False for DESI data and True for mocks
 
-        is_sv: bool
-        True if all the read data belong to SV. False otherwise
-
         Raise
         -----
         DataError if no quasars were found
         """
         grouped_catalogue = self.catalogue.group_by(["HEALPIX", "SURVEY"])
 
-        is_sv = True
         is_mock = False
         forests_by_targetid = {}
 
         arguments = []
         for group in grouped_catalogue.groups:
             healpix, survey = group["HEALPIX", "SURVEY"][0]
-            if survey not in ["sv", "sv1", "sv2", "sv3"]:
-                is_sv = False
 
             filename, is_mock_aux = self.get_filename(survey, healpix)
             if is_mock_aux:
@@ -147,7 +141,7 @@ def read_data(self):
             raise DataError("No quasars found, stopping here")
         self.forests = list(forests_by_targetid.values())
 
-        return is_mock, is_sv
+        return is_mock
 
 
 # Class to read in parallel

py/picca/delta_extraction/data_catalogues/desi_tile.py

@@ -54,20 +54,11 @@ def read_data(self):
         is_mock: bool
         False as DESI data are not mocks
 
-        is_sv: bool
-        True if all the read data belong to SV. False otherwise
-
         Raise
         -----
         DataError if the analysis type is PK 1D and resolution data is not present
         DataError if no quasars were found
         """
-        if np.any((self.catalogue['TILEID'] < 60000) &
-                  (self.catalogue['TILEID'] >= 1000)):
-            is_sv = False
-        else:
-            is_sv = True
-
         coadd_name = "spectra" if self.use_non_coadded_spectra else "coadd"
 
         files_in = sorted(
@@ -77,7 +68,7 @@ def read_data(self):
 
         if "cumulative" in self.input_directory:
             petal_tile_night = [
-                f"{entry['PETAL_LOC']}-{entry['TILEID']}-thru{entry['LAST_NIGHT']}"
+                f"{entry['PETAL_LOC']}-{entry['TILEID']}-thru{entry['LASTNIGHT']}"
                 for entry in self.catalogue
             ]
         else:
@@ -131,7 +122,7 @@ def read_data(self):
 
         self.forests = list(forests_by_targetid.values())
 
-        return False, is_sv
+        return False
 
 
 # Class to read in parallel
@@ -248,7 +239,7 @@ def read_file(self, filename, catalogue):
         if "cumulative" in self.input_directory:
             select = ((catalogue['TILEID'] == tile_spec) &
                       (catalogue['PETAL_LOC'] == petal_spec) &
-                      (catalogue['LAST_NIGHT'] == night_spec))
+                      (catalogue['LASTNIGHT'] == night_spec))
         else:
             select = ((catalogue['TILEID'] == tile_spec) &
                       (catalogue['PETAL_LOC'] == petal_spec) &

py/picca/delta_extraction/quasar_catalogues/desi_quasar_catalogue.py

@@ -126,29 +126,49 @@ def add_healpix(self):
         self.catalogue.sort("HEALPIX")
 
     def read_catalogue(self):
-        """Read the DESI quasar catalogue"""
+        """Read the DESI quasar catalogue
+
+        Raise
+        -----
+        QuasarCatalogueError if the catalogue has missing columns or is
+        empty after the filters are applied
+        """
         self.logger.progress(f'Reading catalogue from {self.filename}')
         extnames = [ext.get_extname() for ext in fitsio.FITS(self.filename)]
         if "QSO_CAT" in extnames:
             extension = "QSO_CAT"
         elif "ZCATALOG" in extnames:
             extension = "ZCATALOG"
         else:
-            raise QuasarCatalogueError(
-                f"Could not find valid quasar catalog extension in fits file: {self.filename}")
+            # TODO: this is a patch that should be removed before merging with master
+            # The extension=1 line should be removed and the raise uncommented
+            extension = 1
+            #raise QuasarCatalogueError(
+            #    f"Could not find valid quasar catalog extension in fits file: {self.filename}")
         catalogue = Table(fitsio.read(self.filename, ext=extension))
 
         if 'TARGET_RA' in catalogue.colnames:
             catalogue.rename_column('TARGET_RA', 'RA')
             catalogue.rename_column('TARGET_DEC', 'DEC')
 
+        # mandatory columns
         keep_columns = ['RA', 'DEC', 'Z', 'TARGETID']
+        for col in keep_columns:
+            if col not in catalogue.colnames:
+                raise QuasarCatalogueError(
+                    f"Missing required column {col} in quasar catalogue")
+
+        # optional columns
         if 'TILEID' in catalogue.colnames:
             keep_columns += ['TILEID', 'PETAL_LOC']
+            if 'PETAL_LOC' not in catalogue.colnames:
+                raise QuasarCatalogueError(
+                    "When TILEID is in the catalogue, PETAL_LOC is also "
+                    "expected to be present but it is not.")
         if 'NIGHT' in catalogue.colnames:
             keep_columns += ['NIGHT']
-        if 'LAST_NIGHT' in catalogue.colnames:
-            keep_columns += ['LAST_NIGHT']
+        if 'LASTNIGHT' in catalogue.colnames:
+            keep_columns += ['LASTNIGHT']
         if 'SURVEY' in catalogue.colnames:
             keep_columns += ['SURVEY']
         if 'DESI_TARGET' in catalogue.colnames:

py/picca/delta_extraction/utils.py

@@ -82,8 +82,10 @@
     "LY10": 919.3514,
 }
 
-ACCEPTED_BLINDING_STRATEGIES = ["none", "minimal", "corr_yshift"]
-
+ACCEPTED_BLINDING_STRATEGIES = [
+    "none", "minimal", "desi_m2", "desi_y1", "desi_y3"]
+# TODO: add tags here when we are allowed to unblind them
+UNBLINDABLE_STRATEGIES = []
 
 def class_from_string(class_name, module_name):
     """Return a class from a string. The class must be saved in a module

py/picca/fitter2/data.py

@@ -245,6 +245,10 @@ def __init__(self,dic_init):
         if (self._blinding == 'minimal') and (('fix_at' not in self.par_fixed.keys()) or (not self.par_fixed['fix_at'])):
             raise ValueError("Running with minimal blinding, please fix at (ap is fixed already)!")
 
+        if self._blinding == 'desi_y3':
+            raise ValueError("Running on DESI Y3 data and fitter2 does not support "
+                             "full-shape blinding. Use Vega instead.")
+
         self.dm_met = {}
         self.rp_met = {}
         self.rt_met = {}

py/picca/tests/delta_extraction/abstract_test.py

@@ -153,8 +153,9 @@ def compare_fits(self, orig_file, new_file):
                     self.assertTrue(key in new_header)
                     if not key in ["CHECKSUM", "DATASUM"]:
                         if (orig_header[key] != new_header[key] and
-                                not np.isclose(orig_header[key],
-                                               new_header[key])):
+                                (isinstance(orig_header[key], str) or not
+                                     np.isclose(orig_header[key],
+                                                new_header[key]))):
                             print(f"\nOriginal file: {orig_file}")
                             print(f"New file: {new_file}")
                             print(f"\n For header {orig_header['EXTNAME']}")