pointprocesses/spatial/
processes.rs

1use super::domains::Domain;
2
3use rand::prelude::*;
4use rand::rngs::SmallRng;
5use rand_distr::{Poisson, Distribution};
6
7use ndarray::stack;
8use ndarray::prelude::*;
9
10static XORSHIFT_ERR: &str = "Unable to create XorShift rng from thread local rng";
11
12/// A higher-dimensional homogeneous Poisson process, for parallepipedal domains.
13pub fn poisson_process(lambda: f64, domain: &Domain) -> Array2<f64>
14{
15    let ref mut rng = thread_rng();
16    let far = &domain.far;
17    let close = &domain.close;
18
19    // dimension of space
20    let d = far.shape()[0];
21    let area = (0..d).fold(1.0, |area, i| {
22        area * (far[i] - close[i])
23    });
24
25    // get number of events to generate
26    // events outside of the set will be rejected
27    let fish = Poisson::new(lambda*area).unwrap();
28    let num_events: u64 = fish.sample(rng);
29    
30    let mut srng = SmallRng::from_rng(rng).expect(XORSHIFT_ERR);
31
32    let events: Vec<Array2<f64>> = (0..num_events).map(|_| {
33        // generate a point inside the bounding box
34        let mut ev: Array1<f64> = Array::zeros((d,));
35
36        for i in 0..d {
37            ev[i] = srng.gen_range(close[i], far[i]);
38        }
39
40        // if it's in, then keep it
41        ev.into_shape((1,d)).unwrap()
42    }).collect();
43
44    let events_ref: Vec<ArrayView2<f64>> = events.iter().map(|ev| {
45        ev.view()
46    }).collect();
47
48    stack(Axis(0), events_ref.as_slice()).unwrap()
49}
50
51/// Poisson process on a d-dimensional region with variable intensity, using a rejection sampling algorithm.
52pub fn variable_poisson<F>(lambda: F, max_lambda: f64, domain: &Domain) -> Array2<f64>
53    where F: Fn(&Array1<f64>) -> f64 + Sync + Send
54{
55    let close = &domain.close;
56    let far = &domain.far;
57
58    let d = close.shape()[0];
59    let area = (0..d).fold(1.0, |area, i| {
60        area * (far[i] - close[i])
61    });
62
63    // get number of events to generate
64    // events outside of the set will be rejected
65    let ref mut rng = thread_rng();
66    let fish = Poisson::new(max_lambda*area).unwrap();
67    let num_events: u64 = fish.sample(rng);
68    
69    let mut srng = rand::rngs::SmallRng::from_rng(rng).expect(XORSHIFT_ERR);
70
71    let events: Vec<Array2<f64>> = (0..num_events).filter_map(|_| {
72        // generate a point inside the bounding box
73        let mut ev: Array1<f64> = Array::zeros((d,));
74        let intens = max_lambda*random::<f64>();
75
76        for i in 0..d {
77            ev[i] = srng.gen_range(close[i], far[i]);
78        }
79
80        // if it's in, then keep it
81        if intens < lambda(&ev) {
82            Some(ev.into_shape((1,d)).unwrap())
83        } else {
84            None
85        }
86    }).collect();
87
88    let events_ref: Vec<ArrayView2<f64>> = events.iter().map(|ev| {
89        ev.view()
90    }).collect();
91
92    stack(Axis(0), &events_ref).unwrap()
93}
pointprocesses/spatial/processes.rs

pointprocesses/spatial/
processes.rs
