Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/tests/cell.rs

20.5 KiB, 1 run

created by r1870400018:44202, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! Closed-form correctness oracle for the cube-sphere quadtree cell grid.
2//!
3//! There is no external reference implementation: these invariants ARE the specification.
4//! Every check is offline and closed-form, so the suite is deterministic and hermetic.
5
6use oxedyne_fe2o3_geom::cell::{Cell, Corner};
7
8use oxedyne_fe2o3_core::prelude::*;
9
10use std::collections::{HashMap, HashSet};
11
12// ---------------------------------------------------------------------------------------------
13// Test-local deterministic RNG (SplitMix64) -- no dev-dependency introduced.
14// ---------------------------------------------------------------------------------------------
15
16struct Rng(u64);
17
18impl Rng {
19 fn new(seed: u64) -> Self { Rng(seed) }
20
21 fn next_u64(&mut self) -> u64 {
22 self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
23 let mut z = self.0;
24 z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
25 z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
26 z ^ (z >> 31)
27 }
28
29 fn next_f64(&mut self) -> f64 { // [0, 1)
30 (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64)
31 }
32
33 fn lat(&mut self) -> f64 { self.next_f64() * 180.0 - 90.0 }
34 fn lon(&mut self) -> f64 { self.next_f64() * 360.0 - 180.0 }
35
36 fn cell_at(&mut self, level: u8) -> Cell {
37 let n: u32 = 1u32 << level;
38 let face = (self.next_u64() % 6) as u8;
39 let i = (self.next_u64() % n as u64) as u32;
40 let j = (self.next_u64() % n as u64) as u32;
41 match Cell::from_face_ij(face, level, i, j) {
42 Ok(c) => c,
43 Err(e) => panic!("from_face_ij failed for a valid cell: {}", e),
44 }
45 }
46}
47
48// ---------------------------------------------------------------------------------------------
49// Independent geometric helpers used as oracles (deliberately NOT the production code path).
50// ---------------------------------------------------------------------------------------------
51
52fn dot(a: &[f64; 3], b: &[f64; 3]) -> f64 { a[0]*b[0] + a[1]*b[1] + a[2]*b[2] }
53
54fn cross(a: &[f64; 3], b: &[f64; 3]) -> [f64; 3] {
55 [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
56}
57
58fn norm(v: &[f64; 3]) -> [f64; 3] {
59 let len = dot(v, v).sqrt();
60 [v[0]/len, v[1]/len, v[2]/len]
61}
62
63fn latlon_to_vec(lat: f64, lon: f64) -> [f64; 3] {
64 let (la, lo) = (lat.to_radians(), lon.to_radians());
65 [la.cos()*lo.cos(), la.cos()*lo.sin(), la.sin()]
66}
67
68/// Independent point-in-cell test using the four great-circle edges of the cell.
69///
70/// Each cell edge is the intersection of a plane through the origin with the sphere, so a
71/// point is inside when it lies on the interior side of all four edge planes. This is a
72/// wholly separate implementation from `Cell::contains`, used only to cross-check it.
73fn naive_contains(cell: &Cell, p: &[f64; 3]) -> bool {
74 let c: [Corner; 4] = cell.corners();
75 let centre = cell.centre_vec();
76 let vs = [&c[0].vec, &c[1].vec, &c[2].vec, &c[3].vec];
77 for e in 0..4 {
78 let a = vs[e];
79 let b = vs[(e + 1) % 4];
80 let n = cross(a, b); // edge-plane normal
81 // Orient the normal so the cell centre is on the positive side.
82 let s = if dot(&n, &centre) >= 0.0 { 1.0 } else { -1.0 };
83 if s * dot(&n, p) < 0.0 {
84 return false;
85 }
86 }
87 true
88}
89
90/// Spherical-excess area of a spherical triangle of unit vectors (Van Oosterom & Strackee).
91fn tri_area(a: &[f64; 3], b: &[f64; 3], c: &[f64; 3]) -> f64 {
92 let triple = dot(a, &cross(b, c)).abs();
93 let den = 1.0 + dot(a, b) + dot(b, c) + dot(c, a);
94 2.0 * triple.atan2(den)
95}
96
97fn cell_area_oracle(cell: &Cell) -> f64 {
98 let c = cell.corners();
99 tri_area(&c[0].vec, &c[1].vec, &c[2].vec) + tri_area(&c[0].vec, &c[2].vec, &c[3].vec)
100}
101
102fn all_cells(level: u8) -> Vec<Cell> {
103 let n: u32 = 1u32 << level;
104 let mut out = Vec::new();
105 for face in 0u8..6 {
106 for i in 0..n {
107 for j in 0..n {
108 match Cell::from_face_ij(face, level, i, j) {
109 Ok(c) => out.push(c),
110 Err(e) => panic!("enumeration built an invalid cell: {}", e),
111 }
112 }
113 }
114 }
115 out
116}
117
118fn same_vec(a: &[f64; 3], b: &[f64; 3]) -> bool { dot(a, b) > 1.0 - 1e-12 }
119
120fn shared_corner_count(a: &Cell, b: &Cell) -> usize {
121 let ca = a.corners();
122 let cb = b.corners();
123 let mut n = 0;
124 for x in &ca {
125 for y in &cb {
126 if same_vec(&x.vec, &y.vec) {
127 n += 1;
128 break;
129 }
130 }
131 }
132 n
133}
134
135// ---------------------------------------------------------------------------------------------
136// Invariant 1: centre round-trip is exact at every level.
137// ---------------------------------------------------------------------------------------------
138
139#[test]
140fn inv01_centre_round_trip() {
141 let mut rng = Rng::new(0x1234_5678);
142 for level in 0u8..=26 {
143 for _ in 0..200 {
144 let cell = rng.cell_at(level);
145 let (lat, lon) = cell.centre();
146 let got = match Cell::at(lat, lon, level) {
147 Ok(c) => c,
148 Err(e) => panic!("at() failed on a cell centre: {}", e),
149 };
150 assert_eq!(got, cell, "centre round-trip broke at level {}", level);
151 }
152 }
153}
154
155// ---------------------------------------------------------------------------------------------
156// Invariant 2: parent is exact -- parent(at(p, b), a) == at(p, a) for all a < b.
157// ---------------------------------------------------------------------------------------------
158
159#[test]
160fn inv02_parent_exact() {
161 let mut rng = Rng::new(0xDEAD_BEEF);
162 for _ in 0..4000 {
163 let lat = rng.lat();
164 let lon = rng.lon();
165 let b = (rng.next_u64() % 27) as u8;
166 if b == 0 { continue; }
167 let a = (rng.next_u64() % b as u64) as u8; // a < b
168 let deep = match Cell::at(lat, lon, b) { Ok(c) => c, Err(_) => continue };
169 let via_parent = match deep.parent(a) { Ok(c) => c, Err(e) => panic!("parent: {}", e) };
170 let direct = match Cell::at(lat, lon, a) { Ok(c) => c, Err(e) => panic!("at: {}", e) };
171 assert_eq!(via_parent, direct, "parent inexact: a={} b={}", a, b);
172 }
173}
174
175// ---------------------------------------------------------------------------------------------
176// Invariant 3: children -- exactly 4, each parents back, (i,j) partition the parent's cell.
177// ---------------------------------------------------------------------------------------------
178
179#[test]
180fn inv03_children_partition() {
181 let mut rng = Rng::new(0x0BADF00D);
182 for _ in 0..2000 {
183 let level = (rng.next_u64() % 26) as u8; // 0..=25 so children exist
184 let parent = rng.cell_at(level);
185 let kids = match parent.children() { Ok(k) => k, Err(e) => panic!("children: {}", e) };
186 let (pi, pj) = parent.coords();
187 let mut seen = HashSet::new();
188 for kid in &kids {
189 assert_eq!(kid.level(), level + 1);
190 assert_eq!(kid.face(), parent.face());
191 let back = match kid.parent(level) { Ok(c) => c, Err(e) => panic!("{}", e) };
192 assert_eq!(back, parent, "child did not parent back");
193 let (ci, cj) = kid.coords();
194 assert!(ci / 2 == pi && cj / 2 == pj, "child outside parent's index block");
195 assert!(seen.insert((ci, cj)), "duplicate child");
196 }
197 assert_eq!(seen.len(), 4);
198 }
199}
200
201// ---------------------------------------------------------------------------------------------
202// Invariant 4: neighbour symmetry and the cube-vertex census (exactly 24 cells with 7).
203// ---------------------------------------------------------------------------------------------
204
205#[test]
206fn inv04_neighbour_symmetry_and_census() {
207 // Level 0: each face touches its four side faces only.
208 for cell in all_cells(0) {
209 let nb = match cell.neighbours() { Ok(n) => n, Err(e) => panic!("{}", e) };
210 assert_eq!(nb.len(), 4, "level-0 face should have four neighbours");
211 }
212
213 for level in 1u8..=6 {
214 let cells = all_cells(level);
215 let mut nmap: HashMap<u64, HashSet<u64>> = HashMap::new();
216 for cell in &cells {
217 let nb = match cell.neighbours() { Ok(n) => n, Err(e) => panic!("{}", e) };
218 let set: HashSet<u64> = nb.iter().map(|c| c.bits()).collect();
219 assert!(!set.contains(&cell.bits()), "cell listed itself as neighbour");
220 nmap.insert(cell.bits(), set);
221 }
222 // Symmetry: b in N(a) iff a in N(b).
223 for (a, set) in &nmap {
224 for b in set {
225 let back = match nmap.get(b) { Some(s) => s, None => panic!("neighbour off census") };
226 assert!(back.contains(a), "neighbour relation not symmetric at level {}", level);
227 }
228 }
229 // Census: exactly 24 cells have 7 neighbours, the rest have 8.
230 let mut sevens = 0usize;
231 for set in nmap.values() {
232 match set.len() {
233 7 => sevens += 1,
234 8 => {},
235 other => panic!("cell has {} neighbours at level {}", other, level),
236 }
237 }
238 assert_eq!(sevens, 24, "vertex census wrong at level {}", level);
239 }
240}
241
242// ---------------------------------------------------------------------------------------------
243// Invariant 5: geometric adjacency -- neighbours equal the cells hit just outside the boundary.
244// ---------------------------------------------------------------------------------------------
245
246#[test]
247fn inv05_geometric_adjacency() {
248 let mut rng = Rng::new(0xFEED_FACE);
249 for _ in 0..600 {
250 let level = 2 + (rng.next_u64() % 5) as u8; // 2..=6
251 let cell = rng.cell_at(level);
252 let neigh: HashSet<u64> = match cell.neighbours() {
253 Ok(n) => n.iter().map(|c| c.bits()).collect(),
254 Err(e) => panic!("{}", e),
255 };
256 // Sample points just outside each edge midpoint and each corner of the boundary.
257 let corners = cell.corners();
258 let centre = norm(&cell.centre_vec());
259 let mut geo: HashSet<u64> = HashSet::new();
260 let boundary_samples = {
261 let mut pts: Vec<[f64; 3]> = Vec::new();
262 // Edge midpoints.
263 for e in 0..4 {
264 let a = &corners[e].vec;
265 let b = &corners[(e + 1) % 4].vec;
266 let mid = norm(&[a[0]+b[0], a[1]+b[1], a[2]+b[2]]);
267 pts.push(nudge_out(&mid, &centre));
268 }
269 // Corners.
270 for e in 0..4 {
271 pts.push(nudge_out(&corners[e].vec, &centre));
272 }
273 pts
274 };
275 for p in &boundary_samples {
276 let (lat, lon) = vec_latlon(p);
277 let hit = match Cell::at(lat, lon, level) { Ok(c) => c, Err(e) => panic!("{}", e) };
278 if hit.bits() != cell.bits() {
279 geo.insert(hit.bits());
280 }
281 }
282 assert_eq!(geo, neigh, "geometric adjacency disagreed with neighbours()");
283 }
284}
285
286fn nudge_out(p: &[f64; 3], centre: &[f64; 3]) -> [f64; 3] {
287 // Move from the boundary point a hair further from the centre.
288 let d = 1e-5;
289 let out = [p[0] - centre[0], p[1] - centre[1], p[2] - centre[2]];
290 let outn = norm(&out);
291 norm(&[p[0] + d*outn[0], p[1] + d*outn[1], p[2] + d*outn[2]])
292}
293
294fn vec_latlon(v: &[f64; 3]) -> (f64, f64) {
295 let u = norm(v);
296 (u[2].clamp(-1.0, 1.0).asin().to_degrees(), u[1].atan2(u[0]).to_degrees())
297}
298
299// ---------------------------------------------------------------------------------------------
300// Invariant 6: edge-neighbours share exactly 2 corners, corner-neighbours exactly 1.
301// ---------------------------------------------------------------------------------------------
302
303#[test]
304fn inv06_shared_corners() {
305 let mut rng = Rng::new(0xC0FF_EE00);
306 let mut tested = 0;
307 while tested < 200 {
308 let level = 3 + (rng.next_u64() % 4) as u8; // 3..=6
309 let n = 1u32 << level;
310 // Choose an interior cell (not touching a cube vertex) so it has 8 neighbours.
311 let i = 1 + (rng.next_u64() % (n as u64 - 2)) as u32;
312 let j = 1 + (rng.next_u64() % (n as u64 - 2)) as u32;
313 let face = (rng.next_u64() % 6) as u8;
314 let cell = match Cell::from_face_ij(face, level, i, j) { Ok(c) => c, Err(e) => panic!("{}", e) };
315 let nb = match cell.neighbours() { Ok(n) => n, Err(e) => panic!("{}", e) };
316 assert_eq!(nb.len(), 8, "interior cell should have 8 neighbours");
317 let mut twos = 0;
318 let mut ones = 0;
319 for other in &nb {
320 match shared_corner_count(&cell, other) {
321 2 => twos += 1,
322 1 => ones += 1,
323 k => panic!("neighbour shares {} corners", k),
324 }
325 }
326 assert_eq!(twos, 4, "expected 4 edge-neighbours");
327 assert_eq!(ones, 4, "expected 4 corner-neighbours");
328 tested += 1;
329 }
330}
331
332// ---------------------------------------------------------------------------------------------
333// Invariant 7: tiling and evenness -- areas sum to 4*pi and lie in [avg/sqrt2, avg*sqrt2].
334// ---------------------------------------------------------------------------------------------
335
336#[test]
337fn inv07_tiling_and_evenness() {
338 let root2 = 2.0f64.sqrt();
339 for level in 1u8..=5 {
340 let cells = all_cells(level);
341 let count = cells.len() as f64;
342 let mut sum = 0.0;
343 let mut min = f64::INFINITY;
344 let mut max = 0.0f64;
345 for cell in &cells {
346 let a = cell.area();
347 // The production area must match the independent oracle.
348 let ao = cell_area_oracle(cell);
349 assert!((a - ao).abs() < 1e-12, "area disagrees with oracle");
350 sum += a;
351 if a < min { min = a; }
352 if a > max { max = a; }
353 }
354 assert!((sum - 4.0*std::f64::consts::PI).abs() < 1e-9,
355 "areas did not tile the sphere at level {}: sum={}", level, sum);
356 let avg = 4.0*std::f64::consts::PI / count;
357 assert!(min >= avg/root2 - 1e-12 && max <= avg*root2 + 1e-12,
358 "evenness band violated at level {}: min/avg={} max/avg={}", level, min/avg, max/avg);
359 }
360}
361
362// ---------------------------------------------------------------------------------------------
363// Invariant 8: contains agrees with at, including at poles, antimeridian, edges and vertices.
364// ---------------------------------------------------------------------------------------------
365
366#[test]
367fn inv08_contains_boundaries() {
368 // Random interior points: contains, at and the independent naive oracle must all agree.
369 let mut rng = Rng::new(0xABCD_1234);
370 for _ in 0..3000 {
371 let level = (rng.next_u64() % 12) as u8;
372 let lat = rng.lat();
373 let lon = rng.lon();
374 let cell = match Cell::at(lat, lon, level) { Ok(c) => c, Err(e) => panic!("{}", e) };
375 assert!(match cell.contains(lat, lon) { Ok(b) => b, Err(e) => panic!("{}", e) },
376 "cell does not contain the point it was built from");
377 let p = latlon_to_vec(lat, lon);
378 // The naive geometric oracle must agree for a clearly-interior point.
379 assert!(naive_contains(&cell, &p), "naive oracle disagreed with at()");
380 }
381
382 // Special boundary points at several levels: resolution is total and deterministic.
383 let mut specials: Vec<(f64, f64)> = vec![
384 (90.0, 0.0), (-90.0, 0.0), // poles
385 (0.0, 180.0), (0.0, -180.0), // antimeridian
386 (0.0, 45.0), (0.0, 135.0), // face edges (|x|=|y|)
387 (45.0, 0.0), // face edge (|x|=|z|)
388 ];
389 let vlat = (1.0f64/3.0).sqrt().asin().to_degrees(); // 35.264 deg
390 for k in 0..4 {
391 let lon = 45.0 + 90.0 * k as f64;
392 specials.push((vlat, lon));
393 specials.push((-vlat, lon));
394 }
395 for (lat, lon) in specials {
396 for level in 0u8..=8 {
397 let c1 = match Cell::at(lat, lon, level) { Ok(c) => c, Err(e) => panic!("{}", e) };
398 let c2 = match Cell::at(lat, lon, level) { Ok(c) => c, Err(e) => panic!("{}", e) };
399 assert_eq!(c1, c2, "at() not deterministic at boundary ({}, {})", lat, lon);
400 assert!(c1.face() < 6 && c1.level() == level);
401 assert!(match c1.contains(lat, lon) { Ok(b) => b, Err(e) => panic!("{}", e) },
402 "boundary point not contained by its own cell ({}, {})", lat, lon);
403 }
404 }
405}
406
407// ---------------------------------------------------------------------------------------------
408// Invariant 9: fixtures -- face centres, string round-trip, malformed ids refused.
409// ---------------------------------------------------------------------------------------------
410
411#[test]
412fn inv09_fixtures() {
413 // Six face centres at level 0.
414 let expect: [(f64, f64); 6] = [
415 (0.0, 0.0), // +x
416 (0.0, 90.0), // +y
417 (90.0, 0.0), // +z (north pole)
418 (0.0, 180.0), // -x
419 (0.0, -90.0), // -y
420 (-90.0, 0.0), // -z (south pole)
421 ];
422 for (face, (elat, elon)) in expect.iter().enumerate() {
423 let c = match Cell::from_face_ij(face as u8, 0, 0, 0) { Ok(c) => c, Err(e) => panic!("{}", e) };
424 let (lat, lon) = c.centre();
425 assert!((lat - elat).abs() < 1e-9, "face {} centre lat {} != {}", face, lat, elat);
426 // Longitude is undefined at the poles, so only check away from them; +/-180 name
427 // the same meridian, so compare modulo a full turn.
428 if elat.abs() < 89.0 {
429 let d = ((lon - elon).rem_euclid(360.0) + 180.0).rem_euclid(360.0) - 180.0;
430 assert!(d.abs() < 1e-9, "face {} centre lon {} != {}", face, lon, elon);
431 }
432 }
433
434 // String round-trip over random cells at all levels.
435 let mut rng = Rng::new(0x5150_5150);
436 for level in 0u8..=26 {
437 for _ in 0..50 {
438 let c = rng.cell_at(level);
439 let s = c.to_string();
440 assert_eq!(s.len(), 16, "id string not 16 hex chars");
441 let back: Cell = match s.parse() { Ok(c) => c, Err(e) => panic!("parse: {}", e) };
442 assert_eq!(back, c, "string round-trip failed");
443 }
444 }
445
446 // Malformed ids are refused.
447 let good = match Cell::from_face_ij(2, 5, 7, 11) { Ok(c) => c, Err(e) => panic!("{}", e) };
448 let bits = good.bits();
449 // Unknown scheme (nibble != 1).
450 assert!(Cell::from_bits((bits & !(0xFu64 << 60)) | (0x2u64 << 60)).is_err(), "bad scheme accepted");
451 assert!(Cell::from_bits(bits & !(0xFu64 << 60)).is_err(), "zero scheme accepted");
452 // Face >= 6.
453 assert!(Cell::from_bits((bits & !(0x7u64 << 57)) | (6u64 << 57)).is_err(), "face 6 accepted");
454 assert!(Cell::from_bits((bits & !(0x7u64 << 57)) | (7u64 << 57)).is_err(), "face 7 accepted");
455 // Level > 26.
456 assert!(Cell::from_bits((bits & !(0x1Fu64 << 52)) | (27u64 << 52)).is_err(), "level 27 accepted");
457 // Non-zero unused low bit (level 5 uses top 10 of the 52-bit field; bit 0 must be zero).
458 assert!(Cell::from_bits(bits | 1u64).is_err(), "dirty low bit accepted");
459
460 // Bad string forms.
461 assert!("".parse::<Cell>().is_err());
462 assert!("xyz".parse::<Cell>().is_err());
463 assert!("deadbeef".parse::<Cell>().is_err(), "short string accepted");
464 assert!("00000000000000001".parse::<Cell>().is_err(), "long string accepted");
465 // A well-formed level-0 +x cell string.
466 let z = match Cell::from_face_ij(0, 0, 0, 0) { Ok(c) => c, Err(e) => panic!("{}", e) };
467 assert_eq!(z.to_string(), "1000000000000000");
468}
469
470// ---------------------------------------------------------------------------------------------
471// Invariant 10: k-ring BFS -- ring(0) is the cell, rings are disjoint, ring(1) == neighbours.
472// ---------------------------------------------------------------------------------------------
473
474#[test]
475fn inv10_ring() {
476 let mut rng = Rng::new(0x9999_1111);
477 for _ in 0..80 {
478 let level = 3 + (rng.next_u64() % 3) as u8; // 3..=5
479 let cell = rng.cell_at(level);
480 let r0 = match cell.ring(0) { Ok(r) => r, Err(e) => panic!("{}", e) };
481 assert_eq!(r0.len(), 1);
482 assert_eq!(r0[0], cell);
483
484 let r1: HashSet<u64> = match cell.ring(1) {
485 Ok(r) => r.iter().map(|c| c.bits()).collect(),
486 Err(e) => panic!("{}", e),
487 };
488 let nb: HashSet<u64> = match cell.neighbours() {
489 Ok(n) => n.iter().map(|c| c.bits()).collect(),
490 Err(e) => panic!("{}", e),
491 };
492 assert_eq!(r1, nb, "ring(1) must equal the neighbour set");
493
494 // Rings 0..=3 are pairwise disjoint.
495 let mut seen: HashSet<u64> = HashSet::new();
496 for k in 0..=3u32 {
497 let rk = match cell.ring(k) { Ok(r) => r, Err(e) => panic!("{}", e) };
498 for c in rk {
499 assert!(seen.insert(c.bits()), "rings overlapped at k={}", k);
500 }
501 }
502 }
503}