Add MOM filtering to return values in a MOC and associated weights

CHANGELOG.md

@@ -1,14 +1,24 @@
 # `moc` Change Log
 
+## 0.16.0
+
+Released 2024-XX-XX
+
+### Added
+
+* MOM filtering to return values in a MOC and associated weights
+
+
+
 ## 0.15.0
 
 Released 2024-06-27
 
-## Fixed
+### Fixed
 
 * Remove spurious coma in empty MOC JSON serialization
 
-## Added
+### Added
 
 * Add methods `all_cells_with_unidirectional_neig`
 * Re-export 'OrdinalMap' and 'OrdinalSet'
@@ -20,7 +30,7 @@ Released 2024-06-27
 
 Released 2024-05-28
 
-## Added
+### Added
 
 * Metohd 'all_cells_with_unidirectional_neigs' for AladinLite
 * Re-export `cdshealpix::compass_point::OrdinalMap` and `cdshealpix::compass_point::OrdinalSet` in `moc::range`

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "moc"
-version = "0.15.0"
+version = "0.16.0"
 authors = [
   "F.-X. Pineau <[email protected]>",
   "Matthieu Baumann <[email protected]>"

src/mom/mod.rs

@@ -11,7 +11,7 @@
 //!     + `(key << 2 + 3, value / 4)`
 
 use std::{
-  f64,
+  f64::{self, consts::FRAC_PI_3},
   marker::PhantomData,
   ops::{AddAssign, Mul},
 };
@@ -67,6 +67,13 @@ pub trait HpxMOMIterator<T: Idx, V: Value<T>>: MOMIterator<T, Hpx<T>, V> {
     }
     sum
   }
+
+  fn retain_values_with_weights_in_hpxmoc(
+    self,
+    moc: &RangeMOC<T, Hpx<T>>,
+  ) -> HpxMOMFilter<T, V, Self> {
+    HpxMOMFilter::new(self, moc)
+  }
 }
 
 pub struct HpxMomIter<T: Idx, Q: MocQty<T>, V: Value<T>, I: Sized + Iterator<Item = (T, V)>> {
@@ -101,3 +108,56 @@ impl<T: Idx, Q: MocQty<T>, V: Value<T>, I: Sized + Iterator<Item = (T, V)>> HpxM
   for HpxMomIter<T, Q, V, I>
 {
 }
+
+/// Filter cell which are in a given MOC and Map to return
+/// a value together with the sky area it covers.
+pub struct HpxMOMFilter<'a, T, V, I>
+where
+  T: Idx,
+  V: Value<T>,
+  I: HpxMOMIterator<T, V>,
+{
+  it: I,
+  moc: &'a RangeMOC<T, Hpx<T>>,
+  _phantom: PhantomData<V>,
+}
+
+impl<'a, T, V, I> HpxMOMFilter<'a, T, V, I>
+where
+  T: Idx,
+  V: Value<T>,
+  I: HpxMOMIterator<T, V>,
+{
+  pub fn new(it: I, moc: &'a RangeMOC<T, Hpx<T>>) -> Self {
+    Self {
+      it,
+      moc,
+      _phantom: PhantomData,
+    }
+  }
+}
+
+impl<'a, T, V, I> Iterator for HpxMOMFilter<'a, T, V, I>
+where
+  T: Idx,
+  V: Value<T>,
+  I: HpxMOMIterator<T, V>,
+{
+  type Item = (V, f64); // (value, weight)
+
+  fn next(&mut self) -> Option<Self::Item> {
+    while let Some((uniq, val)) = self.it.next() {
+      let (depth, ipix) = Hpx::<T>::from_uniq_hpx(uniq);
+      let cell_fraction = self.moc.cell_fraction(depth, ipix);
+      if cell_fraction > 0.0 {
+        let cell_area = FRAC_PI_3 / (1_u64 << (depth << 1) as u32) as f64;
+        return Some((val, cell_area * cell_fraction));
+      }
+    }
+    None
+  }
+
+  fn size_hint(&self) -> (usize, Option<usize>) {
+    (0, self.it.size_hint().1)
+  }
+}

src/storage/u64idx/mod.rs

@@ -45,7 +45,9 @@ use crate::{
     CellMOC2IntoIterator, CellOrCellRangeMOC2IntoIterator, RangeMOC2IntoIterator, RangeMOC2Iterator,
   },
   qty::{Frequency, Hpx, MocQty, Time},
-  storage::u64idx::op1::{op1_mom_sum, op1_mom_sum_from_data, op1_mom_sum_from_path},
+  storage::u64idx::op1::{
+    op1_mom_filter, op1_mom_sum, op1_mom_sum_from_data, op1_mom_sum_from_path,
+  },
 };
 
 pub mod common;
@@ -1829,6 +1831,24 @@ impl U64MocStore {
     op1_mom_sum(index, mom_it)
   }
 
+  /// Filter the value of the given multi-order map to return only the one which are in the given MOC,
+  /// together with the associated sky area (or weight).
+  /// # Params
+  /// * `index`: index pf the S-MOC in the storage
+  /// * `mom_it`: iterator on non-overlapping `(uniq, value)` pairs.
+  /// # Output
+  /// * result made of first the vector of values, and the vector of associated weights.
+  pub fn multiordermap_filter_in_moc<I>(
+    &self,
+    index: usize,
+    mom_it: I,
+  ) -> Result<(Vec<f64>, Vec<f64>), String>
+  where
+    I: Sized + Iterator<Item = (u64, f64)>,
+  {
+    op1_mom_filter(index, mom_it)
+  }
+
   /// Sum the value of the multi-order map in the given path which are in the given MOC.
   /// Remark: we have no information and cannot make any guess on the order if te `UNIQ` cell
   /// in the iterator.

src/storage/u64idx/op1.rs

@@ -295,3 +295,19 @@ pub(crate) fn op1_mom_sum_from_data(index: usize, mom_data: &[u8]) -> Result<f64
     InternalMoc::TimeSpace(_) => Err(String::from("MOM Sum not implemented for ST-MOCs.")),
   })
 }
+
+/// Retuns the `(values, weights)` from the `values` of the input MOM which are in the MOC.
+pub(crate) fn op1_mom_filter<I>(index: usize, it: I) -> Result<(Vec<f64>, Vec<f64>), String>
+where
+  I: Sized + Iterator<Item = (u64, f64)>,
+{
+  store::exec_on_one_readonly_moc(index, move |moc| match moc {
+    InternalMoc::Space(moc) => {
+      let mom_it = HpxMomIter::<u64, Hpx<u64>, f64, _>::new(it);
+      Ok(mom_it.retain_values_with_weights_in_hpxmoc(&moc).unzip())
+    }
+    InternalMoc::Time(_) => Err(String::from("MOM Filter not implemented for T-MOCs.")),
+    InternalMoc::Frequency(_) => Err(String::from("MOM Filter not implemented for F-MOCs.")),
+    InternalMoc::TimeSpace(_) => Err(String::from("MOM Filter not implemented for ST-MOCs.")),
+  })
+}