Add depth_threshold when build chi2 MOMs

Author: fxpineau

crates/cli/src/map.rs

@@ -1,5 +1,4 @@
-use std::marker::PhantomData;
-use std::{error::Error, path::PathBuf};
+use std::{error::Error, path::PathBuf, marker::PhantomData};
 
 use clap::{Args, Subcommand, ValueEnum};
 use log::{error, warn};

crates/cli/src/mom.rs

@@ -18,7 +18,8 @@ use hpxlib::nested::{
     img::{ColorMapFunctionType, PosConversion, Val},
     mom::{
       impls::{bslice::FITSMom, zvec::MomVecImpl},
-      new_chi2_count_merger, Mom as MomTrait, ViewableMom, WritableMom,
+      new_chi2_count_merger_no_depth_threshold, new_chi2_count_merger_with_depth_threshold,
+      Mom as MomTrait, ViewableMom, WritableMom,
     },
     skymap::SkyMapValue,
   },
@@ -64,14 +65,22 @@ pub struct Count {
   /// Completeness of the chi2 distribution of 3 degrees of freedom.
   #[clap(default_value_t = 16.266)]
   threshold: f64,
+  /// Depth threshold
+  #[clap(short = 't', long, value_name = "DEPTH_MIN")]
+  depth_threshold: Option<u8>,
   /// Count MOM destination FITS file.
   output: PathBuf,
   #[clap(subcommand)]
   input: PosOrHash,
 }
 impl Count {
   pub fn exec(self) -> Result<(), Box<dyn Error>> {
-    self.input.exec(self.depth, self.threshold, self.output)
+    self.input.exec(
+      self.depth,
+      self.threshold,
+      self.depth_threshold,
+      self.output,
+    )
   }
 }
 
@@ -89,22 +98,30 @@ enum PosOrHash {
   },
 }
 impl PosOrHash {
-  pub fn exec(self, depth: u8, threshold: f64, output: PathBuf) -> Result<(), Box<dyn Error>> {
+  pub fn exec(
+    self,
+    depth: u8,
+    threshold: f64,
+    depth_threshold: Option<u8>,
+    output: PathBuf,
+  ) -> Result<(), Box<dyn Error>> {
     match self {
-      Self::Pos { input } => build_count_mom_from_pos(depth, threshold, input),
-      Self::Hash { input } => build_count_mom_from_hash(depth, threshold, input),
+      Self::Pos { input } => build_count_mom_from_pos(depth, threshold, depth_threshold, input),
+      Self::Hash { input } => build_count_mom_from_hash(depth, threshold, depth_threshold, input),
     }
     .and_then(|mom| mom.to_fits_file(output, "DENSITY").map_err(|e| e.into()))
   }
 }
 fn build_count_mom_from_pos(
   depth: u8,
   threshold: f64,
+  depth_threshold: Option<u8>,
   input: PosListInput,
 ) -> Result<MomVecImpl<u64, u32>, Box<dyn Error>> {
   struct Op {
     depth: u8,
     threshold: f64,
+    depth_threshold: Option<u8>,
   }
   impl PosItOperation for Op {
     type R = MomVecImpl<u64, u32>;
@@ -123,25 +140,37 @@ fn build_count_mom_from_pos(
         self.depth,
         SortedHashIt::new(self.depth, it),
       );
-      let chi2_merger = new_chi2_count_merger(self.threshold);
-      MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
-        self.depth,
-        hash_count_it,
-        chi2_merger,
-      )
+      match self.depth_threshold {
+        None => MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
+          self.depth,
+          hash_count_it,
+          new_chi2_count_merger_no_depth_threshold(self.threshold),
+        ),
+        Some(depth_threshold) => MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
+          self.depth,
+          hash_count_it,
+          new_chi2_count_merger_with_depth_threshold(self.threshold, depth_threshold),
+        ),
+      }
       .map_err(|e| e.into())
     }
   }
-  input.exec(Op { depth, threshold })
+  input.exec(Op {
+    depth,
+    threshold,
+    depth_threshold,
+  })
 }
 fn build_count_mom_from_hash(
   depth: u8,
   threshold: f64,
+  depth_threshold: Option<u8>,
   input: HashListInput,
 ) -> Result<MomVecImpl<u64, u32>, Box<dyn Error>> {
   struct Op {
     depth: u8,
     threshold: f64,
+    depth_threshold: Option<u8>,
   }
   impl HashItOperation for Op {
     type R = MomVecImpl<u64, u32>;
@@ -157,16 +186,26 @@ fn build_count_mom_from_hash(
         }
       });
       let hash_count_it = SortedHash2HashCountIncludingZeroIt::new_from_infallible(self.depth, it);
-      let chi2_merger = new_chi2_count_merger(self.threshold);
-      MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
-        self.depth,
-        hash_count_it,
-        chi2_merger,
-      )
+      match self.depth_threshold {
+        None => MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
+          self.depth,
+          hash_count_it,
+          new_chi2_count_merger_no_depth_threshold(self.threshold),
+        ),
+        Some(depth_threshold) => MomVecImpl::<u64, u32>::from_hpx_sorted_entries_fallible(
+          self.depth,
+          hash_count_it,
+          new_chi2_count_merger_with_depth_threshold(self.threshold, depth_threshold),
+        ),
+      }
       .map_err(|e| e.into())
     }
   }
-  input.exec(Op { depth, threshold })
+  input.exec(Op {
+    depth,
+    threshold,
+    depth_threshold,
+  })
 }
 
 #[derive(Debug, Subcommand)]
@@ -229,6 +268,9 @@ pub enum OpType {
     /// Completeness of the chi2 distribution of 3 degrees of freedom below which cells are post-merged.
     #[clap(short, long, default_value_t = 16.266)]
     threshold: f64,
+    /// Depth threshold
+    #[clap(short = 't', long, value_name = "DEPTH_MIN")]
+    depth_threshold: Option<u8>,
   },
 }
 
@@ -257,6 +299,7 @@ impl Operation {
         cte,
         value_name,
         threshold,
+        depth_threshold,
       } => match (l, r) {
         (FITSMom::U32U32(_l), FITSMom::U32U32(_r)) => {
           Err(String::from("Method not available for integers").into())
@@ -267,6 +310,7 @@ impl Operation {
             r.get_mom(),
             cte as f32,
             threshold,
+            depth_threshold,
           )
           .to_fits_file(self.output, value_name)
           .map_err(|e| e.into())
@@ -277,6 +321,7 @@ impl Operation {
             r.get_mom(),
             cte,
             threshold,
+            depth_threshold,
           )
           .to_fits_file(self.output, value_name)
           .map_err(|e| e.into())
@@ -290,6 +335,7 @@ impl Operation {
             r.get_mom(),
             cte as f32,
             threshold,
+            depth_threshold,
           )
           .to_fits_file(self.output, value_name)
           .map_err(|e| e.into())
@@ -300,6 +346,7 @@ impl Operation {
             r.get_mom(),
             cte,
             threshold,
+            depth_threshold,
           )
           .to_fits_file(self.output, value_name)
           .map_err(|e| e.into())

src/nested/map/mom/impls/zvec.rs

@@ -767,6 +767,8 @@ where
   /// * `chi2_of_3dof_threshold`: threshold on the value of the chi square distribution with 3
   /// degrees of freedom below which we consider the 4 values of 4 sibling cells as coming
   /// from the same normal distribution which mean and variance comes from a poisson distribution.
+  /// * `depth_threshold`: threshold on `depth` to avoid making to low resolution cells, i.e MOM minimum depth
+  ///
   /// Here a few typical values corresponding the the given completeness:
   /// * Completeness = 90.0% =>  6.251
   /// * Completeness = 95.0% =>  7.815
@@ -778,6 +780,7 @@ where
     rhs: R,
     cte: V,
     chi2_of_3dof_threshold: f64,
+    depth_threshold: Option<u8>,
   ) -> Self
   where
     L: Mom<'s, ZUniqHType = Z, ValueType = V>,
@@ -796,7 +799,7 @@ where
       rhs,
       split,
       op,
-      new_chi2_density_merger(chi2_of_3dof_threshold),
+      new_chi2_density_merger(chi2_of_3dof_threshold, depth_threshold),
     )
   }
 }
@@ -811,7 +814,9 @@ mod tests {
     n_hash,
     nested::map::{
       img::{to_mom_png_file, ColorMapFunctionType, PosConversion},
-      mom::{impls::zvec::MomVecImpl, new_chi2_count_ref_merger, Mom, ZUniqHashT},
+      mom::{
+        impls::zvec::MomVecImpl, new_chi2_count_ref_merger_no_depth_threshold, Mom, ZUniqHashT,
+      },
       skymap::SkyMapEnum,
     },
   };
@@ -882,7 +887,7 @@ mod tests {
       SkyMapEnum::ImplicitU64I32(skymap) => {
         // println!("Skymap size: {}", skymap.len());
 
-        let merger = new_chi2_count_ref_merger(16.266);
+        let merger = new_chi2_count_ref_merger_no_depth_threshold(16.266);
         let mut mom = MomVecImpl::from_skymap_ref(&skymap, merger);
         /*println!("Mom len: {}", mom.entries.len());
         for (z, v) in mom.entries {