Component norm for MC periodogram

Author: hombit

src/features/_periodogram_peaks.rs

@@ -50,13 +50,13 @@ impl PeriodogramPeaks {
             .flat_map(|i| {
                 vec![
                     format!(
-                        "period of the {}{} highest peak of periodogram",
+                        "period of the {}{} highest peak",
                         i + 1,
                         number_ending(i + 1),
                     ),
                     format!(
                         "Spectral density to spectral density standard deviation ratio of \
-                            the {}{} highest peak of periodogram",
+                            the {}{} highest peak",
                         i + 1,
                         number_ending(i + 1)
                     ),

src/features/periodogram.rs

@@ -45,6 +45,10 @@ where
     max_freq_factor: f32,
     nyquist: NyquistFreq,
     pub(crate) feature_extractor: FeatureExtractor<T, F>,
+    // In can be re-defined in MultiColorPeriodogram
+    pub(crate) name_prefix: String,
+    // In can be re-defined in MultiColorPeriodogram
+    pub(crate) description_suffix: String,
     periodogram_algorithm: PeriodogramPower<T>,
     properties: Box<EvaluatorProperties>,
 }
@@ -100,13 +104,13 @@ where
             feature
                 .get_names()
                 .iter()
-                .map(|name| "periodogram_".to_owned() + name),
+                .map(|name| format!("{}_{}", self.name_prefix, name)),
         );
         self.properties.descriptions.extend(
             feature
                 .get_descriptions()
                 .into_iter()
-                .map(|desc| format!("{} of periodogram", desc)),
+                .map(|desc| format!("{} {}", desc, self.description_suffix)),
         );
         self.feature_extractor.add_feature(feature);
         self
@@ -149,41 +153,44 @@ where
 {
     /// New [Periodogram] that finds given number of peaks
     pub fn new(peaks: usize) -> Self {
-        let peaks = PeriodogramPeaks::new(peaks);
-        let peak_names = peaks
-            .get_names()
-            .into_iter()
-            .map(ToOwned::to_owned)
-            .collect();
-        let peak_descriptions = peaks
-            .get_descriptions()
-            .into_iter()
-            .map(ToOwned::to_owned)
-            .collect();
-        let peaks_size_hint = peaks.size_hint();
-        let peaks_min_ts_length = peaks.min_ts_length();
+        Self::with_name_description(
+            peaks,
+            "periodogram",
+            "of periodogram (interpreting frequency as time, power as magnitude)",
+        )
+    }
+
+    pub(crate) fn with_name_description(
+        peaks: usize,
+        name_prefix: impl ToString,
+        description_suffix: impl ToString,
+    ) -> Self {
         let info = EvaluatorInfo {
-            size: peaks_size_hint,
-            min_ts_length: usize::max(peaks_min_ts_length, 2),
+            size: 0,
+            min_ts_length: 2,
             t_required: true,
             m_required: true,
             w_required: false,
             sorting_required: true,
             variability_required: false,
         };
-        Self {
+        let mut slf = Self {
             properties: EvaluatorProperties {
                 info,
-                names: peak_names,
-                descriptions: peak_descriptions,
+                names: vec![],
+                descriptions: vec![],
             }
             .into(),
             resolution: Self::default_resolution(),
+            name_prefix: name_prefix.to_string(),
+            description_suffix: description_suffix.to_string(),
             max_freq_factor: Self::default_max_freq_factor(),
             nyquist: AverageNyquistFreq.into(),
-            feature_extractor: FeatureExtractor::new(vec![peaks.into()]),
+            feature_extractor: FeatureExtractor::new(vec![]),
             periodogram_algorithm: PeriodogramPowerFft::new().into(),
-        }
+        };
+        slf.add_feature(PeriodogramPeaks::new(peaks).into());
+        slf
     }
 }
 
@@ -286,7 +293,7 @@ where
             nyquist,
             feature_extractor,
             periodogram_algorithm,
-            properties: _,
+            ..
         } = f;
 
         let mut features = feature_extractor.into_vec();

src/multicolor/features/mod.rs

@@ -8,4 +8,4 @@ mod color_of_minimum;
 pub use color_of_minimum::ColorOfMinimum;
 
 mod multi_color_periodogram;
-pub use multi_color_periodogram::MultiColorPeriodogram;
+pub use multi_color_periodogram::{MultiColorPeriodogram, MultiColorPeriodogramNormalisation};

src/multicolor/features/multi_color_periodogram.rs

@@ -13,6 +13,20 @@ use crate::periodogram::{self, NyquistFreq, PeriodogramPower};
 use ndarray::Array1;
 use std::fmt::Debug;
 
+/// Normalisation of periodogram across passbands
+#[derive(Clone, Debug, Serialize, Deserialize, JsonSchema)]
+pub enum MultiColorPeriodogramNormalisation {
+    /// Weight individual periodograms by the number of observations in each passband.
+    /// Useful if no weight is given to observations
+    Count,
+    /// Weight individual periodograms by $\chi^2 = \sum \left(\frac{m_i - \bar{m}}{\delta_i}\right)^2$
+    ///
+    /// Be aware that if no weight are given to observations
+    /// (i.e. via [TimeSeries::new_without_weight]) unity weights are assumed and this is NOT
+    /// equivalent to [::Count], but weighting by magnitude variance.
+    Chi2,
+}
+
 /// Multi-passband periodogram
 #[derive(Clone, Debug, Serialize, Deserialize, JsonSchema)]
 #[serde(
@@ -25,14 +39,26 @@ where
 {
     // We use it to not reimplement some internals
     monochrome: Periodogram<T, F>,
-    properties: Box<EvaluatorProperties>,
+    normalization: MultiColorPeriodogramNormalisation,
 }
 
 impl<T, F> MultiColorPeriodogram<T, F>
 where
     T: Float,
-    F: FeatureEvaluator<T>,
+    F: FeatureEvaluator<T> + From<PeriodogramPeaks>,
 {
+    pub fn new(peaks: usize, normalization: MultiColorPeriodogramNormalisation) -> Self {
+        let monochrome = Periodogram::with_name_description(
+            peaks,
+            "multicolor_periodogram",
+            "of multi-color periodogram (interpreting frequency as time, power as magnitude)",
+        );
+        Self {
+            monochrome,
+            normalization,
+        }
+    }
+
     #[inline]
     pub fn default_peaks() -> usize {
         PeriodogramPeaks::default_peaks()
@@ -103,19 +129,53 @@ where
         'a: 'mcts,
         P: PassbandTrait,
     {
-        let unnormed_power = mcts
-            .mapping_mut()
+        let ts_weights = {
+            let mut a: Array1<_> = match self.normalization {
+                MultiColorPeriodogramNormalisation::Count => {
+                    mcts.mapping_mut().values().map(|ts| ts.lenf()).collect()
+                }
+                MultiColorPeriodogramNormalisation::Chi2 => mcts
+                    .mapping_mut()
+                    .values_mut()
+                    .map(|ts| ts.get_m_chi2())
+                    .collect(),
+            };
+            let norm = a.sum();
+            if norm.is_zero() {
+                match self.normalization {
+                    MultiColorPeriodogramNormalisation::Count => {
+                        return Err(MultiColorEvaluatorError::all_time_series_short(
+                            mcts.mapping_mut(),
+                            self.min_ts_length(),
+                        ));
+                    }
+                    MultiColorPeriodogramNormalisation::Chi2 => {
+                        return Err(MultiColorEvaluatorError::AllTimeSeriesAreFlat);
+                    }
+                }
+            }
+            a /= norm;
+            a
+        };
+        mcts.mapping_mut()
             .values_mut()
-            .filter(|ts| self.monochrome.check_ts_length(ts).is_ok())
-            .map(|ts| p.power(ts) * ts.lenf())
-            .reduce(|acc, x| acc + x)
+            .zip(ts_weights.iter())
+            .filter(|(ts, _ts_weight)| self.monochrome.check_ts_length(ts).is_ok())
+            .map(|(ts, &ts_weight)| {
+                let mut power = p.power(ts);
+                power *= ts_weight;
+                power
+            })
+            .reduce(|mut acc, power| {
+                acc += &power;
+                acc
+            })
             .ok_or_else(|| {
                 MultiColorEvaluatorError::all_time_series_short(
                     mcts.mapping_mut(),
-                    self.monochrome.min_ts_length(),
+                    self.min_ts_length(),
                 )
-            })?;
-        Ok(unnormed_power / mcts.total_lenf())
+            })
     }
 
     pub fn power<'slf, 'a, 'mcts, P>(
@@ -174,7 +234,7 @@ where
     <F as TryInto<PeriodogramPeaks>>::Error: Debug,
 {
     fn get_info(&self) -> &EvaluatorInfo {
-        &self.properties.info
+        self.monochrome.get_info()
     }
 }
 
@@ -185,15 +245,11 @@ where
     <F as TryInto<PeriodogramPeaks>>::Error: Debug,
 {
     fn get_names(&self) -> Vec<&str> {
-        self.properties.names.iter().map(String::as_str).collect()
+        self.monochrome.get_names()
     }
 
     fn get_descriptions(&self) -> Vec<&str> {
-        self.properties
-            .descriptions
-            .iter()
-            .map(String::as_str)
-            .collect()
+        self.monochrome.get_descriptions()
     }
 }