math_audio_optimisation/
init_sobol.rs

1/// Generate Halton quasi-random sequence for initialization.
2///
3/// Uses Van der Corput sequences in the first `dimensions` prime bases.
4/// This provides better parameter space coverage than pure random
5/// initialization and avoids the correlated-base bug of the old fallback.
6///
7/// # Arguments
8/// * `dimensions` - Number of dimensions (parameters)
9/// * `num_samples` - Number of samples to generate
10/// * `bounds` - Parameter bounds for scaling
11///
12/// # Returns
13/// Vector of parameter vectors sampled from the quasi-random sequence
14pub fn init_halton(dimensions: usize, num_samples: usize, bounds: &[(f64, f64)]) -> Vec<Vec<f64>> {
15    let bases = halton_bases(dimensions);
16    let mut samples = Vec::with_capacity(num_samples);
17
18    for i in 0..num_samples {
19        let mut sample = Vec::with_capacity(dimensions);
20
21        for (dim, &(lower, upper)) in bounds.iter().enumerate().take(dimensions) {
22            let quasi_random = van_der_corput(i + 1, bases[dim]);
23            sample.push(lower + quasi_random * (upper - lower));
24        }
25
26        samples.push(sample);
27    }
28
29    samples
30}
31
32/// Deprecated alias for [`init_halton`].
33/// The old name claimed a Sobol sequence, but the implementation was a
34/// Halton-style Van der Corput construction.
35#[deprecated(since = "0.5.9", note = "Use init_halton instead")]
36pub fn init_sobol(dimensions: usize, num_samples: usize, bounds: &[(f64, f64)]) -> Vec<Vec<f64>> {
37    init_halton(dimensions, num_samples, bounds)
38}
39
40/// Return the first `n` prime numbers to use as Halton bases.
41fn halton_bases(n: usize) -> Vec<usize> {
42    let mut primes = Vec::with_capacity(n);
43    let mut candidate = 2;
44    while primes.len() < n {
45        if is_prime(candidate) {
46            primes.push(candidate);
47        }
48        candidate += 1;
49    }
50    primes
51}
52
53fn is_prime(n: usize) -> bool {
54    if n < 2 {
55        return false;
56    }
57    if n == 2 {
58        return true;
59    }
60    if n.is_multiple_of(2) {
61        return false;
62    }
63    let mut i = 3;
64    while i * i <= n {
65        if n.is_multiple_of(i) {
66            return false;
67        }
68        i += 2;
69    }
70    true
71}
72
73#[cfg(test)]
74fn gcd(a: usize, b: usize) -> usize {
75    if b == 0 { a } else { gcd(b, a % b) }
76}
77
78/// Van der Corput sequence for quasi-random number generation.
79fn van_der_corput(mut n: usize, base: usize) -> f64 {
80    let mut result = 0.0;
81    let mut f = 1.0 / base as f64;
82
83    while n > 0 {
84        result += (n % base) as f64 * f;
85        n /= base;
86        f /= base as f64;
87    }
88
89    result
90}
91
92#[cfg(test)]
93mod tests {
94    use super::*;
95
96    #[test]
97    fn test_van_der_corput() {
98        let val1 = van_der_corput(1, 2);
99        let val2 = van_der_corput(2, 2);
100        let val3 = van_der_corput(3, 2);
101
102        assert!((0.0..1.0).contains(&val1));
103        assert!((0.0..1.0).contains(&val2));
104        assert!((0.0..1.0).contains(&val3));
105
106        assert_ne!(val1, val2);
107        assert_ne!(val2, val3);
108    }
109
110    #[test]
111    fn test_init_halton() {
112        let bounds = vec![(0.0, 10.0), (0.1, 5.0), (-12.0, 12.0)];
113        let samples = init_halton(3, 5, &bounds);
114
115        assert_eq!(samples.len(), 5);
116        for sample in &samples {
117            assert_eq!(sample.len(), 3);
118            assert!((0.0..=10.0).contains(&sample[0]));
119            assert!((0.1..=5.0).contains(&sample[1]));
120            assert!((-12.0..=12.0).contains(&sample[2]));
121        }
122    }
123
124    #[test]
125    fn test_halton_bases_are_coprime() {
126        // The old fallback used bases like 2 + (dim % 10), which are not coprime
127        // (e.g. dim 8 used base 10, dim 10 used base 2 again). A proper Halton
128        // sequence uses the first n primes, which are pairwise coprime.
129        let bases = halton_bases(20);
130        for i in 0..bases.len() {
131            for j in (i + 1)..bases.len() {
132                assert_eq!(
133                    gcd(bases[i], bases[j]),
134                    1,
135                    "bases {} and {} must be coprime",
136                    bases[i],
137                    bases[j]
138                );
139            }
140        }
141    }
142
143    #[test]
144    fn test_halton_no_cycling() {
145        // With the old dim%10 cycling, dim 10 reused base 2, causing correlation.
146        // Halton uses unique primes for each dimension.
147        let bases = halton_bases(15);
148        let unique: std::collections::HashSet<_> = bases.iter().cloned().collect();
149        assert_eq!(
150            unique.len(),
151            bases.len(),
152            "each dimension must have a unique base"
153        );
154    }
155}
math_audio_optimisation/init_sobol.rs

math_audio_optimisation/
init_sobol.rs
