libitofin 0.4.0

//! Lattice rules for low-discrepancy sequences, ported from
//! `ql/math/randomnumbers/latticerules.{hpp,cpp}` and
//! `ql/math/randomnumbers/latticersg.{hpp,cpp}`.
//!
//! The generating vectors are the Cools, Kuo, Nuyens rules
//! (SIAM J. Sci. Comp., 2006); the data tables live in
//! `lattice_tables.rs`, generated by `scripts/gen_sobol_tables.py`.

use super::lattice_tables;

/// Rank-1 lattice rule generating vectors.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum LatticeRule {
    A,
    B,
    C,
    D,
}

impl LatticeRule {
    /// Returns the generating vector of the rule, like QuantLib's `getRule`.
    ///
    /// # Panics
    ///
    /// Panics if `n` is outside the range QuantLib requires for these rules.
    pub fn get_rule(self, n: i64) -> Vec<f64> {
        assert!(
            n >= 1024 && (n as f64) <= 2.9f64.powi(20),
            "N must be between 2 to 10 and 2 to the 20 for these lattice rules"
        );
        let table = match self {
            LatticeRule::A => lattice_tables::LATTICE_A,
            LatticeRule::B => lattice_tables::LATTICE_B,
            LatticeRule::C => lattice_tables::LATTICE_C,
            LatticeRule::D => lattice_tables::LATTICE_D,
        };
        table[..3600].to_vec()
    }
}

/// Rank-1 lattice rule sequence generator.
pub struct LatticeRsg {
    dimensionality: usize,
    n: usize,
    i: usize,
    z: Vec<f64>,
    sequence: Vec<f64>,
}

impl LatticeRsg {
    /// Creates a generator over `n` points using the generating vector `z`.
    ///
    /// # Panics
    ///
    /// Panics if `z` has fewer than `dimensionality` components; QuantLib
    /// reads out of bounds there.
    pub fn new(dimensionality: usize, z: Vec<f64>, n: usize) -> Self {
        assert!(
            z.len() >= dimensionality,
            "generating vector has {} components but the dimensionality is {dimensionality}",
            z.len()
        );
        Self {
            dimensionality,
            n,
            i: 0,
            z,
            sequence: vec![0.0; dimensionality],
        }
    }

    /// Skips `n` samples in the low-discrepancy sequence.
    pub fn skip_to(&mut self, n: usize) {
        self.i += n;
    }

    /// Returns the next sample in the sequence.
    pub fn next_sequence(&mut self) -> &[f64] {
        for (value, z) in self.sequence.iter_mut().zip(&self.z) {
            let theta = self.i as f64 * z / self.n as f64;
            *value = theta % 1.0;
        }
        self.i += 1;
        &self.sequence
    }

    /// Returns the last sample in the sequence.
    pub fn last_sequence(&self) -> &[f64] {
        &self.sequence
    }

    /// Returns the dimensionality of the sequence.
    pub fn dimension(&self) -> usize {
        self.dimensionality
    }
}

#[cfg(test)]
mod tests {
    use super::{LatticeRsg, LatticeRule};

    #[test]
    fn get_rule_returns_generating_vectors() {
        let cases = [
            (
                LatticeRule::A,
                [1.0, 182667.0, 469891.0, 498753.0, 110745.0],
            ),
            (LatticeRule::B, [1.0, 182667.0, 213731.0, 255351.0, 96013.0]),
            (
                LatticeRule::C,
                [1.0, 433461.0, 103659.0, 481853.0, 186513.0],
            ),
            (
                LatticeRule::D,
                [1.0, 182667.0, 279195.0, 223491.0, 205755.0],
            ),
        ];
        for (rule, head) in cases {
            let z = rule.get_rule(1024);
            assert_eq!(z.len(), 3600, "{rule:?}");
            assert_eq!(&z[..5], &head, "{rule:?}");
        }
    }

    #[test]
    #[should_panic(expected = "N must be between")]
    fn get_rule_rejects_small_n() {
        LatticeRule::A.get_rule(512);
    }

    #[test]
    #[should_panic(expected = "generating vector has")]
    fn new_rejects_short_generating_vector() {
        LatticeRsg::new(5, vec![1.0, 182667.0], 1024);
    }

    #[test]
    fn sequence_is_scaled_lattice_points_modulo_one() {
        let z = LatticeRule::A.get_rule(1024);
        let mut rsg = LatticeRsg::new(3, z, 1024);
        assert_eq!(rsg.dimension(), 3);
        assert_eq!(rsg.next_sequence(), &[0.0, 0.0, 0.0]);
        let second = rsg.next_sequence();
        assert_eq!(second, &[0.0009765625, 0.3857421875, 0.8779296875]);
        assert_eq!(
            rsg.last_sequence(),
            &[0.0009765625, 0.3857421875, 0.8779296875]
        );
    }

    #[test]
    fn skipping_matches_drawing() {
        let z = LatticeRule::D.get_rule(1024);
        let mut rsg1 = LatticeRsg::new(10, z.clone(), 1024);
        for _ in 0..7 {
            rsg1.next_sequence();
        }
        let mut rsg2 = LatticeRsg::new(10, z, 1024);
        rsg2.skip_to(7);
        for _ in 0..20 {
            let s1 = rsg1.next_sequence().to_vec();
            let s2 = rsg2.next_sequence();
            assert_eq!(s1, s2);
        }
    }
}