1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

use generalized::domains::Set;

use ndarray::stack;
use ndarray::prelude::*;

use rand::thread_rng;
use rand::distributions::Uniform;
use rand::distributions::Poisson;
use rand::prelude::*;

/// A higher-dimensional homogeneous Poisson process.
pub fn poisson_process<T>(lambda: f64, domain: &T) -> Array2<f64> 
    where T: Set {
    let bounds = domain.bounding_box();
    let mut area = 1.0;

    let n = bounds.shape()[0];
    let d = bounds.shape()[1];

    for i in 0..n {
        area *= bounds[[1,i]] - bounds[[0,i]];
    }

    // get number of events to generate
    // events outside of the set will be rejected
    let ref mut rng = thread_rng();
    let num_events = Poisson::new(lambda*area).sample(rng) as usize;

    let mut res = unsafe {
        Array::uninitialized((1,d))
    };
    
    for _ in 0..num_events {
        // generate a point inside the bounding box
        let mut ev = Array::zeros((d,));

        for i in 0..d {
            ev[i] = rng.sample(Uniform::new(bounds[[0,i]], bounds[[1,i]]));
        }

        // if it's in, then keep it
        if domain.contains(&ev) {
            res = stack(
                Axis(0), 
                &[res.view(), ev.into_shape((1,d)).unwrap().view()]
                ).unwrap();
        }
    }

    res
}

/// Poisson process on a d-dimensional region with variable intensity, using a rejection sampling algorithm.
pub fn variable_poisson<F, T>(lambda: F, max_lambda: f64, domain: &T) -> Array2<f64>
where F: Fn(&Array1<f64>) -> f64,
      T: Set
{
    let bounds = domain.bounding_box();
    let mut area = 1.0;

    let n = bounds.shape()[0];
    let d = bounds.shape()[1];

    for i in 0..n {
        area *= bounds[[1,i]] - bounds[[0,i]];
    }

    // get number of events to generate
    // events outside of the set will be rejected
    let ref mut rng = thread_rng();
    let num_events = Poisson::new(max_lambda*area).sample(rng) as usize;

    let mut res = unsafe {
        Array::uninitialized((1,d))
    };
    
    for _ in 0..num_events {
        // generate a point inside the bounding box
        let mut ev = Array::zeros((d,));
        let intens = max_lambda*random::<f64>();

        for i in 0..d {
            ev[i] = rng.sample(Uniform::new(bounds[[0,i]], bounds[[1,i]]));
        }

        // if the point lies in the domain and the simulated intensity
        // fits, then add it.
        if domain.contains(&ev) && intens < lambda(&ev) {
            res = stack(
                Axis(0), 
                &[res.view(), ev.into_shape((1,d)).unwrap().view()]
                ).unwrap();
        }
    }

    res
}