somatize_study/sampler/
halton.rs

1//! Spread on purpose, one prime per knob.
2
3use super::Sampler;
4use super::drawing::{draw, splitmix};
5use crate::{Point, Space};
6
7/// Cover the space evenly instead of drawing from it evenly.
8///
9/// What separates it from [`Random`](super::Random) is not what it looks at —
10/// both look at nothing — but what it promises. Random is uniform *in
11/// expectation*; this is uniform *by construction, for every prefix*: of the
12/// first `base²` trials exactly one lands in each cell of a `base²` grid, and no
13/// arrangement of the indices makes it otherwise. For a study handed out of a
14/// folder that is the difference between a collision being unlikely and there
15/// being no way to arrange one.
16///
17/// Knob `d` is read in base the `d`-th prime, which is why the promise thins out
18/// once there are many knobs: the high primes need a long prefix before they
19/// look like anything. [`Sobol`](super::Sobol) has no such seam, at the price of
20/// a table. Its point is a function of the **seed and the index**.
21#[derive(Debug, Clone, PartialEq, Eq, Hash)]
22pub struct Halton {
23    /// The seed, which here permutes the digits rather than drawing them: with
24    /// no scramble a Halton sequence is one fixed sequence, and two studies of
25    /// the same space would search it in exactly the same order.
26    pub seed: u64,
27}
28
29impl Halton {
30    /// The `trial`-th point. It never runs out, and it never looks at what the
31    /// finished trials did.
32    pub fn ask(
33        &self,
34        space: &Space,
35        trial: usize,
36        _seen: &[(Point, Option<f64>)],
37    ) -> Option<Point> {
38        if space.is_empty() {
39            return None;
40        }
41        Some(Point::of(
42            space
43                .dimensions()
44                .iter()
45                .enumerate()
46                .map(|(which, (name, dimension))| {
47                    let u = radical(prime(which), trial as u64, self.seed, which);
48                    (name.clone(), draw(dimension, u))
49                })
50                .collect(),
51        ))
52    }
53}
54
55/// The index written in `base`, read back with its digits reversed and
56/// scrambled — a number in `0.0..1.0`. Reversing is the trick: consecutive
57/// indices differ in their last digit, which becomes the first.
58///
59/// **Every place is scrambled, including the zeroes the index never reaches.**
60/// Scrambling only the digits written down would put a one-digit index and a
61/// two-digit one on different grids, and the first `base²` trials would stop
62/// landing one per cell. It also settles trial zero, which has no digits.
63fn radical(base: u64, index: u64, seed: u64, which: usize) -> f64 {
64    let mut state = seed ^ (which as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
65    // Never zero, so `digit -> multiplier * digit + offset` is a permutation of
66    // the digits and not a mangling of them: the cover survives the seed. And it
67    // is one multiplier per knob, which is what pulls two large primes back out
68    // of step with each other.
69    let multiplier = 1 + splitmix(&mut state) % (base - 1);
70    let (mut left, mut place, mut value) = (index, 1.0 / base as f64, 0.0);
71    for _ in 0..places(base) {
72        // A fresh offset per place, so it is a permutation *of each place* and
73        // not one rotation of the whole number.
74        let offset = splitmix(&mut state) % base;
75        value += ((multiplier * (left % base) + offset) % base) as f64 * place;
76        left /= base;
77        place /= base as f64;
78    }
79    value
80}
81
82/// How many places are worth scrambling: the first count whose grid is finer
83/// than an `f64` tells apart. Past it every digit falls below the rounding.
84fn places(base: u64) -> u32 {
85    let (mut span, mut count) = (1u64, 0);
86    while span < 1 << 53 {
87        span *= base;
88        count += 1;
89    }
90    count
91}
92
93/// The `which`-th prime, counting from zero. Found rather than tabulated: there
94/// is no ceiling to write down, and a space has a handful of knobs.
95fn prime(which: usize) -> u64 {
96    let (mut found, mut candidate) = (0, 1u64);
97    loop {
98        candidate += 1;
99        if (2u64..)
100            .take_while(|by| by * by <= candidate)
101            .all(|by| candidate % by != 0)
102        {
103            if found == which {
104                return candidate;
105            }
106            found += 1;
107        }
108    }
109}
110
111impl From<Halton> for Sampler {
112    fn from(how: Halton) -> Self {
113        Self::Halton(how)
114    }
115}
somatize_study/sampler/halton.rs

somatize_study/sampler/
halton.rs
