diff --git a/reCTBN/Cargo.toml b/reCTBN/Cargo.toml index b0a691b..4749b23 100644 --- a/reCTBN/Cargo.toml +++ b/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" } diff --git a/reCTBN/src/parameter_learning.rs b/reCTBN/src/parameter_learning.rs index 3f505f9..7e45e58 100644 --- a/reCTBN/src/parameter_learning.rs +++ b/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( &self, net: &T, - dataset: &tools::Dataset, + dataset: &Dataset, node: usize, parent_set: Option>, ) -> Params; @@ -19,7 +19,7 @@ pub trait ParameterLearning { pub fn sufficient_statistics( net: &T, - dataset: &tools::Dataset, + dataset: &Dataset, node: usize, parent_set: &BTreeSet, ) -> (Array3, Array2) { @@ -76,7 +76,7 @@ impl ParameterLearning for MLE { fn fit( &self, net: &T, - dataset: &tools::Dataset, + dataset: &Dataset, node: usize, parent_set: Option>, ) -> Params { @@ -123,7 +123,7 @@ impl ParameterLearning for BayesianApproach { fn fit( &self, net: &T, - dataset: &tools::Dataset, + dataset: &Dataset, node: usize, parent_set: Option>, ) -> Params { @@ -170,26 +170,3 @@ impl ParameterLearning for BayesianApproach { return n; } } - -pub struct Cache { - parameter_learning: P, - dataset: tools::Dataset, -} - -impl Cache

{ - pub fn new(parameter_learning: P, dataset: tools::Dataset) -> Cache

{ - Cache { - parameter_learning, - dataset, - } - } - pub fn fit( - &mut self, - net: &T, - node: usize, - parent_set: Option>, - ) -> Params { - self.parameter_learning - .fit(net, &self.dataset, node, parent_set) - } -} diff --git a/reCTBN/src/process.rs b/reCTBN/src/process.rs index dc297bc..45c5e0a 100644 --- a/reCTBN/src/process.rs +++ b/reCTBN/src/process.rs @@ -21,7 +21,7 @@ pub type NetworkProcessState = Vec; /// 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; /// Add an **directed edge** between a two nodes of the network. diff --git a/reCTBN/src/structure_learning.rs b/reCTBN/src/structure_learning.rs index b272e22..a4c6ea1 100644 --- a/reCTBN/src/structure_learning.rs +++ b/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(&self, net: T, dataset: &tools::Dataset) -> T + fn fit_transform(&self, net: T, dataset: &Dataset) -> T where T: process::NetworkProcess; } diff --git a/reCTBN/src/structure_learning/constraint_based_algorithm.rs b/reCTBN/src/structure_learning/constraint_based_algorithm.rs index 670c8ed..f49b194 100644 --- a/reCTBN/src/structure_learning/constraint_based_algorithm.rs +++ b/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>, Params>, + cache_persistent_big: HashMap>, 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( + &mut self, + net: &T, + dataset: &Dataset, + node: usize, + parent_set: Option>, + ) -> 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 { + parameter_learning: P, + Ftest: F, + Chi2test: ChiSquare, +} + +impl CTPC

{ + pub fn new(parameter_learning: P, Ftest: F, Chi2test: ChiSquare) -> CTPC

{ + CTPC { + parameter_learning, + Ftest, + Chi2test, + } + } +} + +impl StructureLearningAlgorithm for CTPC

{ + fn fit_transform(&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)> = 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 = 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 + } +} diff --git a/reCTBN/src/structure_learning/hypothesis_test.rs b/reCTBN/src/structure_learning/hypothesis_test.rs index 344c995..4c02929 100644 --- a/reCTBN/src/structure_learning/hypothesis_test.rs +++ b/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( @@ -15,7 +16,8 @@ pub trait HypothesisTest { child_node: usize, parent_node: usize, separation_set: &BTreeSet, - cache: &mut parameter_learning::Cache

, + dataset: &Dataset, + cache: &mut Cache

, ) -> 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, + cim_1: &Array3, + j: usize, + M2: &Array3, + cim_2: &Array3, + ) -> 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( + &self, + net: &T, + child_node: usize, + parent_node: usize, + separation_set: &BTreeSet, + dataset: &Dataset, + cache: &mut Cache

, + ) -> 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, - cache: &mut parameter_learning::Cache

, + dataset: &Dataset, + cache: &mut Cache

, ) -> 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 diff --git a/reCTBN/src/structure_learning/score_based_algorithm.rs b/reCTBN/src/structure_learning/score_based_algorithm.rs index 16e9056..d65ea88 100644 --- a/reCTBN/src/structure_learning/score_based_algorithm.rs +++ b/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 { score_function: S, @@ -21,7 +21,7 @@ impl HillClimbing { } impl StructureLearningAlgorithm for HillClimbing { - fn fit_transform(&self, net: T, dataset: &tools::Dataset) -> T + fn fit_transform(&self, net: T, dataset: &Dataset) -> T where T: process::NetworkProcess, { diff --git a/reCTBN/src/tools.rs b/reCTBN/src/tools.rs index 38ebd49..1fdb661 100644 --- a/reCTBN/src/tools.rs +++ b/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, events: Array2, @@ -29,6 +30,7 @@ impl Trajectory { } } +#[derive(Clone)] pub struct Dataset { trajectories: Vec, } diff --git a/reCTBN/tests/structure_learning.rs b/reCTBN/tests/structure_learning.rs index 2ec64b2..9a69b45 100644 --- a/reCTBN/tests/structure_learning.rs +++ b/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