Repository: mlondschien/changeforest (language: rust)

Issue:
Add `ModelSelectionResult`
This could have an entry `value`, corresponding to p-value for classifier based and rel. max gain for the parametric setting. This would allow pruning of the tree / result.


Discussion / hints:

Fix patch:
diff --git a/hdcd-python/src/result.rs b/hdcd-python/src/result.rs
index 02d149b..7039b0b 100644
--- a/hdcd-python/src/result.rs
+++ b/hdcd-python/src/result.rs
@@ -3,10 +3,36 @@
 
 use hdcd::gain::GainResult;
 use hdcd::optimizer::OptimizerResult;
-use hdcd::BinarySegmentationResult;
+use hdcd::{BinarySegmentationResult, ModelSelectionResult};
 use numpy::{PyArray1, PyArray2, ToPyArray};
 use pyo3::prelude::*;
 
+#[pyclass]
+#[derive(Clone, Debug)]
+pub struct MyModelSelectionResult {
+    pub result: ModelSelectionResult,
+}
+
+#[pymethods]
+impl MyModelSelectionResult {
+    #[getter]
+    pub fn is_significant(&self) -> bool {
+        self.result.is_significant
+    }
+
+    #[getter]
+    pub fn p_value(&self) -> Option<f64> {
+        self.result.p_value
+    }
+}
+
+#[pyproto]
+impl pyo3::class::basic::PyObjectProtocol for MyModelSelectionResult {
+    fn __repr__(&self) -> PyResult<String> {
+        Ok(format!("{}", self.result))
+    }
+}
+
 #[pyclass]
 #[derive(Clone, Debug)]
 pub struct MyGainResult {
@@ -143,6 +169,13 @@ impl MyBinarySegmentationResult {
         })
     }
 
+    #[getter]
+    fn model_selection_result(&self) -> MyModelSelectionResult {
+        MyModelSelectionResult {
+            result: self.result.model_selection_result.clone(),
+        }
+    }
+
     #[getter]
     fn left(&self) -> Option<Self> {
         // This seems overly complicated. We turn Option<Box<Self>> into
diff --git a/hdcd-r/R/display.R b/hdcd-r/R/display.R
index fde4459..f4b39da 100644
--- a/hdcd-r/R/display.R
+++ b/hdcd-r/R/display.R
@@ -27,7 +27,8 @@ to_data_frame = function(node, parent_indent="", prefix="") {
             name=name,
             best_split=node$best_split,
             max_gain=node$max_gain,
-            is_significant=node$is_significant
+            p_value=node$model_selection_result$p_value,
+            is_significant=node$model_selection_result$is_significant
         ),
         left,
         right
diff --git a/hdcd-r/R/hdcd.R b/hdcd-r/R/hdcd.R
index ec6d769..8d72a0f 100644
--- a/hdcd-r/R/hdcd.R
+++ b/hdcd-r/R/hdcd.R
@@ -5,7 +5,7 @@ BinarySegmentationResult = R6::R6Class(
         stop = NULL,
         best_split = NULL,
         max_gain = NULL,
-        is_significant = NULL,
+        model_selection_result = NULL,
         gain_results = NULL,
         segments = NULL,
         left = NULL,
@@ -16,7 +16,7 @@ BinarySegmentationResult = R6::R6Class(
             stop = NULL,
             best_split = NULL,
             max_gain = NULL,
-            is_significant = NULL,
+            model_selection_result = NULL,
             gain_results = NULL,
             segments = NULL,
             left = NULL,
@@ -26,7 +26,7 @@ BinarySegmentationResult = R6::R6Class(
             self$stop = stop
             self$best_split = best_split
             self$max_gain = max_gain
-            self$is_significant = is_significant
+            self$model_selection_result = model_selection_result
             self$gain_results = gain_results
             self$segments = segments
             self$left = left
@@ -47,7 +47,7 @@ BinarySegmentationResult = R6::R6Class(
                 split_points = append(split_points, self$left$split_points())
             }
 
-            if (!is.null(self$best_split) & self$is_significant) {
+            if (!is.null(self$best_split) & self$model_selection_result$is_significant) {
                 split_points = append(split_points, c(self$best_split))
             }
 
diff --git a/hdcd-r/src/rust/src/result.rs b/hdcd-r/src/rust/src/result.rs
index 7fddf8d..b68882e 100644
--- a/hdcd-r/src/rust/src/result.rs
+++ b/hdcd-r/src/rust/src/result.rs
@@ -2,7 +2,25 @@ use extendr_api::prelude::*;
 //use extendr_api::robj_ndarray::TryFrom;
 use hdcd::gain::GainResult;
 use hdcd::optimizer::OptimizerResult;
-use hdcd::BinarySegmentationResult;
+use hdcd::{BinarySegmentationResult, ModelSelectionResult};
+
+pub struct MyModelSelectionResult {
+    model_selection_result: ModelSelectionResult,
+}
+
+impl From<MyModelSelectionResult> for Robj {
+    fn from(my_model_selection_result: MyModelSelectionResult) -> Self {
+        List::from_values(&[
+            r!(my_model_selection_result
+                .model_selection_result
+                .is_significant),
+            r!(my_model_selection_result.model_selection_result.p_value),
+        ])
+        .into_robj()
+        .set_names(&["is_significant", "p_value"])
+        .expect("From<ModelSelectionResult> failed")
+    }
+}
 
 pub struct MyGainResult {
     gain_result: GainResult,
@@ -116,12 +134,17 @@ impl From<MyBinarySegmentationResult> for Robj {
             None => ().into(),
         };
 
+        let model_selection_result: Robj = MyModelSelectionResult {
+            model_selection_result: my_result.result.model_selection_result,
+        }
+        .into();
+
         List::from_values(&[
             r!(my_result.result.start as i32),
             r!(my_result.result.stop as i32),
             r!(my_result.result.best_split.map(|u| u as i32)),
             r!(my_result.result.max_gain),
-            r!(my_result.result.is_significant),
+            r!(model_selection_result),
             r!(gain_results),
             r!(segments),
             r!(left),
@@ -133,7 +156,7 @@ impl From<MyBinarySegmentationResult> for Robj {
             "stop",
             "best_split",
             "max_gain",
-            "is_significant",
+            "model_selection_result",
             "gain_results",
             "segments",
             "left",
diff --git a/src/binary_segmentation.rs b/src/binary_segmentation.rs
index ad4d757..cb3338c 100644
--- a/src/binary_segmentation.rs
+++ b/src/binary_segmentation.rs
@@ -1,6 +1,6 @@
 use crate::gain::GainResult;
 use crate::optimizer::OptimizerResult;
-use crate::Segmentation;
+use crate::{ModelSelectionResult, Segmentation};
 
 pub struct BinarySegmentationTree {
     pub start: usize,
@@ -8,7 +8,7 @@ pub struct BinarySegmentationTree {
     pub n: usize,
     pub split: Option<usize>,
     pub max_gain: Option<f64>,
-    pub is_significant: bool,
+    pub model_selection_result: ModelSelectionResult,
     pub left: Option<Box<BinarySegmentationTree>>,
     pub right: Option<Box<BinarySegmentationTree>>,
     pub optimizer_result: Option<OptimizerResult>,
@@ -20,12 +20,12 @@ impl BinarySegmentationTree {
             start: 0,
             stop: X.nrows(),
             n: X.nrows(),
-            split: Option::None,
-            max_gain: Option::None,
-            is_significant: false,
-            left: Option::None,
-            right: Option::None,
-            optimizer_result: Option::None,
+            split: None,
+            max_gain: None,
+            model_selection_result: ModelSelectionResult::default(),
+            left: None,
+            right: None,
+            optimizer_result: None,
         }
     }
 
@@ -34,12 +34,12 @@ impl BinarySegmentationTree {
             start: self.start,
             stop: split,
             n: self.n,
-            split: Option::None,
-            max_gain: Option::None,
-            is_significant: false,
-            left: Option::None,
-            right: Option::None,
-            optimizer_result: Option::None,
+            split: None,
+            max_gain: None,
+            model_selection_result: ModelSelectionResult::default(),
+            left: None,
+            right: None,
+            optimizer_result: None,
         })
     }
 
@@ -48,12 +48,12 @@ impl BinarySegmentationTree {
             start: split,
             stop: self.stop,
             n: self.n,
-            split: Option::None,
-            max_gain: Option::None,
-            is_significant: false,
-            left: Option::None,
-            right: Option::None,
-            optimizer_result: Option::None,
+            split: None,
+            max_gain: None,
+            model_selection_result: ModelSelectionResult::default(),
+            left: None,
+            right: None,
+            optimizer_result: None,
         })
     }
 
@@ -68,9 +68,9 @@ impl BinarySegmentationTree {
             self.split = Some(optimizer_result.best_split);
             self.max_gain = Some(optimizer_result.max_gain);
 
-            self.is_significant = segmentation.is_significant(&optimizer_result);
+            self.model_selection_result = segmentation.model_selection(&optimizer_result);
 
-            if self.is_significant {
+            if self.model_selection_result.is_significant {
                 let mut left = self.new_left(optimizer_result.best_split);
                 left.grow(segmentation);
                 self.left = Some(left);
@@ -92,7 +92,7 @@ pub struct BinarySegmentationResult {
     pub stop: usize,
     pub best_split: Option<usize>,
     pub max_gain: Option<f64>,
-    pub is_significant: bool,
+    pub model_selection_result: ModelSelectionResult,
     pub gain_results: Option<Vec<GainResult>>,
     pub left: Option<Box<BinarySegmentationResult>>,
     pub right: Option<Box<BinarySegmentationResult>>,
@@ -116,7 +116,7 @@ impl BinarySegmentationResult {
             stop: tree.stop,
             best_split: tree.split,
             max_gain: tree.max_gain,
-            is_significant: tree.is_significant,
+            model_selection_result: tree.model_selection_result,
             gain_results,
             left,
             right,
@@ -132,7 +132,7 @@ impl BinarySegmentationResult {
         }
 
         if let Some(best_split) = self.best_split {
-            if self.is_significant {
+            if self.model_selection_result.is_significant {
                 split_points.push(best_split);
             }
         }
@@ -201,7 +201,7 @@ mod tests {
         assert_eq!(result.start, 0);
         assert_eq!(result.stop, 100);
         assert_eq!(result.best_split, Some(25));
-        assert!(result.is_significant);
+        assert!(result.model_selection_result.is_significant);
         assert!(result.gain_results.is_some());
 
         let right = result.right.as_ref().unwrap();
@@ -209,7 +209,7 @@ mod tests {
         assert_eq!(right.start, 25);
         assert_eq!(right.stop, 100);
         assert_eq!(right.best_split, Some(40));
-        assert!(right.is_significant);
+        assert!(right.model_selection_result.is_significant);
         assert!(right.gain_results.is_some());
 
         let left = result.left.as_ref().unwrap();
@@ -217,7 +217,7 @@ mod tests {
         assert_eq!(left.start, 0);
         assert_eq!(left.stop, 25);
         assert_eq!(left.best_split, Some(10));
-        assert!(!left.is_significant);
+        assert!(!left.model_selection_result.is_significant);
         assert!(left.gain_results.is_some()); // even though is_significant is false
 
         let result = result.with_segments(segmentation);
diff --git a/src/gain/change_in_mean.rs b/src/gain/change_in_mean.rs
index 0b97ea5..40b51f0 100644
--- a/src/gain/change_in_mean.rs
+++ b/src/gain/change_in_mean.rs
@@ -1,5 +1,6 @@
 use crate::control::Control;
 use crate::gain::{Gain, GainResult};
+use crate::ModelSelectionResult;
 use std::cell::{Ref, RefCell};
 
 pub struct ChangeInMean<'a, 'b> {
@@ -60,8 +61,11 @@ impl<'a, 'b> Gain for ChangeInMean<'a, 'b> {
         result / (s * s_1 * s_2)
     }
 
-    fn is_significant(&self, max_gain: f64, _: &GainResult) -> bool {
-        max_gain > self.control.minimal_gain_to_split * (self.n() as f64)
+    fn model_selection(&self, max_gain: f64, _: &GainResult) -> ModelSelectionResult {
+        ModelSelectionResult {
+            is_significant: max_gain > self.control.minimal_gain_to_split * (self.n() as f64),
+            p_value: None,
+        }
     }
 
     fn control(&self) -> &Control {
diff --git a/src/gain/classifier_gain.rs b/src/gain/classifier_gain.rs
index c054264..0e38356 100644
--- a/src/gain/classifier_gain.rs
+++ b/src/gain/classifier_gain.rs
@@ -1,6 +1,7 @@
 use crate::control::Control;
 use crate::gain::{ApproxGain, ApproxGainResult, Gain, GainResult};
 use crate::Classifier;
+use crate::ModelSelectionResult;
 use ndarray::{s, Array1, Array2, Axis};
 use rand::{rngs::StdRng, SeedableRng};
 
@@ -25,7 +26,7 @@ where
             .single_likelihood(&predictions, start, stop, split)
     }
 
-    fn is_significant(&self, _: f64, gain_result: &GainResult) -> bool {
+    fn model_selection(&self, _: f64, gain_result: &GainResult) -> ModelSelectionResult {
         let likelihoods: &Array2<f64>;
         let start: usize;
         let stop: usize;
@@ -53,8 +54,6 @@ where
             }
         }
 
-        // assert_eq!(max_gain, max_gain_);
-
         let mut p_value: u32 = 1;
 
         for _ in 0..n_permutations {
@@ -67,7 +66,13 @@ where
                 }
             }
         }
-        (p_value as f64 / (n_permutations + 1) as f64) < self.control().model_selection_alpha
+
+        let p_value = p_value as f64 / (n_permutations + 1) as f64;
+        let is_significant = p_value < self.control().model_selection_alpha;
+        ModelSelectionResult {
+            is_significant,
+            p_value: Some(p_value),
+        }
     }
 
     fn control(&self) -> &Control {
diff --git a/src/gain/gain.rs b/src/gain/gain.rs
index c1d3236..ec5507d 100644
--- a/src/gain/gain.rs
+++ b/src/gain/gain.rs
@@ -1,5 +1,6 @@
 use crate::control::Control;
 use crate::gain::{ApproxGainResult, FullGainResult, GainResult};
+use crate::ModelSelectionResult;
 
 pub trait Gain {
     /// Total number of observations.
@@ -34,7 +35,7 @@ pub trait Gain {
     }
 
     /// Does a certain split corresponds to a true change point?
-    fn is_significant(&self, max_gain: f64, gain_result: &GainResult) -> bool;
+    fn model_selection(&self, max_gain: f64, gain_result: &GainResult) -> ModelSelectionResult;
 
     /// Hyperparameters.
     fn control(&self) -> &Control;
diff --git a/src/lib.rs b/src/lib.rs
index 630a3aa..a0bf7c2 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -8,6 +8,7 @@ mod binary_segmentation;
 pub mod classifier;
 mod control;
 pub mod gain;
+mod model_selection_result;
 pub mod optimizer;
 mod segmentation;
 
@@ -15,6 +16,7 @@ pub use binary_segmentation::{BinarySegmentationResult, BinarySegmentationTree};
 pub use classifier::Classifier;
 pub use control::Control;
 pub use gain::{ClassifierGain, Gain};
+pub use model_selection_result::ModelSelectionResult;
 pub use optimizer::Optimizer;
 pub use segmentation::{Segmentation, SegmentationType};
 pub mod utils;
diff --git a/src/model_selection_result.rs b/src/model_selection_result.rs
new file mode 100644
index 0000000..9c355df
--- /dev/null
+++ b/src/model_selection_result.rs
@@ -0,0 +1,28 @@
+use std::fmt;
+
+#[derive(Clone, Debug)]
+pub struct ModelSelectionResult {
+    pub is_significant: bool,
+    pub p_value: Option<f64>,
+}
+
+// https://doc.rust-lang.org/rust-by-example/hello/print/print_display.html
+impl fmt::Display for ModelSelectionResult {
+    // This trait requires `fmt` with this exact signature.
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(
+            f,
+            "ModelSelectionResult(is_significant={}, p_value={:?})",
+            self.is_significant, self.p_value
+        )
+    }
+}
+
+impl Default for ModelSelectionResult {
+    fn default() -> Self {
+        ModelSelectionResult {
+            is_significant: false,
+            p_value: None,
+        }
+    }
+}
diff --git a/src/optimizer/grid_search.rs b/src/optimizer/grid_search.rs
index f688802..9fb0251 100644
--- a/src/optimizer/grid_search.rs
+++ b/src/optimizer/grid_search.rs
@@ -1,6 +1,6 @@
 use crate::gain::GainResult;
 use crate::optimizer::OptimizerResult;
-use crate::{Control, Gain, Optimizer};
+use crate::{Control, Gain, ModelSelectionResult, Optimizer};
 
 pub struct GridSearch<T: Gain> {
     pub gain: T,
@@ -46,10 +46,10 @@ where
         })
     }
 
-    fn is_significant(&self, optimizer_result: &OptimizerResult) -> bool {
+    fn model_selection(&self, optimizer_result: &OptimizerResult) -> ModelSelectionResult {
         let gain_result = optimizer_result.gain_results.last().unwrap();
         self.gain
-            .is_significant(optimizer_result.max_gain, gain_result)
+            .model_selection(optimizer_result.max_gain, gain_result)
     }
 }
 
diff --git a/src/optimizer/optimizer.rs b/src/optimizer/optimizer.rs
index 35a0ff9..1235a3d 100644
--- a/src/optimizer/optimizer.rs
+++ b/src/optimizer/optimizer.rs
@@ -1,5 +1,5 @@
 use crate::optimizer::OptimizerResult;
-use crate::Control;
+use crate::{Control, ModelSelectionResult};
 
 pub trait Optimizer {
     /// Find the element of `split_candidates` to split segment `[start, stop)`.
@@ -8,7 +8,7 @@ pub trait Optimizer {
     fn find_best_split(&self, start: usize, stop: usize) -> Result<OptimizerResult, &str>;
 
     /// Does a certain split corresponds to a true change point?
-    fn is_significant(&self, optimizer_result: &OptimizerResult) -> bool;
+    fn model_selection(&self, optimizer_result: &OptimizerResult) -> ModelSelectionResult;
 
     /// Total number of observations.
     fn n(&self) -> usize;
diff --git a/src/optimizer/two_step_search.rs b/src/optimizer/two_step_search.rs
index 80e937c..c09b02d 100644
--- a/src/optimizer/two_step_search.rs
+++ b/src/optimizer/two_step_search.rs
@@ -1,6 +1,6 @@
 use crate::gain::{ApproxGain, GainResult};
 use crate::optimizer::OptimizerResult;
-use crate::{Control, Gain, Optimizer};
+use crate::{Control, Gain, ModelSelectionResult, Optimizer};
 
 pub struct TwoStepSearch<T: Gain> {
     pub gain: T,
@@ -62,10 +62,10 @@ where
         })
     }
 
-    fn is_significant(&self, optimizer_result: &OptimizerResult) -> bool {
+    fn model_selection(&self, optimizer_result: &OptimizerResult) -> ModelSelectionResult {
         let gain_result = optimizer_result.gain_results.first().unwrap();
         self.gain
-            .is_significant(optimizer_result.max_gain, gain_result)
+            .model_selection(optimizer_result.max_gain, gain_result)
     }
 }
 
diff --git a/src/segmentation.rs b/src/segmentation.rs
index 290a3bd..c031a6a 100644
--- a/src/segmentation.rs
+++ b/src/segmentation.rs
@@ -1,11 +1,11 @@
 use crate::optimizer::OptimizerResult;
+use crate::ModelSelectionResult;
 use crate::Optimizer;
 use rand::{
     distributions::{Distribution, Uniform},
     rngs::StdRng,
     SeedableRng,
 };
-
 pub enum SegmentationType {
     BS,
     WBS,
@@ -112,8 +112,8 @@ impl<'a> Segmentation<'a> {
         }
     }
 
-    pub fn is_significant(&self, optimizer_result: &OptimizerResult) -> bool {
-        self.optimizer.is_significant(optimizer_result)
+    pub fn model_selection(&self, optimizer_result: &OptimizerResult) -> ModelSelectionResult {
+        self.optimizer.model_selection(optimizer_result)
     }
 }
 


Test patch:
diff --git a/hdcd-r/tests/testthat/test-display.R b/hdcd-r/tests/testthat/test-display.R
index a7fdf78..0458ab0 100644
--- a/hdcd-r/tests/testthat/test-display.R
+++ b/hdcd-r/tests/testthat/test-display.R
@@ -12,18 +12,18 @@ test_that("display", {
     output = capture.output(print(result))
 
     expected = c(
-        "                  name best_split  max_gain is_significant",
-        "1  (0, 10]                      5  5.379459           TRUE",
-        "2   ¦--(0, 5]                   2  1.383814           TRUE",
-        "3   ¦    ¦--(0, 2]             NA        NA          FALSE",
-        "4   ¦    °--(2, 5]              3  0.000000           TRUE",
-        "5   ¦        ¦--(2, 3]         NA        NA          FALSE",
-        "6   ¦        °--(3, 5]         NA        NA          FALSE",
-        "7   °--(5, 10]                  7 -4.616186           TRUE",
-        "8       ¦--(5, 7]              NA        NA          FALSE",
-        "9       °--(7, 10]              8  0.000000           TRUE",
-        "10          ¦--(7, 8]          NA        NA          FALSE",
-        "11          °--(8, 10]         NA        NA          FALSE"
+        "                  name best_split  max_gain p_value is_significant",
+        "1  (0, 10]                      5  5.379459    0.03           TRUE",
+        "2   ¦--(0, 5]                   2  1.383814    0.01           TRUE",
+        "3   ¦    ¦--(0, 2]             NA        NA      NA          FALSE",
+        "4   ¦    °--(2, 5]              3  0.000000    0.01           TRUE",
+        "5   ¦        ¦--(2, 3]         NA        NA      NA          FALSE",
+        "6   ¦        °--(3, 5]         NA        NA      NA          FALSE",
+        "7   °--(5, 10]                  7 -4.616186    0.01           TRUE",
+        "8       ¦--(5, 7]              NA        NA      NA          FALSE",
+        "9       °--(7, 10]              8  0.000000    0.01           TRUE",
+        "10          ¦--(7, 8]          NA        NA      NA          FALSE",
+        "11          °--(8, 10]         NA        NA      NA          FALSE"
     )
 
     expect_lists_equal(output, expected)
diff --git a/src/testing.rs b/src/testing.rs
index 038bc8e..5e90af4 100644
--- a/src/testing.rs
+++ b/src/testing.rs
@@ -1,6 +1,6 @@
 use crate::gain::{gain_from_likelihoods, ApproxGain, ApproxGainResult, Gain, GainResult};
 use crate::optimizer::OptimizerResult;
-use crate::{Control, Optimizer};
+use crate::{Control, ModelSelectionResult, Optimizer};
 use ndarray::{s, stack, Array, Array1, Array2, ArrayView2, Axis};
 use ndarray_rand::rand_distr::Uniform;
 use ndarray_rand::RandomExt;
@@ -40,8 +40,11 @@ impl<'a> Gain for ChangeInMean<'a> {
             - slice.sum_axis(Axis(0)).mapv(|a| a.powi(2)).sum() / n_slice
     }
 
-    fn is_significant(&self, max_gain: f64, _: &GainResult) -> bool {
-        max_gain > self.control.minimal_gain_to_split * (self.n() as f64)
+    fn model_selection(&self, max_gain: f64, _: &GainResult) -> ModelSelectionResult {
+        ModelSelectionResult {
+            is_significant: max_gain > self.control.minimal_gain_to_split * (self.n() as f64),
+            p_value: None,
+        }
     }
 
     fn control(&self) -> &Control {
@@ -117,9 +120,11 @@ impl<'a> Optimizer for TrivialOptimizer<'a> {
         })
     }
 
-    #[allow(unused_variables)]
-    fn is_significant(&self, optimizer_result: &OptimizerResult) -> bool {
-        optimizer_result.stop <= 50
+    fn model_selection(&self, optimizer_result: &OptimizerResult) -> ModelSelectionResult {
+        ModelSelectionResult {
+            is_significant: optimizer_result.stop <= 50,
+            p_value: None,
+        }
     }
 
     fn control(&self) -> &Control {