pointprocesses/temporal/
poisson.rs

1//! Poisson processes.
2use super::traits::*;
3use rand::prelude::*;
4use rand_distr::{Uniform, Poisson, Distribution};
5
6use ndarray::array;
7use ndarray::prelude::*;
8use ndarray_parallel::prelude::*;
9
10use rayon::prelude::*;
11
12
13/// Homogeneous, constant intensity Poisson process.
14/// The intensity of the process is given by the average
15/// number of events between two instants:
16/// $$
17/// \lambda = \lim_{h\to 0}
18/// \frac{\mathbb E[N_{t+h} - N_t]}{h}
19/// $$
20#[derive(Debug)]
21pub struct PoissonProcess {
22    /// Process intensity.
23    lambda: f64
24}
25
26impl PoissonProcess {
27    pub fn new(lambda: f64) -> Self {
28        PoissonProcess {
29            lambda
30        }
31    }
32}
33
34impl DeterministicIntensity for PoissonProcess {
35    fn intensity(&self, _t: f64) -> f64 {
36        self.lambda
37    }
38}
39
40/// Poisson process with variable intensity.
41/// The average number of events between $t$ and $t+dt$ is
42/// $$ \mathbb{E}[ dN_t ] = \lambda(t) dt $$
43/// where $\lambda(t)$ is a deterministic process.
44#[derive(Debug)]
45pub struct VariablePoissonProcess<F>
46where F: Fn(f64) -> f64 + Send + Sync
47{
48    /// Upper bound on the intensity function of the process.
49    max_lambda: f64,
50    /// Process intensity function.
51    func: F
52}
53
54impl<F> VariablePoissonProcess<F>
55where F: Fn(f64) -> f64 + Send + Sync
56{
57    pub fn new(func: F, max_lambda: f64) -> Self {
58        VariablePoissonProcess {
59            max_lambda,
60            func
61        }
62    }
63
64    /// Get the max_lambda upper bound on the intensity.
65    pub fn get_max_lambda(&self) -> f64 {
66        self.max_lambda
67    }
68}
69
70impl<F> DeterministicIntensity for VariablePoissonProcess<F>
71where F: Fn(f64) -> f64 + Send + Sync
72{
73    fn intensity(&self, t: f64) -> f64 {
74        (self.func)(t)
75    }
76}
77
78
79impl TemporalProcess for PoissonProcess {
80    fn sample(&self, tmax: f64) -> TimeProcessResult {
81        let lambda = self.lambda;
82        let mut rng = thread_rng();
83        let fish = Poisson::new(tmax * lambda).unwrap();
84        let num_events: u64 = fish.sample(&mut rng);
85        let num_events = num_events as usize;
86        
87        let mut events_vec: Vec<_> = (0..num_events).into_par_iter()
88            .map(|_| {
89                // get reference to local thread rng
90                let mut rng = thread_rng();
91                let u = Uniform::new(0.0, tmax);
92                u.sample(&mut rng)
93        }).collect();
94        events_vec.sort_by(|a, b| a.partial_cmp(b).unwrap());
95        let timestamps = Array1::<f64>::from_vec(events_vec);
96        let mut intensities = Array1::<f64>::zeros(num_events as usize);
97        for i in 0..num_events as usize {
98            intensities[i] = lambda;
99        }
100        TimeProcessResult {
101            timestamps, intensities
102        }
103    }
104}
105
106impl<F> TemporalProcess for VariablePoissonProcess<F>
107where F: Fn(f64) -> f64 + Send + Sync
108{
109    fn sample(&self, tmax: f64) -> TimeProcessResult {
110        // Parallelized, multithreaded algorithm for sampling
111        // from the process.
112
113        let mut rng = thread_rng();
114        let max_lambda = self.max_lambda;
115        let lambda = &self.func;
116        let fish = Poisson::new(tmax*max_lambda).unwrap();
117        let num_events: u64 = fish.sample(&mut rng);
118        let num_events = num_events as usize;
119        let lambda = std::sync::Arc::from(lambda);
120        
121        // Get timestamp and intensity values of events distributed
122        // according to a homogeneous Poisson process
123        // and keep those who are under the intensity curve
124        let mut events: Vec<Array1<f64>> = (0..num_events)
125                .into_par_iter().filter_map(|_| {
126            let mut rng = thread_rng();
127            let timestamp = rng.gen::<f64>()*tmax;
128            let lambda_val = rng.gen::<f64>()*max_lambda;
129
130            if lambda_val < lambda(timestamp) {
131                Some(array![timestamp, lambda(timestamp)])
132            } else {
133                None
134            }
135        }).collect();
136        events.sort_by(|a, b| a[0].partial_cmp(&b[0]).unwrap());
137
138        let num_events = events.len();
139
140        let mut timestamps = Array1::<f64>::zeros(num_events);
141        let mut intensities = Array1::<f64>::zeros(num_events);
142        if num_events > 0 {
143            for i in 0..num_events {
144                timestamps[i] = events[i][0];
145                intensities[i] = events[i][1];
146
147            }
148        }
149        TimeProcessResult { timestamps, intensities }
150    }
151}
pointprocesses/temporal/poisson.rs

pointprocesses/temporal/
poisson.rs
