diff --git a/Cargo.toml b/Cargo.toml index 553e294..53c74f6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,18 +1,5 @@ -[package] -name = "reCTBN" -version = "0.1.0" -edition = "2021" +[workspace] -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -ndarray = {version="*", features=["approx"]} -thiserror = "*" -rand = "*" -bimap = "*" -enum_dispatch = "*" -statrs = "*" -rand_chacha = "*" - -[dev-dependencies] -approx = "*" +members = [ + "reCTBN", +] diff --git a/reCTBN/Cargo.toml b/reCTBN/Cargo.toml new file mode 100644 index 0000000..547a8b8 --- /dev/null +++ b/reCTBN/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "reCTBN" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +ndarray = {version="*", features=["approx-0_5"]} +thiserror = "*" +rand = "*" +bimap = "*" +enum_dispatch = "*" +statrs = "*" +rand_chacha = "*" + +[dev-dependencies] +approx = { package = "approx", version = "0.5" } diff --git a/src/ctbn.rs b/reCTBN/src/ctbn.rs similarity index 100% rename from src/ctbn.rs rename to reCTBN/src/ctbn.rs diff --git a/src/lib.rs b/reCTBN/src/lib.rs similarity index 90% rename from src/lib.rs rename to reCTBN/src/lib.rs index 8c57af2..280bd21 100644 --- a/src/lib.rs +++ b/reCTBN/src/lib.rs @@ -6,5 +6,6 @@ pub mod ctbn; pub mod network; pub mod parameter_learning; pub mod params; +pub mod sampling; pub mod structure_learning; pub mod tools; diff --git a/src/network.rs b/reCTBN/src/network.rs similarity index 100% rename from src/network.rs rename to reCTBN/src/network.rs diff --git a/src/parameter_learning.rs b/reCTBN/src/parameter_learning.rs similarity index 100% rename from src/parameter_learning.rs rename to reCTBN/src/parameter_learning.rs diff --git a/src/params.rs b/reCTBN/src/params.rs similarity index 100% rename from src/params.rs rename to reCTBN/src/params.rs diff --git a/reCTBN/src/sampling.rs b/reCTBN/src/sampling.rs new file mode 100644 index 0000000..0660939 --- /dev/null +++ b/reCTBN/src/sampling.rs @@ -0,0 +1,111 @@ +use crate::{ + network::Network, + params::{self, ParamsTrait}, +}; +use rand::SeedableRng; +use rand_chacha::ChaCha8Rng; + +pub trait Sampler: Iterator { + fn reset(&mut self); +} + +pub struct ForwardSampler<'a, T> +where + T: Network, +{ + net: &'a T, + rng: ChaCha8Rng, + current_time: f64, + current_state: Vec, + next_transitions: Vec>, +} + +impl<'a, T: Network> ForwardSampler<'a, T> { + pub fn new(net: &'a T, seed: Option) -> ForwardSampler<'a, T> { + let rng: ChaCha8Rng = match seed { + //If a seed is present use it to initialize the random generator. + Some(seed) => SeedableRng::seed_from_u64(seed), + //Otherwise create a new random generator using the method `from_entropy` + None => SeedableRng::from_entropy(), + }; + let mut fs = ForwardSampler { + net: net, + rng: rng, + current_time: 0.0, + current_state: vec![], + next_transitions: vec![], + }; + fs.reset(); + return fs; + } +} + +impl<'a, T: Network> Iterator for ForwardSampler<'a, T> { + type Item = (f64, Vec); + + fn next(&mut self) -> Option { + let ret_time = self.current_time.clone(); + let ret_state = self.current_state.clone(); + + for (idx, val) in self.next_transitions.iter_mut().enumerate() { + if let None = val { + *val = Some( + self.net + .get_node(idx) + .get_random_residence_time( + self.net + .get_node(idx) + .state_to_index(&self.current_state[idx]), + self.net.get_param_index_network(idx, &self.current_state), + &mut self.rng, + ) + .unwrap() + + self.current_time, + ); + } + } + + let next_node_transition = self + .next_transitions + .iter() + .enumerate() + .min_by(|x, y| x.1.unwrap().partial_cmp(&y.1.unwrap()).unwrap()) + .unwrap() + .0; + + self.current_time = self.next_transitions[next_node_transition].unwrap().clone(); + + self.current_state[next_node_transition] = self + .net + .get_node(next_node_transition) + .get_random_state( + self.net + .get_node(next_node_transition) + .state_to_index(&self.current_state[next_node_transition]), + self.net + .get_param_index_network(next_node_transition, &self.current_state), + &mut self.rng, + ) + .unwrap(); + + self.next_transitions[next_node_transition] = None; + + for child in self.net.get_children_set(next_node_transition) { + self.next_transitions[child] = None; + } + + Some((ret_time, ret_state)) + } +} + +impl<'a, T: Network> Sampler for ForwardSampler<'a, T> { + fn reset(&mut self) { + self.current_time = 0.0; + self.current_state = self + .net + .get_node_indices() + .map(|x| self.net.get_node(x).get_random_state_uniform(&mut self.rng)) + .collect(); + self.next_transitions = self.net.get_node_indices().map(|_| Option::None).collect(); + } +} diff --git a/src/structure_learning.rs b/reCTBN/src/structure_learning.rs similarity index 100% rename from src/structure_learning.rs rename to reCTBN/src/structure_learning.rs diff --git a/src/structure_learning/constraint_based_algorithm.rs b/reCTBN/src/structure_learning/constraint_based_algorithm.rs similarity index 100% rename from src/structure_learning/constraint_based_algorithm.rs rename to reCTBN/src/structure_learning/constraint_based_algorithm.rs diff --git a/src/structure_learning/hypothesis_test.rs b/reCTBN/src/structure_learning/hypothesis_test.rs similarity index 100% rename from src/structure_learning/hypothesis_test.rs rename to reCTBN/src/structure_learning/hypothesis_test.rs diff --git a/src/structure_learning/score_based_algorithm.rs b/reCTBN/src/structure_learning/score_based_algorithm.rs similarity index 100% rename from src/structure_learning/score_based_algorithm.rs rename to reCTBN/src/structure_learning/score_based_algorithm.rs diff --git a/src/structure_learning/score_function.rs b/reCTBN/src/structure_learning/score_function.rs similarity index 100% rename from src/structure_learning/score_function.rs rename to reCTBN/src/structure_learning/score_function.rs diff --git a/reCTBN/src/tools.rs b/reCTBN/src/tools.rs new file mode 100644 index 0000000..70bbf76 --- /dev/null +++ b/reCTBN/src/tools.rs @@ -0,0 +1,106 @@ +use ndarray::prelude::*; + +use crate::sampling::{ForwardSampler, Sampler}; +use crate::{network, params}; + +pub struct Trajectory { + time: Array1, + events: Array2, +} + +impl Trajectory { + pub fn new(time: Array1, events: Array2) -> Trajectory { + //Events and time are two part of the same trajectory. For this reason they must have the + //same number of sample. + if time.shape()[0] != events.shape()[0] { + panic!("time.shape[0] must be equal to events.shape[0]"); + } + Trajectory { time, events } + } + + pub fn get_time(&self) -> &Array1 { + &self.time + } + + pub fn get_events(&self) -> &Array2 { + &self.events + } +} + +pub struct Dataset { + trajectories: Vec, +} + +impl Dataset { + pub fn new(trajectories: Vec) -> Dataset { + //All the trajectories in the same dataset must represent the same process. For this reason + //each trajectory must represent the same number of variables. + if trajectories + .iter() + .any(|x| trajectories[0].get_events().shape()[1] != x.get_events().shape()[1]) + { + panic!("All the trajectories mus represents the same number of variables"); + } + Dataset { trajectories } + } + + pub fn get_trajectories(&self) -> &Vec { + &self.trajectories + } +} + +pub fn trajectory_generator( + net: &T, + n_trajectories: u64, + t_end: f64, + seed: Option, +) -> Dataset { + //Tmp growing vector containing generated trajectories. + let mut trajectories: Vec = Vec::new(); + + //Random Generator object + + let mut sampler = ForwardSampler::new(net, seed); + //Each iteration generate one trajectory + for _ in 0..n_trajectories { + //History of all the moments in which something changed + let mut time: Vec = Vec::new(); + //Configuration of the process variables at time t initialized with an uniform + //distribution. + let mut events: Vec> = Vec::new(); + + //Current Time and Current State + let (mut t, mut current_state) = sampler.next().unwrap(); + //Generate new samples until ending time is reached. + while t < t_end { + time.push(t); + events.push(current_state); + (t, current_state) = sampler.next().unwrap(); + } + + current_state = events.last().unwrap().clone(); + events.push(current_state); + + //Add t_end as last time. + time.push(t_end.clone()); + + //Add the sampled trajectory to trajectories. + trajectories.push(Trajectory::new( + Array::from_vec(time), + Array2::from_shape_vec( + (events.len(), events.last().unwrap().len()), + events + .iter() + .flatten() + .map(|x| match x { + params::StateType::Discrete(x) => x.clone(), + }) + .collect(), + ) + .unwrap(), + )); + sampler.reset(); + } + //Return a dataset object with the sampled trajectories. + Dataset::new(trajectories) +} diff --git a/tests/ctbn.rs b/reCTBN/tests/ctbn.rs similarity index 100% rename from tests/ctbn.rs rename to reCTBN/tests/ctbn.rs diff --git a/tests/parameter_learning.rs b/reCTBN/tests/parameter_learning.rs similarity index 99% rename from tests/parameter_learning.rs rename to reCTBN/tests/parameter_learning.rs index 5de02d7..7d09b07 100644 --- a/tests/parameter_learning.rs +++ b/reCTBN/tests/parameter_learning.rs @@ -10,6 +10,7 @@ use reCTBN::tools::*; use utils::*; extern crate approx; +use crate::approx::AbsDiffEq; fn learn_binary_cim(pl: T) { let mut net = CtbnNetwork::new(); @@ -427,7 +428,7 @@ fn learn_mixed_discrete_cim(pl: T) { [0.8, 0.6, 0.2, -1.6] ], ]), - 0.1 + 0.2 )); } diff --git a/tests/params.rs b/reCTBN/tests/params.rs similarity index 100% rename from tests/params.rs rename to reCTBN/tests/params.rs diff --git a/tests/structure_learning.rs b/reCTBN/tests/structure_learning.rs similarity index 100% rename from tests/structure_learning.rs rename to reCTBN/tests/structure_learning.rs diff --git a/tests/tools.rs b/reCTBN/tests/tools.rs similarity index 100% rename from tests/tools.rs rename to reCTBN/tests/tools.rs diff --git a/tests/utils.rs b/reCTBN/tests/utils.rs similarity index 100% rename from tests/utils.rs rename to reCTBN/tests/utils.rs diff --git a/src/tools.rs b/src/tools.rs deleted file mode 100644 index 448b26f..0000000 --- a/src/tools.rs +++ /dev/null @@ -1,187 +0,0 @@ -use ndarray::prelude::*; -use rand_chacha::rand_core::SeedableRng; -use rand_chacha::ChaCha8Rng; - -use crate::params::ParamsTrait; -use crate::{network, params}; - -pub struct Trajectory { - time: Array1, - events: Array2, -} - -impl Trajectory { - pub fn new(time: Array1, events: Array2) -> Trajectory { - //Events and time are two part of the same trajectory. For this reason they must have the - //same number of sample. - if time.shape()[0] != events.shape()[0] { - panic!("time.shape[0] must be equal to events.shape[0]"); - } - Trajectory { time, events } - } - - pub fn get_time(&self) -> &Array1 { - &self.time - } - - pub fn get_events(&self) -> &Array2 { - &self.events - } -} - -pub struct Dataset { - trajectories: Vec, -} - -impl Dataset { - pub fn new(trajectories: Vec) -> Dataset { - //All the trajectories in the same dataset must represent the same process. For this reason - //each trajectory must represent the same number of variables. - if trajectories - .iter() - .any(|x| trajectories[0].get_events().shape()[1] != x.get_events().shape()[1]) - { - panic!("All the trajectories mus represents the same number of variables"); - } - Dataset { trajectories } - } - - pub fn get_trajectories(&self) -> &Vec { - &self.trajectories - } -} - -pub fn trajectory_generator( - net: &T, - n_trajectories: u64, - t_end: f64, - seed: Option, -) -> Dataset { - //Tmp growing vector containing generated trajectories. - let mut trajectories: Vec = Vec::new(); - - //Random Generator object - let mut rng: ChaCha8Rng = match seed { - //If a seed is present use it to initialize the random generator. - Some(seed) => SeedableRng::seed_from_u64(seed), - //Otherwise create a new random generator using the method `from_entropy` - None => SeedableRng::from_entropy(), - }; - - //Each iteration generate one trajectory - for _ in 0..n_trajectories { - //Current time of the sampling process - let mut t = 0.0; - //History of all the moments in which something changed - let mut time: Vec = Vec::new(); - //Configuration of the process variables at time t initialized with an uniform - //distribution. - let mut current_state: Vec = net - .get_node_indices() - .map(|x| net.get_node(x).get_random_state_uniform(&mut rng)) - .collect(); - //History of all the configurations of the process variables. - let mut events: Vec> = Vec::new(); - //Vector containing to time to the next transition for each variable. - let mut next_transitions: Vec> = - net.get_node_indices().map(|_| Option::None).collect(); - - //Add the starting time for the trajectory. - time.push(t.clone()); - //Add the starting configuration of the trajectory. - events.push( - current_state - .iter() - .map(|x| match x { - params::StateType::Discrete(state) => state.clone(), - }) - .collect(), - ); - //Generate new samples until ending time is reached. - while t < t_end { - //Generate the next transition time for each uninitialized variable. - for (idx, val) in next_transitions.iter_mut().enumerate() { - if let None = val { - *val = Some( - net.get_node(idx) - .get_random_residence_time( - net.get_node(idx).state_to_index(¤t_state[idx]), - net.get_param_index_network(idx, ¤t_state), - &mut rng, - ) - .unwrap() - + t, - ); - } - } - - //Get the variable with the smallest transition time. - let next_node_transition = next_transitions - .iter() - .enumerate() - .min_by(|x, y| x.1.unwrap().partial_cmp(&y.1.unwrap()).unwrap()) - .unwrap() - .0; - //Check if the next transition take place after the ending time. - if next_transitions[next_node_transition].unwrap() > t_end { - break; - } - //Get the time in which the next transition occurs. - t = next_transitions[next_node_transition].unwrap().clone(); - //Add the transition time to next - time.push(t.clone()); - - //Compute the new state of the transitioning variable. - current_state[next_node_transition] = net - .get_node(next_node_transition) - .get_random_state( - net.get_node(next_node_transition) - .state_to_index(¤t_state[next_node_transition]), - net.get_param_index_network(next_node_transition, ¤t_state), - &mut rng, - ) - .unwrap(); - - //Add the new state to events - events.push(Array::from_vec( - current_state - .iter() - .map(|x| match x { - params::StateType::Discrete(state) => state.clone(), - }) - .collect(), - )); - //Reset the next transition time for the transitioning node. - next_transitions[next_node_transition] = None; - - //Reset the next transition time for each child of the transitioning node. - for child in net.get_children_set(next_node_transition) { - next_transitions[child] = None - } - } - - //Add current_state as last state. - events.push( - current_state - .iter() - .map(|x| match x { - params::StateType::Discrete(state) => state.clone(), - }) - .collect(), - ); - //Add t_end as last time. - time.push(t_end.clone()); - - //Add the sampled trajectory to trajectories. - trajectories.push(Trajectory::new( - Array::from_vec(time), - Array2::from_shape_vec( - (events.len(), current_state.len()), - events.iter().flatten().cloned().collect(), - ) - .unwrap(), - )); - } - //Return a dataset object with the sampled trajectories. - Dataset::new(trajectories) -}