prodef 0.1.0

A simple Rust crate for handling probability distributions.
Documentation 
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

//! A module that implements a particle based probability density functions.

mod normal;

use nalgebra::{
    DefaultAllocator, Dim, Dyn, OMatrix, OVector, RealField, Scalar, allocator::Allocator,
};
use rand::Rng;
use rand_distr::{Uniform, uniform::SampleUniform};
use serde::{Deserialize, Serialize};

/// A population of particles representing a probability density function.
#[derive(Clone, Debug, Deserialize, Serialize)]
#[serde(bound(
    serialize = "D: Serialize, B: Serialize, K: Serialize, OMatrix<T, D, Dyn>: Serialize, OVector<T, Dyn>: Serialize"
))]
#[serde(bound(
    deserialize = "D: Deserialize<'de>, B: Deserialize<'de>, K: Deserialize<'de>, OMatrix<T, D, Dyn>: Deserialize<'de>, OVector<T, Dyn>: Deserialize<'de>"
))]
pub struct ParticleDensity<T, D, B, K>
where
    T: Scalar,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<D, Dyn>,
{
    particles: OMatrix<T, D, Dyn>,
    opt_weights: Option<OVector<T, Dyn>>,
    domain: B,
    kernel: K,
}

impl<T, D, B, K> ParticleDensity<T, D, B, K>
where
    T: RealField + SampleUniform,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<D, Dyn>,
{
    /// Get a reference to the kernel density.
    pub fn kernel(&self) -> &K {
        &self.kernel
    }

    /// Returns the number of particles
    #[allow(clippy::len_without_is_empty)]
    pub fn len(&self) -> usize {
        self.particles.ncols()
    }

    /// Get a reference to the particles.
    pub fn particles(&self) -> &OMatrix<T, D, Dyn> {
        &self.particles
    }

    /// Draw a random sample particle.
    pub fn sample_particle(&self, rng: &mut impl Rng) -> OVector<T, D> {
        let pdx = {
            match &self.opt_weights {
                Some(weights) =>
                // Here we abuse try_fold to return particle index early wrapped within Err().
                {
                    let uniform = Uniform::new(T::zero(), T::one()).unwrap();

                    // Select particle index by weight.
                    let wdx = rng.sample(uniform);

                    match weights
                        .iter()
                        .enumerate()
                        .try_fold(T::zero(), |acc, (idx, weight)| {
                            let next_weight = acc + weight.clone();
                            if wdx < next_weight {
                                Err(idx)
                            } else {
                                Ok(next_weight)
                            }
                        }) {
                        Ok(_) => weights.len() - 1,
                        Err(idx) => idx,
                    }
                }
                None => {
                    let uniform = Uniform::new(0, self.particles.ncols()).unwrap();

                    // Select particle index by weight.
                    rng.sample(uniform)
                }
            }
        };

        self.particles.column(pdx).clone_owned()
    }

    /// Get a reference to the particle weights.
    pub fn weights(&self) -> &Option<OVector<T, Dyn>> {
        &self.opt_weights
    }
}