imgal 0.3.0

A fast and open-source scientific image processing and algorithm library.
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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

use std::ops::{Bound, RangeBounds};

use crate::prelude::*;

const INCREMENT: u64 = 1442695040888963407;
const MULTIPLIER: u64 = 6364136223846793005;

/// A 32-bit Permuted Congruential Generator (PCG) pseudo-random number
/// generator.
///
/// The PCG stores a 64-bit internal state and uses an xorshift and bit rotation
/// to produce u32 values. The generator is deterministic when given the same
/// seed.
#[derive(Clone, Debug)]
pub struct Pcg {
    /// The PCG state.
    state: u64,
}

impl Pcg {
    /// Fork the PCG with a random seed value.
    ///
    /// # Description
    ///
    /// Forks the current PCG with a random seed value, returning a new PCG.
    ///
    /// # Returns
    ///
    /// * `Pcg`: A new randomly seeded PCG.
    pub fn fork(&mut self) -> Self {
        Self::new(self.next_u32() as u64)
    }

    /// Create a new seeded PCG.
    ///
    /// # Description
    ///
    /// Creates a new PCG with the given `seed` value.
    ///
    /// # Arguments
    ///
    /// * `seed`: The seed value for the PCG.
    ///
    /// # Returns
    ///
    /// * `Pcg`: A seeded PCG.
    pub fn new(seed: u64) -> Self {
        let mut rng = Self { state: 0 };
        // adding the increment and seed gives the PCG something to work with,
        // this "warms" up the RNG
        rng.state = rng.state.wrapping_add(INCREMENT);
        rng.next_u32();
        rng.state = rng.state.wrapping_add(seed);
        rng.next_u32();
        rng
    }

    /// Return a pseudo-random f32 value.
    ///
    /// # Description
    ///
    /// Returns a pseudo-random f32 value in the half-open interval [0, 1).
    ///
    /// # Returns
    ///
    /// * `f32`: A pseudo-random f32 value in the half-open interval [0, 1).
    pub fn next_f32(&mut self) -> f32 {
        (self.next_u32() >> 8) as f32 / (1u32 << 24) as f32
    }

    /// Return a pseudo-random u32 value.
    ///
    /// # Description
    ///
    /// Returns a pseudo-random u32 value.
    ///
    /// # Returns
    ///
    /// * `u32`: A pseudo-random u32 value.
    pub fn next_u32(&mut self) -> u32 {
        let old_state = self.state;
        self.state = old_state.wrapping_mul(MULTIPLIER).wrapping_add(INCREMENT);
        let xor_shifted = (((old_state >> 18) ^ old_state) >> 27) as u32;
        let rot = (old_state >> 59) as u32;
        xor_shifted.rotate_right(rot)
    }

    /// Return a pseudo-random u32 value within a given range.
    ///
    /// # Description
    ///
    /// Returns a pseudo-random u32 value within a given range. The start value
    /// of the range must be smaller than the end value.
    ///
    /// # Arguments
    ///
    /// * `r`: The positive range to limit valid pseudo-random u32 values.
    /// # Returns
    ///
    /// * `Ok(u32)`: A pseudo-random u32 value within the given range, `r`.
    /// * `Err(ImgalError)`: If the range start value is larger than the range
    ///   end.
    pub fn next_u32_range<R>(&mut self, r: R) -> Result<u32, ImgalError>
    where
        R: RangeBounds<u32>,
    {
        let start = match r.start_bound() {
            Bound::Included(&s) => s,
            Bound::Excluded(&s) => s.saturating_add(1),
            Bound::Unbounded => 0,
        };
        let end = match r.end_bound() {
            Bound::Included(&e) => e.saturating_add(1),
            Bound::Excluded(&e) => e,
            Bound::Unbounded => u32::MAX,
        };
        if start >= end {
            return Err(ImgalError::InvalidPositiveRange {
                start: start as usize,
                end: end as usize,
            });
        }
        let diff = end - start;
        // this threshold value is used to avoid "modulo bias" when 2^32 (i.e. u32)
        // can't be evenly divided by the range diff
        let threshold = diff.wrapping_neg() % diff;
        loop {
            let v = self.next_u32();
            if v >= threshold {
                return Ok(start + (v % diff));
            }
        }
    }
}