sim_lib_interference_solve/
analysis.rs

1//! Deterministic scalar-field statistics and fringe candidates.
2
3use std::fmt;
4
5use sim_lib_interference_core::SamplingCertificate;
6
7use crate::{Observable, ProjectionIdentity, ScalarProjection, ScalarSample};
8
9/// Aggregate statistics for every cell in one scalar field.
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub struct FieldStats {
12    /// Number of analyzed cells.
13    pub count: usize,
14    /// Smallest scalar sample.
15    pub minimum: f64,
16    /// Largest scalar sample.
17    pub maximum: f64,
18    /// Row-major online mean.
19    pub mean: f64,
20    /// Population variance computed in deterministic row-major order.
21    pub population_variance: f64,
22}
23
24/// Physical interpretation of a strict local scalar extremum.
25#[derive(Clone, Copy, Debug, PartialEq, Eq)]
26pub enum ExtremumKind {
27    /// Strict local minimum and therefore a possible destructive node.
28    NodeCandidate,
29    /// Strict local maximum and therefore a possible constructive antinode.
30    AntinodeCandidate,
31}
32
33/// One strict local extremum in target-grid coordinates.
34#[derive(Clone, Copy, Debug, PartialEq)]
35pub struct Extremum {
36    /// Zero-based target row.
37    pub row: usize,
38    /// Zero-based target column.
39    pub column: usize,
40    /// Projected scalar value at this cell.
41    pub value: f64,
42    /// Node or antinode interpretation.
43    pub kind: ExtremumKind,
44}
45
46/// Statistics and fringe candidates with inseparable physical provenance.
47#[derive(Clone, Debug, PartialEq)]
48pub struct FringeReport {
49    /// Sampling evidence inherited unchanged from the source solve.
50    pub sampling: SamplingCertificate,
51    /// Exact projection parameters of the analyzed scalar field.
52    pub projection: ProjectionIdentity,
53    /// Whole-field scalar statistics.
54    pub stats: FieldStats,
55    /// Strict local extrema in row-major order.
56    pub extrema: Vec<Extremum>,
57    /// `(maximum - minimum) / (maximum + minimum)`.
58    ///
59    /// This is `None` when the field's maximum amplitude is at or below the
60    /// caller's amplitude floor.
61    pub michelson_contrast: Option<f64>,
62}
63
64/// A scalar projection could not be honestly analyzed as a fringe field.
65#[derive(Clone, Debug, PartialEq)]
66pub enum AnalysisError {
67    /// The amplitude floor was negative or non-finite.
68    InvalidAmplitudeFloor {
69        /// Rejected floor.
70        value: f64,
71    },
72    /// Fringe analysis requires a non-negative amplitude-like observable.
73    IncompatibleObservable {
74        /// Rejected projection observable.
75        observable: Observable,
76    },
77    /// An amplitude-like projection unexpectedly contained a masked cell.
78    MaskedSample {
79        /// Zero-based target row.
80        row: usize,
81        /// Zero-based target column.
82        column: usize,
83    },
84    /// A derived statistic was not finite.
85    NonFiniteStatistic {
86        /// Stable statistic name.
87        name: &'static str,
88        /// Rejected value.
89        value: f64,
90    },
91    /// Extrema storage could not be reserved.
92    AllocationFailed {
93        /// Maximum number of extrema requested.
94        cells: usize,
95    },
96}
97
98impl fmt::Display for AnalysisError {
99    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
100        match self {
101            Self::InvalidAmplitudeFloor { value } => write!(
102                formatter,
103                "analysis amplitude floor must be finite and non-negative: {value:?}"
104            ),
105            Self::IncompatibleObservable { observable } => write!(
106                formatter,
107                "fringe analysis requires amplitude or magnitude-squared, not {observable:?}"
108            ),
109            Self::MaskedSample { row, column } => {
110                write!(formatter, "analysis input is masked at ({row}, {column})")
111            }
112            Self::NonFiniteStatistic { name, value } => {
113                write!(
114                    formatter,
115                    "analysis statistic `{name}` is non-finite: {value:?}"
116                )
117            }
118            Self::AllocationFailed { cells } => {
119                write!(
120                    formatter,
121                    "could not reserve extrema storage for {cells} cells"
122                )
123            }
124        }
125    }
126}
127
128impl std::error::Error for AnalysisError {}
129
130/// Analyzes an amplitude or normalized squared-magnitude projection.
131///
132/// Extrema use the fixed eight-cell Moore neighbourhood, clipped at field
133/// edges. A candidate must be strictly less than or strictly greater than
134/// every available neighbour, so plateaus never produce traversal-dependent
135/// representatives. Results are emitted in row-major order. Michelson
136/// contrast is omitted when the maximum represented amplitude is at or below
137/// `amplitude_floor`.
138pub fn analyze_fringes(
139    projection: &ScalarProjection,
140    amplitude_floor: f64,
141) -> Result<FringeReport, AnalysisError> {
142    if !amplitude_floor.is_finite() || amplitude_floor < 0.0 {
143        return Err(AnalysisError::InvalidAmplitudeFloor {
144            value: amplitude_floor,
145        });
146    }
147    let certificate = projection.certificate();
148    let observable = certificate.observable();
149    if !matches!(
150        observable,
151        Observable::Amplitude | Observable::MagnitudeSquared
152    ) {
153        return Err(AnalysisError::IncompatibleObservable { observable });
154    }
155
156    let stats = field_stats(projection)?;
157    let extrema = find_extrema(projection)?;
158    let maximum_amplitude = match observable {
159        Observable::Amplitude => stats.maximum,
160        Observable::MagnitudeSquared => stats.maximum.sqrt(),
161        _ => unreachable!("observable was admitted above"),
162    };
163    let michelson_contrast = if maximum_amplitude <= amplitude_floor {
164        None
165    } else {
166        let minimum_to_maximum = stats.minimum / stats.maximum;
167        let contrast = (1.0 - minimum_to_maximum) / (1.0 + minimum_to_maximum);
168        Some(require_finite("michelson-contrast", contrast)?)
169    };
170
171    Ok(FringeReport {
172        sampling: certificate.source_sampling_certificate(),
173        projection: certificate.identity(),
174        stats,
175        extrema,
176        michelson_contrast,
177    })
178}
179
180fn field_stats(projection: &ScalarProjection) -> Result<FieldStats, AnalysisError> {
181    let mut minimum = scalar_value(projection, 0)?;
182    let mut maximum = minimum;
183    let mut mean = 0.0;
184    let mut squared_deviation = 0.0;
185    for index in 0..projection.samples().len() {
186        let value = scalar_value(projection, index)?;
187        minimum = minimum.min(value);
188        maximum = maximum.max(value);
189        let count = (index + 1) as f64;
190        let delta = value - mean;
191        mean += delta / count;
192        squared_deviation += delta * (value - mean);
193        require_finite("mean", mean)?;
194        require_finite("squared-deviation", squared_deviation)?;
195    }
196    let population_variance = require_finite(
197        "population-variance",
198        squared_deviation / projection.samples().len() as f64,
199    )?;
200    Ok(FieldStats {
201        count: projection.samples().len(),
202        minimum,
203        maximum,
204        mean,
205        population_variance,
206    })
207}
208
209fn find_extrema(projection: &ScalarProjection) -> Result<Vec<Extremum>, AnalysisError> {
210    let rows = projection.rows();
211    let columns = projection.columns();
212    let mut extrema = Vec::new();
213    extrema
214        .try_reserve(projection.samples().len())
215        .map_err(|_| AnalysisError::AllocationFailed {
216            cells: projection.samples().len(),
217        })?;
218    for row in 0..rows {
219        for column in 0..columns {
220            let value = scalar_value(projection, row * columns + column)?;
221            let mut has_neighbour = false;
222            let mut below_all = true;
223            let mut above_all = true;
224            for neighbour_row in row.saturating_sub(1)..=(row + 1).min(rows - 1) {
225                for neighbour_column in column.saturating_sub(1)..=(column + 1).min(columns - 1) {
226                    if neighbour_row == row && neighbour_column == column {
227                        continue;
228                    }
229                    has_neighbour = true;
230                    let neighbour =
231                        scalar_value(projection, neighbour_row * columns + neighbour_column)?;
232                    below_all &= value < neighbour;
233                    above_all &= value > neighbour;
234                }
235            }
236            let kind = match (has_neighbour && below_all, has_neighbour && above_all) {
237                (true, false) => Some(ExtremumKind::NodeCandidate),
238                (false, true) => Some(ExtremumKind::AntinodeCandidate),
239                _ => None,
240            };
241            if let Some(kind) = kind {
242                extrema.push(Extremum {
243                    row,
244                    column,
245                    value,
246                    kind,
247                });
248            }
249        }
250    }
251    Ok(extrema)
252}
253
254fn scalar_value(projection: &ScalarProjection, index: usize) -> Result<f64, AnalysisError> {
255    match projection.samples()[index] {
256        ScalarSample::Value(value) => Ok(value),
257        ScalarSample::Masked => Err(AnalysisError::MaskedSample {
258            row: index / projection.columns(),
259            column: index % projection.columns(),
260        }),
261    }
262}
263
264fn require_finite(name: &'static str, value: f64) -> Result<f64, AnalysisError> {
265    value
266        .is_finite()
267        .then_some(value)
268        .ok_or(AnalysisError::NonFiniteStatistic { name, value })
269}
270
271#[cfg(test)]
272#[path = "analysis_tests.rs"]
273mod tests;
sim_lib_interference_solve/analysis.rs

sim_lib_interference_solve/
analysis.rs
