sim_lib_interference_compute/
diff.rs

1//! Fixed-tolerance differential reporting against the deterministic f64 oracle.
2
3use std::{f64::consts::PI, fmt};
4
5use sim_lib_interference_solve::HostPhasorField;
6
7use crate::DenseF32Field;
8
9/// A quantity checked by the portable f32 differential contract.
10#[derive(Clone, Copy, Debug, Eq, PartialEq)]
11pub enum ConformanceMetric {
12    /// Real Cartesian component.
13    Real,
14    /// Imaginary Cartesian component.
15    Imaginary,
16    /// Complex amplitude.
17    Amplitude,
18    /// Wrapped phase in radians, for cells above the amplitude floor.
19    Phase,
20    /// Squared complex magnitude.
21    MagnitudeSquared,
22}
23
24/// Absolute and relative tolerances for one scalar quantity.
25#[derive(Clone, Copy, Debug, PartialEq)]
26pub struct ScalarTolerance {
27    /// Fixed absolute error allowance.
28    pub absolute: f64,
29    /// Fixed relative error allowance.
30    pub relative: f64,
31}
32
33/// Published fixed tolerances for dense f32 interference conformance.
34#[derive(Clone, Copy, Debug, PartialEq)]
35pub struct DifferentialTolerances {
36    /// Real and imaginary component tolerance.
37    pub component: ScalarTolerance,
38    /// Complex amplitude tolerance.
39    pub amplitude: ScalarTolerance,
40    /// Wrapped phase tolerance in radians.
41    pub phase: ScalarTolerance,
42    /// Squared-magnitude tolerance.
43    pub magnitude_squared: ScalarTolerance,
44    /// Reference amplitude below which phase is undefined and not compared.
45    pub phase_amplitude_floor: f64,
46}
47
48impl Default for DifferentialTolerances {
49    fn default() -> Self {
50        Self {
51            component: ScalarTolerance {
52                absolute: 2.0e-5,
53                relative: 2.0e-4,
54            },
55            amplitude: ScalarTolerance {
56                absolute: 2.0e-5,
57                relative: 2.0e-4,
58            },
59            phase: ScalarTolerance {
60                absolute: 3.0e-4,
61                relative: 1.0e-4,
62            },
63            magnitude_squared: ScalarTolerance {
64                absolute: 4.0e-5,
65                relative: 4.0e-4,
66            },
67            phase_amplitude_floor: 1.0e-5,
68        }
69    }
70}
71
72/// Maximum observed error and its fixed tolerance at one cell.
73#[derive(Clone, Copy, Debug, PartialEq)]
74pub struct DifferentialMaximum {
75    /// Absolute error.
76    pub error: f64,
77    /// Absolute-plus-relative limit at this cell.
78    pub limit: f64,
79    /// Row containing the maximum.
80    pub row: usize,
81    /// Column containing the maximum.
82    pub column: usize,
83}
84
85impl DifferentialMaximum {
86    /// Returns whether the maximum satisfies its fixed tolerance.
87    pub fn passed(self) -> bool {
88        self.error <= self.limit
89    }
90
91    fn ratio(self) -> f64 {
92        if self.limit == 0.0 {
93            if self.error == 0.0 {
94                0.0
95            } else {
96                f64::INFINITY
97            }
98        } else {
99            self.error / self.limit
100        }
101    }
102}
103
104/// Complete comparison report shared by portable and accelerated providers.
105#[derive(Clone, Debug, PartialEq)]
106pub struct DifferentialReport {
107    /// Maximum real-component error.
108    pub real: DifferentialMaximum,
109    /// Maximum imaginary-component error.
110    pub imaginary: DifferentialMaximum,
111    /// Maximum amplitude error.
112    pub amplitude: DifferentialMaximum,
113    /// Maximum wrapped-phase error, absent when every cell is below the floor.
114    pub phase: Option<DifferentialMaximum>,
115    /// Maximum squared-magnitude error.
116    pub magnitude_squared: DifferentialMaximum,
117    /// Number of cells whose phase was compared.
118    pub phase_cells: usize,
119    /// Metric with the largest error-to-limit ratio.
120    pub worst_metric: ConformanceMetric,
121    /// Row containing the worst normalized error.
122    pub worst_row: usize,
123    /// Column containing the worst normalized error.
124    pub worst_column: usize,
125}
126
127impl DifferentialReport {
128    /// Returns true when every compared quantity is within tolerance.
129    pub fn passed(&self) -> bool {
130        self.real.passed()
131            && self.imaginary.passed()
132            && self.amplitude.passed()
133            && self.phase.is_none_or(DifferentialMaximum::passed)
134            && self.magnitude_squared.passed()
135    }
136
137    /// Returns the largest absolute Cartesian-component error.
138    pub fn max_component_absolute_error(&self) -> f64 {
139        self.real.error.max(self.imaginary.error)
140    }
141
142    /// Returns the largest compared absolute wrapped-phase error.
143    ///
144    /// A field containing only below-floor cancellation cells returns zero
145    /// because phase is deliberately undefined there.
146    pub fn max_phase_absolute_error(&self) -> f64 {
147        self.phase.map_or(0.0, |maximum| maximum.error)
148    }
149}
150
151/// Structural or non-finite input rejected before differential reporting.
152#[derive(Clone, Debug, PartialEq)]
153pub struct DifferentialError {
154    detail: String,
155}
156
157impl DifferentialError {
158    fn new(detail: impl Into<String>) -> Self {
159        Self {
160            detail: detail.into(),
161        }
162    }
163
164    /// Returns the stable diagnostic.
165    pub fn detail(&self) -> &str {
166        &self.detail
167    }
168}
169
170impl fmt::Display for DifferentialError {
171    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
172        formatter.write_str(&self.detail)
173    }
174}
175
176impl std::error::Error for DifferentialError {}
177
178/// Compares one dense f32 result with the deterministic f64 reference field.
179pub fn compare_dense_to_reference(
180    reference: &HostPhasorField,
181    candidate: &DenseF32Field,
182    tolerances: DifferentialTolerances,
183) -> Result<DifferentialReport, DifferentialError> {
184    compare_candidate(
185        reference,
186        CandidateField {
187            name: "dense",
188            rows: candidate.rows(),
189            columns: candidate.columns(),
190            real: CandidateComponent::F32(candidate.real()),
191            imaginary: CandidateComponent::F32(candidate.imaginary()),
192        },
193        tolerances,
194    )
195}
196
197/// Compares one materialized provider result with the deterministic f64 oracle.
198///
199/// Provider results use the same component, amplitude, phase-floor, wrapped
200/// phase, and squared-magnitude reporting path as [`compare_dense_to_reference`].
201/// The candidate is expected to be a materialized f32 Tensor field represented
202/// by the runtime's host phasor container.
203pub fn compare_materialized_to_reference(
204    reference: &HostPhasorField,
205    candidate: &HostPhasorField,
206    tolerances: DifferentialTolerances,
207) -> Result<DifferentialReport, DifferentialError> {
208    compare_candidate(
209        reference,
210        CandidateField {
211            name: "materialized",
212            rows: candidate.rows(),
213            columns: candidate.columns(),
214            real: CandidateComponent::F64(candidate.real()),
215            imaginary: CandidateComponent::F64(candidate.imaginary()),
216        },
217        tolerances,
218    )
219}
220
221fn compare_candidate(
222    reference: &HostPhasorField,
223    candidate: CandidateField<'_>,
224    tolerances: DifferentialTolerances,
225) -> Result<DifferentialReport, DifferentialError> {
226    validate_inputs(reference, candidate, tolerances)?;
227    let mut maxima = Maxima::default();
228    let columns = reference.columns();
229    for index in 0..reference.len() {
230        let row = index / columns;
231        let column = index % columns;
232        let rr = reference.real()[index];
233        let ri = reference.imaginary()[index];
234        let cr = candidate.real.get(index);
235        let ci = candidate.imaginary.get(index);
236        let reference_amplitude = rr.hypot(ri);
237        let candidate_amplitude = cr.hypot(ci);
238        maxima.observe(
239            ConformanceMetric::Real,
240            (rr - cr).abs(),
241            limit(tolerances.component, rr, cr),
242            row,
243            column,
244        );
245        maxima.observe(
246            ConformanceMetric::Imaginary,
247            (ri - ci).abs(),
248            limit(tolerances.component, ri, ci),
249            row,
250            column,
251        );
252        maxima.observe(
253            ConformanceMetric::Amplitude,
254            (reference_amplitude - candidate_amplitude).abs(),
255            limit(
256                tolerances.amplitude,
257                reference_amplitude,
258                candidate_amplitude,
259            ),
260            row,
261            column,
262        );
263        if reference_amplitude >= tolerances.phase_amplitude_floor {
264            maxima.phase_cells += 1;
265            let reference_phase = ri.atan2(rr);
266            let candidate_phase = ci.atan2(cr);
267            maxima.observe(
268                ConformanceMetric::Phase,
269                wrapped_phase_error(reference_phase, candidate_phase),
270                limit(tolerances.phase, reference_phase, candidate_phase),
271                row,
272                column,
273            );
274        }
275        let reference_squared = rr.mul_add(rr, ri * ri);
276        let candidate_squared = cr.mul_add(cr, ci * ci);
277        maxima.observe(
278            ConformanceMetric::MagnitudeSquared,
279            (reference_squared - candidate_squared).abs(),
280            limit(
281                tolerances.magnitude_squared,
282                reference_squared,
283                candidate_squared,
284            ),
285            row,
286            column,
287        );
288    }
289    Ok(maxima.finish())
290}
291
292fn validate_inputs(
293    reference: &HostPhasorField,
294    candidate: CandidateField<'_>,
295    tolerances: DifferentialTolerances,
296) -> Result<(), DifferentialError> {
297    if reference.rows() != candidate.rows || reference.columns() != candidate.columns {
298        return Err(DifferentialError::new(format!(
299            "reference shape [{}, {}] differs from {} shape [{}, {}]",
300            reference.rows(),
301            reference.columns(),
302            candidate.name,
303            candidate.rows,
304            candidate.columns
305        )));
306    }
307    if candidate.real.len() != reference.len() || candidate.imaginary.len() != reference.len() {
308        return Err(DifferentialError::new(format!(
309            "{} component lengths [{}, {}] differ from reference length {}",
310            candidate.name,
311            candidate.real.len(),
312            candidate.imaginary.len(),
313            reference.len()
314        )));
315    }
316    for (name, tolerance) in [
317        ("component", tolerances.component),
318        ("amplitude", tolerances.amplitude),
319        ("phase", tolerances.phase),
320        ("magnitude-squared", tolerances.magnitude_squared),
321    ] {
322        if !tolerance.absolute.is_finite()
323            || tolerance.absolute < 0.0
324            || !tolerance.relative.is_finite()
325            || tolerance.relative < 0.0
326        {
327            return Err(DifferentialError::new(format!(
328                "{name} tolerances must be finite and non-negative"
329            )));
330        }
331    }
332    if !tolerances.phase_amplitude_floor.is_finite() || tolerances.phase_amplitude_floor <= 0.0 {
333        return Err(DifferentialError::new(
334            "phase amplitude floor must be finite and positive",
335        ));
336    }
337    for (name, values) in [
338        ("reference real", reference.real()),
339        ("reference imaginary", reference.imaginary()),
340    ] {
341        if let Some((index, value)) = values
342            .iter()
343            .copied()
344            .enumerate()
345            .find(|(_, value)| !value.is_finite())
346        {
347            return Err(DifferentialError::new(format!(
348                "{name} cell {index} is non-finite: {value}"
349            )));
350        }
351    }
352    for (component, values) in [("real", candidate.real), ("imaginary", candidate.imaginary)] {
353        if let Some((index, value)) = values.first_non_finite() {
354            return Err(DifferentialError::new(format!(
355                "{} {component} cell {index} is non-finite: {value}",
356                candidate.name
357            )));
358        }
359    }
360    Ok(())
361}
362
363#[derive(Clone, Copy)]
364enum CandidateComponent<'a> {
365    F32(&'a [f32]),
366    F64(&'a [f64]),
367}
368
369impl CandidateComponent<'_> {
370    fn len(self) -> usize {
371        match self {
372            Self::F32(values) => values.len(),
373            Self::F64(values) => values.len(),
374        }
375    }
376
377    fn get(self, index: usize) -> f64 {
378        match self {
379            Self::F32(values) => f64::from(values[index]),
380            Self::F64(values) => values[index],
381        }
382    }
383
384    fn first_non_finite(self) -> Option<(usize, f64)> {
385        (0..self.len())
386            .map(|index| (index, self.get(index)))
387            .find(|(_, value)| !value.is_finite())
388    }
389}
390
391#[derive(Clone, Copy)]
392struct CandidateField<'a> {
393    name: &'static str,
394    rows: usize,
395    columns: usize,
396    real: CandidateComponent<'a>,
397    imaginary: CandidateComponent<'a>,
398}
399
400fn limit(tolerance: ScalarTolerance, reference: f64, candidate: f64) -> f64 {
401    tolerance.absolute + tolerance.relative * reference.abs().max(candidate.abs())
402}
403
404fn wrapped_phase_error(left: f64, right: f64) -> f64 {
405    let difference = (left - right).abs().rem_euclid(2.0 * PI);
406    difference.min(2.0 * PI - difference)
407}
408
409#[derive(Clone, Copy, Debug)]
410struct ObservedMaximum {
411    value: DifferentialMaximum,
412    initialized: bool,
413}
414
415impl Default for ObservedMaximum {
416    fn default() -> Self {
417        Self {
418            value: DifferentialMaximum {
419                error: 0.0,
420                limit: 0.0,
421                row: 0,
422                column: 0,
423            },
424            initialized: false,
425        }
426    }
427}
428
429impl ObservedMaximum {
430    fn observe(&mut self, error: f64, limit: f64, row: usize, column: usize) {
431        let candidate = DifferentialMaximum {
432            error,
433            limit,
434            row,
435            column,
436        };
437        if !self.initialized
438            || candidate.error > self.value.error
439            || (candidate.error == self.value.error && candidate.ratio() > self.value.ratio())
440        {
441            self.value = candidate;
442            self.initialized = true;
443        }
444    }
445}
446
447#[derive(Default)]
448struct Maxima {
449    real: ObservedMaximum,
450    imaginary: ObservedMaximum,
451    amplitude: ObservedMaximum,
452    phase: ObservedMaximum,
453    magnitude_squared: ObservedMaximum,
454    phase_cells: usize,
455    worst: Option<(ConformanceMetric, DifferentialMaximum)>,
456}
457
458impl Maxima {
459    fn observe(
460        &mut self,
461        metric: ConformanceMetric,
462        error: f64,
463        limit: f64,
464        row: usize,
465        column: usize,
466    ) {
467        let candidate = DifferentialMaximum {
468            error,
469            limit,
470            row,
471            column,
472        };
473        match metric {
474            ConformanceMetric::Real => self.real.observe(error, limit, row, column),
475            ConformanceMetric::Imaginary => self.imaginary.observe(error, limit, row, column),
476            ConformanceMetric::Amplitude => self.amplitude.observe(error, limit, row, column),
477            ConformanceMetric::Phase => self.phase.observe(error, limit, row, column),
478            ConformanceMetric::MagnitudeSquared => {
479                self.magnitude_squared.observe(error, limit, row, column);
480            }
481        }
482        if self
483            .worst
484            .is_none_or(|(_, current)| candidate.ratio() > current.ratio())
485        {
486            self.worst = Some((metric, candidate));
487        }
488    }
489
490    fn finish(self) -> DifferentialReport {
491        let phase = (self.phase_cells != 0).then_some(self.phase.value);
492        let (worst_metric, worst) = self
493            .worst
494            .expect("a valid interference field always has at least one cell");
495        DifferentialReport {
496            real: self.real.value,
497            imaginary: self.imaginary.value,
498            amplitude: self.amplitude.value,
499            phase,
500            magnitude_squared: self.magnitude_squared.value,
501            phase_cells: self.phase_cells,
502            worst_metric,
503            worst_row: worst.row,
504            worst_column: worst.column,
505        }
506    }
507}
sim_lib_interference_compute/diff.rs

sim_lib_interference_compute/
diff.rs
