qfall_math/rational/q/sample/
gauss.rs

1// Copyright 2024 Marvin Beckmann
2//
3// This file is part of qFALL-math.
4//
5// qFALL-math is free software: you can redistribute it and/or modify it under
6// the terms of the Mozilla Public License Version 2.0 as published by the
7// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
8
9//! This module contains sampling algorithms for gaussian distributions over [`Q`].
10
11use crate::{
12    error::MathError, macros::seeded::seedable_function, rational::Q,
13    utils::sample::uniform::SamplerRng,
14};
15use probability::{
16    distribution::{Gaussian, Sample},
17    source,
18};
19use rand::Rng;
20
21impl Q {
22    seedable_function!(
23        /// Chooses a [`Q`] instance according to the continuous Gaussian distribution.
24        ///
25        /// Parameters:
26        /// - `center`: specifies the position of the center
27        /// - `sigma`: specifies the standard deviation
28        #[seeded]
29        /// - `seed`: specifies the 256-bit seed of the PRNG used for sampling
30        #[optionally_seeded]
31        /// - `seed`: specifies an optional 256-bit seed of the PRNG used for sampling.
32        #[optionally_seeded]
33        ///   If `None` is provided, a fresh [`ThreadRng`](rand::rngs::ThreadRng) is used instead.
34        ///
35        /// Returns new [`Q`] sample chosen according to the specified continuous Gaussian
36        /// distribution or a [`MathError`] if the specified parameters were not chosen
37        /// appropriately (`sigma > 0`).
38        ///
39        /// # Examples
40        /// ```
41        /// use qfall_math::rational::Q;
42        ///
43        #[unseeded]
44        /// let sample = Q::sample_gauss(0, 1).unwrap();
45        #[seeded]
46        /// let sample = Q::sample_gauss_seeded(0, 1, [42; 32]).unwrap();
47        /// ```
48        ///
49        /// # Errors and Failures
50        /// - Returns a [`MathError`] of type [`NonPositive`](MathError::NonPositive)
51        ///   if `sigma <= 0`.
52        pub(crate) fn sample_gauss(
53            center: impl Into<Q>,
54            sigma: impl Into<f64>,
55            seed: Option<[u8; 32]>,
56        ) -> Result<Q, MathError> {
57            let center = center.into();
58            let sigma = sigma.into();
59            if sigma <= 0.0 {
60                return Err(MathError::NonPositive(format!(
61                    "The sigma has to be positive and not zero, but the provided value is {sigma}."
62                )));
63            }
64            let mut rng = SamplerRng::new(seed);
65            let mut source = source::default(rng.next_u64());
66
67            // Instead of sampling with a center of c, we sample with center 0 and add the
68            // center later. These are equivalent and this way we can sample in larger ranges
69            let sampler = Gaussian::new(0.0, sigma);
70            let sample = center + Q::from(sampler.sample(&mut source));
71
72            Ok(sample)
73        }
74    );
75}
76
77#[cfg(test)]
78mod test_sample_gauss {
79    use crate::rational::Q;
80    use crate::utils::sample::test_seed;
81
82    /// Test correct distribution with a confidence level of 99.7% -> 3 standard
83    /// deviations.
84    #[test]
85    fn in_concentration_bound() {
86        for seeded in [false, true] {
87            let range = 3;
88            for (mu, sigma) in [(i64::MAX, 1), (0, 20), (i64::MIN, 100)] {
89                assert!(
90                    range * sigma
91                        >= (Q::from(mu)
92                            - Q::sample_gauss_optionally_seeded(mu, sigma, seeded.then(test_seed))
93                                .unwrap())
94                        .abs()
95                )
96            }
97        }
98    }
99
100    /// Ensure that an error is returned if `sigma` is not positive
101    #[test]
102    fn non_positive_sigma() {
103        for seeded in [false, true] {
104            for (mu, sigma) in [(0, 0), (0, -1)] {
105                assert!(
106                    Q::sample_gauss_optionally_seeded(mu, sigma, seeded.then(test_seed)).is_err()
107                )
108            }
109        }
110    }
111
112    /// Checks whether the same seed results in the same sample.
113    #[test]
114    fn same_seed_same_sample() {
115        use crate::rational::Q;
116        let sample_0 = Q::sample_gauss_seeded(0, 1, [42; 32]).unwrap();
117        let sample_1 = Q::sample_gauss_seeded(0, 1, [42; 32]).unwrap();
118
119        assert_eq!(sample_0, sample_1);
120    }
121
122    /// Checks whether the seeded variant outputs the known sample for a fixed seed.
123    #[test]
124    fn known_answer() {
125        use crate::rational::Q;
126        let sample = Q::sample_gauss_seeded(0, 1, [42; 32]).unwrap();
127
128        assert_eq!("2896371026169371/2251799813685248", sample.to_string());
129    }
130
131    /// Checks whether the unseeded and the seeded function are available.
132    #[test]
133    fn availability() {
134        use crate::rational::Q;
135
136        // explicitly checks the existence of the unseeded and the seeded function
137        let _ = Q::sample_gauss(0, 1);
138        let _ = Q::sample_gauss_seeded(0, 1, test_seed());
139    }
140}
qfall_math/rational/q/sample/gauss.rs

qfall_math/rational/q/sample/
gauss.rs
