Merge pull request #1087 from igmhub/new_tile_petal_night

add reading of exposure info from coadds

Author: julienguy

py/picca/delta_extraction/astronomical_objects/desi_forest.py

@@ -37,7 +37,7 @@ class DesiForest(Forest):
     targetid: int
     Targetid of the object
 
-    tile: list of int
+    tileid: list of int
     Identifier of the tile used in the observation. None for no info
     """
     def __init__(self, **kwargs):
@@ -68,10 +68,20 @@ def __init__(self, **kwargs):
                                           "Missing variable 'targetid'")
         del kwargs["targetid"]
 
-        self.tile = []
-        if kwargs.get("tile") is not None:
-            self.tile.append(kwargs.get("tile"))
-            del kwargs["tile"]
+        self.tileid = []
+        if kwargs.get("tileid") is not None:
+            self.tileid.append(kwargs.get("tileid"))
+            del kwargs["tileid"]
+
+        self.fiber = []
+        if kwargs.get("fiber") is not None:
+            self.fiber.append(kwargs.get("fiber"))
+            del kwargs["fiber"]
+
+        self.expid = []
+        if kwargs.get("expid") is not None:
+            self.expid.append(kwargs.get("expid"))
+            del kwargs["expid"]
 
         # call parent constructor
         kwargs["los_id"] = self.targetid
@@ -98,7 +108,9 @@ def coadd(self, other):
                                           f"{type(other).__name__}")
         self.night += other.night
         self.petal += other.petal
-        self.tile += other.tile
+        self.tileid += other.tileid
+        self.expid += other.expid
+        self.fiber += other.fiber
         super().coadd(other)
 
     def get_header(self):
@@ -120,19 +132,29 @@ def get_header(self):
             },
             {
                 'name': 'NIGHT',
-                'value': "-".join(str(night) for night in self.night),
+                'value': ",".join(str(night) for night in self.night),
                 'comment': "Observation night(s)"
             },
             {
                 'name': 'PETAL',
-                'value': "-".join(str(petal) for petal in self.petal),
+                'value': ",".join(str(petal) for petal in self.petal),
                 'comment': 'Observation petal(s)'
             },
             {
-                'name': 'TILE',
-                'value': "-".join(str(tile) for tile in self.tile),
+                'name': 'TILEID',
+                'value': ",".join(str(tileid) for tileid in self.tileid),
                 'comment': 'Observation tile(s)'
             },
+            {
+                'name': 'EXPID',
+                'value': ",".join(str(expid) for expid in self.expid),
+                'comment': 'Observation expid(s)'
+            },
+            {
+                'name': 'FIBER',
+                'value': ",".join(str(fiber) for fiber in self.fiber),
+                'comment': 'Observation fiber(s)'
+            },
         ]
 
         return header
@@ -150,9 +172,11 @@ def get_metadata(self):
         metadata = super().get_metadata()
         metadata += [
             self.targetid,
-            "-".join(str(night) for night in self.night),
-            "-".join(str(petal) for petal in self.petal),
-            "-".join(str(tile) for tile in self.tile),
+            ",".join(str(n) for night in self.night for n in night),
+            ",".join(str(p) for petal in self.petal for p in petal),
+            ",".join(str(t) for tileid in self.tileid for t in tileid),
+            ",".join(str(e) for expid in self.expid for e in expid),
+            ",".join(str(f) for fiber in self.fiber for f in fiber),
         ]
         return metadata
 
@@ -168,7 +192,8 @@ def get_metadata_dtype(cls):
         data
         """
         dtype = super().get_metadata_dtype()
-        dtype += [('TARGETID', int), ('NIGHT', 'S12'), ('PETAL', 'S12'), ('TILE', 'S12')]
+        dtype += [('TARGETID', int), ('NIGHT', 'S150'), ('PETAL', 'S150'),
+                  ('TILEID', 'S150'), ('EXPID', 'S150'),('FIBER', 'S150')]
         return dtype
 
     @classmethod
@@ -182,5 +207,5 @@ def get_metadata_units(cls):
         A list with the units of the line-of-sight data
         """
         units = super().get_metadata_units()
-        units += ["", "", "", ""]
+        units += ["", "", "", "", "", ""]
         return units

py/picca/delta_extraction/data_catalogues/desi_data.py

@@ -24,7 +24,7 @@
 accepted_options = update_accepted_options(accepted_options,
     ["use non-coadded spectra",
      "uniquify night targetid",
-     "keep single exposures", 
+     "keep single exposures",
      "wave solution"])
 
 defaults = update_default_options(defaults, {
@@ -38,7 +38,7 @@
 defaults = update_default_options(defaults, defaults_quasar_catalogue)
 
 def verify_exposures_shape(forests_by_targetid):
-    """Verify that the exposures have the same shape. 
+    """Verify that the exposures have the same shape.
     If not, it removes them from the dictionnary of forests by targetid.
     Only works for use_non_coadded_spectra and keep_single_exposures options.
 
@@ -134,7 +134,7 @@ class DesiData(Data):
     If True, remove the quasars taken on the same night.
 
     keep_single_exposures: bool
-    If True, the date loadded from non-coadded spectra are not coadded. 
+    If True, the date loadded from non-coadded spectra are not coadded.
     Otherwise, coadd the spectra here.
     """
 
@@ -295,7 +295,7 @@ class DesiDataFileHandler():
     If True, remove the quasars taken on the same night.
 
     keep_single_exposures: bool
-    If True, the date loadded from non-coadded spectra are not coadded. 
+    If True, the date loadded from non-coadded spectra are not coadded.
     Otherwise, coadd the spectra here.
     """
 
@@ -314,7 +314,7 @@ def __init__(self,
         for details
 
         keep_single_exposures: bool
-        If True, the date loadded from non-coadded spectra are not coadded. 
+        If True, the date loadded from non-coadded spectra are not coadded.
         Otherwise, coadd the spectra here.
 
         uniquify_night_targetid: bool
@@ -352,7 +352,8 @@ def format_data(self,
                     catalogue,
                     spectrographs_data,
                     targetid_spec,
-                    reso_from_truth=False):
+                    reso_from_truth=False,
+                    metadata_dict=None):
         """After data has been read, format it into DesiForest instances
 
         Instances will be DesiForest or DesiPk1dForest depending on analysis_type
@@ -399,6 +400,51 @@ def format_data(self,
                         f"for {targetid}")
             else:
                 w_t = w_t[0]
+            if metadata_dict is not None and not self.use_non_coadded_spectra:
+                exp_w_t = np.where(metadata_dict["EXP_TARGETID"] == targetid)[0]
+                expid = metadata_dict["EXP_EXPID"][exp_w_t]
+                night = metadata_dict["EXP_NIGHT"][exp_w_t]
+                petal = metadata_dict["EXP_PETAL"][exp_w_t]
+                tileid = metadata_dict["EXP_TILEID"][exp_w_t]
+                fiber = metadata_dict["EXP_FIBER"][exp_w_t]
+                metadata_dict_targetid = {'expid': expid,
+                                          'night': night,
+                                          'petal': petal,
+                                          'fiber': fiber,
+                                          'tileid': tileid}
+            elif metadata_dict is not None and self.use_non_coadded_spectra:
+                try:
+                    len(metadata_dict["EXPID"][w_t])
+                    expid = metadata_dict["EXPID"][w_t]
+                except TypeError:
+                    expid = [metadata_dict["EXPID"][w_t]]
+                try:
+                    len(metadata_dict["NIGHT"][w_t])
+                    night = metadata_dict["NIGHT"][w_t]
+                except TypeError:
+                    night = [metadata_dict["NIGHT"][w_t]]
+                try:
+                    len(metadata_dict["PETAL"][w_t])
+                    petal = metadata_dict["PETAL"][w_t]
+                except TypeError:
+                    petal = [metadata_dict["PETAL"][w_t]]
+                try:
+                    len(metadata_dict["FIBER"][w_t])
+                    fiber = metadata_dict["FIBER"][w_t]
+                except TypeError:
+                    fiber = [metadata_dict["FIBER"][w_t]]
+                try:
+                    len(metadata_dict["TILEID"][w_t])
+                    tileid = metadata_dict["TILEID"][w_t]
+                except TypeError:
+                    tileid = [metadata_dict["TILEID"][w_t]]
+                metadata_dict_targetid = {'expid': expid,
+                                          'night': night,
+                                          'petal': petal,
+                                          'fiber': fiber,
+                                          'tileid': tileid}
+            else:
+                metadata_dict_targetid = None
             # Construct DesiForest instance
             # Fluxes from the different spectrographs will be coadded
             for spec in spectrographs_data.values():
@@ -431,7 +477,8 @@ def format_data(self,
                                 ivar_i,
                                 w_t,
                                 reso_from_truth,
-                                num_data)
+                                num_data,
+                                metadata_dict=metadata_dict_targetid)
                 else:
                     forests_by_targetid, num_data  = self.update_forest_dictionary(
                             forests_by_targetid,
@@ -443,7 +490,8 @@ def format_data(self,
                             ivar,
                             w_t,
                             reso_from_truth,
-                            num_data)
+                            num_data,
+                            metadata_dict=metadata_dict_targetid)
         return forests_by_targetid, num_data
 
     def update_forest_dictionary(self,
@@ -456,7 +504,8 @@ def update_forest_dictionary(self,
                                  ivar,
                                  w_t,
                                  reso_from_truth,
-                                 num_data):
+                                 num_data,
+                                 metadata_dict=None):
         """Add new forests to the current forest dictonary
 
         Arguments
@@ -491,7 +540,7 @@ def update_forest_dictionary(self,
 
         num_data: int
         The number of instances loaded
-        
+
         Return
         ------
         forests_by_targetid: dict
@@ -508,6 +557,8 @@ def update_forest_dictionary(self,
             "dec": row['DEC'],
             "z": row['Z'],
         }
+        if metadata_dict is not None:
+            args.update(metadata_dict)
         args["log_lambda"] = np.log10(spec['WAVELENGTH'])
 
 
@@ -582,3 +633,51 @@ def read_file(self, filename, catalogue):
         """
         raise DataError(
             "Function 'read_data' was not overloaded by child class")
+
+    def get_metadata_dict(self,fibermap,exp_fibermap,index_unique) :
+        """
+        Constructs a dictionary containing metadata extracted from the provided fibermap 
+        and exp_fibermap.
+
+        Parameters:
+        - fibermap (numpy.ndarray): The primary fibermap data structure containing metadata.
+        - exp_fibermap (numpy.ndarray or None): The exposure-specific fibermap data structure 
+          containing metadata. If None, the function uses the `fibermap` and `index_unique` 
+          to extract metadata.
+        - index_unique (numpy.ndarray): An array of indices used to select unique entries from 
+          the `fibermap`.
+
+        Returns:
+        - metadata_dict (dict): A dictionary where keys are metadata field names and values 
+          are numpy arrays containing the corresponding metadata values. The keys are prefixed 
+          with 'EXP_' if `use_non_coadded_spectra` is False and `exp_fibermap` is not None.
+
+        Notes:
+        - The function checks if certain keys (`TARGETID`, `NIGHT`, `EXPID`, `PETAL_LOC`, 
+          `FIBER`, `TILEID`) exist in the `input_fibermap`. If a key exists, it extracts the 
+          corresponding data; otherwise, it fills the entry with zeros.
+        - The `use_non_coadded_spectra` attribute of the class instance determines whether to 
+          use the `exp_fibermap` or the `fibermap` combined with `index_unique` for extracting 
+          metadata.
+        """
+        input_fibermap = fibermap
+        ikeys=["TARGETID","NIGHT","EXPID","PETAL_LOC","FIBER","TILEID"]
+        if not self.use_non_coadded_spectra :
+            if exp_fibermap is not None :
+                input_fibermap = exp_fibermap
+            okeys=["EXP_TARGETID","EXP_NIGHT","EXP_EXPID","EXP_PETAL","EXP_FIBER","EXP_TILEID"]
+        else :
+            okeys=["TARGETID","NIGHT","EXPID","PETAL","FIBER","TILEID"]
+
+        metadata_dict = {}
+
+        for ikey,okey in zip(ikeys,okeys) :
+            if ikey in input_fibermap.dtype.names :
+                if exp_fibermap is not None :
+                    metadata_dict[okey] = exp_fibermap[ikey]
+                else :
+                    metadata_dict[okey] = fibermap[ikey][index_unique]
+            else :
+                metadata_dict[okey] = np.zeros(index_unique.size,dtype=int)
+
+        return  metadata_dict

py/picca/delta_extraction/data_catalogues/desi_healpix.py

@@ -208,6 +208,10 @@ def read_file(self, filename, catalogue):
             return {}, 0
         # Read targetid from fibermap to match to catalogue later
         fibermap = hdul['FIBERMAP'].read()
+        exp_fibermap = None
+        if not self.use_non_coadded_spectra :
+            if 'EXP_FIBERMAP' in hdul  :
+                exp_fibermap = hdul['EXP_FIBERMAP'].read()
 
         index_unique = np.full(fibermap.shape,True)
         if self.uniquify_night_targetid:
@@ -275,10 +279,13 @@ def _read_resolution(color, indices):
         if hdul_truth is not None:
             hdul_truth.close()
 
+        metadata_dict=self.get_metadata_dict(fibermap,exp_fibermap,index_unique)
+
         forests_by_targetid, num_data = self.format_data(
             catalogue,
             spectrographs_data,
             fibermap["TARGETID"][index_unique],
-            reso_from_truth=reso_from_truth)
+            reso_from_truth=reso_from_truth,
+            metadata_dict=metadata_dict)
 
         return forests_by_targetid, num_data

py/picca/delta_extraction/data_catalogues/desi_healpix_fast.py

@@ -122,6 +122,9 @@ def read_data(self, is_mock=False):
                 self.forests = combine_results(imap_it)
                 t1 = time.time()
                 self.logger.progress(f"Time spent meerging threads: {t1-t0}")
+        else:
+            raise NotImplementedError('fast healpix reading is not implemented'
+                                      'for analyses with "num processors=1"')
 
         if len(self.forests) == 0:
             raise DataError("No quasars found, stopping here")
@@ -198,7 +201,10 @@ def read_file(self, filename, catalogue):
             return {}, 0
         # Read targetid from fibermap to match to catalogue later
         fibermap = hdul['FIBERMAP'].read()
-
+        if 'EXP_FIBERMAP' in hdul :
+            exp_fibermap = hdul['EXP_FIBERMAP'].read()
+        else :
+            exp_fibermap = fibermap
         colors = ["B", "R"]
         if "Z_FLUX" in hdul:
             colors.append("Z")
@@ -222,6 +228,7 @@ def read_file(self, filename, catalogue):
             # It should be there by construction
             targetid = row["TARGETID"]
             w_t = np.where(fibermap["TARGETID"] == targetid)[0]
+            w_t_exp = np.where(exp_fibermap["TARGETID"] == targetid)[0]
             if len(w_t) == 0:
                 self.logger.warning(
                     f"Error reading {targetid}. Ignoring object")
@@ -240,8 +247,16 @@ def read_file(self, filename, catalogue):
                 "ra": row['RA'],
                 "dec": row['DEC'],
                 "z": row['Z'],
-                "log_lambda": log_lambda,
+                "log_lambda": log_lambda
             }
+            ikeys=["NIGHT","PETAL_LOC","FIBER","TILEID","EXPID"]
+            okeys=["night","petal","fiber","tileid","expid"]
+            for ikey,okey in zip(ikeys,okeys) :
+                if ikey in exp_fibermap.dtype.names :
+                    args[okey] = exp_fibermap[ikey][w_t_exp]
+                else :
+                    args[okey] = np.zeros(w_t_exp.size,dtype=int)
+
             forest = DesiForest(**args)
             forest.rebin()
             forests.append(forest)