sim_lib_interference_solve/
reference.rs

1//! Configuration, traversal, and evidence for the CPU reference solver.
2
3use sim_lib_interference_core::{
4    Emitter, InterferenceError, InterferenceProblem, Point3M, RequestPreflight, SamplingPlane,
5    SamplingPolicy, SamplingThresholds, WorkBudget, contribution_at,
6};
7
8use crate::{
9    HostPhasorField, ReferenceSolveError,
10    complex::{CompensatedSum, Complex64},
11};
12
13/// Immutable evidence retained for one complete reference solve.
14#[derive(Clone, Debug, PartialEq)]
15pub struct SolveEvidence {
16    preflight: RequestPreflight,
17    completed_cells: u64,
18    completed_emitter_evaluations: u64,
19}
20
21impl SolveEvidence {
22    /// Returns the sampling policy, certificate, and work admission record.
23    pub fn preflight(&self) -> RequestPreflight {
24        self.preflight
25    }
26
27    /// Returns the number of row-major output cells completed.
28    pub fn completed_cells(&self) -> u64 {
29        self.completed_cells
30    }
31
32    /// Returns the number of canonical source evaluations completed.
33    pub fn completed_emitter_evaluations(&self) -> u64 {
34        self.completed_emitter_evaluations
35    }
36}
37
38/// Configuration for deterministic CPU `f64` reference solving.
39#[derive(Clone, Copy, Debug, PartialEq)]
40pub struct ReferencePhasorSolver {
41    sampling_policy: SamplingPolicy,
42    sampling_thresholds: SamplingThresholds,
43    work_budget: WorkBudget,
44}
45
46impl ReferencePhasorSolver {
47    /// Constructs a solver from explicit admission policy.
48    pub fn new(
49        sampling_policy: SamplingPolicy,
50        sampling_thresholds: SamplingThresholds,
51        work_budget: WorkBudget,
52    ) -> Self {
53        Self {
54            sampling_policy,
55            sampling_thresholds,
56            work_budget,
57        }
58    }
59
60    /// Returns the configured sampling policy.
61    pub fn sampling_policy(self) -> SamplingPolicy {
62        self.sampling_policy
63    }
64
65    /// Returns the configured sampling thresholds.
66    pub fn sampling_thresholds(self) -> SamplingThresholds {
67        self.sampling_thresholds
68    }
69
70    /// Returns the configured work budget.
71    pub fn work_budget(self) -> WorkBudget {
72        self.work_budget
73    }
74
75    /// Solves one complete field or returns no field.
76    ///
77    /// Work and sampling admission run first, followed by an allocation-free
78    /// row-major geometry pass. Component storage is reserved only after those
79    /// checks pass.
80    pub fn solve(
81        self,
82        problem: &InterferenceProblem,
83        plane: &SamplingPlane,
84    ) -> Result<(HostPhasorField, SolveEvidence), ReferenceSolveError> {
85        let preflight = RequestPreflight::admit(
86            problem,
87            plane,
88            self.sampling_policy,
89            self.sampling_thresholds,
90            self.work_budget,
91        )
92        .map_err(ReferenceSolveError::from)?;
93        preflight_geometry(problem, plane)?;
94        let field = evaluate_field(problem, plane)?;
95        let evidence = SolveEvidence {
96            preflight,
97            completed_cells: preflight.work_estimate.cells,
98            completed_emitter_evaluations: preflight.work_estimate.emitter_evaluations,
99        };
100        Ok((field, evidence))
101    }
102}
103
104impl Default for ReferencePhasorSolver {
105    fn default() -> Self {
106        Self::new(
107            SamplingPolicy::Strict,
108            SamplingThresholds::default(),
109            WorkBudget::default(),
110        )
111    }
112}
113
114fn preflight_geometry(
115    problem: &InterferenceProblem,
116    plane: &SamplingPlane,
117) -> Result<(), ReferenceSolveError> {
118    for row in 0..plane.rows() {
119        for column in 0..plane.columns() {
120            let at =
121                plane
122                    .point_at(row, column)
123                    .map_err(|cause| ReferenceSolveError::CellGeometry {
124                        row,
125                        column,
126                        cause: Box::new(cause),
127                    })?;
128            for source in &problem.sources {
129                validate_source_geometry(problem, source, at).map_err(|cause| {
130                    ReferenceSolveError::SourceAtCell {
131                        source_id: source.id().to_owned(),
132                        row,
133                        column,
134                        cause: Box::new(cause),
135                    }
136                })?;
137            }
138        }
139    }
140    Ok(())
141}
142
143fn validate_source_geometry(
144    problem: &InterferenceProblem,
145    source: &Emitter,
146    at: Point3M,
147) -> Result<(), InterferenceError> {
148    match source {
149        Emitter::Point { id, position, .. } => {
150            let distance = at.distance_to(*position);
151            if !distance.is_finite() {
152                return Err(InterferenceError::NonFinitePropagation {
153                    source_id: id.clone(),
154                    name: "point-distance-metres",
155                    value: distance,
156                });
157            }
158            if distance <= problem.singularity_radius.get() {
159                return Err(InterferenceError::SingularPointSample {
160                    source_id: id.clone(),
161                    distance_metres: distance,
162                    singularity_radius_metres: problem.singularity_radius.get(),
163                });
164            }
165        }
166        Emitter::ForwardPlane {
167            id,
168            through,
169            direction,
170            ..
171        } => {
172            let signed_distance = direction.signed_distance_metres(*through, at);
173            if !signed_distance.is_finite() {
174                return Err(InterferenceError::NonFinitePropagation {
175                    source_id: id.clone(),
176                    name: "plane-signed-distance-metres",
177                    value: signed_distance,
178                });
179            }
180            if signed_distance < 0.0 {
181                return Err(InterferenceError::BehindForwardPlane {
182                    source_id: id.clone(),
183                    signed_distance_metres: signed_distance,
184                });
185            }
186        }
187    }
188    Ok(())
189}
190
191fn evaluate_field(
192    problem: &InterferenceProblem,
193    plane: &SamplingPlane,
194) -> Result<HostPhasorField, ReferenceSolveError> {
195    let mut field = HostPhasorField::try_zeroed(plane.rows(), plane.columns())?;
196    for row in 0..plane.rows() {
197        for column in 0..plane.columns() {
198            let at =
199                plane
200                    .point_at(row, column)
201                    .map_err(|cause| ReferenceSolveError::CellGeometry {
202                        row,
203                        column,
204                        cause: Box::new(cause),
205                    })?;
206            let value = accumulate_cell(problem, at, row, column)?;
207            let index = row * plane.columns() + column;
208            field.set_index(index, value);
209        }
210    }
211    Ok(field)
212}
213
214fn accumulate_cell(
215    problem: &InterferenceProblem,
216    at: Point3M,
217    row: usize,
218    column: usize,
219) -> Result<Complex64, ReferenceSolveError> {
220    let mut real_sum = CompensatedSum::default();
221    let mut imaginary_sum = CompensatedSum::default();
222
223    for source in &problem.sources {
224        let source_id = source.id();
225        let (real, imaginary) = contribution_at(problem, source, at).map_err(|cause| {
226            ReferenceSolveError::SourceAtCell {
227                source_id: source_id.to_owned(),
228                row,
229                column,
230                cause: Box::new(cause),
231            }
232        })?;
233        real_sum.add(real);
234        require_finite_accumulation(source_id, row, column, "real", real_sum)?;
235        imaginary_sum.add(imaginary);
236        require_finite_accumulation(source_id, row, column, "imaginary", imaginary_sum)?;
237    }
238
239    Ok(Complex64::new(real_sum.total(), imaginary_sum.total()))
240}
241
242fn require_finite_accumulation(
243    source_id: &str,
244    row: usize,
245    column: usize,
246    component: &'static str,
247    accumulation: CompensatedSum,
248) -> Result<(), ReferenceSolveError> {
249    let value = accumulation.total();
250    if accumulation.is_finite() {
251        Ok(())
252    } else {
253        Err(ReferenceSolveError::NonFiniteAccumulation {
254            source_id: source_id.to_owned(),
255            row,
256            column,
257            component,
258            value,
259        })
260    }
261}
262
263#[cfg(test)]
264mod tests {
265    use sim_lib_interference_core::{
266        Emitter, FieldAmplitude, Hertz, InterferenceProblem, MetresPerSecond, NepersPerMetre,
267        Point3M, PositiveMetres, Radians, SamplingPlane, SamplingPolicy, SamplingThresholds,
268        ScalarMedium, SourceSet, UnitVector3, WorkBudget,
269    };
270
271    use super::{ReferencePhasorSolver, evaluate_field};
272
273    fn point(x: f64, y: f64, z: f64) -> Point3M {
274        Point3M::from_metres(x, y, z).unwrap()
275    }
276
277    fn plane(rows: usize, columns: usize) -> SamplingPlane {
278        SamplingPlane::new(
279            point(0.0, 0.0, 0.0),
280            UnitVector3::new(1.0, 0.0, 0.0).unwrap(),
281            UnitVector3::new(0.0, 1.0, 0.0).unwrap(),
282            PositiveMetres::new(1.0).unwrap(),
283            PositiveMetres::new(1.0).unwrap(),
284            rows,
285            columns,
286        )
287        .unwrap()
288    }
289
290    fn problem(sources: Vec<Emitter>) -> InterferenceProblem {
291        InterferenceProblem::new(
292            Hertz::new(1.0).unwrap(),
293            ScalarMedium::new(
294                MetresPerSecond::new(100.0).unwrap(),
295                NepersPerMetre::new(0.0).unwrap(),
296            ),
297            SourceSet::new(sources).unwrap(),
298            PositiveMetres::new(0.001).unwrap(),
299        )
300    }
301
302    fn plane_source(id: &str, amplitude: f64) -> Emitter {
303        Emitter::ForwardPlane {
304            id: id.to_owned(),
305            through: point(0.0, 0.0, 0.0),
306            direction: UnitVector3::new(0.0, 0.0, 1.0).unwrap(),
307            amplitude: FieldAmplitude::new(amplitude).unwrap(),
308            phase: Radians::new(0.0).unwrap(),
309        }
310    }
311
312    #[test]
313    fn solver_configuration_is_explicit_and_immutable() {
314        let thresholds = SamplingThresholds::new(12.0, 6.0, 0.025, 0.075).unwrap();
315        let budget = WorkBudget {
316            max_cells: 100,
317            max_emitter_evaluations: 200,
318            max_host_bytes: 3_200,
319            max_result_bytes: 1_600,
320            max_certificate_stencil_work: 700,
321        };
322        let solver = ReferencePhasorSolver::new(SamplingPolicy::Annotate, thresholds, budget);
323
324        assert_eq!(solver.sampling_policy(), SamplingPolicy::Annotate);
325        assert_eq!(solver.sampling_thresholds(), thresholds);
326        assert_eq!(solver.work_budget(), budget);
327    }
328
329    #[test]
330    fn row_major_cells_use_canonical_sources_and_neumaier_components() {
331        let problem = problem(vec![
332            plane_source("c-small", 1.0),
333            plane_source("a-large", 1.0e16),
334            plane_source("b-small", 1.0),
335        ]);
336        assert_eq!(
337            problem.sources.iter().map(Emitter::id).collect::<Vec<_>>(),
338            ["a-large", "b-small", "c-small"]
339        );
340
341        let field = evaluate_field(&problem, &plane(2, 3)).unwrap();
342        assert_eq!(field.real(), &[1.0e16 + 2.0; 6]);
343        assert_eq!(field.imaginary(), &[0.0; 6]);
344        assert_eq!(field.cell(0, 2), Some((field.real()[2], 0.0)));
345        assert_eq!(field.cell(1, 0), Some((field.real()[3], 0.0)));
346    }
347
348    #[test]
349    fn solve_retains_complete_sampling_and_work_preflight() {
350        let problem = problem(vec![plane_source("plane", 2.0)]);
351        let (field, evidence) = ReferencePhasorSolver::default()
352            .solve(&problem, &plane(2, 3))
353            .unwrap();
354
355        assert_eq!(field.len(), 6);
356        assert_eq!(evidence.completed_cells(), 6);
357        assert_eq!(evidence.completed_emitter_evaluations(), 6);
358        assert_eq!(evidence.preflight().work_estimate.cells, 6);
359        assert_eq!(evidence.preflight().work_estimate.host_bytes, 96);
360        assert_eq!(evidence.preflight().sampling_policy, SamplingPolicy::Strict);
361    }
362
363    #[test]
364    fn request_and_cell_geometry_fail_before_field_construction() {
365        let problem = problem(vec![plane_source("plane", 1.0)]);
366        let no_cells = WorkBudget {
367            max_cells: 0,
368            ..WorkBudget::default()
369        };
370        let budget_error = ReferencePhasorSolver::new(
371            SamplingPolicy::Annotate,
372            SamplingThresholds::default(),
373            no_cells,
374        )
375        .solve(&problem, &plane(1, 1))
376        .unwrap_err();
377        assert!(matches!(
378            budget_error,
379            crate::ReferenceSolveError::Request { cause }
380                if matches!(*cause, sim_lib_interference_core::InterferenceError::WorkBudgetExceeded { .. })
381        ));
382
383        let extreme_plane = SamplingPlane::new(
384            point(f64::MAX, 0.0, 0.0),
385            UnitVector3::new(1.0, 0.0, 0.0).unwrap(),
386            UnitVector3::new(0.0, 1.0, 0.0).unwrap(),
387            PositiveMetres::new(f64::MAX).unwrap(),
388            PositiveMetres::new(1.0).unwrap(),
389            1,
390            1,
391        )
392        .unwrap();
393        let geometry_error = ReferencePhasorSolver::new(
394            SamplingPolicy::Annotate,
395            SamplingThresholds::default(),
396            WorkBudget::default(),
397        )
398        .solve(&problem, &extreme_plane)
399        .unwrap_err();
400        assert!(matches!(
401            geometry_error,
402            crate::ReferenceSolveError::CellGeometry {
403                row: 0,
404                column: 0,
405                ..
406            }
407        ));
408    }
409
410    #[test]
411    fn singular_and_behind_samples_name_source_and_cell() {
412        let singular_problem = InterferenceProblem::new(
413            Hertz::new(1.0).unwrap(),
414            ScalarMedium::new(
415                MetresPerSecond::new(100.0).unwrap(),
416                NepersPerMetre::new(0.0).unwrap(),
417            ),
418            SourceSet::new(vec![Emitter::Point {
419                id: "near-point".to_owned(),
420                position: point(0.5, 0.5, 0.01),
421                amplitude_at_reference: FieldAmplitude::new(1.0).unwrap(),
422                phase: Radians::new(0.0).unwrap(),
423            }])
424            .unwrap(),
425            PositiveMetres::new(0.02).unwrap(),
426        );
427        let solver = ReferencePhasorSolver::new(
428            SamplingPolicy::Annotate,
429            SamplingThresholds::default(),
430            WorkBudget::default(),
431        );
432        let singular_error = solver.solve(&singular_problem, &plane(1, 1)).unwrap_err();
433        assert!(matches!(
434            singular_error,
435            crate::ReferenceSolveError::SourceAtCell {
436                source_id,
437                row: 0,
438                column: 0,
439                cause,
440            } if source_id == "near-point"
441                && matches!(*cause, sim_lib_interference_core::InterferenceError::SingularPointSample { .. })
442        ));
443
444        let behind_problem = problem(vec![Emitter::ForwardPlane {
445            id: "forward-only".to_owned(),
446            through: point(0.0, 0.0, 1.0),
447            direction: UnitVector3::new(0.0, 0.0, 1.0).unwrap(),
448            amplitude: FieldAmplitude::new(1.0).unwrap(),
449            phase: Radians::new(0.0).unwrap(),
450        }]);
451        let behind_error = solver.solve(&behind_problem, &plane(2, 2)).unwrap_err();
452        assert!(matches!(
453            behind_error,
454            crate::ReferenceSolveError::SourceAtCell {
455                source_id,
456                row: 0,
457                column: 0,
458                cause,
459            } if source_id == "forward-only"
460                && matches!(*cause, sim_lib_interference_core::InterferenceError::BehindForwardPlane { .. })
461        ));
462    }
463}
sim_lib_interference_solve/reference.rs

sim_lib_interference_solve/
reference.rs
