reCTBN/Cargo.toml

@@ -13,6 +13,8 @@ bimap = "~0.6"
 enum_dispatch = "~0.3"
 statrs = "~0.16"
 rand_chacha = "~0.3"
+itertools = "~0.10"
+rayon = "~1.6"
 
 [dev-dependencies]
 approx = { package = "approx", version = "~0.5" }

reCTBN/src/parameter_learning.rs

@@ -5,13 +5,13 @@ use std::collections::BTreeSet;
 use ndarray::prelude::*;
 
 use crate::params::*;
-use crate::{process, tools};
+use crate::{process, tools::Dataset};
 
-pub trait ParameterLearning {
+pub trait ParameterLearning: Sync {
     fn fit<T: process::NetworkProcess>(
         &self,
         net: &T,
-        dataset: &tools::Dataset,
+        dataset: &Dataset,
         node: usize,
         parent_set: Option<BTreeSet<usize>>,
     ) -> Params;
@@ -19,7 +19,7 @@ pub trait ParameterLearning {
 
 pub fn sufficient_statistics<T: process::NetworkProcess>(
     net: &T,
-    dataset: &tools::Dataset,
+    dataset: &Dataset,
     node: usize,
     parent_set: &BTreeSet<usize>,
 ) -> (Array3<usize>, Array2<f64>) {
@@ -76,7 +76,7 @@ impl ParameterLearning for MLE {
     fn fit<T: process::NetworkProcess>(
         &self,
         net: &T,
-        dataset: &tools::Dataset,
+        dataset: &Dataset,
         node: usize,
         parent_set: Option<BTreeSet<usize>>,
     ) -> Params {
@@ -123,7 +123,7 @@ impl ParameterLearning for BayesianApproach {
     fn fit<T: process::NetworkProcess>(
         &self,
         net: &T,
-        dataset: &tools::Dataset,
+        dataset: &Dataset,
         node: usize,
         parent_set: Option<BTreeSet<usize>>,
     ) -> Params {
@@ -164,26 +164,3 @@ impl ParameterLearning for BayesianApproach {
         return n;
     }
 }
-
-pub struct Cache<P: ParameterLearning> {
-    parameter_learning: P,
-    dataset: tools::Dataset,
-}
-
-impl<P: ParameterLearning> Cache<P> {
-    pub fn new(parameter_learning: P, dataset: tools::Dataset) -> Cache<P> {
-        Cache {
-            parameter_learning,
-            dataset,
-        }
-    }
-    pub fn fit<T: process::NetworkProcess>(
-        &mut self,
-        net: &T,
-        node: usize,
-        parent_set: Option<BTreeSet<usize>>,
-    ) -> Params {
-        self.parameter_learning
-            .fit(net, &self.dataset, node, parent_set)
-    }
-}

reCTBN/src/process.rs

@@ -21,7 +21,7 @@ pub type NetworkProcessState = Vec<params::StateType>;
 
 /// It defines the required methods for a structure used as a Probabilistic Graphical Models (such
 /// as a CTBN).
-pub trait NetworkProcess {
+pub trait NetworkProcess: Sync {
     fn initialize_adj_matrix(&mut self);
     fn add_node(&mut self, n: params::Params) -> Result<usize, NetworkError>;
     /// Add an **directed edge** between a two nodes of the network.

reCTBN/src/structure_learning.rs

@@ -4,10 +4,10 @@ pub mod constraint_based_algorithm;
 pub mod hypothesis_test;
 pub mod score_based_algorithm;
 pub mod score_function;
-use crate::{process, tools};
+use crate::{process, tools::Dataset};
 
 pub trait StructureLearningAlgorithm {
-    fn fit_transform<T>(&self, net: T, dataset: &tools::Dataset) -> T
+    fn fit_transform<T>(&self, net: T, dataset: &Dataset) -> T
     where
         T: process::NetworkProcess;
 }

reCTBN/src/structure_learning/constraint_based_algorithm.rs

@@ -1,5 +1,164 @@
 //! Module containing constraint based algorithms like CTPC and Hiton.
 
-//pub struct CTPC {
+use crate::params::Params;
-//
+use itertools::Itertools;
-//}
+use rayon::iter::{IntoParallelIterator, ParallelIterator};
+use rayon::prelude::ParallelExtend;
+use std::collections::{BTreeSet, HashMap};
+use std::mem;
+use std::usize;
+
+use super::hypothesis_test::*;
+use crate::parameter_learning::ParameterLearning;
+use crate::process;
+use crate::structure_learning::StructureLearningAlgorithm;
+use crate::tools::Dataset;
+
+pub struct Cache<'a, P: ParameterLearning> {
+    parameter_learning: &'a P,
+    cache_persistent_small: HashMap<Option<BTreeSet<usize>>, Params>,
+    cache_persistent_big: HashMap<Option<BTreeSet<usize>>, Params>,
+    parent_set_size_small: usize,
+}
+
+impl<'a, P: ParameterLearning> Cache<'a, P> {
+    pub fn new(parameter_learning: &'a P) -> Cache<'a, P> {
+        Cache {
+            parameter_learning,
+            cache_persistent_small: HashMap::new(),
+            cache_persistent_big: HashMap::new(),
+            parent_set_size_small: 0,
+        }
+    }
+    pub fn fit<T: process::NetworkProcess>(
+        &mut self,
+        net: &T,
+        dataset: &Dataset,
+        node: usize,
+        parent_set: Option<BTreeSet<usize>>,
+    ) -> Params {
+        let parent_set_len = parent_set.as_ref().unwrap().len();
+        if parent_set_len > self.parent_set_size_small + 1 {
+            //self.cache_persistent_small = self.cache_persistent_big;
+            mem::swap(
+                &mut self.cache_persistent_small,
+                &mut self.cache_persistent_big,
+            );
+            self.cache_persistent_big = HashMap::new();
+            self.parent_set_size_small += 1;
+        }
+
+        if parent_set_len > self.parent_set_size_small {
+            match self.cache_persistent_big.get(&parent_set) {
+                // TODO: Better not clone `params`, useless clock cycles, RAM use and I/O
+                // not cloning requires a minor and reasoned refactoring across the library
+                Some(params) => params.clone(),
+                None => {
+                    let params =
+                        self.parameter_learning
+                            .fit(net, dataset, node, parent_set.clone());
+                    self.cache_persistent_big.insert(parent_set, params.clone());
+                    params
+                }
+            }
+        } else {
+            match self.cache_persistent_small.get(&parent_set) {
+                // TODO: Better not clone `params`, useless clock cycles, RAM use and I/O
+                // not cloning requires a minor and reasoned refactoring across the library
+                Some(params) => params.clone(),
+                None => {
+                    let params =
+                        self.parameter_learning
+                            .fit(net, dataset, node, parent_set.clone());
+                    self.cache_persistent_small
+                        .insert(parent_set, params.clone());
+                    params
+                }
+            }
+        }
+    }
+}
+
+pub struct CTPC<P: ParameterLearning> {
+    parameter_learning: P,
+    Ftest: F,
+    Chi2test: ChiSquare,
+}
+
+impl<P: ParameterLearning> CTPC<P> {
+    pub fn new(parameter_learning: P, Ftest: F, Chi2test: ChiSquare) -> CTPC<P> {
+        CTPC {
+            parameter_learning,
+            Ftest,
+            Chi2test,
+        }
+    }
+}
+
+impl<P: ParameterLearning> StructureLearningAlgorithm for CTPC<P> {
+    fn fit_transform<T>(&self, net: T, dataset: &Dataset) -> T
+    where
+        T: process::NetworkProcess,
+    {
+        //Check the coherence between dataset and network
+        if net.get_number_of_nodes() != dataset.get_trajectories()[0].get_events().shape()[1] {
+            panic!("Dataset and Network must have the same number of variables.")
+        }
+
+        //Make the network mutable.
+        let mut net = net;
+
+        net.initialize_adj_matrix();
+
+        let mut learned_parent_sets: Vec<(usize, BTreeSet<usize>)> = vec![];
+        learned_parent_sets.par_extend(net.get_node_indices().into_par_iter().map(|child_node| {
+            let mut cache = Cache::new(&self.parameter_learning);
+            let mut candidate_parent_set: BTreeSet<usize> = net
+                .get_node_indices()
+                .into_iter()
+                .filter(|x| x != &child_node)
+                .collect();
+            let mut separation_set_size = 0;
+            while separation_set_size < candidate_parent_set.len() {
+                let mut candidate_parent_set_TMP = candidate_parent_set.clone();
+                for parent_node in candidate_parent_set.iter() {
+                    for separation_set in candidate_parent_set
+                        .iter()
+                        .filter(|x| x != &parent_node)
+                        .map(|x| *x)
+                        .combinations(separation_set_size)
+                    {
+                        let separation_set = separation_set.into_iter().collect();
+                        if self.Ftest.call(
+                            &net,
+                            child_node,
+                            *parent_node,
+                            &separation_set,
+                            dataset,
+                            &mut cache,
+                        ) && self.Chi2test.call(
+                            &net,
+                            child_node,
+                            *parent_node,
+                            &separation_set,
+                            dataset,
+                            &mut cache,
+                        ) {
+                            candidate_parent_set_TMP.remove(parent_node);
+                            break;
+                        }
+                    }
+                }
+                candidate_parent_set = candidate_parent_set_TMP;
+                separation_set_size += 1;
+            }
+            (child_node, candidate_parent_set)
+        }));
+        for (child_node, candidate_parent_set) in learned_parent_sets {
+            for parent_node in candidate_parent_set.iter() {
+                net.add_edge(*parent_node, child_node);
+            }
+        }
+        net
+    }
+}

reCTBN/src/structure_learning/hypothesis_test.rs

@@ -3,10 +3,11 @@
 use std::collections::BTreeSet;
 
 use ndarray::{Array3, Axis};
-use statrs::distribution::{ChiSquared, ContinuousCDF};
+use statrs::distribution::{ChiSquared, ContinuousCDF, FisherSnedecor};
 
 use crate::params::*;
-use crate::{parameter_learning, process};
+use crate::structure_learning::constraint_based_algorithm::Cache;
+use crate::{parameter_learning, process, tools::Dataset};
 
 pub trait HypothesisTest {
     fn call<T, P>(
@@ -15,7 +16,8 @@ pub trait HypothesisTest {
         child_node: usize,
         parent_node: usize,
         separation_set: &BTreeSet<usize>,
-        cache: &mut parameter_learning::Cache<P>,
+        dataset: &Dataset,
+        cache: &mut Cache<P>,
     ) -> bool
     where
         T: process::NetworkProcess,
@@ -37,7 +39,99 @@ pub struct ChiSquare {
     alpha: f64,
 }
 
-pub struct F {}
+pub struct F {
+    alpha: f64,
+}
+
+impl F {
+    pub fn new(alpha: f64) -> F {
+        F { alpha }
+    }
+
+    pub fn compare_matrices(
+        &self,
+        i: usize,
+        M1: &Array3<usize>,
+        cim_1: &Array3<f64>,
+        j: usize,
+        M2: &Array3<usize>,
+        cim_2: &Array3<f64>,
+    ) -> bool {
+        let M1 = M1.index_axis(Axis(0), i).mapv(|x| x as f64);
+        let M2 = M2.index_axis(Axis(0), j).mapv(|x| x as f64);
+        let cim_1 = cim_1.index_axis(Axis(0), i);
+        let cim_2 = cim_2.index_axis(Axis(0), j);
+        let r1 = M1.sum_axis(Axis(1));
+        let r2 = M2.sum_axis(Axis(1));
+        let q1 = cim_1.diag();
+        let q2 = cim_2.diag();
+        for idx in 0..r1.shape()[0] {
+            let s = q2[idx] / q1[idx];
+            let F = FisherSnedecor::new(r1[idx], r2[idx]).unwrap();
+            let s = F.cdf(s);
+            let lim_sx = self.alpha / 2.0;
+            let lim_dx = 1.0 - (self.alpha / 2.0);
+            if s < lim_sx || s > lim_dx {
+                return false;
+            }
+        }
+        true
+    }
+}
+
+impl HypothesisTest for F {
+    fn call<T, P>(
+        &self,
+        net: &T,
+        child_node: usize,
+        parent_node: usize,
+        separation_set: &BTreeSet<usize>,
+        dataset: &Dataset,
+        cache: &mut Cache<P>,
+    ) -> bool
+    where
+        T: process::NetworkProcess,
+        P: parameter_learning::ParameterLearning,
+    {
+        let P_small = match cache.fit(net, &dataset, child_node, Some(separation_set.clone())) {
+            Params::DiscreteStatesContinousTime(node) => node,
+        };
+        let mut extended_separation_set = separation_set.clone();
+        extended_separation_set.insert(parent_node);
+
+        let P_big = match cache.fit(
+            net,
+            &dataset,
+            child_node,
+            Some(extended_separation_set.clone()),
+        ) {
+            Params::DiscreteStatesContinousTime(node) => node,
+        };
+        let partial_cardinality_product: usize = extended_separation_set
+            .iter()
+            .take_while(|x| **x != parent_node)
+            .map(|x| net.get_node(*x).get_reserved_space_as_parent())
+            .product();
+        for idx_M_big in 0..P_big.get_transitions().as_ref().unwrap().shape()[0] {
+            let idx_M_small: usize = idx_M_big % partial_cardinality_product
+                + (idx_M_big
+                    / (partial_cardinality_product
+                        * net.get_node(parent_node).get_reserved_space_as_parent()))
+                    * partial_cardinality_product;
+            if !self.compare_matrices(
+                idx_M_small,
+                P_small.get_transitions().as_ref().unwrap(),
+                P_small.get_cim().as_ref().unwrap(),
+                idx_M_big,
+                P_big.get_transitions().as_ref().unwrap(),
+                P_big.get_cim().as_ref().unwrap(),
+            ) {
+                return false;
+            }
+        }
+        return true;
+    }
+}
 
 impl ChiSquare {
     pub fn new(alpha: f64) -> ChiSquare {
@@ -132,7 +226,8 @@ impl HypothesisTest for ChiSquare {
         child_node: usize,
         parent_node: usize,
         separation_set: &BTreeSet<usize>,
-        cache: &mut parameter_learning::Cache<P>,
+        dataset: &Dataset,
+        cache: &mut Cache<P>,
     ) -> bool
     where
         T: process::NetworkProcess,
@@ -141,14 +236,19 @@ impl HypothesisTest for ChiSquare {
         // Prendo dalla cache l'apprendimento dei parametri, che sarebbe una CIM
         // di dimensione nxn
         //  (CIM, M, T)
-        let P_small = match cache.fit(net, child_node, Some(separation_set.clone())) {
+        let P_small = match cache.fit(net, &dataset, child_node, Some(separation_set.clone())) {
             Params::DiscreteStatesContinousTime(node) => node,
         };
         //
         let mut extended_separation_set = separation_set.clone();
         extended_separation_set.insert(parent_node);
 
-        let P_big = match cache.fit(net, child_node, Some(extended_separation_set.clone())) {
+        let P_big = match cache.fit(
+            net,
+            &dataset,
+            child_node,
+            Some(extended_separation_set.clone()),
+        ) {
             Params::DiscreteStatesContinousTime(node) => node,
         };
         // Commentare qui

reCTBN/src/structure_learning/score_based_algorithm.rs

@@ -4,7 +4,7 @@ use std::collections::BTreeSet;
 
 use crate::structure_learning::score_function::ScoreFunction;
 use crate::structure_learning::StructureLearningAlgorithm;
-use crate::{process, tools};
+use crate::{process, tools::Dataset};
 
 pub struct HillClimbing<S: ScoreFunction> {
     score_function: S,
@@ -21,7 +21,7 @@ impl<S: ScoreFunction> HillClimbing<S> {
 }
 
 impl<S: ScoreFunction> StructureLearningAlgorithm for HillClimbing<S> {
-    fn fit_transform<T>(&self, net: T, dataset: &tools::Dataset) -> T
+    fn fit_transform<T>(&self, net: T, dataset: &Dataset) -> T
     where
         T: process::NetworkProcess,
     {

reCTBN/src/tools.rs

@@ -5,6 +5,7 @@ use ndarray::prelude::*;
 use crate::sampling::{ForwardSampler, Sampler};
 use crate::{params, process};
 
+#[derive(Clone)]
 pub struct Trajectory {
     time: Array1<f64>,
     events: Array2<usize>,
@@ -29,6 +30,7 @@ impl Trajectory {
     }
 }
 
+#[derive(Clone)]
 pub struct Dataset {
     trajectories: Vec<Trajectory>,
 }

reCTBN/tests/structure_learning.rs

@@ -7,9 +7,9 @@ use ndarray::{arr1, arr2, arr3};
 use reCTBN::process::ctbn::*;
 use reCTBN::process::NetworkProcess;
 use reCTBN::parameter_learning::BayesianApproach;
-use reCTBN::parameter_learning::Cache;
 use reCTBN::params;
 use reCTBN::structure_learning::hypothesis_test::*;
+use reCTBN::structure_learning::constraint_based_algorithm::*;
 use reCTBN::structure_learning::score_based_algorithm::*;
 use reCTBN::structure_learning::score_function::*;
 use reCTBN::structure_learning::StructureLearningAlgorithm;
@@ -108,7 +108,7 @@ fn check_compatibility_between_dataset_and_network<T: StructureLearningAlgorithm
         }
     }
 
-    let data = trajectory_generator(&net, 100, 20.0, Some(6347747169756259));
+    let data = trajectory_generator(&net, 100, 30.0, Some(6347747169756259));
 
     let mut net = CtbnNetwork::new();
     let _n1 = net
@@ -392,7 +392,7 @@ pub fn chi_square_compare_matrices() {
             [  700, 800,  0]
         ],
     ]);
-    let chi_sq = ChiSquare::new(0.1);
+    let chi_sq = ChiSquare::new(1e-4);
     assert!(!chi_sq.compare_matrices(i, &M1, j, &M2));
 }
 
@@ -422,7 +422,7 @@ pub fn chi_square_compare_matrices_2() {
         [ 400,  0,  600],
         [  700, 800,  0]]
     ]);
-    let chi_sq = ChiSquare::new(0.1);
+    let chi_sq = ChiSquare::new(1e-4);
     assert!(chi_sq.compare_matrices(i, &M1, j, &M2));
 }
 
@@ -454,7 +454,7 @@ pub fn chi_square_compare_matrices_3() {
             [  700, 800,  0]
         ],
     ]);
-    let chi_sq = ChiSquare::new(0.1);
+    let chi_sq = ChiSquare::new(1e-4);
     assert!(chi_sq.compare_matrices(i, &M1, j, &M2));
 }
 
@@ -466,13 +466,56 @@ pub fn chi_square_call() {
     let N3: usize = 2;
     let N2: usize = 1;
     let N1: usize = 0;
-    let separation_set = BTreeSet::new();
+    let mut separation_set = BTreeSet::new();
+    let parameter_learning = BayesianApproach { alpha: 1, tau:1.0 };
+    let mut cache = Cache::new(&parameter_learning);
+    let chi_sq = ChiSquare::new(1e-4);
+
+    assert!(chi_sq.call(&net, N1, N3, &separation_set, &data, &mut cache));
+    let mut cache = Cache::new(&parameter_learning);
+    assert!(!chi_sq.call(&net, N3, N1, &separation_set, &data, &mut cache));
+    assert!(!chi_sq.call(&net, N3, N2, &separation_set, &data, &mut cache));
+    separation_set.insert(N1);
+    let mut cache = Cache::new(&parameter_learning);
+    assert!(chi_sq.call(&net, N2, N3, &separation_set, &data, &mut cache));
+}
+
+#[test]
+pub fn f_call() {
+
+    let (net, data) = get_mixed_discrete_net_3_nodes_with_data();
+    let N3: usize = 2;
+    let N2: usize = 1;
+    let N1: usize = 0;
+    let mut separation_set = BTreeSet::new();
     let parameter_learning = BayesianApproach { alpha: 1, tau:1.0 };
-    let mut cache = Cache::new(parameter_learning, data);
+    let mut cache = Cache::new(&parameter_learning);
-    let chi_sq = ChiSquare::new(0.0001);
+    let f = F::new(1e-6);
 
-    assert!(chi_sq.call(&net, N1, N3, &separation_set, &mut cache));
+
-    assert!(!chi_sq.call(&net, N3, N1, &separation_set, &mut cache));
+    assert!(f.call(&net, N1, N3, &separation_set, &data, &mut cache));
-    assert!(!chi_sq.call(&net, N3, N2, &separation_set, &mut cache));
+    let mut cache = Cache::new(&parameter_learning);
-    assert!(chi_sq.call(&net, N2, N3, &separation_set, &mut cache));
+    assert!(!f.call(&net, N3, N1, &separation_set, &data, &mut cache));
+    assert!(!f.call(&net, N3, N2, &separation_set, &data, &mut cache));
+    separation_set.insert(N1);
+    let mut cache = Cache::new(&parameter_learning);
+    assert!(f.call(&net, N2, N3, &separation_set, &data, &mut cache));
+}
+
+#[test]
+pub fn learn_ternary_net_2_nodes_ctpc() {
+    let f = F::new(1e-6);
+    let chi_sq = ChiSquare::new(1e-4);
+    let parameter_learning = BayesianApproach { alpha: 1, tau:1.0 };
+    let ctpc = CTPC::new(parameter_learning, f, chi_sq);
+    learn_ternary_net_2_nodes(ctpc);
+}
+
+#[test]
+fn learn_mixed_discrete_net_3_nodes_ctpc() {
+    let f = F::new(1e-6);
+    let chi_sq = ChiSquare::new(1e-4);
+    let parameter_learning = BayesianApproach { alpha: 1, tau:1.0 };
+    let ctpc = CTPC::new(parameter_learning, f, chi_sq);
+    learn_mixed_discrete_net_3_nodes(ctpc);
 }