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

extern crate num;
extern crate rand;
extern crate utils;

use utils::map_range;
use rand::{Rand, random};
use std::cell::RefCell;
use std::fmt::Debug;
use num::{Float, FromPrimitive};

/// A thread local global for holding onto the second generated value without requiring the user
/// stores an instance of some type.
thread_local!(static MAYBE_NEXT_VALUE: RefCell<Option<f64>> = RefCell::new(None));

#[inline]
fn two<F>() -> F where F: Float { let one: F = F::one(); one + one }

/// Gen raw gaussian value with dist. at 0.
#[inline]
pub fn gen_raw<F>() -> F where F: Float + FromPrimitive + Rand {
    MAYBE_NEXT_VALUE.with(|maybe_next_value| {
        let maybe_next = *maybe_next_value.borrow();
        if let Some(next_value) = maybe_next {
            *maybe_next_value.borrow_mut() = None;
            FromPrimitive::from_f64(next_value).unwrap()
        } else {
            let (zero, one, two): (F, F, F) = (F::zero(), F::one(), two::<F>());
            let (mut va, mut vb, mut s): (F, F, F) = (zero, zero, zero);
            while s >= one || s == zero {
                va = two * random::<F>() - one;
                vb = two * random::<F>() - one;
                s = va * vb + va * vb
            };
            let multi = ((-two) * s.abs().ln() / s).abs().sqrt();
            *maybe_next_value.borrow_mut() = (vb * multi).to_f64();
            va * multi
        }
    })
}

/// Gen gaussian value with dist. at 'n' with rand randomness.
/// Generated value will always be 0.0 <= value < 1.0.
#[inline]
pub fn gen<F>(n: F, randomness: f32) -> F
where F: Float + Rand + FromPrimitive + Debug {
    let (zero, one): (F, F) = (F::zero(), F::one());
    assert!(n >= zero && n <= one, "Gaussian::gen : given `n` ({:?}) must \
            be a percentage between 0 and 1.", n);

    // If one was given with no randomness, return it exactly as is.
    if n == one && randomness == 0.0 {
        return one - F::epsilon();
    }

    let mut ans = gen_raw::<F>()
                * FromPrimitive::from_f32(randomness.powf(2.0)).unwrap()
                + (n * two::<F>() - one);
    ans = map_range(ans, -one, one, zero, one);
    if ans >= one || ans < zero {
        gen::<F>(n, randomness)
    } else {
        ans
    }
}

/// Gen gaussian value mapped to a range.
#[inline]
pub fn gen_map<F>(n: F, randomness: f32, min_range: F, max_range: F) -> F
where F: Float + Rand + FromPrimitive + Debug {
    let (zero, one): (F, F) = (F::zero(), F::one());
    let perc = map_range(n, min_range, max_range, zero, one);
    map_range(gen(perc, randomness), zero, one, min_range, max_range)
}