rustyqlib/core/montecarlo/
halton.rs

1//! Halton low-discrepancy sequences: prime radical inverses with a
2//! seeded Cranley-Patterson rotation. Any number of dimensions, so this
3//! is the quasi-random workhorse when the problem's dimension exceeds
4//! the embedded Sobol table ([`SobolSequence`](super::sobol::SobolSequence)).
5
6use crate::core::utils::inv_norm_cdf;
7
8use super::rng::splitmix64;
9
10fn first_primes(n: usize) -> Vec<u64> {
11    let mut primes: Vec<u64> = Vec::with_capacity(n);
12    let mut candidate = 2u64;
13    while primes.len() < n {
14        if primes.iter().take_while(|&&p| p * p <= candidate).all(|&p| candidate % p != 0) {
15            primes.push(candidate);
16        }
17        candidate += 1;
18    }
19    primes
20}
21
22fn radical_inverse(mut i: u64, base: u64) -> f64 {
23    let inv_base = 1.0 / base as f64;
24    let mut f = inv_base;
25    let mut x = 0.0;
26    while i > 0 {
27        x += f * (i % base) as f64;
28        i /= base;
29        f *= inv_base;
30    }
31    x
32}
33
34/// Multi-dimensional low-discrepancy sequence: Halton with prime bases and
35/// a deterministic Cranley-Patterson rotation per dimension (derived from
36/// the seed), mapped to standard normals through the inverse CDF.
37///
38/// Combined with [`BrownianBridge`](super::brownian_bridge::BrownianBridge)
39/// ordering, the well-distributed leading dimensions carry the coarse
40/// structure of each path.
41pub struct QmcSequence {
42    bases: Vec<u64>,
43    shifts: Vec<f64>,
44}
45
46impl QmcSequence {
47    pub fn new(dims: usize, seed: u64) -> Self {
48        let bases = first_primes(dims);
49        let shifts = (0..dims)
50            .map(|d| (splitmix64(seed ^ splitmix64(0xC0FFEE ^ d as u64)) >> 11) as f64
51                / (1u64 << 53) as f64)
52            .collect();
53        QmcSequence { bases, shifts }
54    }
55
56    /// Fill `out` with the standard normals of point `index` (1-based).
57    pub fn normals(&self, index: u64, out: &mut [f64]) {
58        for (d, z) in out.iter_mut().enumerate() {
59            let mut u = radical_inverse(index, self.bases[d]) + self.shifts[d];
60            if u >= 1.0 {
61                u -= 1.0;
62            }
63            *z = inv_norm_cdf(u.clamp(1e-15, 1.0 - 1e-15));
64        }
65    }
66}
67
68#[cfg(test)]
69mod tests {
70    use super::*;
71
72    #[test]
73    fn qmc_sequence_normals_are_standardized() {
74        let steps = 32;
75        let qmc = QmcSequence::new(steps, 42);
76        let mut z = vec![0.0; steps];
77        let n = 8192;
78        let mut mean = vec![0.0; steps];
79        let mut var = vec![0.0; steps];
80        for i in 1..=n {
81            qmc.normals(i as u64, &mut z);
82            for d in 0..steps {
83                mean[d] += z[d];
84                var[d] += z[d] * z[d];
85            }
86        }
87        for d in 0..steps {
88            let m = mean[d] / n as f64;
89            let v = var[d] / n as f64;
90            assert!(m.abs() < 0.05, "dim {d}: mean {m}");
91            assert!((v - 1.0).abs() < 0.1, "dim {d}: var {v}");
92        }
93    }
94}
rustyqlib/core/montecarlo/halton.rs

rustyqlib/core/montecarlo/
halton.rs
