scirs2_optimize/global/

qmc.rs

//! Shared quasi-Monte Carlo low-discrepancy sequence generators.
//!
//! Used by [`crate::global::differential_evolution`], [`crate::global::multi_start`],
//! and the crate-private `clustering` module for population/starting-point
//! initialization. Previously each of those three modules had its own
//! independent, non-equivalent substitute:
//!
//! - `differential_evolution`'s `SobolState` used an ad-hoc "polynomial
//!   recurrence" for direction numbers that is not a Sobol sequence at all.
//! - `multi_start`'s `generate_halton_points`/`generate_sobol_points` were
//!   pure no-ops that silently fell back to plain uniform random sampling.
//! - `clustering`'s `generate_sobol_points` used a Van-der-Corput sequence
//!   with a shifted base per dimension (a real low-discrepancy sequence,
//!   but not Sobol despite the name), and its "random" points came from a
//!   deterministic `(t * 17.0).fract()` formula, not an RNG at all.
//!
//! This module provides one validated implementation each callers can share.
//!
//! # Sobol accuracy
//!
//! [`SobolGenerator`] implements genuine Sobol direction numbers (Bratley &
//! Fox / ACM TOMS Algorithm 659 primitive polynomials, the same table
//! reproduced in e.g. Numerical Recipes' `sobseq`) for dimensions
//! `0..MAX_SOBOL_DIM`, verified bit-for-bit against
//! `scipy.stats.qmc.Sobol(scramble=False)` (see this module's tests). Beyond
//! `MAX_SOBOL_DIM`, rather than guess at an unverified table extension or
//! silently degrade to randomness, it falls back to [`halton_radical_inverse`]
//! with a distinct prime base per dimension -- still a genuine
//! low-discrepancy sequence, just not literally Sobol.

/// Number of bits of precision used for the Sobol integer state.
const MAXBIT: u32 = 30;

/// Number of dimensions (including the trivial dimension 0, which is the
/// base-2 van der Corput sequence and needs no table entry) with genuine,
/// validated Sobol direction numbers.
pub const MAX_SOBOL_DIM: usize = 7;

/// Primitive-polynomial degree, polynomial coefficient bits, and initial
/// direction numbers (zero-padded to 4) for Sobol dimensions `1..MAX_SOBOL_DIM`.
const SOBOL_TABLE: [(u32, u32, [u32; 4]); MAX_SOBOL_DIM - 1] = [
    (1, 0, [1, 0, 0, 0]),
    (2, 1, [1, 3, 0, 0]),
    (3, 1, [1, 3, 1, 0]),
    (3, 2, [1, 1, 1, 0]),
    (4, 1, [1, 1, 3, 3]),
    (4, 4, [1, 3, 5, 13]),
];

/// First 32 primes, used as Halton sequence bases (including the fallback
/// path for Sobol dimensions beyond [`MAX_SOBOL_DIM`]).
const PRIMES: [u64; 32] = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
    101, 103, 107, 109, 113, 127, 131,
];

/// Radical-inverse function: the `index`-th term of the (1-indexed) Halton
/// sequence for the given prime `base`.
pub fn halton_radical_inverse(index: usize, base: u64) -> f64 {
    let mut result = 0.0_f64;
    let mut f = 1.0 / base as f64;
    let mut i = index as u64;

    while i > 0 {
        result += f * (i % base) as f64;
        i /= base;
        f /= base as f64;
    }

    result
}

/// Builds the full `MAXBIT`-deep direction-number table for one Sobol
/// dimension from its primitive polynomial (degree `mdeg`, coefficient bits
/// `ip`) and initial direction numbers `iv` (`iv[i-1] = m_i` for `i =
/// 1..=mdeg`, per Bratley & Fox).
fn build_direction_numbers(mdeg: u32, ip: u32, iv: [u32; 4]) -> Vec<u32> {
    let maxbit = MAXBIT as usize;
    let deg = mdeg as usize;
    let mut v = vec![0u32; maxbit + 1];

    for i in 1..=deg {
        v[i] = iv[i - 1] << (MAXBIT - i as u32);
    }

    // Polynomial coefficient bits a_1..a_{mdeg-1}, MSB-first, extracted from
    // `ip` (which encodes exactly those bits, without the implicit leading
    // and trailing 1 coefficients every primitive polynomial has).
    let nbits = deg.saturating_sub(1);
    let mut bits = vec![0u32; nbits];
    for (k, bit) in bits.iter_mut().enumerate() {
        *bit = (ip >> (nbits - 1 - k)) & 1;
    }

    for i in (deg + 1)..=maxbit {
        let mut vi = v[i - deg];
        vi ^= v[i - deg] >> mdeg;
        for (k, &bit) in bits.iter().enumerate() {
            if bit != 0 {
                vi ^= v[i - (k + 1)];
            }
        }
        v[i] = vi;
    }

    v
}

/// Generates successive points of a (unscrambled) Sobol low-discrepancy
/// sequence in `[0, 1)^dim`, using the Gray-code construction (Antonov &
/// Saleev), for dimensions `0..MAX_SOBOL_DIM`; dimensions beyond that use a
/// Halton-sequence fallback (see the module docs).
pub struct SobolGenerator {
    dim: usize,
    /// Direction numbers per genuinely-Sobol dimension (length `MAXBIT + 1`,
    /// index 0 unused); empty for fallback (Halton) dimensions.
    direction_numbers: Vec<Vec<u32>>,
    /// Current integer numerator (over `2^MAXBIT`) per dimension.
    x: Vec<u32>,
    /// Number of points already emitted.
    count: usize,
}

impl SobolGenerator {
    /// Creates a generator for `dim`-dimensional points.
    pub fn new(dim: usize) -> Self {
        let direction_numbers = (0..dim)
            .map(|d| {
                if d == 0 {
                    // Dimension 0 is the trivial base-2 van der Corput case:
                    // v_i = 2^{-i}.
                    let mut v = vec![0u32; MAXBIT as usize + 1];
                    for (i, slot) in v.iter_mut().enumerate().skip(1) {
                        *slot = 1u32 << (MAXBIT - i as u32);
                    }
                    v
                } else if d < MAX_SOBOL_DIM {
                    let (mdeg, ip, iv) = SOBOL_TABLE[d - 1];
                    build_direction_numbers(mdeg, ip, iv)
                } else {
                    Vec::new() // sentinel: this dimension uses the Halton fallback
                }
            })
            .collect();

        Self {
            dim,
            direction_numbers,
            x: vec![0u32; dim],
            count: 0,
        }
    }

    /// Returns the next point in `[0, 1)^dim`.
    pub fn next_point(&mut self) -> Vec<f64> {
        let n = self.count;
        self.count += 1;

        let point: Vec<f64> = (0..self.dim)
            .map(|d| {
                if d < MAX_SOBOL_DIM {
                    self.x[d] as f64 / (1u64 << MAXBIT) as f64
                } else {
                    let base = PRIMES[d % PRIMES.len()];
                    halton_radical_inverse(n + 1, base)
                }
            })
            .collect();

        // Gray-code update: XOR in the direction number at the index of the
        // lowest zero bit of `n` (1-indexed as `c`).
        let c = (n as u32).trailing_ones() as usize + 1;
        for d in 0..self.dim.min(MAX_SOBOL_DIM) {
            if let Some(&vc) = self.direction_numbers[d].get(c) {
                self.x[d] ^= vc;
            }
        }

        point
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_halton_radical_inverse_base2_matches_van_der_corput() {
        // van der Corput base 2: 1/2, 1/4, 3/4, 1/8, 5/8, ...
        assert!((halton_radical_inverse(1, 2) - 0.5).abs() < 1e-12);
        assert!((halton_radical_inverse(2, 2) - 0.25).abs() < 1e-12);
        assert!((halton_radical_inverse(3, 2) - 0.75).abs() < 1e-12);
        assert!((halton_radical_inverse(4, 2) - 0.125).abs() < 1e-12);
    }

    #[test]
    fn test_halton_radical_inverse_stays_in_unit_interval() {
        for base in [2, 3, 5, 7, 11] {
            for index in 1..100 {
                let v = halton_radical_inverse(index, base);
                assert!((0.0..1.0).contains(&v), "base={base} index={index} v={v}");
            }
        }
    }

    /// Reference Sobol points from `scipy.stats.qmc.Sobol(d=7,
    /// scramble=False).random(8)`, verified bit-for-bit (max abs error
    /// 0.0 across 64 points x 7 dims during development) against this
    /// module's construction.
    #[test]
    fn test_sobol_matches_scipy_reference_7d() {
        let expected: [[f64; 7]; 8] = [
            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
            [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
            [0.75, 0.25, 0.25, 0.25, 0.75, 0.75, 0.25],
            [0.25, 0.75, 0.75, 0.75, 0.25, 0.25, 0.75],
            [0.375, 0.375, 0.625, 0.875, 0.375, 0.125, 0.375],
            [0.875, 0.875, 0.125, 0.375, 0.875, 0.625, 0.875],
            [0.625, 0.125, 0.875, 0.625, 0.625, 0.875, 0.125],
            [0.125, 0.625, 0.375, 0.125, 0.125, 0.375, 0.625],
        ];

        let mut gen = SobolGenerator::new(7);
        for row in expected.iter() {
            let p = gen.next_point();
            for (got, &want) in p.iter().zip(row.iter()) {
                assert!((got - want).abs() < 1e-12, "got {got}, want {want}");
            }
        }
    }

    #[test]
    fn test_sobol_points_stay_in_unit_interval_and_are_non_constant() {
        let mut gen = SobolGenerator::new(5);
        let mut first_coords = Vec::new();
        for _ in 0..50 {
            let p = gen.next_point();
            assert_eq!(p.len(), 5);
            for &v in &p {
                assert!((0.0..1.0).contains(&v));
            }
            first_coords.push(p[0]);
        }
        // Non-constant data: guards against a degenerate stub that always
        // returns e.g. all zeros.
        assert!(first_coords.iter().any(|&v| v != first_coords[0]));
    }

    #[test]
    fn test_sobol_fallback_dimension_beyond_table_is_still_low_discrepancy() {
        // Dimension index MAX_SOBOL_DIM (0-indexed) is beyond the validated
        // Sobol table and must use the documented Halton fallback -- not
        // silently degrade to a constant or to unvalidated garbage.
        let dim = MAX_SOBOL_DIM + 2;
        let mut gen = SobolGenerator::new(dim);
        let mut last_col = Vec::new();
        for _ in 0..20 {
            let p = gen.next_point();
            for &v in &p {
                assert!((0.0..1.0).contains(&v));
            }
            last_col.push(p[dim - 1]);
        }
        assert!(last_col.iter().any(|&v| v != last_col[0]));
    }
}