1use std::marker::PhantomData;
4
5use rand::Rng;
6
7use super::*;
8
9pub struct Lattice<A, P = ()>
20where
21 A: Axis,
22{
23 pairs: Box<[(f64, P)]>,
24 axis: PhantomData<A>,
25}
26
27impl<A, P> std::fmt::Debug for Lattice<A, P>
28where
29 A: Axis,
30 P: std::fmt::Debug,
31{
32 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
33 f.debug_struct("Lattice").field("pairs", &self.pairs).finish()
34 }
35}
36
37impl<A, P> Clone for Lattice<A, P>
38where
39 A: Axis,
40 P: Clone,
41{
42 fn clone(&self) -> Self {
43 Self {
44 pairs: self.pairs.clone(),
45 axis: PhantomData,
46 }
47 }
48}
49
50impl<A, P> FromIterator<(f64, P)> for Lattice<A, P>
51where
52 A: Axis,
53{
54 fn from_iter<I>(iter: I) -> Self
58 where
59 I: IntoIterator<Item = (f64, P)>,
60 {
61 let mut pairs = iter.into_iter().collect::<Vec<_>>();
62 assert!(!pairs.is_empty(), "Lattice must be non-empty");
63 assert!(pairs.iter().all(|(s, _)| s.is_finite()), "Lattice scalars must be finite");
64 pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).expect("finite scalar keys"));
65 assert!(
66 pairs.windows(2).all(|w| w[0].0 < w[1].0),
67 "Lattice scalars must be strictly ascending (no duplicates)",
68 );
69 Self {
70 pairs: pairs.into_boxed_slice(),
71 axis: PhantomData,
72 }
73 }
74}
75
76impl<A> FromIterator<f64> for Lattice<A, ()>
77where
78 A: Axis,
79{
80 fn from_iter<I>(iter: I) -> Self
82 where
83 I: IntoIterator<Item = f64>,
84 {
85 iter.into_iter().map(|s| (s, ())).collect()
86 }
87}
88
89impl<A, P> Lattice<A, P>
90where
91 A: Axis,
92{
93 pub fn len(&self) -> usize {
95 self.pairs.len()
96 }
97
98 pub fn is_empty(&self) -> bool {
101 self.pairs.is_empty()
102 }
103
104 pub fn scalars(&self) -> impl ExactSizeIterator<Item = f64> + '_ {
106 self.pairs.iter().map(|(s, _)| *s)
107 }
108
109 pub fn payload(&self, anchor: Anchor) -> &P {
111 &self.pairs[anchor.idx()].1
112 }
113
114 pub fn bracket(&self, observed: Scalar<A>) -> Bracket {
119 let x = observed.value();
120 let i = self.pairs.len() - 1;
121 if x <= self.pairs[0].0 {
122 Bracket::new(Anchor::new(0), Anchor::new(0))
123 } else if x >= self.pairs[i].0 {
124 Bracket::new(Anchor::new(i), Anchor::new(i))
125 } else {
126 Bracket::from(self.pairs.partition_point(|(s, _)| *s < x))
127 }
128 }
129}
130
131impl<A, P> Lattice<A, P>
138where
139 A: Axis,
140{
141 pub fn snap(&self, observed: Scalar<A>) -> Anchor {
142 let x = observed.value();
143 let i = self
144 .scalars()
145 .enumerate()
146 .min_by(|(_, a), (_, b)| (a - x).abs().partial_cmp(&(b - x).abs()).expect("finite anchors"))
147 .map(|(i, _)| i)
148 .expect("non-empty lattice");
149 Anchor::new(i)
150 }
151
152 pub fn pharmonic(&self, bracket: Bracket, observed: Scalar<A>) -> f64 {
160 if bracket.is_clamped() {
161 unreachable!("pharmonic requires distinct bracketing anchors")
162 } else {
163 let a = self.pairs[bracket.lo().idx()].0;
164 let b = self.pairs[bracket.hi().idx()].0;
165 let x = observed.value();
166 ((b - x) * (1.0 + a)) / ((b - a) * (1.0 + x))
167 }
168 }
169
170 pub fn harmonic<R>(&self, observed: Scalar<A>, rng: &mut R) -> Anchor
171 where
172 R: Rng + ?Sized,
173 {
174 let bracket = self.bracket(observed);
175 if bracket.is_clamped() || rng.random::<f64>() < self.pharmonic(bracket, observed) {
176 bracket.lo()
177 } else {
178 bracket.hi()
179 }
180 }
181
182 pub fn phargmax(&self, observed: Scalar<A>) -> Anchor {
183 let bracket = self.bracket(observed);
184 if bracket.is_clamped() || self.pharmonic(bracket, observed) >= 0.5 {
185 bracket.lo()
186 } else {
187 bracket.hi()
188 }
189 }
190}
191
192#[cfg(test)]
193mod tests {
194 use super::*;
195
196 struct T;
197 impl Axis for T {}
198
199 fn obs(x: f64) -> Scalar<T> {
200 Scalar::new(x)
201 }
202
203 fn lat(xs: impl IntoIterator<Item = f64>) -> Lattice<T> {
204 xs.into_iter().collect()
205 }
206
207 #[test]
208 fn from_scalars_constructs_unit_payload() {
209 let l = lat([0.5, 1.0, 2.0]);
210 assert_eq!(l.scalars().collect::<Vec<_>>(), vec![0.5, 1.0, 2.0]);
211 assert_eq!(l.len(), 3);
212 }
213
214 #[test]
215 fn from_iter_sorts_unsorted_pairs() {
216 let l = [(2.0, "two"), (0.5, "half"), (1.0, "one")]
217 .into_iter()
218 .collect::<Lattice<T, _>>();
219 assert_eq!(l.scalars().collect::<Vec<_>>(), vec![0.5, 1.0, 2.0]);
220 assert_eq!(*l.payload(Anchor::new(0)), "half");
221 assert_eq!(*l.payload(Anchor::new(2)), "two");
222 }
223
224 #[test]
225 #[should_panic(expected = "non-empty")]
226 fn rejects_empty() {
227 let _: Lattice<T> = std::iter::empty::<f64>().collect();
228 }
229
230 #[test]
231 #[should_panic(expected = "finite")]
232 fn rejects_non_finite() {
233 let _: Lattice<T, _> = [(0.0, ()), (f64::NAN, ()), (1.0, ())].into_iter().collect();
234 }
235
236 #[test]
237 #[should_panic(expected = "ascending")]
238 fn rejects_duplicate_keys() {
239 let _: Lattice<T, _> = [(1.0, 'a'), (1.0, 'b')].into_iter().collect();
240 }
241
242 #[test]
243 fn bracket_clamps_below_extreme() {
244 let l = lat([0.5, 1.0, 2.0]);
245 let b = l.bracket(obs(0.1));
246 assert!(b.is_clamped());
247 assert_eq!(b.lo(), Anchor::new(0));
248 let b = l.bracket(obs(0.5));
249 assert!(b.is_clamped());
250 assert_eq!(b.lo(), Anchor::new(0));
251 }
252
253 #[test]
254 fn bracket_clamps_above_extreme() {
255 let l = lat([0.5, 1.0, 2.0]);
256 let b = l.bracket(obs(2.0));
257 assert!(b.is_clamped());
258 assert_eq!(b.lo(), Anchor::new(2));
259 let b = l.bracket(obs(3.0));
260 assert!(b.is_clamped());
261 assert_eq!(b.lo(), Anchor::new(2));
262 }
263
264 #[test]
265 fn bracket_inside_returns_distinct_anchors() {
266 let l = lat([0.5, 1.0, 2.0]);
267 let b = l.bracket(obs(0.75));
268 assert!(!b.is_clamped());
269 assert_eq!(b.lo(), Anchor::new(0));
270 assert_eq!(b.hi(), Anchor::new(1));
271 }
272
273 #[test]
274 fn pharmonic_formula_exact() {
275 let l = lat([0.5, 1.0]);
276 let bracket = l.bracket(obs(0.75));
277 let expected = (1.0 - 0.75) * (1.0 + 0.5) / ((1.0 - 0.5) * (1.0 + 0.75));
278 let p = l.pharmonic(bracket, obs(0.75));
279 assert!((p - expected).abs() < 1e-12);
280 }
281}