Merge dbb4970 into 1d92ac4

bin/picca_cf.py

@@ -53,10 +53,6 @@ def main(cmdargs):
                         required=True,
                         help='Directory to delta files')
 
-    parser.add_argument('--from-image',
-                        action="store_true",
-                        help='Read delta from image format')
-
     parser.add_argument('--in-dir2',
                         type=str,
                         default=None,
@@ -261,7 +257,6 @@ def main(cmdargs):
                                                   cosmo,
                                                   max_num_spec=args.nspec,
                                                   no_project=args.no_project,
-                                                  from_image=args.from_image,
                                                   nproc=args.nproc,
                                                   rebin_factor=args.rebin_factor)
     del z_max
@@ -292,7 +287,6 @@ def main(cmdargs):
             cosmo,
             max_num_spec=args.nspec,
             no_project=args.no_project,
-            from_image=args.from_image,
             nproc=args.nproc,
             rebin_factor=args.rebin_factor)
         del z_max2

bin/picca_cf_angl.py

@@ -53,13 +53,6 @@ def main(cmdargs):
                         required=True,
                         help='Directory to delta files')
 
-    parser.add_argument('--from-image',
-                        type=str,
-                        default=None,
-                        required=False,
-                        help='Read delta from image format',
-                        nargs='*')
-
     parser.add_argument('--in-dir2',
                         type=str,
                         default=None,

bin/picca_xcf.py

@@ -57,13 +57,6 @@ def main(cmdargs):
                         required=True,
                         help='Directory to delta files')
 
-    parser.add_argument('--from-image',
-                        type=str,
-                        default=None,
-                        required=False,
-                        help='Read delta from image format',
-                        nargs='*')
-
     parser.add_argument('--drq',
                         type=str,
                         default=None,
@@ -299,7 +292,6 @@ def main(cmdargs):
                                                   cosmo=cosmo,
                                                   max_num_spec=args.nspec,
                                                   no_project=args.no_project,
-                                                  from_image=args.from_image,
                                                   nproc=args.nproc,
                                                   rebin_factor=args.rebin_factor)
     xcf.data = data

bin/picca_xwick.py

@@ -55,13 +55,6 @@ def main(cmdargs):
                         required=True,
                         help='Directory to delta files')
 
-    parser.add_argument('--from-image',
-                        type=str,
-                        default=None,
-                        required=False,
-                        help='Read delta from image format',
-                        nargs='*')
-
     parser.add_argument('--drq',
                         type=str,
                         default=None,

py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py

@@ -195,42 +195,45 @@ def rebin(self):
         w2: array of bool
         Masking array for the rebinned ivar solution
 
+        bins: array of int
+        Bins of log_lambda with respect to Forest.log_lambda_grid
+
         Raise
         -----
         AstronomicalObjectError if Forest.wave_solution is not 'lin' or 'log'
         """
-        rebin_ivar, orig_ivar, w1, w2, wslice_inner = super().rebin()
+        rebin_ivar, orig_ivar, w1, w2, wslice_inner, bins = super().rebin()
         if len(rebin_ivar) == 0 or np.sum(w2) == 0:
             self.resolution_matrix = np.array([[]])
-            return [], [], [], np.array([]), np.array([])
+            return [], [], [], np.array([]), np.array([]), np.array([])
 
         # apply mask due to cuts in bin
         self.resolution_matrix = self.resolution_matrix[:, w1]
 
         # Find non-empty bins
-        binned_arr_size = self.log_lambda_index.max() + 1
+        binned_arr_size = bins.max() + 1
 
         # Do a simple average when ivar=0
         orig_ivar_2 = orig_ivar[w1]
         w__ = orig_ivar_2>0
         orig_ivar_2[~w__] = 1
-        rebin_reso_ivar = np.bincount(self.log_lambda_index,
+        rebin_reso_ivar = np.bincount(bins,
                                       weights=orig_ivar_2, minlength=binned_arr_size)
 
         # rebin resolution_matrix
         rebin_reso_matrix_aux = np.zeros(
             (self.resolution_matrix.shape[0], binned_arr_size))
         for index, reso_matrix_col in enumerate(self.resolution_matrix):
             rebin_reso_matrix_aux[index, :] = np.bincount(
-                self.log_lambda_index, weights=orig_ivar_2 * reso_matrix_col)
+                bins, weights=orig_ivar_2 * reso_matrix_col)
 
         # apply mask due to rebinned inverse vairane
         self.resolution_matrix = rebin_reso_matrix_aux[:, wslice_inner] / rebin_reso_ivar[
             np.newaxis, wslice_inner]
 
         # return weights and binning solution to be used by child classes if
         # required
-        return rebin_ivar, orig_ivar, w1, w2, wslice_inner
+        return rebin_ivar, orig_ivar, w1, w2, wslice_inner, bins
 
     @classmethod
     def update_class_variables(cls):

py/picca/delta_extraction/astronomical_objects/forest.py

@@ -11,7 +11,13 @@
 from picca.delta_extraction.utils import find_bins
 
 defaults = {
-    "mask fields": ["flux", "ivar", "transmission_correction", "log_lambda"],
+    "mask fields": [
+        "flux",
+        "ivar",
+        "transmission_correction",
+        "log_lambda",
+        "log_lambda_index"
+    ],
 }
 
 @njit
@@ -173,7 +179,7 @@ def rebin(log_lambda, flux, ivar, transmission_correction, z, wave_solution,
         w1 = np.zeros(0, dtype=bool_)
         w2 = np.zeros(0, dtype=bool_)
         return (log_lambda, flux, ivar, transmission_correction, mean_snr, bins,
-                rebin_ivar, orig_ivar, w1, w2, w2)
+                bins, rebin_ivar, orig_ivar, w1, w2, w2)
 
     log_lambda = log_lambda[w1]
     flux = flux[w1]
@@ -218,14 +224,15 @@ def rebin(log_lambda, flux, ivar, transmission_correction, z, wave_solution,
     else:  # we have already checked that it will always be "lin" at this point
         log_lambda = np.log10(10**log_lambda_grid[0] + pixel_step *
                     np.arange(binned_arr_size)[wslice_inner])
+    rebin_bins = find_bins(log_lambda, log_lambda_grid, wave_solution)
 
     # finally update control variables
     snr = flux * np.sqrt(ivar)
     mean_snr = np.mean(snr[w2])
 
     # return weights and binning solution to be used by child classes if
     # required
-    return (log_lambda, flux, ivar, transmission_correction, mean_snr, bins,
+    return (log_lambda, flux, ivar, transmission_correction, mean_snr, bins, rebin_bins,
             rebin_ivar, orig_ivar, w1, w2, wslice_inner)
 
 class Forest(AstronomicalObject):
@@ -656,20 +663,23 @@ def rebin(self):
         w2: array of bool
         Masking array for the rebinned ivar solution
 
+        bins: array of int
+        Bins of log lambda with respect to Forest.log_lambda_grid
+
         Raise
         -----
         AstronomicalObjectError if Forest.wave_solution is not 'lin' or 'log'
         AstronomicalObjectError if ivar only has zeros
         """
         (self.log_lambda, self.flux, self.ivar, self.transmission_correction,
-         self.mean_snr, self.log_lambda_index, rebin_ivar, orig_ivar, w1,
+         self.mean_snr, bins, self.log_lambda_index, rebin_ivar, orig_ivar, w1,
          w2, wslice_inner) = rebin(self.log_lambda, self.flux, self.ivar,
                      self.transmission_correction, self.z, Forest.wave_solution,
                      Forest.log_lambda_grid, Forest.log_lambda_rest_frame_grid)
 
         # return weights and binning solution to be used by child classes if
         # required
-        return rebin_ivar, orig_ivar, w1, w2, wslice_inner
+        return rebin_ivar, orig_ivar, w1, w2, wslice_inner, bins
 
     @classmethod
     def set_class_variables(cls, lambda_min, lambda_max, lambda_min_rest_frame,

py/picca/delta_extraction/astronomical_objects/pk1d_forest.py

@@ -285,32 +285,35 @@ def rebin(self):
         w2: array of bool
         Masking array for the rebinned ivar solution
 
+        bins: array of int
+        Bins of log_lambda with respect to Forest.log_lambda_grid
+
         Raise
         -----
         AstronomicalObjectError if Forest.wave_solution is not 'lin' or 'log'
         """
-        rebin_ivar, orig_ivar, w1, w2, wslice_inner = super().rebin()
+        rebin_ivar, orig_ivar, w1, w2, wslice_inner, bins = super().rebin()
         if len(rebin_ivar) == 0 or np.sum(w2) == 0:
             self.exposures_diff = np.array([])
             self.reso = np.array([])
             self.reso_pix = np.array([])
-            return [], [], [], np.array([]), np.array([])
+            return [], [], [], np.array([]), np.array([]), np.array([])
 
         # apply mask due to cuts in bin
         self.exposures_diff = self.exposures_diff[w1]
         self.reso = self.reso[w1]
         self.reso_pix = self.reso_pix[w1]
 
         # Find non-empty bins
-        binned_arr_size = self.log_lambda_index.max() + 1
+        binned_arr_size = bins.max() + 1
         final_arr_size = np.sum(wslice_inner)
 
         # rebin exposures_diff and reso
-        rebin_exposures_diff = np.bincount(self.log_lambda_index,
+        rebin_exposures_diff = np.bincount(bins,
             weights=orig_ivar[w1] * self.exposures_diff, minlength=binned_arr_size)
-        rebin_reso = np.bincount(self.log_lambda_index,
+        rebin_reso = np.bincount(bins,
                                  weights=orig_ivar[w1] * self.reso, minlength=binned_arr_size)
-        rebin_reso_pix = np.bincount(self.log_lambda_index, weights=orig_ivar[w1] * self.reso_pix,
+        rebin_reso_pix = np.bincount(bins, weights=orig_ivar[w1] * self.reso_pix,
                                      minlength=binned_arr_size)
 
         # Remove empty bins but not ivar
@@ -338,4 +341,4 @@ def rebin(self):
 
         # return weights and binning solution to be used by child classes if
         # required
-        return rebin_ivar, orig_ivar, w1, w2, wslice_inner
+        return rebin_ivar, orig_ivar, w1, w2, wslice_inner, bins

py/picca/io.py

@@ -1415,13 +1415,11 @@ def read_blinding(in_dir):
     return blinding
 
 
-def read_delta_file(filename, from_image=False, rebin_factor=None):
+def read_delta_file(filename, rebin_factor=None):
     """Extracts deltas from a single file.
     Args:
         filename: str
             Path to the file to read        
-        from_image: bool
-            If True read deltas in ImageHDU format.
         rebin_factor: int - default: None
             Factor to rebin the lambda grid by. If None, no rebinning is done.
     Returns:
@@ -1431,7 +1429,8 @@ def read_delta_file(filename, from_image=False, rebin_factor=None):
     """
 
     hdul = fitsio.FITS(filename)
-    if from_image:
+    # If there is an extension called lambda format is image
+    if 'LAMBDA' in hdul:
         deltas = Delta.from_image(hdul)
     else:
         deltas = [Delta.from_fitsio(hdu) for hdu in hdul[1:]]
@@ -1454,7 +1453,6 @@ def read_deltas(in_dir,
                 cosmo,
                 max_num_spec=None,
                 no_project=False,
-                from_image=False,
                 nproc=None,
                 rebin_factor=None):
     """Reads deltas and computes their redshifts.
@@ -1483,8 +1481,6 @@ def read_deltas(in_dir,
         no_project: bool - default: False
             If True, project the deltas (see equation 5 of du Mas des Bourboux
             et al. 2020)
-        from_image: bool - default: False
-            If True, read deltas in ImageHDU format. 
         nproc: int - default: None
             Number of cpus for parallelization. If None, uses all available.
         rebin_factor: int - default: None
@@ -1516,7 +1512,7 @@ def read_deltas(in_dir,
     if rebin_factor is not None:
         userprint(f"Rebinning deltas by a factor of {rebin_factor}\n")
 
-    arguments = [(f, from_image, rebin_factor) for f in files]
+    arguments = [(f, rebin_factor) for f in files]
     pool = Pool(processes=nproc)
     results = pool.starmap(read_delta_file, arguments)
     pool.close()

py/picca/tests/delta_extraction/astronomical_object_tests.py

@@ -527,11 +527,11 @@ def assert_forest_object(self, test_obj, kwargs):
         self.assertTrue(np.allclose(test_obj.ivar, ivar))
 
         if isinstance(test_obj, DesiPk1dForest):
-            self.assertTrue(len(Forest.mask_fields) == 8)
+            self.assertTrue(len(Forest.mask_fields) == 9)
         elif isinstance(test_obj, Pk1dForest):
-            self.assertTrue(len(Forest.mask_fields) == 7)
+            self.assertTrue(len(Forest.mask_fields) == 8)
         else:
-            self.assertTrue(len(Forest.mask_fields) == 4)
+            self.assertTrue(len(Forest.mask_fields) == 5)
         self.assertTrue(Forest.mask_fields[0] == "flux")
         self.assertTrue(Forest.mask_fields[1] == "ivar")
         self.assertTrue(Forest.mask_fields[2] == "transmission_correction")

py/picca/tests/test_3_cor.py

@@ -173,7 +173,6 @@ def test_cf_image_data(self):
         cmd += " --np 15"
         cmd += " --nt 15"
         cmd += " --nproc 1"
-        cmd += " --from-image"
         print(repr(cmd))
         picca.bin.picca_cf.main(cmd.split()[1:])