axiolid_predicates/scene.rs
1//! Degeneracy-controlled scene generation for benchmarking and testing.
2//!
3//! The question a predicate suite must answer is not "how fast is it on random
4//! data" -- random points in a box are never degenerate, so that measures only
5//! the filter's fast path. The useful question is how throughput and the
6//! escalation rate move as the fraction of degenerate inputs rises.
7//!
8//! Scenes are deterministic given a seed, so a measurement is reproducible.
9
10use axiolid_core::{Point2, Point3};
11
12/// Fraction of inputs that are exactly degenerate.
13///
14/// The tiers span the range that matters: clean data, the incidental
15/// degeneracies of authored CAD models, the systematic ones of grid-aligned
16/// architecture, and a stress tier well past anything real.
17#[derive(Debug, Clone, Copy, PartialEq, Eq)]
18pub enum DegeneracyRate {
19 /// No degenerate inputs.
20 None,
21 /// One in ten thousand.
22 Rare,
23 /// One in a hundred.
24 Occasional,
25 /// One in ten.
26 Frequent,
27}
28
29impl DegeneracyRate {
30 /// Every tier, ascending.
31 pub const ALL: [Self; 4] = [Self::None, Self::Rare, Self::Occasional, Self::Frequent];
32
33 /// The fraction as a ratio.
34 #[must_use]
35 pub const fn fraction(self) -> f64 {
36 match self {
37 Self::None => 0.0,
38 Self::Rare => 0.0001,
39 Self::Occasional => 0.01,
40 Self::Frequent => 0.1,
41 }
42 }
43
44 /// Human-readable label for reports.
45 #[must_use]
46 pub const fn label(self) -> &'static str {
47 match self {
48 Self::None => "0%",
49 Self::Rare => "0.01%",
50 Self::Occasional => "1%",
51 Self::Frequent => "10%",
52 }
53 }
54
55 /// Whether the sample at `index` should be degenerate.
56 ///
57 /// Deterministic striping rather than a random draw, so a given index is
58 /// degenerate in every run and a regression is reproducible.
59 #[must_use]
60 pub fn is_degenerate(self, index: usize) -> bool {
61 match self {
62 Self::None => false,
63 Self::Rare => index % 10_000 == 0,
64 Self::Occasional => index % 100 == 0,
65 Self::Frequent => index % 10 == 0,
66 }
67 }
68}
69
70/// Deterministic xorshift, so scenes are reproducible across runs and machines.
71#[derive(Debug, Clone)]
72pub struct SceneRng {
73 state: u64,
74}
75
76impl SceneRng {
77 /// Seed the generator. Zero is remapped: xorshift is stuck at zero.
78 #[must_use]
79 pub const fn new(seed: u64) -> Self {
80 Self {
81 state: if seed == 0 {
82 0x2545_F491_4F6C_DD1D
83 } else {
84 seed
85 },
86 }
87 }
88
89 /// Next raw value.
90 pub fn next_u64(&mut self) -> u64 {
91 self.state ^= self.state << 13;
92 self.state ^= self.state >> 7;
93 self.state ^= self.state << 17;
94 self.state
95 }
96
97 /// Integer coordinate in `[-bound, bound)`.
98 ///
99 /// Integers, not arbitrary floats: they keep the degenerate constructions
100 /// below exactly degenerate, which is the entire point of the scene.
101 pub fn coordinate(&mut self, bound: i64) -> f64 {
102 (((self.next_u64() >> 33) as i64) % (2 * bound) - bound) as f64
103 }
104}
105
106/// A 2D orientation query: three points.
107pub type Orient2Case = [Point2; 3];
108
109/// A 3D orientation query: four points.
110pub type Orient3Case = [Point3; 4];
111
112/// Generate `count` orientation queries with the given degeneracy rate.
113///
114/// Degenerate cases are exactly collinear by construction (`c` is an integer
115/// affine combination of `a` and `b`), so they are not merely close to the
116/// boundary -- they are on it, which is what forces the exact path.
117#[must_use]
118pub fn orient2_scene(count: usize, rate: DegeneracyRate, seed: u64) -> Vec<Orient2Case> {
119 let mut rng = SceneRng::new(seed);
120 (0..count)
121 .map(|index| {
122 let a = Point2::new(rng.coordinate(1_000), rng.coordinate(1_000));
123 let b = Point2::new(rng.coordinate(1_000), rng.coordinate(1_000));
124 let c = if rate.is_degenerate(index) {
125 let k = ((rng.next_u64() >> 40) as i64 % 7) - 3;
126 Point2::new(
127 a.x + (k as f64) * (b.x - a.x),
128 a.y + (k as f64) * (b.y - a.y),
129 )
130 } else {
131 Point2::new(rng.coordinate(1_000), rng.coordinate(1_000))
132 };
133 [a, b, c]
134 })
135 .collect()
136}
137
138/// Generate `count` 3D orientation queries with the given degeneracy rate.
139///
140/// Degenerate cases are exactly coplanar: `d = a + i*(b-a) + j*(c-a)` for
141/// integers `i, j`, so no rounding can move the point off the plane.
142#[must_use]
143pub fn orient3_scene(count: usize, rate: DegeneracyRate, seed: u64) -> Vec<Orient3Case> {
144 let mut rng = SceneRng::new(seed);
145 (0..count)
146 .map(|index| {
147 let p = |r: &mut SceneRng| {
148 Point3::new(
149 r.coordinate(1_000),
150 r.coordinate(1_000),
151 r.coordinate(1_000),
152 )
153 };
154 let (a, b, c) = (p(&mut rng), p(&mut rng), p(&mut rng));
155 let d = if rate.is_degenerate(index) {
156 let i = (((rng.next_u64() >> 40) as i64 % 5) - 2) as f64;
157 let j = (((rng.next_u64() >> 40) as i64 % 5) - 2) as f64;
158 Point3::new(
159 a.x + i * (b.x - a.x) + j * (c.x - a.x),
160 a.y + i * (b.y - a.y) + j * (c.y - a.y),
161 a.z + i * (b.z - a.z) + j * (c.z - a.z),
162 )
163 } else {
164 p(&mut rng)
165 };
166 [a, b, c, d]
167 })
168 .collect()
169}
170
171/// Generate 3D orientation queries that are *nearly* coplanar.
172///
173/// Distinct from [`orient3_scene`]'s degenerate cases, which lie exactly on the
174/// plane and are therefore answered zero. These sit one unit-in-the-last-place
175/// off it: the filter must defer, and the exact path must then return a
176/// non-zero sign. A suite containing only exact degeneracies cannot tell a
177/// correct exact path from one that always answers zero.
178#[must_use]
179pub fn near_coplanar_scene(count: usize, seed: u64) -> Vec<Orient3Case> {
180 let mut rng = SceneRng::new(seed);
181 (0..count)
182 .map(|_| {
183 let p = |r: &mut SceneRng| {
184 Point3::new(
185 r.coordinate(1_000),
186 r.coordinate(1_000),
187 r.coordinate(1_000),
188 )
189 };
190 let (a, b, c) = (p(&mut rng), p(&mut rng), p(&mut rng));
191 let i = (((rng.next_u64() >> 40) as i64 % 5) - 2) as f64;
192 let j = (((rng.next_u64() >> 40) as i64 % 5) - 2) as f64;
193 let z = a.z + i * (b.z - a.z) + j * (c.z - a.z);
194 [
195 a,
196 b,
197 c,
198 Point3::new(
199 a.x + i * (b.x - a.x) + j * (c.x - a.x),
200 a.y + i * (b.y - a.y) + j * (c.y - a.y),
201 // One ULP off the plane: far too small for the filter to
202 // resolve, but a definite side that exact arithmetic must
203 // recover.
204 next_ulp(z),
205 ),
206 ]
207 })
208 .collect()
209}
210
211/// The next representable f64 above `value`.
212///
213/// Bit increment rather than `next_after`, which core does not expose. The
214/// step is exactly one ULP, the smallest perturbation that still moves the
215/// point off the plane.
216#[must_use]
217fn next_ulp(value: f64) -> f64 {
218 if value.is_nan() || value == f64::INFINITY {
219 return value;
220 }
221 let bits = value.to_bits();
222 f64::from_bits(if value >= 0.0 { bits + 1 } else { bits - 1 })
223}