Skip to main content

polydat_core/iteration/comprehension/strategies/
shells.rs

1// Copyright 2024-2026 Jonathan Shook
2// SPDX-License-Identifier: Apache-2.0
3
4//! `Shells` strategy — spec §3.6.
5//!
6//! Emits multi-indices stratified by concentric shells around
7//! the lattice center, outermost first. A "shell" is the set
8//! of multi-indices at Chebyshev distance `d` from the center
9//! (max-norm). Discrete `Lattice` is the native shape;
10//! continuous rejected per spec §3.6 ("ill-defined without
11//! discretization parameter").
12//!
13//! Emission within a shell uses Lex order as tiebreak so the
14//! walk is fully deterministic.
15//!
16//! ## References
17//!
18//! - The shell metric is the Chebyshev / L∞ (max-norm) distance,
19//!   named for P. L. Chebyshev; see e.g. M. M. Deza & E. Deza,
20//!   *Encyclopedia of Distances*, 4th ed., Springer (2016), §1.1.
21//!   A "shell" is the set of points at a fixed L∞ distance from the
22//!   centre — the square (hyper-cube) ring at radius `d`. This
23//!   differs from [`super::extrema`], whose strata are by *interior
24//!   count* (k-faces), not a distance. Outermost-first ordering and
25//!   the shell metric are cross-checked in
26//!   `tests::shells_are_chebyshev_strata_outermost_first`.
27
28use super::{
29    EvaluatedInput, MultiIndex, Strategy, Tuple, index_fn_size, index_fn_supports_lookup,
30    lex::lex_multi_indices, multi_index_to_flat,
31};
32use crate::iteration::comprehension::metadata::IndexFn;
33use crate::iteration::comprehension::strategy::StrategyName;
34
35/// Concentric L-infinity shells from a chosen origin.
36pub struct Shells;
37
38impl Strategy for Shells {
39    fn name(&self) -> StrategyName {
40        StrategyName::Shells
41    }
42
43    fn accepts_input(&self, idx: Option<&IndexFn>) -> bool {
44        match idx {
45            None => false,
46            Some(i) => !i.has_continuous_axis(),
47        }
48    }
49
50    fn has_closed_form_for(&self, idx: &IndexFn) -> bool {
51        matches!(idx, IndexFn::Lattice { .. })
52    }
53
54    fn apply(&self, input: &EvaluatedInput, truncation: Option<u64>) -> Vec<Tuple> {
55        if index_fn_supports_lookup(&input.index_fn) {
56            let mis = shells_multi_indices(&input.index_fn, truncation);
57            mis.into_iter()
58                .filter_map(|mi| multi_index_to_flat(&input.index_fn, &mi))
59                .filter_map(|flat| input.tuples.get(flat).cloned())
60                .collect()
61        } else {
62            naive_lex_prefix(&input.tuples, truncation)
63        }
64    }
65}
66
67fn naive_lex_prefix(input: &[Tuple], truncation: Option<u64>) -> Vec<Tuple> {
68    match truncation {
69        Some(n) => input.iter().take(n as usize).cloned().collect(),
70        None => input.to_vec(),
71    }
72}
73
74pub(crate) fn shells_multi_indices(idx: &IndexFn, truncation: Option<u64>) -> Vec<MultiIndex> {
75    let total = index_fn_size(idx);
76    let n = match truncation {
77        Some(t) => t.min(total),
78        None => total,
79    };
80
81    let axis_sizes = match idx {
82        IndexFn::Lattice { axis_sizes } => axis_sizes.clone(),
83        _ => return lex_multi_indices(idx, truncation),
84    };
85
86    if axis_sizes.is_empty() {
87        return Vec::new();
88    }
89
90    let centers: Vec<f64> = axis_sizes.iter().map(|s| (*s as f64 - 1.0) / 2.0).collect();
91
92    // Bucket by *-2 + rounded-int Chebyshev to avoid fp issues
93    // with half-integer centers.
94    let mut buckets: std::collections::BTreeMap<i64, Vec<MultiIndex>> =
95        std::collections::BTreeMap::new();
96
97    enumerate_all(
98        &axis_sizes,
99        &mut Vec::with_capacity(axis_sizes.len()),
100        &mut |mi| {
101            let r = chebyshev_distance(mi, &centers);
102            let bucket_key = (r * 2.0).round() as i64;
103            buckets.entry(bucket_key).or_default().push(mi.clone());
104        },
105    );
106
107    let mut out: Vec<MultiIndex> = Vec::with_capacity(n as usize);
108    for (_key, mut shell) in buckets.into_iter().rev() {
109        shell.sort();
110        for mi in shell {
111            if out.len() as u64 >= n {
112                return out;
113            }
114            out.push(mi);
115        }
116    }
117    out
118}
119
120fn chebyshev_distance(mi: &[u64], center: &[f64]) -> f64 {
121    mi.iter()
122        .zip(center.iter())
123        .map(|(c, ctr)| ((*c as f64) - ctr).abs())
124        .fold(0.0f64, f64::max)
125}
126
127fn enumerate_all(
128    axis_sizes: &[u64],
129    current: &mut Vec<u64>,
130    callback: &mut dyn FnMut(&MultiIndex),
131) {
132    if current.len() == axis_sizes.len() {
133        callback(current);
134        return;
135    }
136    let size = axis_sizes[current.len()];
137    for v in 0..size {
138        current.push(v);
139        enumerate_all(axis_sizes, current, callback);
140        current.pop();
141    }
142}
143
144#[cfg(test)]
145mod tests {
146    use super::*;
147
148    #[test]
149    fn shells_3x3() {
150        let idx = IndexFn::Lattice {
151            axis_sizes: vec![3, 3],
152        };
153        let out = shells_multi_indices(&idx, None);
154        assert_eq!(out.len(), 9);
155        // Center (1,1) is innermost.
156        assert_eq!(out[8], vec![1, 1]);
157    }
158
159    #[test]
160    fn shells_outermost_first_for_5x5() {
161        let idx = IndexFn::Lattice {
162            axis_sizes: vec![5, 5],
163        };
164        let out = shells_multi_indices(&idx, Some(4));
165        for mi in &out {
166            let chebyshev = chebyshev_distance(mi, &[2.0, 2.0]);
167            assert!(
168                (chebyshev - 2.0).abs() < 1e-9,
169                "expected distance 2.0, got {chebyshev:?}"
170            );
171        }
172    }
173
174    #[test]
175    fn shells_are_chebyshev_strata_outermost_first() {
176        // 5×5 about centre (2,2): three L∞ shells — radius 2 (the 16
177        // boundary points), radius 1 (the 8-point inner ring), radius
178        // 0 (the centre). The full walk visits them outermost-first,
179        // so Chebyshev distance is monotonically non-increasing.
180        let idx = IndexFn::Lattice {
181            axis_sizes: vec![5, 5],
182        };
183        let out = shells_multi_indices(&idx, None);
184        assert_eq!(out.len(), 25);
185        let dists: Vec<f64> = out
186            .iter()
187            .map(|mi| chebyshev_distance(mi, &[2.0, 2.0]))
188            .collect();
189        assert!(
190            dists.windows(2).all(|w| w[0] >= w[1] - 1e-9),
191            "shell distances not outermost-first: {dists:?}"
192        );
193        // Exactly three distinct radii {2,1,0} with the documented
194        // populations 16 / 8 / 1.
195        let r2 = dists.iter().filter(|d| (**d - 2.0).abs() < 1e-9).count();
196        let r1 = dists.iter().filter(|d| (**d - 1.0).abs() < 1e-9).count();
197        let r0 = dists.iter().filter(|d| **d < 1e-9).count();
198        assert_eq!((r2, r1, r0), (16, 8, 1));
199    }
200
201    #[test]
202    fn rejects_continuous() {
203        use crate::iteration::comprehension::cardinality::{Interval, ProductMeasure};
204        let cont = IndexFn::Continuous {
205            intervals: vec![Interval::closed(0.0, 1.0)],
206            measure: ProductMeasure::Uniform,
207        };
208        assert!(!Shells.accepts_input(Some(&cont)));
209    }
210}