sim_lib_interference_compute/
tile.rs1use sim_lib_interference_core::{Point3M, SamplingPlane};
4use sim_lib_numbers_tensor::bounded_element_count;
5
6use crate::{LoweringError, PhaseBudget, PreflightCheck, TileProfile};
7
8#[derive(Clone, Copy, Debug, PartialEq)]
10pub struct PlaneTile {
11 row_start: usize,
12 column_start: usize,
13 rows: usize,
14 columns: usize,
15 center: Point3M,
16 max_offset_metres: f64,
17 segments_per_tensor: usize,
18 tensor_bytes: u64,
19}
20
21impl PlaneTile {
22 pub fn row_start(self) -> usize {
24 self.row_start
25 }
26
27 pub fn column_start(self) -> usize {
29 self.column_start
30 }
31
32 pub fn rows(self) -> usize {
34 self.rows
35 }
36
37 pub fn columns(self) -> usize {
39 self.columns
40 }
41
42 pub fn shape(self) -> [usize; 2] {
44 [self.rows, self.columns]
45 }
46
47 pub fn center(self) -> Point3M {
49 self.center
50 }
51
52 pub fn max_offset_metres(self) -> f64 {
54 self.max_offset_metres
55 }
56
57 pub fn segments_per_tensor(self) -> usize {
59 self.segments_per_tensor
60 }
61
62 pub fn tensor_bytes(self) -> u64 {
64 self.tensor_bytes
65 }
66}
67
68#[derive(Clone, Debug, PartialEq)]
70pub struct TilePlan {
71 tiles: Vec<PlaneTile>,
72 tile_rows: usize,
73 tile_columns: usize,
74 conservative_max_abs_phase_rad: f64,
75}
76
77impl TilePlan {
78 pub fn new(
80 plane: SamplingPlane,
81 wavenumber_rad_per_metre: f64,
82 budget: PhaseBudget,
83 profile: TileProfile,
84 ) -> Result<Self, LoweringError> {
85 budget.validate()?;
86 let element_limit = profile.effective_element_limit()?;
87 if !wavenumber_rad_per_metre.is_finite() || wavenumber_rad_per_metre <= 0.0 {
88 return Err(LoweringError::new(
89 PreflightCheck::Constant,
90 "real wavenumber must be finite and positive",
91 ));
92 }
93 if !(wavenumber_rad_per_metre as f32).is_finite() {
94 return Err(LoweringError::new(
95 PreflightCheck::Constant,
96 "real wavenumber is not representable as finite f32",
97 ));
98 }
99 admit_result_allocation(plane, profile)?;
100
101 let phase_radius = budget.max_abs_residual_phase_rad / wavenumber_rad_per_metre;
102 let (mut tile_rows, mut tile_columns) = phase_bounded_shape(plane, phase_radius);
103 tile_columns = tile_columns.min(element_limit).max(1);
104 tile_rows = tile_rows
105 .min(element_limit.checked_div(tile_columns).unwrap_or(0))
106 .max(1);
107
108 let row_tiles = plane.rows().div_ceil(tile_rows);
109 let column_tiles = plane.columns().div_ceil(tile_columns);
110 let tile_count = row_tiles.checked_mul(column_tiles).ok_or_else(|| {
111 LoweringError::new(PreflightCheck::Shape, "tile count overflowed usize")
112 })?;
113 let mut tiles = Vec::with_capacity(tile_count);
114 let mut conservative_max_abs_phase_rad = 0.0_f64;
115
116 for row_start in (0..plane.rows()).step_by(tile_rows) {
117 let rows = tile_rows.min(plane.rows() - row_start);
118 for column_start in (0..plane.columns()).step_by(tile_columns) {
119 let columns = tile_columns.min(plane.columns() - column_start);
120 let tile = build_tile(plane, profile, row_start, column_start, rows, columns)?;
121 let phase = wavenumber_rad_per_metre * tile.max_offset_metres;
122 if !phase.is_finite()
123 || phase > budget.max_abs_residual_phase_rad * (1.0 + 8.0 * f64::EPSILON)
124 {
125 return Err(LoweringError::new(
126 PreflightCheck::PhaseBudget,
127 format!(
128 "tile ({row_start},{column_start}) predicts residual phase {phase} rad above {}",
129 budget.max_abs_residual_phase_rad
130 ),
131 ));
132 }
133 conservative_max_abs_phase_rad = conservative_max_abs_phase_rad.max(phase);
134 tiles.push(tile);
135 }
136 }
137
138 Ok(Self {
139 tiles,
140 tile_rows,
141 tile_columns,
142 conservative_max_abs_phase_rad,
143 })
144 }
145
146 pub fn tiles(&self) -> &[PlaneTile] {
148 &self.tiles
149 }
150
151 pub fn tile_rows(&self) -> usize {
153 self.tile_rows
154 }
155
156 pub fn tile_columns(&self) -> usize {
158 self.tile_columns
159 }
160
161 pub fn conservative_max_abs_phase_rad(&self) -> f64 {
163 self.conservative_max_abs_phase_rad
164 }
165}
166
167fn admit_result_allocation(
168 plane: SamplingPlane,
169 profile: TileProfile,
170) -> Result<(), LoweringError> {
171 let cells = u64::try_from(plane.cell_count()).map_err(|_| {
172 LoweringError::new(
173 PreflightCheck::ResultAllocation,
174 "plane cell count does not fit u64",
175 )
176 })?;
177 let bytes = cells.checked_mul(8).ok_or_else(|| {
178 LoweringError::new(
179 PreflightCheck::ResultAllocation,
180 "two-component f32 result byte count overflowed",
181 )
182 })?;
183 if bytes > profile.max_result_bytes {
184 return Err(LoweringError::new(
185 PreflightCheck::ResultAllocation,
186 format!(
187 "two-component f32 result requires {bytes} bytes above {}",
188 profile.max_result_bytes
189 ),
190 ));
191 }
192 Ok(())
193}
194
195fn phase_bounded_shape(plane: SamplingPlane, radius: f64) -> (usize, usize) {
196 let rows = plane.rows();
197 let columns = plane.columns();
198 if rows == 1 {
199 return (
200 1,
201 axis_cells_within_radius(radius, plane.cell_size_u_m(), columns),
202 );
203 }
204 if columns == 1 {
205 return (
206 axis_cells_within_radius(radius, plane.cell_size_v_m(), rows),
207 1,
208 );
209 }
210 let axis_radius = radius / 2.0_f64.sqrt();
211 (
212 axis_cells_within_radius(axis_radius, plane.cell_size_v_m(), rows),
213 axis_cells_within_radius(axis_radius, plane.cell_size_u_m(), columns),
214 )
215}
216
217fn axis_cells_within_radius(radius: f64, spacing: f64, extent: usize) -> usize {
218 if radius >= (extent.saturating_sub(1) as f64) * spacing * 0.5 {
219 return extent;
220 }
221 let span_cells = (2.0 * radius / spacing).floor();
222 if span_cells >= usize::MAX as f64 {
223 extent
224 } else {
225 (span_cells as usize).saturating_add(1).clamp(1, extent)
226 }
227}
228
229fn build_tile(
230 plane: SamplingPlane,
231 profile: TileProfile,
232 row_start: usize,
233 column_start: usize,
234 rows: usize,
235 columns: usize,
236) -> Result<PlaneTile, LoweringError> {
237 let shape = [rows, columns];
238 let elements = bounded_element_count(&shape).map_err(|error| {
239 LoweringError::new(
240 PreflightCheck::Shape,
241 format!("tile shape {shape:?} is invalid: {error}"),
242 )
243 })?;
244 let bytes = u64::try_from(elements)
245 .ok()
246 .and_then(|count| count.checked_mul(4))
247 .ok_or_else(|| {
248 LoweringError::new(
249 PreflightCheck::ResultAllocation,
250 "tile Tensor byte count overflowed",
251 )
252 })?;
253 if elements > profile.max_elements_per_tile || bytes > profile.max_tensor_bytes {
254 return Err(LoweringError::new(
255 PreflightCheck::ResultAllocation,
256 format!("tile shape {shape:?} exceeds its element or byte limit"),
257 ));
258 }
259 let segments_u64 = bytes.div_ceil(profile.max_segment_bytes);
260 let segments = usize::try_from(segments_u64).map_err(|_| {
261 LoweringError::new(
262 PreflightCheck::ResultAllocation,
263 "tile segment count does not fit usize",
264 )
265 })?;
266 if segments > profile.max_segments_per_tensor {
267 return Err(LoweringError::new(
268 PreflightCheck::ResultAllocation,
269 format!(
270 "tile shape {shape:?} needs {segments} segments above {}",
271 profile.max_segments_per_tensor
272 ),
273 ));
274 }
275
276 let center = tile_center(plane, row_start, column_start, rows, columns)?;
277 let half_u = (columns.saturating_sub(1) as f64) * plane.cell_size_u_m() * 0.5;
278 let half_v = (rows.saturating_sub(1) as f64) * plane.cell_size_v_m() * 0.5;
279 let max_offset_metres = half_u.hypot(half_v);
280 Ok(PlaneTile {
281 row_start,
282 column_start,
283 rows,
284 columns,
285 center,
286 max_offset_metres,
287 segments_per_tensor: segments,
288 tensor_bytes: bytes,
289 })
290}
291
292fn tile_center(
293 plane: SamplingPlane,
294 row_start: usize,
295 column_start: usize,
296 rows: usize,
297 columns: usize,
298) -> Result<Point3M, LoweringError> {
299 let offset_u = (column_start as f64 + columns as f64 * 0.5) * plane.cell_size_u_m();
300 let offset_v = (row_start as f64 + rows as f64 * 0.5) * plane.cell_size_v_m();
301 let [origin_x, origin_y, origin_z] = plane.origin().coordinates_metres();
302 let [ux, uy, uz] = plane.u_axis().components();
303 let [vx, vy, vz] = plane.v_axis().components();
304 Point3M::from_metres(
305 origin_x + offset_u * ux + offset_v * vx,
306 origin_y + offset_u * uy + offset_v * vy,
307 origin_z + offset_u * uz + offset_v * vz,
308 )
309 .map_err(|error| {
310 LoweringError::new(
311 PreflightCheck::Shape,
312 format!("tile center is not finite: {error}"),
313 )
314 })
315}