fcmaes-core 0.1.3

//! Random number generation for the optimizers.
//!
//! Wraps `rand_pcg::Pcg64` behind the small set of scalar and vector
//! distributions required by the optimizers.
//!
//! # Examples
//!
//! Every optimizer takes an [`Rng`], so seeding it is what makes a run
//! reproducible:
//!
//! ```
//! use fcmaes_core::Rng;
//!
//! let mut first = Rng::new(42);
//! let mut second = Rng::new(42);
//! assert_eq!(first.uniform01(), second.uniform01());
//!
//! // A different seed gives an independent stream.
//! let mut other = Rng::new(43);
//! assert_ne!(Rng::new(42).uniform01(), other.uniform01());
//! ```
//!
//! # Reference
//!
//! M. E. O'Neill, [“PCG: A Family of Simple Fast Space-Efficient Statistically
//! Good Algorithms for Random Number
//! Generation”](https://www.pcg-random.org/paper.html) (2014).

use rand::{Rng as _, SeedableRng};
use rand_distr::StandardNormal;
use rand_pcg::Pcg64;

/// Deterministic, seedable RNG shared by all fcmaes optimizers.
#[derive(Clone, Debug)]
pub struct Rng {
    inner: Pcg64,
}

impl Rng {
    /// Seed from a single 64-bit value (convenience for `seed`/`runid` pairs).
    pub fn new(seed: u64) -> Self {
        Self {
            inner: Pcg64::seed_from_u64(seed),
        }
    }

    /// Seed from a full 128-bit state + stream selector.
    pub fn from_state_stream(state: u128, stream: u128) -> Self {
        Self {
            inner: Pcg64::new(state, stream),
        }
    }

    /// Draw a uniform value in `[0, 1)`.
    #[inline]
    pub fn uniform01(&mut self) -> f64 {
        self.inner.r#gen::<f64>()
    }

    /// Draw a full-width seed for an independent child operation.
    #[inline]
    pub fn next_u64(&mut self) -> u64 {
        self.inner.r#gen::<u64>()
    }

    /// Standard normal sample `N(0, 1)`.
    #[inline]
    pub fn gaussian(&mut self) -> f64 {
        self.inner.sample(StandardNormal)
    }

    /// Draw a normal value with mean `mu` and standard deviation `sdev`.
    #[inline]
    pub fn normreal(&mut self, mu: f64, sdev: f64) -> f64 {
        self.gaussian() * sdev + mu
    }

    /// Draw an integer in `[0, max)` by scaling a uniform value.
    ///
    /// Returns zero when `max <= 0`.
    #[inline]
    pub fn int_below(&mut self, max: i64) -> i64 {
        if max <= 0 {
            return 0;
        }
        (max as f64 * self.uniform01()) as i64
    }

    /// Vector of `dim` uniform `[0, 1)` samples.
    pub fn uniform_vec(&mut self, dim: usize) -> Vec<f64> {
        (0..dim).map(|_| self.uniform01()).collect()
    }

    /// Vector of `dim` standard-normal samples.
    pub fn normal_vec(&mut self, dim: usize) -> Vec<f64> {
        (0..dim).map(|_| self.gaussian()).collect()
    }
}

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

    #[test]
    fn same_seed_reproduces_stream() {
        let mut a = Rng::new(12345);
        let mut b = Rng::new(12345);
        for _ in 0..100 {
            assert_eq!(a.uniform01(), b.uniform01());
        }
    }

    #[test]
    fn different_seed_diverges() {
        let mut a = Rng::new(1);
        let mut b = Rng::new(2);
        // Extremely unlikely the first draws coincide.
        assert_ne!(a.uniform01(), b.uniform01());
    }

    #[test]
    fn full_width_seeds_reproduce_and_advance() {
        let mut a = Rng::new(42);
        let mut b = Rng::new(42);
        let first = a.next_u64();
        assert_eq!(first, b.next_u64());
        assert_ne!(first, a.next_u64());
    }

    #[test]
    fn uniform01_in_range() {
        let mut rng = Rng::new(7);
        for _ in 0..10_000 {
            let u = rng.uniform01();
            assert!((0.0..1.0).contains(&u));
        }
    }

    #[test]
    fn int_below_in_range() {
        let mut rng = Rng::new(7);
        for _ in 0..10_000 {
            let v = rng.int_below(5);
            assert!((0..5).contains(&v));
        }
        assert_eq!(rng.int_below(0), 0);
        assert_eq!(rng.int_below(-3), 0);
    }

    #[test]
    fn gaussian_statistics_are_sane() {
        let mut rng = Rng::new(2024);
        let n = 200_000;
        let mut sum = 0.0;
        let mut sumsq = 0.0;
        for _ in 0..n {
            let g = rng.gaussian();
            sum += g;
            sumsq += g * g;
        }
        let mean = sum / n as f64;
        let var = sumsq / n as f64 - mean * mean;
        assert!(mean.abs() < 0.02, "mean={mean}");
        assert!((var - 1.0).abs() < 0.03, "var={var}");
    }

    #[test]
    fn state_stream_normal_and_vector_helpers() {
        let mut first = Rng::from_state_stream(123, 7);
        let mut second = Rng::from_state_stream(123, 7);
        assert_eq!(first.uniform_vec(8), second.uniform_vec(8));
        let normal = first.normal_vec(16);
        assert_eq!(normal.len(), 16);
        assert!(normal.iter().all(|value| value.is_finite()));
        let sample = first.normreal(10.0, 0.0);
        assert_eq!(sample, 10.0);
        assert!(first.uniform_vec(0).is_empty());
        assert!(first.normal_vec(0).is_empty());
    }
}