sciforge_lib/maths/probability/

sampling.rs

pub fn latin_hypercube(n_samples: usize, n_dims: usize, seed: u64) -> Vec<Vec<f64>> {
    let mut rng = SampRng::new(seed);
    let mut result = vec![vec![0.0; n_dims]; n_samples];
    for d in 0..n_dims {
        let mut perm: Vec<usize> = (0..n_samples).collect();
        for i in (1..n_samples).rev() {
            let j = rng.next_usize(i + 1);
            perm.swap(i, j);
        }
        for (ri, &pi) in result.iter_mut().zip(&perm) {
            ri[d] = (pi as f64 + rng.next_f64()) / n_samples as f64;
        }
    }
    result
}

pub fn stratified_sampling(n_strata: usize, seed: u64) -> Vec<f64> {
    let mut rng = SampRng::new(seed);
    (0..n_strata)
        .map(|i| (i as f64 + rng.next_f64()) / n_strata as f64)
        .collect()
}

pub fn sobol_sequence_1d(n: usize) -> Vec<f64> {
    let mut result = vec![0.0; n];
    for i in 1..n {
        let mut c = 0u32;
        let mut v = i;
        while v & 1 == 0 {
            c += 1;
            v >>= 1;
        }
        let direction = 1u64 << (63 - c);
        let prev_bits = (result[i - 1] * (1u64 << 63) as f64) as u64;
        result[i] = (prev_bits ^ direction) as f64 / (1u64 << 63) as f64;
    }
    result
}

pub fn halton_sequence(n: usize, base: u64) -> Vec<f64> {
    (0..n).map(|i| halton_element(i as u64, base)).collect()
}

fn halton_element(mut index: u64, base: u64) -> f64 {
    let mut result = 0.0;
    let mut f = 1.0 / base as f64;
    index += 1;
    while index > 0 {
        result += f * (index % base) as f64;
        index /= base;
        f /= base as f64;
    }
    result
}

pub fn inverse_transform_exponential(u: f64, lambda: f64) -> f64 {
    -(1.0 - u).ln() / lambda
}

pub fn box_muller(u1: f64, u2: f64) -> (f64, f64) {
    let r = (-2.0 * u1.ln()).sqrt();
    let theta = 2.0 * std::f64::consts::PI * u2;
    (r * theta.cos(), r * theta.sin())
}

pub fn reservoir_sampling(stream: &[f64], k: usize, seed: u64) -> Vec<f64> {
    let mut rng = SampRng::new(seed);
    let mut reservoir = stream[..k.min(stream.len())].to_vec();
    for (i, &si) in stream.iter().enumerate().skip(k) {
        let j = rng.next_usize(i + 1);
        if j < k {
            reservoir[j] = si;
        }
    }
    reservoir
}

struct SampRng {
    state: u64,
}
impl SampRng {
    fn new(seed: u64) -> Self {
        Self { state: seed }
    }
    fn next_u64(&mut self) -> u64 {
        self.state = self
            .state
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        self.state
    }
    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }
    fn next_usize(&mut self, max: usize) -> usize {
        (self.next_f64() * max as f64) as usize % max
    }
}

pub fn alias_table_build(probs: &[f64]) -> (Vec<f64>, Vec<usize>) {
    let n = probs.len();
    let avg = 1.0 / n as f64;
    let mut prob = vec![0.0; n];
    let mut alias = vec![0; n];
    let mut small = Vec::new();
    let mut large = Vec::new();
    let mut scaled: Vec<f64> = probs.iter().map(|p| p * n as f64).collect();
    for (i, &si) in scaled.iter().enumerate() {
        if si < 1.0 {
            small.push(i);
        } else {
            large.push(i);
        }
    }
    while let (Some(s), Some(l)) = (small.pop(), large.pop()) {
        prob[s] = scaled[s];
        alias[s] = l;
        scaled[l] = scaled[l] + scaled[s] - 1.0;
        if scaled[l] < 1.0 {
            small.push(l);
        } else {
            large.push(l);
        }
    }
    for &l in &large {
        prob[l] = 1.0;
    }
    for &s in &small {
        prob[s] = 1.0;
    }
    let _ = avg;
    (prob, alias)
}

pub fn alias_sample(prob: &[f64], alias: &[usize], seed: u64) -> usize {
    let mut rng = SampRng::new(seed);
    let n = prob.len();
    let i = rng.next_usize(n);
    if rng.next_f64() < prob[i] {
        i
    } else {
        alias[i]
    }
}

pub fn systematic_sampling(n_samples: usize, seed: u64) -> Vec<f64> {
    let mut rng = SampRng::new(seed);
    let u0 = rng.next_f64() / n_samples as f64;
    (0..n_samples)
        .map(|i| u0 + i as f64 / n_samples as f64)
        .collect()
}

pub fn importance_resampling(weights: &[f64], n_samples: usize, seed: u64) -> Vec<usize> {
    let mut rng = SampRng::new(seed);
    let total: f64 = weights.iter().sum();
    let normalized: Vec<f64> = weights.iter().map(|w| w / total).collect();
    let mut cumul = vec![0.0; normalized.len()];
    cumul[0] = normalized[0];
    for i in 1..normalized.len() {
        cumul[i] = cumul[i - 1] + normalized[i];
    }
    let mut indices = Vec::with_capacity(n_samples);
    for _ in 0..n_samples {
        let u = rng.next_f64();
        let idx = cumul.partition_point(|&c| c < u);
        indices.push(idx.min(weights.len() - 1));
    }
    indices
}

pub fn van_der_corput(n: usize, base: u64) -> Vec<f64> {
    (0..n)
        .map(|i| {
            let mut result = 0.0;
            let mut f = 1.0 / base as f64;
            let mut idx = (i + 1) as u64;
            while idx > 0 {
                result += f * (idx % base) as f64;
                idx /= base;
                f /= base as f64;
            }
            result
        })
        .collect()
}

pub fn hammersley_sequence(n: usize, base: u64) -> Vec<(f64, f64)> {
    let vdc = van_der_corput(n, base);
    (0..n).map(|i| (i as f64 / n as f64, vdc[i])).collect()
}

pub fn weighted_sampling(items: &[f64], weights: &[f64], n_samples: usize, seed: u64) -> Vec<f64> {
    let mut rng = SampRng::new(seed);
    let total: f64 = weights.iter().sum();
    let mut cumul = Vec::with_capacity(weights.len());
    let mut s = 0.0;
    for w in weights {
        s += w / total;
        cumul.push(s);
    }
    let mut samples = Vec::with_capacity(n_samples);
    for _ in 0..n_samples {
        let u = rng.next_f64();
        let idx = cumul.partition_point(|&c| c < u).min(items.len() - 1);
        samples.push(items[idx]);
    }
    samples
}

pub fn poisson_disk_sampling_1d(
    domain_min: f64,
    domain_max: f64,
    min_dist: f64,
    seed: u64,
) -> Vec<f64> {
    let mut rng = SampRng::new(seed);
    let mut points = Vec::new();
    let mut candidate = domain_min + rng.next_f64() * (domain_max - domain_min);
    points.push(candidate);
    let max_attempts = 30;
    let mut active = vec![0usize];
    while !active.is_empty() {
        let idx = rng.next_usize(active.len());
        let center = points[active[idx]];
        let mut found = false;
        for _ in 0..max_attempts {
            let offset = min_dist + rng.next_f64() * min_dist;
            candidate = if rng.next_f64() < 0.5 {
                center + offset
            } else {
                center - offset
            };
            if candidate < domain_min || candidate > domain_max {
                continue;
            }
            let ok = points.iter().all(|&p| (p - candidate).abs() >= min_dist);
            if ok {
                active.push(points.len());
                points.push(candidate);
                found = true;
                break;
            }
        }
        if !found {
            active.swap_remove(idx);
        }
    }
    points
}