Merge branch 'dev' into '8-feature-constraint-based-structure-learning-algorithm-for-ctbn'

Synchronized with dev

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",
+]

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" }

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;

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<params::StateType>,
+    next_transitions: Vec<Option<f64>>,
+}
+
+impl<'a, T: Network> ForwardSampler<'a, T> {
+    pub fn new(net: &'a T, seed: Option<u64>) -> 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<params::StateType>);
+
+    fn next(&mut self) -> Option<Self::Item> {
+        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();
+    }
+}

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<f64>,
+    events: Array2<usize>,
+}
+
+impl Trajectory {
+    pub fn new(time: Array1<f64>, events: Array2<usize>) -> 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<f64> {
+        &self.time
+    }
+
+    pub fn get_events(&self) -> &Array2<usize> {
+        &self.events
+    }
+}
+
+pub struct Dataset {
+    trajectories: Vec<Trajectory>,
+}
+
+impl Dataset {
+    pub fn new(trajectories: Vec<Trajectory>) -> 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<Trajectory> {
+        &self.trajectories
+    }
+}
+
+pub fn trajectory_generator<T: network::Network>(
+    net: &T,
+    n_trajectories: u64,
+    t_end: f64,
+    seed: Option<u64>,
+) -> Dataset {
+    //Tmp growing vector containing generated trajectories.
+    let mut trajectories: Vec<Trajectory> = 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<f64> = Vec::new();
+        //Configuration of the process variables at time t initialized with an uniform
+        //distribution.
+        let mut events: Vec<Vec<params::StateType>> = 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)
+}

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<T: ParameterLearning>(pl: T) {
     let mut net = CtbnNetwork::new();
@@ -427,7 +428,7 @@ fn learn_mixed_discrete_cim<T: ParameterLearning>(pl: T) {
                 [0.8, 0.6, 0.2, -1.6]
             ],
         ]),
-        0.1
+        0.2
     ));
 }
 

src/tools.rs

@@ -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<f64>,
-    events: Array2<usize>,
-}
-
-impl Trajectory {
-    pub fn new(time: Array1<f64>, events: Array2<usize>) -> 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<f64> {
-        &self.time
-    }
-
-    pub fn get_events(&self) -> &Array2<usize> {
-        &self.events
-    }
-}
-
-pub struct Dataset {
-    trajectories: Vec<Trajectory>,
-}
-
-impl Dataset {
-    pub fn new(trajectories: Vec<Trajectory>) -> 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<Trajectory> {
-        &self.trajectories
-    }
-}
-
-pub fn trajectory_generator<T: network::Network>(
-    net: &T,
-    n_trajectories: u64,
-    t_end: f64,
-    seed: Option<u64>,
-) -> Dataset {
-    //Tmp growing vector containing generated trajectories.
-    let mut trajectories: Vec<Trajectory> = 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<f64> = Vec::new();
-        //Configuration of the process variables at time t initialized with an uniform
-        //distribution.
-        let mut current_state: Vec<params::StateType> = 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<Array1<usize>> = Vec::new();
-        //Vector containing to time to the next transition for each variable.
-        let mut next_transitions: Vec<Option<f64>> =
-            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(&current_state[idx]),
-                                net.get_param_index_network(idx, &current_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(&current_state[next_node_transition]),
-                    net.get_param_index_network(next_node_transition, &current_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)
-}