From 6784bd168d381ab48082284073ec9130c5f61796 Mon Sep 17 00:00:00 2001 From: Jiekang Tian Date: Thu, 27 Aug 2026 18:32:09 +0800 Subject: [PATCH] feat(craft): add validated modeling and amplitude reports --- crates/analysis/src/craft.rs | 84 +++ crates/analysis/src/craft_tests.rs | 50 ++ crates/app/src/shot/craft_shot.rs | 3 +- crates/app/src/ui/canvas/craft_results.rs | 142 +++++ crates/app/src/ui/commands_craft_tests.rs | 2 +- crates/app/src/ui/tools/craft.rs | 28 +- crates/app/src/ui/tools/craft/results.rs | 290 ++++++++-- .../src/ui/tools/craft/results_diagnostics.rs | 32 ++ crates/app/src/ui/tools/craft/setup.rs | 87 +-- crates/core/src/project/craft_tests.rs | 161 +++++- crates/core/src/project/mod.rs | 19 + crates/core/src/state/app_impl.rs | 2 + crates/core/src/state/app_impl_compute.rs | 43 ++ .../core/src/state/app_impl_compute_tests.rs | 14 +- crates/core/src/state/craft.rs | 94 +++- crates/core/src/state/mod.rs | 2 + crates/core/src/state/reports.rs | 191 +++++++ crates/core/src/state/ui_state.rs | 4 + crates/core/src/state/ui_state_craft.rs | 1 + crates/processing/src/craft.rs | 499 ++++++------------ crates/processing/src/craft/diagnostics.rs | 70 ++- crates/processing/src/craft/fitting.rs | 496 +++++++++++++++++ crates/processing/src/craft/preflight.rs | 44 +- crates/processing/src/craft/regions.rs | 126 ++++- crates/processing/src/craft/report.rs | 222 ++++++++ crates/processing/src/craft/resolution.rs | 237 +++++---- crates/processing/src/craft/stability.rs | 193 +++++++ crates/processing/src/craft_tests.rs | 362 ++++++++++--- docs/src/content/docs/guides/craft.md | 47 +- docs/src/content/docs/zh-cn/guides/craft.md | 39 +- 30 files changed, 2846 insertions(+), 738 deletions(-) create mode 100644 crates/app/src/ui/tools/craft/results_diagnostics.rs create mode 100644 crates/core/src/state/reports.rs create mode 100644 crates/processing/src/craft/fitting.rs create mode 100644 crates/processing/src/craft/report.rs create mode 100644 crates/processing/src/craft/stability.rs diff --git a/crates/analysis/src/craft.rs b/crates/analysis/src/craft.rs index ef170e3..772620b 100644 --- a/crates/analysis/src/craft.rs +++ b/crates/analysis/src/craft.rs @@ -94,6 +94,90 @@ pub enum CraftFitError { Singular, } +/// Replace the leading samples of a complex record by backward linear +/// prediction from the immediately following observed samples. +/// +/// CRAFT uses this after digital filtering because a finite FIR must invent a +/// prehistory at the acquisition boundary. The prediction is fitted in reverse +/// time, so the supplied autoregressive order has the same meaning as a +/// conventional forward linear predictor. +pub fn backward_linear_predict( + samples: &mut [Complex64], + predicted_count: usize, + training_count: usize, + order: usize, +) -> Result<(), CraftFitError> { + if predicted_count == 0 { + return Ok(()); + } + if order == 0 + || training_count <= order + || predicted_count + .checked_add(training_count) + .is_none_or(|required| required > samples.len()) + || samples + .iter() + .take(predicted_count + training_count) + .any(|value| !value.re.is_finite() || !value.im.is_finite()) + { + return Err(CraftFitError::InvalidInput); + } + + let training = &samples[predicted_count..predicted_count + training_count]; + let reversed = training.iter().rev().copied().collect::>(); + let scale = reversed + .iter() + .map(|value| value.norm()) + .fold(0.0_f64, f64::max); + if scale <= f64::MIN_POSITIVE { + return Err(CraftFitError::Singular); + } + let equation_count = reversed.len() - order; + let mut design = DMatrix::::zeros(equation_count * 2, order * 2); + let mut observed = DVector::::zeros(equation_count * 2); + for row in 0..equation_count { + let target = reversed[row + order]; + observed[row * 2] = target.re / scale; + observed[row * 2 + 1] = target.im / scale; + for lag in 0..order { + let basis = reversed[row + order - lag - 1] / scale; + design[(row * 2, lag * 2)] = basis.re; + design[(row * 2, lag * 2 + 1)] = -basis.im; + design[(row * 2 + 1, lag * 2)] = basis.im; + design[(row * 2 + 1, lag * 2 + 1)] = basis.re; + } + } + // Scale the singular-value cutoff by the design energy so rank detection + // remains stable across differently normalized input records. + let rank_tolerance = (5e-14 * design.norm_squared()).sqrt().max(1e-12); + let solution = design + .svd(true, true) + .solve(&observed, rank_tolerance) + .map_err(|_| CraftFitError::Singular)?; + let coefficients = solution + .as_slice() + .as_chunks::<2>() + .0 + .iter() + .map(|pair| Complex64::new(pair[0], pair[1])) + .collect::>(); + let mut history = reversed; + for index in 0..predicted_count { + let predicted = coefficients + .iter() + .enumerate() + .fold(Complex64::new(0.0, 0.0), |sum, (lag, coefficient)| { + sum + coefficient * history[history.len() - lag - 1] + }); + if !predicted.re.is_finite() || !predicted.im.is_finite() { + return Err(CraftFitError::Singular); + } + samples[predicted_count - index - 1] = predicted; + history.push(predicted); + } + Ok(()) +} + /// Fit a fixed set of initial component frequencies. Model-order selection and /// residual candidate discovery live in `plotx-processing`, beside its FFT. pub fn fit_damped_sinusoids_cancellable( diff --git a/crates/analysis/src/craft_tests.rs b/crates/analysis/src/craft_tests.rs index 4d24cfc..ffcfa52 100644 --- a/crates/analysis/src/craft_tests.rs +++ b/crates/analysis/src/craft_tests.rs @@ -23,6 +23,56 @@ fn synthetic( (times, samples) } +#[test] +fn backward_prediction_restores_filtered_record_leading_points() { + let components = [ + (13.0, 4.0, 0.3, 0.8), + (-21.0, 2.5, -0.4, 1.7), + (37.0, 1.2, 1.1, 2.4), + ]; + let (_, expected) = synthetic(&components, 300, 500.0); + let mut samples = expected.clone(); + samples[..5].fill(Complex64::new(100.0, -50.0)); + + backward_linear_predict(&mut samples, 5, 256, 32).unwrap(); + + for index in 0..5 { + assert!( + (samples[index] - expected[index]).norm() < 1e-7, + "index={index} predicted={:?} expected={:?}", + samples[index], + expected[index] + ); + } +} + +#[test] +fn backward_prediction_restores_a_short_single_exponential() { + let (_, expected) = synthetic(&[(0.0, 3.0, 0.4, 5.0)], 192, 4_000.0); + let mut samples = expected.clone(); + samples[..5].fill(Complex64::new(100.0, -50.0)); + + backward_linear_predict(&mut samples, 5, 187, 16).unwrap(); + + for index in 0..5 { + assert!( + (samples[index] - expected[index]).norm() < 1e-7, + "index={index} predicted={:?} expected={:?}", + samples[index], + expected[index] + ); + } +} + +#[test] +fn backward_prediction_rejects_an_underspecified_fit() { + let mut samples = vec![Complex64::new(1.0, 0.0); 12]; + assert_eq!( + backward_linear_predict(&mut samples, 5, 7, 7), + Err(CraftFitError::InvalidInput) + ); +} + #[test] fn recovers_single_damped_sinusoid() { let (times, samples) = synthetic(&[(123.4, 7.5, 0.37, 2.2)], 2048, 2000.0); diff --git a/crates/app/src/shot/craft_shot.rs b/crates/app/src/shot/craft_shot.rs index e1ce8dc..d152913 100644 --- a/crates/app/src/shot/craft_shot.rs +++ b/crates/app/src/shot/craft_shot.rs @@ -12,8 +12,7 @@ pub(super) fn setup(app: &mut PlotxApp, ctx: &egui::Context) -> Result<(), Strin .data .clone(); let mut params = CraftParams::conventional(); - params.max_fit_window_width_hz = data.spectral_width_hz; - params.max_components_per_fit_window = 8; + params.maximum_model_order = 8; let invocation = CraftInvocation::acquisition(&data, params); let result = process_craft_cancellable(&data, &invocation, &|| false) .map_err(|error| format!("CRAFT screenshot analysis failed: {error}"))?; diff --git a/crates/app/src/ui/canvas/craft_results.rs b/crates/app/src/ui/canvas/craft_results.rs index 5e54a47..8cf2b3c 100644 --- a/crates/app/src/ui/canvas/craft_results.rs +++ b/crates/app/src/ui/canvas/craft_results.rs @@ -45,6 +45,17 @@ pub(crate) fn handle_and_paint_craft_result( else { return; }; + paint_craft_ranges(CraftRangePaintContext { + app, + dataset, + run, + stored, + nmr, + plot, + figure, + painter, + ui, + }); if let Some(selected) = app.session.ui.craft_selected_component && let Some(component) = stored .components @@ -102,3 +113,134 @@ pub(crate) fn handle_and_paint_craft_result( .open_task_tab(plotx_core::state::TaskDockTab::Craft); } } + +struct CraftRangePaintContext<'a> { + app: &'a PlotxApp, + dataset: plotx_core::state::DatasetId, + run: plotx_core::state::CraftRunId, + stored: &'a plotx_core::state::StoredCraftRun, + nmr: &'a plotx_core::state::NmrDataset, + plot: PlotRect, + figure: &'a plotx_figure::Figure, + painter: &'a egui::Painter, + ui: &'a Ui, +} + +fn paint_craft_ranges(context: CraftRangePaintContext<'_>) { + let CraftRangePaintContext { + app, + dataset, + run, + stored, + nmr, + plot, + figure, + painter, + ui, + } = context; + let carrier = stored + .provenance + .invocation + .reference + .effective_carrier_ppm(); + let observe = nmr.data.observe_freq_mhz; + let modeling = stored + .diagnostics + .modeling_windows + .iter() + .map(|window| { + ( + carrier + window.modeling_band_hz.0 / observe, + carrier + window.modeling_band_hz.1 / observe, + ) + }) + .collect::>(); + let regions = stored + .region_summaries + .iter() + .map(|region| (region.start_ppm, region.end_ppm)) + .collect::>(); + let report_segments = app + .session + .ui + .craft_selected_report + .and_then(|id| app.doc.report(id)) + .filter(|record| { + record.source + == plotx_core::state::ReportSource { + dataset, + craft_run: run, + } + }) + .and_then(|record| { + serde_json::from_value::( + record.snapshot.clone(), + ) + .ok() + }) + .map(|report| { + report + .segments + .into_iter() + .map(|segment| { + ( + carrier + segment.start_hz / observe, + carrier + segment.end_hz / observe, + ) + }) + .collect::>() + }) + .unwrap_or_default(); + + paint_range_track( + &modeling, + plot.top + 2.0, + plot, + figure, + painter, + ui.visuals().weak_text_color().linear_multiply(0.45), + ); + paint_range_track( + ®ions, + plot.top + 7.0, + plot, + figure, + painter, + ui.visuals().selection.stroke.color.linear_multiply(0.75), + ); + paint_range_track( + &report_segments, + plot.top + 12.0, + plot, + figure, + painter, + ui.visuals().warn_fg_color.linear_multiply(0.75), + ); +} + +fn paint_range_track( + ranges: &[(f64, f64)], + y: f32, + plot: PlotRect, + figure: &plotx_figure::Figure, + painter: &egui::Painter, + color: egui::Color32, +) { + for &(left, right) in ranges { + let first = x_to_screen(left, plot, figure.x.min, figure.x.span(), figure.x.reversed); + let second = x_to_screen( + right, + plot, + figure.x.min, + figure.x.span(), + figure.x.reversed, + ); + let rect = egui::Rect::from_min_max( + Pos2::new(first.min(second).max(plot.left), y), + Pos2::new(first.max(second).min(plot.right()), y + 3.0), + ); + if rect.is_positive() { + painter.rect_filled(rect, 0.0, color); + } + } +} diff --git a/crates/app/src/ui/commands_craft_tests.rs b/crates/app/src/ui/commands_craft_tests.rs index fdd5de1..d54442a 100644 --- a/crates/app/src/ui/commands_craft_tests.rs +++ b/crates/app/src/ui/commands_craft_tests.rs @@ -2,7 +2,7 @@ use super::tests::app_with_nmr; use super::*; fn use_short_fixture_filter(app: &mut PlotxApp) { - app.session.ui.craft_overrides.filter_taps = Some(31); + app.session.ui.craft_overrides.fir_filter_taps = Some(31); } #[test] diff --git a/crates/app/src/ui/tools/craft.rs b/crates/app/src/ui/tools/craft.rs index 1cdc797..0f61ac9 100644 --- a/crates/app/src/ui/tools/craft.rs +++ b/crates/app/src/ui/tools/craft.rs @@ -8,6 +8,7 @@ use super::task_card::{self, TaskCardGeometry}; use crate::ui::commands::{self, CommandId}; mod results; +mod results_diagnostics; mod setup; mod spectrum; @@ -304,6 +305,18 @@ fn command_button(app: &mut PlotxApp, command: CommandId, label: &str, primary: #[cfg(test)] mod tests { use super::*; + + fn preview_sample_indices(point_count: usize, sample_count: usize) -> Vec { + let count = point_count.min(sample_count.max(2)); + if count == 0 { + return Vec::new(); + } + if count == 1 { + return vec![0]; + } + let last = point_count - 1; + (0..count).map(|index| index * last / (count - 1)).collect() + } use num_complex::Complex64; use plotx_core::state::{CraftRunId, NmrDataset, StoredCraftRun}; use plotx_io::{Domain, NmrData}; @@ -340,8 +353,9 @@ mod tests { residual_rss: 1.0, normalized_residual: 1.0, maximum_condition_number: Some(1.0), - fit_windows: Vec::new(), + modeling_windows: Vec::new(), warnings: Vec::new(), + stability: Default::default(), }, synthetic_fid: Vec::new(), residual_fid: Vec::new(), @@ -354,7 +368,7 @@ mod tests { let first = NmrDataset::load(time_domain_data("first")); let mut second = NmrDataset::load(time_domain_data("second")); let mut provenance_params = CraftParams::ssfp(); - provenance_params.min_amplitude_to_noise = 8.5; + provenance_params.minimum_amplitude_to_noise = 8.5; second .craft_runs .push(stored_run(&second.data, provenance_params.clone())); @@ -364,10 +378,10 @@ mod tests { app.set_active_dataset(Some(0)); open_for_active(&mut app); - app.session.ui.craft_overrides.min_amplitude_to_noise = Some(6.0); + app.session.ui.craft_overrides.minimum_amplitude_to_noise = Some(6.0); open_for_active(&mut app); assert_eq!( - app.session.ui.craft_overrides.min_amplitude_to_noise, + app.session.ui.craft_overrides.minimum_amplitude_to_noise, Some(6.0) ); @@ -382,7 +396,7 @@ mod tests { #[test] fn preview_indices_cover_endpoints_with_a_bounded_sample_count() { - let indices = results::preview_sample_indices(65_536, 310); + let indices = preview_sample_indices(65_536, 310); assert_eq!(indices.len(), 310); assert_eq!(indices.first(), Some(&0)); @@ -408,7 +422,7 @@ mod tests { &dataset.data, dataset.craft_reference(), &plotx_processing::craft::CraftParamOverrides { - filter_taps: Some(31), + fir_filter_taps: Some(31), ..Default::default() }, None, @@ -427,7 +441,7 @@ mod tests { } #[test] - fn detected_signal_width_is_independent_of_internal_fit_window_width() { + fn detected_signal_width_is_independent_of_modeling_bandwidth() { assert!((45.0 / 600.0_f64 - 0.075).abs() < f64::EPSILON); } } diff --git a/crates/app/src/ui/tools/craft/results.rs b/crates/app/src/ui/tools/craft/results.rs index f8247c2..a693e6f 100644 --- a/crates/app/src/ui/tools/craft/results.rs +++ b/crates/app/src/ui/tools/craft/results.rs @@ -3,6 +3,7 @@ use plotx_core::state::{ CraftAnalysisIntent, CraftComponentSort, CraftResultTab, CraftTaskPage, PlotxApp, StoredCraftRun, }; +use plotx_processing::craft::{CraftAmplitudeReport, CraftReportDefinition}; use plotx_processing::craft::{CraftComponent, CraftProfile, CraftRegionId, CraftRunStatus}; use crate::ui::commands::CommandId; @@ -123,6 +124,11 @@ pub(super) fn show(app: &mut PlotxApp, index: usize, ui: &mut Ui) { CraftResultTab::Diagnostics, "Diagnostics", ); + ui.selectable_value( + &mut app.session.ui.craft_result_tab, + CraftResultTab::Reports, + "Reports", + ); }); ui.separator(); @@ -130,9 +136,223 @@ pub(super) fn show(app: &mut PlotxApp, index: usize, ui: &mut Ui) { CraftResultTab::Overview => overview(app, &nmr, &run, ui), CraftResultTab::Components => components(app, &nmr, &run, ui), CraftResultTab::Diagnostics => diagnostics(app, &nmr, &run, ui), + CraftResultTab::Reports => reports(app, index, &nmr, &run, ui), } } +fn reports( + app: &mut PlotxApp, + index: usize, + nmr: &plotx_core::state::NmrDataset, + run: &StoredCraftRun, + ui: &mut Ui, +) { + let source = plotx_core::state::ReportSource { + dataset: nmr.resource_id, + craft_run: run.id, + }; + let report_ids = app + .doc + .reports_for_source(source) + .map(|r| r.id) + .collect::>(); + let quantitative_ready = run.diagnostics.status == CraftRunStatus::Complete + && run.diagnostics.stability.passed + && !run.is_stale_for(&nmr.data, nmr.craft_reference()); + ui.horizontal_wrapped(|ui| { + if ui + .add_enabled(quantitative_ready, egui::Button::new("New report")) + .on_disabled_hover_text( + "A stable, current CRAFT run is required for a quantitative amplitude report.", + ) + .clicked() + { + let definition = CraftReportDefinition { + threshold_an: run.provenance.invocation.params.minimum_amplitude_to_noise, + segment_width_hz: 1.0, + regions: Vec::new(), + }; + if let Ok(snapshot) = run.amplitude_report(definition.clone()) + && let (Ok(definition), Ok(snapshot)) = ( + serde_json::to_value(definition), + serde_json::to_value(snapshot), + ) + { + let id = app.doc.create_report(plotx_core::state::NewAnalysisReport { + name: format!("CRAFT report {}", report_ids.len() + 1), + kind: plotx_core::state::ReportKindId::new("craft_amplitude"), + source, + definition, + snapshot, + source_fingerprint: run.provenance.input_sha256.clone(), + schema_version: 1, + }); + app.session.ui.craft_selected_report = Some(id); + } + } + if !report_ids.is_empty() { + let mut selected = app + .session + .ui + .craft_selected_report + .filter(|id| report_ids.contains(id)) + .or_else(|| report_ids.first().copied()); + egui::ComboBox::from_id_salt(("craft_report", run.id.0)) + .selected_text( + selected + .map(|id| format!("Report {}", id.0 + 1)) + .unwrap_or_default(), + ) + .show_ui(ui, |ui| { + for id in &report_ids { + ui.selectable_value( + &mut selected, + Some(*id), + format!("Report {}", id.0 + 1), + ); + } + }); + app.session.ui.craft_selected_report = selected; + if let Some(id) = selected { + if ui.button("Delete").clicked() { + app.doc.delete_report(id); + app.session.ui.craft_selected_report = None; + return; + } + if ui.button("Copy").clicked() + && let Ok(copy) = app.doc.copy_report(id, None) + { + app.session.ui.craft_selected_report = Some(copy); + } + if ui.button("Rename").clicked() { + let _ = app + .doc + .rename_report(id, format!("CRAFT report {}", id.0 + 1)); + } + } + } + }); + let Some(id) = app.session.ui.craft_selected_report else { + ui.weak("Create a report to summarize trusted CRAFT components."); + return; + }; + let Some(record) = app.doc.report(id).cloned() else { + return; + }; + match record.status(&app.doc) { + plotx_core::state::ReportStatus::Unavailable => { + ui.colored_label( + ui.visuals().error_fg_color, + "Source CRAFT run is unavailable.", + ); + return; + } + plotx_core::state::ReportStatus::NeedsReview => { + ui.colored_label( + ui.visuals().warn_fg_color, + "Source CRAFT run changed. Review or recreate this report.", + ); + } + plotx_core::state::ReportStatus::Available => {} + } + let mut definition: CraftReportDefinition = + serde_json::from_value(record.definition.clone()).unwrap_or_default(); + let mut changed = false; + ui.horizontal(|ui| { + ui.label("Report threshold A/N"); + changed |= ui + .add( + egui::DragValue::new(&mut definition.threshold_an) + .speed(0.1) + .range(0.001..=1_000.0), + ) + .changed(); + ui.label("Segment width"); + changed |= ui + .add( + egui::DragValue::new(&mut definition.segment_width_hz) + .speed(0.1) + .range(0.001..=1_000_000.0), + ) + .changed(); + ui.label(format!( + "Hz ({:.5} ppm)", + definition.segment_width_hz / nmr.data.observe_freq_mhz + )); + }); + let mut snapshot: CraftAmplitudeReport = serde_json::from_value(record.snapshot.clone()) + .unwrap_or(CraftAmplitudeReport { + schema_version: 1, + definition: definition.clone(), + segments: Vec::new(), + }); + if changed && let Ok(generated_snapshot) = run.amplitude_report(definition.clone()) { + snapshot = generated_snapshot.clone(); + let mut updated = record.clone(); + updated.definition = + serde_json::to_value(&definition).unwrap_or_else(|_| record.definition.clone()); + updated.snapshot = + serde_json::to_value(generated_snapshot).unwrap_or_else(|_| record.snapshot.clone()); + let _ = app.doc.update_report(updated); + } + let component_count: usize = snapshot.segments.iter().map(|s| s.component_count).sum(); + let scalar: f64 = snapshot + .segments + .iter() + .map(|s| s.scalar_amplitude_sum_t0) + .sum(); + let coherent: f64 = snapshot + .segments + .iter() + .map(|s| s.coherent_amplitude_t0) + .sum(); + ui.small(format!( + "{} segment(s) · {} component(s) · scalar {:.5} · coherent {:.5}", + snapshot.segments.len(), + component_count, + scalar, + coherent + )); + if ui + .add_enabled(quantitative_ready, egui::Button::new("Export report…")) + .on_disabled_hover_text( + "Quantitative export is unavailable until CRAFT stability checks pass.", + ) + .clicked() + { + match app.materialize_craft_report_table(index, id) { + Ok(table) => app.open_data_export(table), + Err(message) => app.session.status = message, + } + } + egui::Grid::new(("craft_report_table", id.0)) + .striped(true) + .show(ui, |ui| { + ui.strong("Segment"); + ui.strong("Bounds (Hz)"); + ui.strong("Components"); + ui.strong("Scalar"); + ui.strong("Coherent"); + ui.end_row(); + for (i, segment) in snapshot.segments.iter().enumerate() { + ui.label((i + 1).to_string()); + ui.label(format!("{:.4} .. {:.4}", segment.start_hz, segment.end_hz)); + ui.label( + segment + .component_ids + .iter() + .map(|id| id.0.to_string()) + .collect::>() + .join(", "), + ); + ui.label(format!("{:.5}", segment.scalar_amplitude_sum_t0)); + ui.label(format!("{:.5}", segment.coherent_amplitude_t0)); + ui.end_row(); + } + }); + let _ = index; +} + fn prepare_rerun(app: &mut PlotxApp, run: &StoredCraftRun) { app.session.ui.craft_base_run = Some(run.id); app.session.ui.craft_overrides = Default::default(); @@ -201,6 +421,21 @@ fn overview( run.components.len(), run.diagnostics.normalized_residual )); + ui.small(format!( + "Fixed protocol: {:.0} Hz modeling bandwidth · boundary dispersion: {:.2}%", + run.provenance + .invocation + .params + .profile + .modeling_bandwidth_hz(), + run.diagnostics + .stability + .regions + .iter() + .map(|region| region.metric.relative_dispersion) + .fold(0.0, f64::max) + * 100.0, + )); ui.small(format!( "Chemical-shift reference {:+.5} ppm · effective carrier {:.5} ppm", run.provenance.invocation.reference.offset_ppm, @@ -514,22 +749,22 @@ fn diagnostics( )); ui.small(format!( "A/N {:.2} ({:?}) · model limit {} ({:?}) · linewidth {:.3}–{:.3} Hz ({:?})", - invocation.params.min_amplitude_to_noise, - invocation.sources.min_amplitude_to_noise, - invocation.params.max_components_per_fit_window, - invocation.sources.max_components_per_fit_window, - invocation.params.linewidth_hz.0, - invocation.params.linewidth_hz.1, - invocation.sources.linewidth_hz, + invocation.params.minimum_amplitude_to_noise, + invocation.sources.minimum_amplitude_to_noise, + invocation.params.maximum_model_order, + invocation.sources.maximum_model_order, + invocation.params.component_linewidth_bounds_hz.0, + invocation.params.component_linewidth_bounds_hz.1, + invocation.sources.component_linewidth_bounds_hz, )); ui.small(format!( - "Skip {} points ({:?}) · FIR {} taps · {} available · {} reconstructed · {} fit window(s)", + "Skip {} points ({:?}) · FIR {} taps · {} available · {} reconstructed · {} modeling window(s)", invocation.derived_plan.effective_skip_points, invocation.derived_plan.effective_skip_source, - invocation.derived_plan.actual_filter_taps, + invocation.derived_plan.effective_fir_filter_taps, invocation.derived_plan.available_points, invocation.derived_plan.reconstruction_points, - invocation.derived_plan.fit_windows.len(), + invocation.derived_plan.modeling_windows.len(), )); for issue in &invocation.assessment.issues { ui.colored_label( @@ -541,27 +776,7 @@ fn diagnostics( ); } }); - if !run.diagnostics.fit_windows.is_empty() { - ui.collapsing("Fit-window BIC", |ui| { - for (window, diagnostic) in run.diagnostics.fit_windows.iter().enumerate() { - ui.small(format!( - "Window {} · Region {} · order {}/{} · decimation {} · {} samples · BIC {} · condition {}", - window + 1, - region_number(run, diagnostic.region), - diagnostic.selected_model_order, - diagnostic.evaluated_model_orders, - diagnostic.actual_decimation, - diagnostic.retained_samples, - diagnostic - .bic - .map_or_else(|| "unavailable".into(), |value| format!("{value:.4}")), - diagnostic - .condition_number - .map_or_else(|| "unavailable".into(), |value| format!("{value:.3e}")), - )); - } - }); - } + super::results_diagnostics::show_modeling_windows(run, ui); } fn region_number(run: &StoredCraftRun, id: CraftRegionId) -> usize { @@ -570,16 +785,3 @@ fn region_number(run: &StoredCraftRun, id: CraftRegionId) -> usize { .position(|summary| summary.region == id) .map_or(0, |position| position + 1) } - -#[cfg(test)] -pub(super) fn preview_sample_indices(point_count: usize, sample_count: usize) -> Vec { - let count = point_count.min(sample_count.max(2)); - if count == 0 { - return Vec::new(); - } - if count == 1 { - return vec![0]; - } - let last = point_count - 1; - (0..count).map(|index| index * last / (count - 1)).collect() -} diff --git a/crates/app/src/ui/tools/craft/results_diagnostics.rs b/crates/app/src/ui/tools/craft/results_diagnostics.rs new file mode 100644 index 0000000..dccba2d --- /dev/null +++ b/crates/app/src/ui/tools/craft/results_diagnostics.rs @@ -0,0 +1,32 @@ +use egui::Ui; +use plotx_core::state::StoredCraftRun; + +pub(super) fn show_modeling_windows(run: &StoredCraftRun, ui: &mut Ui) { + if run.diagnostics.modeling_windows.is_empty() { + return; + } + ui.collapsing("Modeling-window validation", |ui| { + for (window, diagnostic) in run.diagnostics.modeling_windows.iter().enumerate() { + let training_bic = diagnostic + .training_bic + .map_or_else(|| "unavailable".into(), |value| format!("{value:.4}")); + let condition = diagnostic + .condition_number + .map_or_else(|| "unavailable".into(), |value| format!("{value:.3e}")); + ui.small(format!( + "Window {} · retain {:.1}..{:.1} Hz · model {:.1}..{:.1} Hz · order {}/{} · decimation {} · {} samples · training residual {:.3} · validation residual {:.3} · training BIC {training_bic} · condition {condition}", + window + 1, + diagnostic.retention_band_hz.0, + diagnostic.retention_band_hz.1, + diagnostic.modeling_band_hz.0, + diagnostic.modeling_band_hz.1, + diagnostic.selected_model_order, + diagnostic.evaluated_model_orders, + diagnostic.decimation_factor, + diagnostic.modeled_sample_count, + diagnostic.training_normalized_residual, + diagnostic.validation_normalized_residual, + )); + } + }); +} diff --git a/crates/app/src/ui/tools/craft/setup.rs b/crates/app/src/ui/tools/craft/setup.rs index 07ea9ca..fc08d60 100644 --- a/crates/app/src/ui/tools/craft/setup.rs +++ b/crates/app/src/ui/tools/craft/setup.rs @@ -69,10 +69,10 @@ fn readiness(intent: CraftAnalysisIntent, invocation: &CraftInvocation, ui: &mut }; ui.colored_label(color, crate::typography::headline(label)); ui.small(format!( - "{} points ({} usable) · {duration} · {} fit window(s) · {} clear signal(s)", + "{} points ({} usable) · {duration} · {} modeling window(s) · {} clear signal(s)", assessment.point_count, assessment.effective_point_count, - assessment.fit_window_count, + assessment.modeling_window_count, assessment.clear_signals.len(), )); for issue in &assessment.issues { @@ -179,43 +179,43 @@ fn settings(app: &mut PlotxApp, index: usize, invocation: &CraftInvocation, ui: ui.add_space(8.0); ui.label(crate::typography::headline( - "3. Confirm acquisition and fit settings", + "3. Confirm acquisition and component settings", )); - ui.collapsing("Advanced fit settings", |ui| { - let mut value = invocation.params.min_amplitude_to_noise; + ui.collapsing("Advanced component settings", |ui| { + let mut value = invocation.params.minimum_amplitude_to_noise; if setting_row( ui, "Minimum A/N", - invocation.sources.min_amplitude_to_noise, + invocation.sources.minimum_amplitude_to_noise, &nmr, - &mut overrides.min_amplitude_to_noise, + &mut overrides.minimum_amplitude_to_noise, |ui| { ui.add(DragValue::new(&mut value).range(0.1..=100.0).speed(0.1)) .changed() }, ) { - overrides.min_amplitude_to_noise = Some(value); + overrides.minimum_amplitude_to_noise = Some(value); } - let mut value = invocation.params.max_components_per_fit_window; + let mut value = invocation.params.maximum_model_order; if setting_row( ui, - "Max components / fit window", - invocation.sources.max_components_per_fit_window, + "Maximum model order", + invocation.sources.maximum_model_order, &nmr, - &mut overrides.max_components_per_fit_window, + &mut overrides.maximum_model_order, |ui| ui.add(DragValue::new(&mut value).range(1..=64)).changed(), ) { - overrides.max_components_per_fit_window = Some(value); + overrides.maximum_model_order = Some(value); } - let mut value = invocation.params.linewidth_hz; + let mut value = invocation.params.component_linewidth_bounds_hz; if setting_row( ui, - "Linewidth range (Hz)", - invocation.sources.linewidth_hz, + "Component linewidth range (Hz)", + invocation.sources.component_linewidth_bounds_hz, &nmr, - &mut overrides.linewidth_hz, + &mut overrides.component_linewidth_bounds_hz, |ui| { let first = ui .add(DragValue::new(&mut value.0).range(0.001..=1_000.0)) @@ -226,52 +226,7 @@ fn settings(app: &mut PlotxApp, index: usize, invocation: &CraftInvocation, ui: .changed() }, ) { - overrides.linewidth_hz = Some(value); - } - - let mut value = invocation.params.max_fit_window_width_hz; - if setting_row( - ui, - "Fit window width (Hz)", - invocation.sources.max_fit_window_width_hz, - &nmr, - &mut overrides.max_fit_window_width_hz, - |ui| { - ui.add(DragValue::new(&mut value).range(10.0..=10_000.0)) - .changed() - }, - ) { - overrides.max_fit_window_width_hz = Some(value); - } - - let mut value = invocation.params.filter_taps; - if setting_row( - ui, - "FIR taps", - invocation.sources.filter_taps, - &nmr, - &mut overrides.filter_taps, - |ui| { - ui.add(DragValue::new(&mut value).range(3..=4_095)) - .changed() - }, - ) { - overrides.filter_taps = Some(value | 1); - } - - let mut value = invocation.params.max_downsampled_points; - if setting_row( - ui, - "Max downsampled points", - invocation.sources.max_downsampled_points, - &nmr, - &mut overrides.max_downsampled_points, - |ui| { - ui.add(DragValue::new(&mut value).range(64..=65_536)) - .changed() - }, - ) { - overrides.max_downsampled_points = Some(value); + overrides.component_linewidth_bounds_hz = Some(value); } if invocation.params.profile == CraftProfile::Ssfp { @@ -311,9 +266,11 @@ fn settings(app: &mut PlotxApp, index: usize, invocation: &CraftInvocation, ui: } } ui.weak(format!( - "Derived plan: skip {} points · {} actual taps · {} reconstructed points", + "Fixed protocol: {:.0} Hz modeling bandwidth · {:.2} s modeling duration · skip {} points · {} FIR taps · {} reconstructed points", + invocation.params.profile.modeling_bandwidth_hz(), + invocation.params.profile.modeling_duration_s(), invocation.derived_plan.effective_skip_points, - invocation.derived_plan.actual_filter_taps, + invocation.derived_plan.effective_fir_filter_taps, invocation.derived_plan.reconstruction_points, )); }); diff --git a/crates/core/src/project/craft_tests.rs b/crates/core/src/project/craft_tests.rs index df8a39a..57bdbad 100644 --- a/crates/core/src/project/craft_tests.rs +++ b/crates/core/src/project/craft_tests.rs @@ -1,10 +1,15 @@ use super::tests::synthetic_1d; use super::*; -use crate::state::{CraftRunId, FieldPayload, StoredCraftRun}; +use crate::state::{ + CraftRunId, FieldPayload, NewAnalysisReport, ReportKindId, ReportSource, ReportStatus, + StoredCraftRun, +}; use plotx_processing::craft::{ - CraftComponent, CraftComponentId, CraftDiagnostics, CraftFitWindowDiagnostic, + CraftComponent, CraftComponentId, CraftDiagnostics, CraftModelingWindowDiagnostic, CraftParamOverrides, CraftParams, CraftReference, CraftRegionId, CraftRegionRatio, - CraftRegionSummary, CraftResult, CraftRunStatus, resolve_craft_invocation, + CraftRegionSummary, CraftReportDefinition, CraftResult, CraftRunStatus, + CraftStabilityDiagnostics, CraftStabilityMetric, CraftStabilityRegion, + resolve_craft_invocation, }; fn sample_run(data: &NmrData) -> StoredCraftRun { @@ -80,31 +85,54 @@ fn sample_run(data: &NmrData) -> StoredCraftRun { residual_rss: 1.0, normalized_residual: 0.02, maximum_condition_number: Some(4.0), - fit_windows: vec![ - CraftFitWindowDiagnostic { - region: CraftRegionId(0), - core_hz: (-1300.0, -1100.0), - padded_hz: (-1320.0, -1080.0), - actual_decimation: 4, - retained_samples: 256, + modeling_windows: vec![ + CraftModelingWindowDiagnostic { + retention_band_hz: (-1300.0, -1100.0), + modeling_band_hz: (-1320.0, -1080.0), + decimation_factor: 4, + modeled_sample_count: 256, evaluated_model_orders: 7, selected_model_order: 1, - bic: Some(-25.0), + training_bic: Some(-25.0), condition_number: Some(3.0), + modeled_duration_s: 1.0, + training_normalized_residual: 0.01, + validation_normalized_residual: 0.02, }, - CraftFitWindowDiagnostic { - region: CraftRegionId(1), - core_hz: (1100.0, 1300.0), - padded_hz: (1080.0, 1320.0), - actual_decimation: 4, - retained_samples: 256, + CraftModelingWindowDiagnostic { + retention_band_hz: (1100.0, 1300.0), + modeling_band_hz: (1080.0, 1320.0), + decimation_factor: 4, + modeled_sample_count: 256, evaluated_model_orders: 7, selected_model_order: 1, - bic: Some(-20.0), + training_bic: Some(-20.0), condition_number: Some(4.0), + modeled_duration_s: 1.0, + training_normalized_residual: 0.01, + validation_normalized_residual: 0.02, }, ], warnings: Vec::new(), + stability: CraftStabilityDiagnostics { + delta_ppm: 0.016, + regions: vec![CraftStabilityRegion { + region: CraftRegionId(0), + metric: CraftStabilityMetric { + median: 1.0, + minimum: 0.999, + maximum: 1.001, + relative_dispersion: 0.002, + }, + component_count_min: 1, + component_count_max: 1, + model_order_min: 1, + model_order_max: 1, + }], + ratio: None, + passed: true, + skipped: Vec::new(), + }, }, synthetic_fid: Vec::new(), residual_fid: Vec::new(), @@ -181,7 +209,7 @@ fn unavailable_craft_diagnostics_survive_project_roundtrip() { run.components[0].linewidth_std_hz = None; run.components[0].phase_std_rad = None; run.diagnostics.maximum_condition_number = None; - run.diagnostics.fit_windows[0].bic = None; + run.diagnostics.modeling_windows[0].training_bic = None; let mut dataset = NmrDataset::load(data); dataset.craft_runs.push(run.clone()); dataset.reconcile_craft_fields(); @@ -200,6 +228,101 @@ fn unavailable_craft_diagnostics_survive_project_roundtrip() { assert_eq!(dataset.craft_runs, vec![run]); } +#[test] +fn only_stable_complete_runs_create_quantitative_reports() { + let data = synthetic_1d(); + let definition = CraftReportDefinition::default(); + let stable = sample_run(&data); + assert!(stable.amplitude_report(definition.clone()).is_ok()); + + let mut needs_review = stable.clone(); + needs_review.diagnostics.status = CraftRunStatus::Partial; + needs_review.diagnostics.stability.passed = false; + assert!( + needs_review + .amplitude_report(definition) + .unwrap_err() + .contains("NeedsReview") + ); + assert!(crate::state::craft_component_table(&needs_review).is_ok()); +} + +#[test] +fn report_status_tracks_stability_and_source_availability() { + let data = synthetic_1d(); + let mut dataset = NmrDataset::load(data.clone()); + let mut run = sample_run(&data); + run.provenance.invocation.reference = dataset.craft_reference(); + let definition = CraftReportDefinition::default(); + let snapshot = run.amplitude_report(definition.clone()).unwrap(); + let source = ReportSource { + dataset: dataset.resource_id, + craft_run: run.id, + }; + dataset.craft_runs.push(run); + dataset.reconcile_craft_fields(); + let mut app = crate::state::PlotxApp::new(); + app.doc.datasets.push(Dataset::Nmr(Box::new(dataset))); + let report = app.doc.create_report(NewAnalysisReport { + name: "CRAFT stability status".to_owned(), + kind: ReportKindId::new("craft_amplitude"), + source, + definition: serde_json::to_value(definition).unwrap(), + snapshot: serde_json::to_value(snapshot).unwrap(), + source_fingerprint: crate::state::craft_input_sha256(&data), + schema_version: 1, + }); + + assert_eq!( + app.doc.report(report).unwrap().status(&app.doc), + ReportStatus::Available + ); + app.doc.datasets[0].as_nmr_mut().unwrap().craft_runs[0] + .diagnostics + .stability + .passed = false; + assert_eq!( + app.doc.report(report).unwrap().status(&app.doc), + ReportStatus::NeedsReview + ); + app.doc.datasets[0].as_nmr_mut().unwrap().craft_runs.clear(); + assert_eq!( + app.doc.report(report).unwrap().status(&app.doc), + ReportStatus::Unavailable + ); +} + +#[test] +fn stability_snapshot_survives_project_roundtrip() { + let data = synthetic_1d(); + let mut run = sample_run(&data); + run.diagnostics.stability.skipped = vec!["contract: overlapping regions".to_owned()]; + run.diagnostics.stability.ratio = Some(CraftStabilityMetric { + median: 0.5, + minimum: 0.498, + maximum: 0.502, + relative_dispersion: 0.008, + }); + let mut dataset = NmrDataset::load(data); + dataset.craft_runs.push(run.clone()); + dataset.reconcile_craft_fields(); + let mut app = crate::state::PlotxApp::new(); + app.doc.datasets.push(Dataset::Nmr(Box::new(dataset))); + let path = super::tests::temp_project("craft-stability"); + let _ = std::fs::remove_file(&path); + + save_project(&app, &path, false).unwrap(); + let loaded = load_project(&path).unwrap(); + let _ = std::fs::remove_file(&path); + + assert_eq!( + loaded.doc.datasets[0].as_nmr().unwrap().craft_runs[0] + .diagnostics + .stability, + run.diagnostics.stability + ); +} + #[test] fn craft_component_table_link_and_board_visibility_survive_roundtrip() { let data = synthetic_1d(); diff --git a/crates/core/src/project/mod.rs b/crates/core/src/project/mod.rs index cbd8e80..ffe0b3e 100644 --- a/crates/core/src/project/mod.rs +++ b/crates/core/src/project/mod.rs @@ -403,6 +403,16 @@ fn save_project_impl( }); } + for report in &doc.reports { + let path = format!("reports/{}.json", report.id.0); + write_json(&mut zip, options, &path, report)?; + manifest.runs.push(Entry { + id: report.id.0.to_string(), + role: "report".to_owned(), + path, + }); + } + let workspace = Workspace { dataset_order: bindings, view_order, @@ -469,6 +479,7 @@ pub fn load_project(path: &Path) -> Result { app.doc.datasets.clear(); app.doc.canvases.clear(); app.doc.assets.clear(); + app.doc.reports.clear(); app.session.project_load_warnings = asset_codec::load_assets(&mut zip, &manifest, &mut app)?; app.doc.project_path = Some(path.to_owned()); // Restore before the canvases below are built: figures stamp the document @@ -532,8 +543,16 @@ pub fn load_project(path: &Path) -> Result { app.doc.automation_runs = manifest .runs .iter() + .filter(|entry| entry.role == "run") .map(|entry| read_json(&mut zip, &entry.path)) .collect::>>()?; + app.doc.reports = manifest + .runs + .iter() + .filter(|entry| entry.role == "report") + .map(|entry| read_json(&mut zip, &entry.path)) + .collect::>>()?; + app.doc.repair_report_allocator(); app.doc.automation_revision = workspace.automation_revision; asset_codec::append_undeclared_image_warnings( &app.doc, diff --git a/crates/core/src/state/app_impl.rs b/crates/core/src/state/app_impl.rs index 070f6a9..2b828f8 100644 --- a/crates/core/src/state/app_impl.rs +++ b/crates/core/src/state/app_impl.rs @@ -45,6 +45,8 @@ impl PlotxApp { project_revision: None, automation_revision: 0, automation_runs: Vec::new(), + reports: Vec::new(), + next_report_id: 0, edit_generation: 0, dirty: false, save_include_view_snapshots: settings.export.include_view_snapshots, diff --git a/crates/core/src/state/app_impl_compute.rs b/crates/core/src/state/app_impl_compute.rs index 01ea236..f986d20 100644 --- a/crates/core/src/state/app_impl_compute.rs +++ b/crates/core/src/state/app_impl_compute.rs @@ -279,6 +279,49 @@ impl PlotxApp { Ok(table_index) } + pub fn materialize_craft_report_table( + &mut self, + dataset: usize, + report_id: crate::state::ReportId, + ) -> Result { + let record = self + .doc + .report(report_id) + .cloned() + .ok_or_else(|| "The CRAFT report is no longer available.".to_owned())?; + match record.status(&self.doc) { + crate::state::ReportStatus::Available => {} + crate::state::ReportStatus::NeedsReview => { + return Err("The CRAFT report needs review and cannot be exported as reliable quantitative data.".to_owned()); + } + crate::state::ReportStatus::Unavailable => { + return Err("The CRAFT report source is unavailable.".to_owned()); + } + } + let report: plotx_processing::craft::CraftAmplitudeReport = + serde_json::from_value(record.snapshot) + .map_err(|error| format!("Could not decode CRAFT report snapshot: {error}"))?; + let mut table = craft_amplitude_report_table(&report)?; + let source_id = self + .doc + .datasets + .get(dataset) + .map(Dataset::resource_id) + .ok_or_else(|| "The CRAFT source dataset is no longer available.".to_owned())?; + table.lineage = Some(DatasetLineage::new( + DerivationKind::CraftComponentTable, + [source_id], + )); + table.name = Some(format!("CRAFT report {}", report_id.0 + 1)); + let sheet = table.board_rect_pt(); + table.board_pos = crate::state::next_board_frame_pos(self, [sheet.width, sheet.height]); + table.board_sheet_visible = false; + let table_index = self.doc.datasets.len(); + self.doc.datasets.push(Dataset::Table(Box::new(table))); + self.mark_document_dirty(); + Ok(table_index) + } + pub fn show_craft_component_table_on_board(&mut self, table: usize) -> Result<(), String> { let table = self .doc diff --git a/crates/core/src/state/app_impl_compute_tests.rs b/crates/core/src/state/app_impl_compute_tests.rs index be1798b..2bf74bd 100644 --- a/crates/core/src/state/app_impl_compute_tests.rs +++ b/crates/core/src/state/app_impl_compute_tests.rs @@ -74,10 +74,9 @@ fn craft_result_is_installed_with_provenance_by_dataset_identity() { let target = app.doc.datasets[0].resource_id(); app.session.ui.craft_task_dataset = Some(target); let mut params = plotx_processing::craft::CraftParams::conventional(); - params.filter_taps = 31; - params.max_fit_window_width_hz = 2_000.0; - params.max_downsampled_points = 512; - params.max_components_per_fit_window = 2; + params.fir_filter_taps = 31; + params.maximum_modeled_sample_count = 512; + params.maximum_model_order = 2; assert!(app.request_craft_analysis( 0, @@ -173,10 +172,9 @@ fn craft_rerun_keeps_requested_parent_without_hijacking_another_task() { let target = app.doc.datasets[0].resource_id(); app.session.ui.craft_task_dataset = Some(target); let mut params = plotx_processing::craft::CraftParams::conventional(); - params.filter_taps = 31; - params.max_fit_window_width_hz = 2_000.0; - params.max_downsampled_points = 512; - params.max_components_per_fit_window = 2; + params.fir_filter_taps = 31; + params.maximum_modeled_sample_count = 512; + params.maximum_model_order = 2; assert!(app.request_craft_analysis( 0, plotx_processing::craft::CraftParamOverrides::from_params(params), diff --git a/crates/core/src/state/craft.rs b/crates/core/src/state/craft.rs index b9915a0..a2da678 100644 --- a/crates/core/src/state/craft.rs +++ b/crates/core/src/state/craft.rs @@ -1,8 +1,9 @@ use super::{FloatSeries, NmrDataset, TableDataset, materialized_float_series_table}; use plotx_io::NmrData; use plotx_processing::craft::{ - CRAFT_ALGORITHM, CRAFT_ALGORITHM_VERSION, CraftComponent, CraftDiagnostics, CraftInvocation, - CraftReference, CraftRegionRatio, CraftRegionSummary, CraftResult, + CRAFT_ALGORITHM, CRAFT_ALGORITHM_VERSION, CraftAmplitudeReport, CraftComponent, + CraftDiagnostics, CraftInvocation, CraftReference, CraftRegionRatio, CraftRegionSummary, + CraftReportDefinition, CraftResult, calculate_craft_report, }; use serde::{Deserialize, Serialize}; use sha2::{Digest, Sha256}; @@ -32,6 +33,17 @@ pub struct StoredCraftRun { } impl StoredCraftRun { + pub fn amplitude_report( + &self, + definition: CraftReportDefinition, + ) -> Result { + if self.diagnostics.status != plotx_processing::craft::CraftRunStatus::Complete + || !self.diagnostics.stability.passed + { + return Err("CRAFT run is marked NeedsReview; quantitative amplitude reports are unavailable until stability checks pass.".to_owned()); + } + calculate_craft_report(&self.components, definition).map_err(|e| e.to_string()) + } pub fn from_result( id: CraftRunId, data: &NmrData, @@ -241,3 +253,81 @@ pub fn craft_component_table(run: &StoredCraftRun) -> Result Result { + let rows = report.segments.len(); + let values = |read: fn(&plotx_processing::craft::CraftReportSegment) -> f64| { + report + .segments + .iter() + .map(|segment| Some(read(segment))) + .collect() + }; + materialized_float_series_table( + ( + "segment".into(), + "".into(), + (1..=rows).map(|i| Some(i as f64)).collect(), + ), + vec![ + FloatSeries { + name: "report threshold A/N".into(), + unit: "".into(), + values: vec![Some(report.definition.threshold_an); rows], + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "segment width".into(), + unit: "Hz".into(), + values: vec![Some(report.definition.segment_width_hz); rows], + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "center".into(), + unit: "Hz".into(), + values: values(|s| s.center_hz), + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "start".into(), + unit: "Hz".into(), + values: values(|s| s.start_hz), + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "end".into(), + unit: "Hz".into(), + values: values(|s| s.end_hz), + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "component count".into(), + unit: "".into(), + values: values(|s| s.component_count as f64), + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "scalar amplitude sum t0".into(), + unit: "".into(), + values: values(|s| s.scalar_amplitude_sum_t0), + uncertainty: None, + fit: None, + }, + FloatSeries { + name: "coherent amplitude t0".into(), + unit: "".into(), + values: values(|s| s.coherent_amplitude_t0), + uncertainty: None, + fit: None, + }, + ], + "plotx.analysis.craft-amplitude-report.v1", + ) + .map_err(|error| error.to_string()) +} diff --git a/crates/core/src/state/mod.rs b/crates/core/src/state/mod.rs index 3639194..3c80360 100644 --- a/crates/core/src/state/mod.rs +++ b/crates/core/src/state/mod.rs @@ -96,6 +96,7 @@ mod plot_interaction; mod plot_object; mod pseudo_map_field; mod region; +mod reports; mod scientific_summary; mod selection; mod series_binding; @@ -188,6 +189,7 @@ pub use plot_interaction::*; pub use plot_object::*; pub(crate) use pseudo_map_field::{DOSY_GRID_COLS, DOSY_GRID_ROWS, dosy_scalar_grid}; pub use region::*; +pub use reports::*; pub use scientific_summary::*; pub use selection::*; pub use series_binding::*; diff --git a/crates/core/src/state/reports.rs b/crates/core/src/state/reports.rs new file mode 100644 index 0000000..449ad45 --- /dev/null +++ b/crates/core/src/state/reports.rs @@ -0,0 +1,191 @@ +use super::{CraftRunId, DatasetId, Document}; +use serde::{Deserialize, Serialize}; + +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)] +pub struct ReportId(pub u64); + +#[derive(Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)] +pub struct ReportKindId(pub String); + +impl ReportKindId { + pub fn new(value: impl Into) -> Self { + Self(value.into()) + } +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)] +pub struct ReportSource { + pub dataset: DatasetId, + pub craft_run: CraftRunId, +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ReportStatus { + Available, + NeedsReview, + Unavailable, +} + +impl AnalysisReportRecord { + pub fn status(&self, document: &Document) -> ReportStatus { + let Some(dataset) = document + .datasets + .iter() + .find(|d| d.resource_id() == self.source.dataset) + else { + return ReportStatus::Unavailable; + }; + let Some(nmr) = dataset.as_nmr() else { + return ReportStatus::Unavailable; + }; + let Some(run) = nmr.craft_run(self.source.craft_run) else { + return ReportStatus::Unavailable; + }; + if self.kind.0 == "craft_amplitude" + && (run.diagnostics.status != plotx_processing::craft::CraftRunStatus::Complete + || !run.diagnostics.stability.passed) + { + return ReportStatus::NeedsReview; + } + if self.kind.0 == "craft_amplitude" { + let valid_definition = serde_json::from_value::< + plotx_processing::craft::CraftReportDefinition, + >(self.definition.clone()) + .ok() + .is_some_and(|definition| definition.validate().is_ok()); + if !valid_definition { + return ReportStatus::NeedsReview; + } + } + if self.kind.0 == "craft_amplitude" + && serde_json::from_value::( + self.snapshot.clone(), + ) + .ok() + .is_none_or(|snapshot| snapshot.validate_against(&run.components).is_err()) + { + return ReportStatus::NeedsReview; + } + if run.provenance.input_sha256 != self.source_fingerprint + || run.is_stale_for(&nmr.data, nmr.craft_reference()) + { + ReportStatus::NeedsReview + } else { + ReportStatus::Available + } + } +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct AnalysisReportRecord { + pub id: ReportId, + pub name: String, + pub kind: ReportKindId, + pub source: ReportSource, + /// Tagged, domain-owned definition and resolved result snapshot. + pub definition: serde_json::Value, + pub snapshot: serde_json::Value, + pub source_fingerprint: String, + pub schema_version: u32, +} + +#[derive(Clone, Debug, PartialEq)] +pub struct NewAnalysisReport { + pub name: String, + pub kind: ReportKindId, + pub source: ReportSource, + pub definition: serde_json::Value, + pub snapshot: serde_json::Value, + pub source_fingerprint: String, + pub schema_version: u32, +} + +impl Document { + pub fn allocate_report_id(&mut self) -> ReportId { + let id = ReportId(self.next_report_id); + self.next_report_id = self + .next_report_id + .checked_add(1) + .expect("report id overflow"); + id + } + + pub fn repair_report_allocator(&mut self) { + let required = self + .reports + .iter() + .map(|r| r.id.0.saturating_add(1)) + .max() + .unwrap_or(0); + self.next_report_id = self.next_report_id.max(required); + } + + pub fn report(&self, id: ReportId) -> Option<&AnalysisReportRecord> { + self.reports.iter().find(|r| r.id == id) + } + pub fn reports_for_source( + &self, + source: ReportSource, + ) -> impl Iterator { + self.reports.iter().filter(move |r| r.source == source) + } + pub fn create_report(&mut self, report: NewAnalysisReport) -> ReportId { + let id = self.allocate_report_id(); + self.reports.push(AnalysisReportRecord { + id, + name: report.name, + kind: report.kind, + source: report.source, + definition: report.definition, + snapshot: report.snapshot, + source_fingerprint: report.source_fingerprint, + schema_version: report.schema_version, + }); + self.mark_dirty(); + id + } + pub fn rename_report(&mut self, id: ReportId, name: String) -> Result<(), String> { + let report = self + .report_mut(id) + .ok_or_else(|| "Report not found".to_owned())?; + report.name = name; + self.mark_dirty(); + Ok(()) + } + pub fn copy_report(&mut self, id: ReportId, name: Option) -> Result { + let source = self + .report(id) + .cloned() + .ok_or_else(|| "Report not found".to_owned())?; + let new_id = self.allocate_report_id(); + self.reports.push(AnalysisReportRecord { + id: new_id, + name: name.unwrap_or_else(|| format!("{} copy", source.name)), + ..source + }); + self.mark_dirty(); + Ok(new_id) + } + pub fn update_report(&mut self, record: AnalysisReportRecord) -> Result<(), String> { + let slot = self + .reports + .iter_mut() + .find(|r| r.id == record.id) + .ok_or_else(|| "Report not found".to_owned())?; + *slot = record; + self.mark_dirty(); + Ok(()) + } + pub fn delete_report(&mut self, id: ReportId) -> bool { + let before = self.reports.len(); + self.reports.retain(|r| r.id != id); + let changed = before != self.reports.len(); + if changed { + self.mark_dirty(); + } + changed + } + fn report_mut(&mut self, id: ReportId) -> Option<&mut AnalysisReportRecord> { + self.reports.iter_mut().find(|r| r.id == id) + } +} diff --git a/crates/core/src/state/ui_state.rs b/crates/core/src/state/ui_state.rs index 845f8f8..52b8481 100644 --- a/crates/core/src/state/ui_state.rs +++ b/crates/core/src/state/ui_state.rs @@ -362,6 +362,7 @@ pub struct UiState { pub craft_selected_run: Option, pub craft_task_page: CraftTaskPage, pub craft_result_tab: CraftResultTab, + pub craft_selected_report: Option, pub craft_component_sort: CraftComponentSort, pub craft_component_region: Option, pub craft_selected_component: Option, @@ -559,6 +560,7 @@ impl Default for UiState { craft_selected_run: None, craft_task_page: CraftTaskPage::Setup, craft_result_tab: CraftResultTab::Overview, + craft_selected_report: None, craft_component_sort: CraftComponentSort::ChemicalShift, craft_component_region: None, craft_selected_component: None, @@ -619,6 +621,8 @@ pub struct Document { pub project_revision: Option, pub automation_revision: u64, pub automation_runs: Vec, + pub reports: Vec, + pub next_report_id: u64, /// Incremented for every persisted edit. Background save completion uses /// this token instead of clearing `dirty` unconditionally. pub edit_generation: u64, diff --git a/crates/core/src/state/ui_state_craft.rs b/crates/core/src/state/ui_state_craft.rs index 6ed5a56..7676934 100644 --- a/crates/core/src/state/ui_state_craft.rs +++ b/crates/core/src/state/ui_state_craft.rs @@ -11,6 +11,7 @@ pub enum CraftResultTab { Overview, Components, Diagnostics, + Reports, } #[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] diff --git a/crates/processing/src/craft.rs b/crates/processing/src/craft.rs index fb0f326..ff48866 100644 --- a/crates/processing/src/craft.rs +++ b/crates/processing/src/craft.rs @@ -1,36 +1,42 @@ //! Complete Reduction to Amplitude Frequency Table (CRAFT) for one-dimensional FIDs. use num_complex::Complex64; -use plotx_analysis::craft::{ - CraftFitBounds, CraftFitError, DampedSinusoid, evaluate_damped_sinusoids_cancellable, - matrix_pencil_estimates, -}; +use plotx_analysis::craft::CraftFitError; use plotx_io::{Domain, NmrData}; use serde::{Deserialize, Serialize}; -use std::f64::consts::{PI, TAU}; mod diagnostics; +mod fitting; mod preflight; mod reconstruction; mod regions; +mod report; mod resolution; +mod stability; pub use diagnostics::{ - CraftDiagnostics, CraftFitWindowDiagnostic, CraftRegionRatio, CraftRegionSummary, - CraftRunStatus, CraftWarning, CraftWarningKind, + CraftDiagnostics, CraftModelingWindowDiagnostic, CraftRegionRatio, CraftRegionSummary, + CraftRunStatus, CraftStabilityDiagnostics, CraftStabilityMetric, CraftStabilityRegion, + CraftWarning, CraftWarningKind, }; +use fitting::{CraftModelingContext, fit_modeling_window}; pub use preflight::{ CraftAssessmentIssue, CraftInputAssessment, CraftIssueAction, CraftIssueCode, CraftIssueSeverity, CraftSignalSuggestion, }; use reconstruction::model_at; pub use reconstruction::{synthesize_craft_fid, synthesize_craft_samples}; -use regions::{HzRegion, build_regions, region_ratio, selections_are_valid, summarize_regions}; +use regions::{build_modeling_windows, region_ratio, selections_are_valid, summarize_regions}; +pub use report::{ + CraftAmplitudeReport, CraftReportDefinition, CraftReportError, CraftReportSegment, + calculate_craft_report, +}; pub use resolution::{ - CraftDerivedPlan, CraftDerivedWindow, CraftParamOverrides, CraftParamSource, - CraftParameterSources, resolve_craft_invocation, + CraftDerivedModelingWindow, CraftDerivedPlan, CraftModelingPolicy, CraftParamOverrides, + CraftParamSource, CraftParameterSources, resolve_craft_invocation, }; +use stability::{components_for_regions, stability_diagnostics}; -pub const CRAFT_ALGORITHM: &str = "plotx-craft-matrix-pencil-bic"; +pub const CRAFT_ALGORITHM: &str = "plotx-craft-matrix-pencil-validation"; pub const CRAFT_ALGORITHM_VERSION: u32 = 1; #[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)] @@ -41,6 +47,22 @@ pub enum CraftProfile { Ssfp, } +impl CraftProfile { + pub const fn modeling_bandwidth_hz(self) -> f64 { + match self { + Self::Conventional => 250.0, + Self::Ssfp => 2_000.0, + } + } + + pub const fn modeling_duration_s(self) -> f64 { + match self { + Self::Conventional => 1.0, + Self::Ssfp => 1.2, + } + } +} + #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] pub struct CraftRegion { pub id: CraftRegionId, @@ -114,13 +136,11 @@ pub struct CraftParams { pub profile: CraftProfile, /// Empty means the complete acquired spectral width. pub regions: Vec, - pub max_components_per_fit_window: usize, - pub min_amplitude_to_noise: f64, - pub linewidth_hz: (f64, f64), - pub filter_taps: usize, - pub padding_fraction: f64, - pub max_fit_window_width_hz: f64, - pub max_downsampled_points: usize, + pub maximum_model_order: usize, + pub minimum_amplitude_to_noise: f64, + pub component_linewidth_bounds_hz: (f64, f64), + pub fir_filter_taps: usize, + pub maximum_modeled_sample_count: usize, pub skip_duration_s: f64, pub reconstruction_duration_s: Option, } @@ -130,13 +150,11 @@ impl CraftParams { Self { profile: CraftProfile::Conventional, regions: Vec::new(), - max_components_per_fit_window: 15, - min_amplitude_to_noise: 3.3, - linewidth_hz: (0.05, 10.0), - filter_taps: 499, - padding_fraction: 0.2, - max_fit_window_width_hz: 500.0, - max_downsampled_points: 8192, + maximum_model_order: 15, + minimum_amplitude_to_noise: 3.3, + component_linewidth_bounds_hz: (0.05, 20.0), + fir_filter_taps: 499, + maximum_modeled_sample_count: 8192, skip_duration_s: 0.0, reconstruction_duration_s: None, } @@ -147,34 +165,29 @@ impl CraftParams { profile: CraftProfile::Ssfp, skip_duration_s: 0.0005, reconstruction_duration_s: Some(1.2), - max_fit_window_width_hz: 2_000.0, ..Self::conventional() } } pub fn discovery() -> Self { Self { - min_amplitude_to_noise: 2.5, + minimum_amplitude_to_noise: 2.5, ..Self::conventional() } } pub fn validate(&self) -> Result<(), CraftError> { - if self.max_components_per_fit_window == 0 - || self.max_components_per_fit_window > 64 - || !self.min_amplitude_to_noise.is_finite() - || self.min_amplitude_to_noise <= 0.0 - || !self.linewidth_hz.0.is_finite() - || !self.linewidth_hz.1.is_finite() - || self.linewidth_hz.0 <= 0.0 - || self.linewidth_hz.0 >= self.linewidth_hz.1 - || self.filter_taps < 3 - || self.filter_taps.is_multiple_of(2) - || !self.padding_fraction.is_finite() - || !(0.0..=1.0).contains(&self.padding_fraction) - || !self.max_fit_window_width_hz.is_finite() - || self.max_fit_window_width_hz <= 0.0 - || self.max_downsampled_points < 64 + if self.maximum_model_order == 0 + || self.maximum_model_order > 64 + || !self.minimum_amplitude_to_noise.is_finite() + || self.minimum_amplitude_to_noise <= 0.0 + || !self.component_linewidth_bounds_hz.0.is_finite() + || !self.component_linewidth_bounds_hz.1.is_finite() + || self.component_linewidth_bounds_hz.0 <= 0.0 + || self.component_linewidth_bounds_hz.0 >= self.component_linewidth_bounds_hz.1 + || self.fir_filter_taps < 3 + || self.fir_filter_taps.is_multiple_of(2) + || self.maximum_modeled_sample_count < 64 || !self.skip_duration_s.is_finite() || self.skip_duration_s < 0.0 || self @@ -206,6 +219,7 @@ pub struct CraftInvocation { pub reference: CraftReference, pub derived_plan: CraftDerivedPlan, pub assessment: CraftInputAssessment, + pub modeling_policy: CraftModelingPolicy, } impl CraftInvocation { @@ -217,6 +231,9 @@ impl CraftInvocation { pub fn validate(&self, data: &NmrData) -> Result<(), CraftError> { self.params.validate()?; self.reference.validate(data)?; + if self.modeling_policy != CraftModelingPolicy::for_params(&self.params) { + return Err(CraftError::InvalidParameters); + } if self.assessment.can_run() { Ok(()) } else { @@ -272,21 +289,10 @@ pub enum CraftError { InvalidReference, #[error("CRAFT analysis was cancelled")] Cancelled, - #[error("CRAFT could not fit a requested region: {0}")] + #[error("CRAFT could not fit a modeling window: {0}")] Fit(#[from] CraftFitError), } -struct RegionResult { - components: Vec, - center_hz: f64, - bic: Option, - condition_number: f64, - decimation: usize, - retained_samples: usize, - evaluated_model_orders: usize, - warning: Option<(CraftWarningKind, String)>, -} - pub fn process_craft_cancellable( data: &NmrData, invocation: &CraftInvocation, @@ -318,59 +324,94 @@ pub fn process_craft_cancellable( return Err(CraftError::InvalidInput); } let input = &data.points[skip..]; + let modeling_context = CraftModelingContext { + input, + skipped_points: skip, + group_delay_points: data.group_delay, + spectral_width_hz: sw, + params, + policy: invocation.modeling_policy, + }; let noise_sigma = estimate_complex_noise(input).max(f64::MIN_POSITIVE); - let regions = build_regions(data, params, reference)?; + let modeling_windows = build_modeling_windows( + data, + params, + reference, + &invocation.assessment.clear_signals, + )?; let mut fitted = Vec::new(); let mut warnings = Vec::new(); - let mut fit_windows = Vec::with_capacity(regions.len()); + let mut window_diagnostics = Vec::with_capacity(modeling_windows.len()); let mut max_condition = 1.0_f64; + let selected_frequency_bands = params + .regions + .iter() + .map(|region| { + let region = region.normalized(); + ( + (region.start_ppm - reference.effective_carrier_ppm()) * data.observe_freq_mhz, + (region.end_ppm - reference.effective_carrier_ppm()) * data.observe_freq_mhz, + ) + }) + .collect::>(); - for (index, region) in regions.iter().copied().enumerate() { + for (index, window) in modeling_windows.iter().copied().enumerate() { if cancelled() { return Err(CraftError::Cancelled); } - let result = fit_region(input, skip, data.group_delay, sw, region, params, cancelled)?; - if let Some((kind, message)) = result.warning { + let result = fit_modeling_window(&modeling_context, window, cancelled)?; + let contributes_to_selection = selected_frequency_bands.is_empty() + || selected_frequency_bands.iter().any(|&(start, end)| { + window.retention_band_hz.0 <= end && window.retention_band_hz.1 >= start + }); + if let Some((kind, message)) = result.warning + && contributes_to_selection + { warnings.push(CraftWarning { kind, - region: Some(region.selection.id), - fit_window: Some(index), - message: format!("Fit window {}: {message}", index + 1), + region: None, + modeling_window: Some(index), + message: format!("Modeling window {}: {message}", index + 1), }); } - fit_windows.push(CraftFitWindowDiagnostic { - region: region.selection.id, - core_hz: region.core, - padded_hz: region.padded, - actual_decimation: result.decimation, - retained_samples: result.retained_samples, + window_diagnostics.push(CraftModelingWindowDiagnostic { + retention_band_hz: window.retention_band_hz, + modeling_band_hz: window.modeling_band_hz, + decimation_factor: result.decimation, + modeled_sample_count: result.modeled_sample_count, evaluated_model_orders: result.evaluated_model_orders, selected_model_order: result.components.len(), - bic: result.bic, + training_bic: result.training_bic, condition_number: result .condition_number .is_finite() .then_some(result.condition_number), + modeled_duration_s: result.modeled_duration_s, + training_normalized_residual: result.training_normalized_residual, + validation_normalized_residual: result.validation_normalized_residual, }); max_condition = max_condition.max(result.condition_number); for component in result.components { let frequency_hz = component.frequency_hz + result.center_hz; - if frequency_hz >= region.core.0 - && frequency_hz < region.core.1 - && component.amplitude / noise_sigma >= params.min_amplitude_to_noise + let is_last_window = index + 1 == modeling_windows.len(); + if frequency_hz >= window.retention_band_hz.0 + && (frequency_hz < window.retention_band_hz.1 + || (is_last_window && frequency_hz <= window.retention_band_hz.1)) { - fitted.push((region.selection.id, frequency_hz, component)); + // Padded modeling bands may overlap. Retention bands assign a + // model to exactly one window before the sub-tables are joined. + fitted.push((frequency_hz, component)); } } } - fitted.sort_by(|left, right| left.1.total_cmp(&right.1)); - let components: Vec = fitted + fitted.sort_by(|left, right| left.0.total_cmp(&right.0)); + let all_components: Vec = fitted .into_iter() .enumerate() - .map(|(id, (region, frequency_hz, component))| CraftComponent { + .map(|(id, (frequency_hz, component))| CraftComponent { id: CraftComponentId(id as u64), - region, + region: CraftRegionId(0), frequency_hz, chemical_shift_ppm: reference.effective_carrier_ppm() + frequency_hz / data.observe_freq_mhz, @@ -378,18 +419,32 @@ pub fn process_craft_cancellable( phase_rad: component.phase_rad, decay_rate_s_inv: component.decay_rate_s_inv, linewidth_hz: component.linewidth_hz, - amplitude_to_noise: component.amplitude / noise_sigma, + amplitude_to_noise: component + .amplitude_std + .filter(|value| *value > 0.0) + .map_or(0.0, |value| component.amplitude / value), amplitude_std: component.amplitude_std, frequency_std_hz: component.frequency_std_hz, linewidth_std_hz: component.linewidth_std_hz, phase_std_rad: component.phase_std_rad, }) .collect(); - if params.min_amplitude_to_noise < 3.3 { + let selections = if params.regions.is_empty() { + let half_width_ppm = sw / (2.0 * data.observe_freq_mhz); + vec![CraftRegion::new( + CraftRegionId(0), + reference.effective_carrier_ppm() - half_width_ppm, + reference.effective_carrier_ppm() + half_width_ppm, + )] + } else { + params.regions.clone() + }; + let components = components_for_regions(&all_components, &selections); + if params.minimum_amplitude_to_noise < 3.3 { warnings.push(CraftWarning { kind: CraftWarningKind::LowAmplitudeThreshold, region: None, - fit_window: None, + modeling_window: None, message: "Discovery threshold is below the strict 3.3 amplitude/noise threshold; confirm weak components independently." .to_owned(), }); @@ -398,7 +453,7 @@ pub fn process_craft_cancellable( warnings.push(CraftWarning { kind: CraftWarningKind::SsfpQuantitation, region: None, - fit_window: None, + modeling_window: None, message: "SSFP response is relaxation-dependent; use this result for screening or relative comparison, not absolute qNMR." .to_owned(), }); @@ -412,7 +467,7 @@ pub fn process_craft_cancellable( .map(|issue| CraftWarning { kind: CraftWarningKind::InputAssessment, region: issue.region, - fit_window: None, + modeling_window: None, message: issue.message.clone(), }), ); @@ -420,18 +475,18 @@ pub fn process_craft_cancellable( warnings.push(CraftWarning { kind: CraftWarningKind::IllConditionedFit, region: None, - fit_window: None, + modeling_window: None, message: "One or more fits are ill-conditioned; inspect overlapping components and uncertainties.".to_owned(), }); } if components.iter().any(|component| { - (component.linewidth_hz - params.linewidth_hz.0).abs() < 1e-6 - || (component.linewidth_hz - params.linewidth_hz.1).abs() < 1e-6 + (component.linewidth_hz - params.component_linewidth_bounds_hz.0).abs() < 1e-6 + || (component.linewidth_hz - params.component_linewidth_bounds_hz.1).abs() < 1e-6 }) { warnings.push(CraftWarning { kind: CraftWarningKind::LinewidthAtBound, region: None, - fit_window: None, + modeling_window: None, message: "One or more linewidths reached a configured fit bound.".to_owned(), }); } @@ -444,7 +499,7 @@ pub fn process_craft_cancellable( warnings.push(CraftWarning { kind: CraftWarningKind::UnboundedUncertainty, region: None, - fit_window: None, + modeling_window: None, message: "One or more components have unbounded uncertainties.".to_owned(), }); } @@ -463,28 +518,13 @@ pub fn process_craft_cancellable( let residual_rss: f64 = residual_fid[skip..].iter().map(Complex64::norm_sqr).sum(); let input_rss: f64 = input.iter().map(Complex64::norm_sqr).sum(); let normalized_residual = (residual_rss / input_rss.max(f64::MIN_POSITIVE)).sqrt(); - let selections = if params.regions.is_empty() { - vec![regions[0].selection] - } else { - params - .regions - .iter() - .map(|requested| { - regions - .iter() - .find(|window| window.selection.id == requested.id) - .map(|window| window.selection) - .expect("validated CRAFT region has at least one fit window") - }) - .collect() - }; let region_summaries = summarize_regions(&components, &selections); for (position, summary) in region_summaries.iter().enumerate() { if summary.component_count == 0 { warnings.push(CraftWarning { kind: CraftWarningKind::EmptyRegion, region: Some(summary.region), - fit_window: None, + modeling_window: None, message: format!( "Region {} contains no retained signal components.", position + 1 @@ -492,6 +532,22 @@ pub fn process_craft_cancellable( }); } } + let stability = stability_diagnostics( + &all_components, + &selections, + &window_diagnostics, + invocation.modeling_policy, + reference, + data, + ); + if !stability.passed { + warnings.push(CraftWarning { + kind: CraftWarningKind::StabilityFailure, + region: None, + modeling_window: None, + message: "Boundary perturbation exceeded the 1% stability tolerance; retain the full fit for review, but do not use it for quantitative reporting.".to_owned(), + }); + } let status = if invocation.assessment.has_warnings() || warnings.iter().any(CraftWarning::blocks_quantitation) { @@ -510,244 +566,15 @@ pub fn process_craft_cancellable( residual_rss, normalized_residual, maximum_condition_number: max_condition.is_finite().then_some(max_condition), - fit_windows, + modeling_windows: window_diagnostics, warnings, + stability, }, synthetic_fid, residual_fid, }) } -fn fit_region( - input: &[Complex64], - skipped_points: usize, - group_delay_points: f64, - sw: f64, - region: HzRegion, - params: &CraftParams, - cancelled: &impl Fn() -> bool, -) -> Result { - let center_hz = (region.padded.0 + region.padded.1) * 0.5; - let padded_width = region.padded.1 - region.padded.0; - // Only the early signal-bearing record is modeled. Include one filter - // length of guard samples so the centered FIR has a fully observed window - // for every retained point. - let filter_input_len = input - .len() - .min(6_000_usize.saturating_add(params.filter_taps)); - let mixed: Vec = input[..filter_input_len] - .iter() - .enumerate() - .map(|(index, &value)| { - // Bruker points after the digital-filter transient are already - // samples of the FID starting at `ceil(delay) - delay`. Keeping the - // raw point number here would reintroduce the removed delay as an - // amplitude extrapolation and a first-order phase ramp. - let time = (skipped_points as f64 + index as f64 - group_delay_points) / sw; - value * Complex64::from_polar(1.0, -TAU * center_hz * time) - }) - .collect(); - let filtered = low_pass_fir( - &mixed, - sw, - padded_width * 0.5, - params.filter_taps, - cancelled, - )?; - // Retain two samples per padded-bandwidth interval so the FIR transition - // band remains below the downsampled Nyquist limit. - let mut decimation = (sw / (2.0 * padded_width).max(f64::MIN_POSITIVE)) - .floor() - .max(1.0) as usize; - decimation = decimation.max(filtered.len().div_ceil(params.max_downsampled_points)); - let filter_half = effective_filter_taps(params.filter_taps, mixed.len()) / 2; - let valid_end = filtered.len().saturating_sub(filter_half); - let phase_search_end = valid_end.min(filter_half.saturating_add(params.filter_taps)); - let phase_start = filtered[filter_half..phase_search_end] - .iter() - .enumerate() - .max_by(|left, right| left.1.norm_sqr().total_cmp(&right.1.norm_sqr())) - .map(|(index, _)| filter_half + index) - .unwrap_or(filter_half); - let useful_end = phase_start.saturating_add(6_000).min(valid_end); - let samples: Vec = filtered[phase_start..useful_end] - .iter() - .step_by(decimation) - .copied() - .collect(); - let times: Vec = (0..samples.len()) - .map(|index| { - (skipped_points as f64 + phase_start as f64 + (index * decimation) as f64 - - group_delay_points) - / sw - }) - .collect(); - if samples.len() < 16 { - return Ok(RegionResult { - components: Vec::new(), - center_hz, - bic: None, - condition_number: 1.0, - decimation, - retained_samples: samples.len(), - evaluated_model_orders: 0, - warning: Some(( - CraftWarningKind::FitWindowFailure, - "too few samples remained after filtering".to_owned(), - )), - }); - } - let relative_bounds = (region.padded.0 - center_hz, region.padded.1 - center_hz); - let n_observations = (samples.len() * 2) as f64; - let initial_rss = samples.iter().map(Complex64::norm_sqr).sum(); - let mut best_bic = bic(initial_rss, n_observations, 1.0); - let mut best_components = Vec::new(); - let mut best_condition = 1.0; - let mut warning = None; - - let fit_bounds = CraftFitBounds { - frequency_hz: relative_bounds, - linewidth_hz: params.linewidth_hz, - }; - let dwell_s = decimation as f64 / sw; - let merge_hz = sw / input.len() as f64; - let max_order = params - .max_components_per_fit_window - .min(samples.len() / 2 - 1); - let mut evaluated_model_orders = 0; - // Matrix-pencil cost grows cubically with its Hankel dimension. The early - // 256 uniformly sampled points contain the same frequency/decay poles and - // keep full-width, long acquisitions bounded; the final LM still uses all - // retained samples. - let pencil_samples = &samples[..samples.len().min(256)]; - for order in 1..=max_order { - evaluated_model_orders += 1; - let Ok(candidate) = matrix_pencil_estimates(pencil_samples, dwell_s, order, fit_bounds) - else { - continue; - }; - if candidate.components.len() != order { - continue; - } - let fit = match evaluate_damped_sinusoids_cancellable( - &samples, - ×, - &candidate.components, - fit_bounds, - cancelled, - ) { - Ok(fit) => Some(fit), - Err(CraftFitError::Cancelled) => return Err(CraftError::Cancelled), - Err(error) => { - warning = Some((CraftWarningKind::FitWindowFailure, error.to_string())); - None - } - }; - if let Some(fit) = fit { - let candidate_bic = bic( - fit.rss, - n_observations, - (fit.components.len() * 4 + 1) as f64, - ); - let separated = fit - .components - .windows(2) - .all(|pair| (pair[1].frequency_hz - pair[0].frequency_hz).abs() >= merge_hz); - if candidate_bic < best_bic && separated && fit.condition_number <= 1e8 { - best_bic = candidate_bic; - best_condition = fit.condition_number; - best_components = fit.components; - } - } - } - if !best_components.is_empty() - && warning - .as_ref() - .is_some_and(|(kind, _)| *kind == CraftWarningKind::FitWindowFailure) - { - warning = None; - } - if best_components.len() == max_order { - warning = Some(( - CraftWarningKind::ModelOrderLimit, - "model order reached the fit-window limit; inspect the residual before quantitation" - .to_owned(), - )); - } - Ok(RegionResult { - components: best_components, - center_hz, - bic: Some(best_bic), - condition_number: best_condition, - decimation, - retained_samples: samples.len(), - evaluated_model_orders, - warning, - }) -} - -fn low_pass_fir( - input: &[Complex64], - sample_rate_hz: f64, - cutoff_hz: f64, - requested_taps: usize, - cancelled: &impl Fn() -> bool, -) -> Result, CraftError> { - if cutoff_hz * 2.0 >= sample_rate_hz * 0.999 { - return Ok(input.to_vec()); - } - let taps = effective_filter_taps(requested_taps, input.len()); - if taps < 3 { - return Ok(input.to_vec()); - } - let half = taps / 2; - let normalized = cutoff_hz / sample_rate_hz; - let mut kernel = Vec::with_capacity(taps); - for index in 0..taps { - let x = index as isize - half as isize; - let sinc = if x == 0 { - 2.0 * normalized - } else { - (TAU * normalized * x as f64).sin() / (PI * x as f64) - }; - let window = 0.42 - 0.5 * (TAU * index as f64 / (taps - 1) as f64).cos() - + 0.08 * (2.0 * TAU * index as f64 / (taps - 1) as f64).cos(); - kernel.push(sinc * window); - } - let sum: f64 = kernel.iter().sum(); - for coefficient in &mut kernel { - *coefficient /= sum; - } - let mut output = vec![Complex64::new(0.0, 0.0); input.len()]; - for (center, filtered) in output - .iter_mut() - .enumerate() - .take(input.len().saturating_sub(half)) - .skip(half) - { - if center % 64 == 0 && cancelled() { - return Err(CraftError::Cancelled); - } - let start = center - half; - *filtered = input[start..start + taps] - .iter() - .zip(&kernel) - .fold(Complex64::new(0.0, 0.0), |sum, (&sample, &coefficient)| { - sum + sample * coefficient - }); - } - Ok(output) -} - -fn effective_filter_taps(requested_taps: usize, input_len: usize) -> usize { - let taps = requested_taps.min(input_len.saturating_sub(1)); - if taps.is_multiple_of(2) { - taps.saturating_sub(1) - } else { - taps - } -} - fn estimate_complex_noise(values: &[Complex64]) -> f64 { let start = values.len().saturating_sub((values.len() / 4).max(64)); let tail = &values[start..]; @@ -777,10 +604,6 @@ fn median(values: &mut [f64]) -> f64 { } } -fn bic(rss: f64, observations: f64, parameters: f64) -> f64 { - observations * (rss.max(f64::MIN_POSITIVE) / observations).ln() + parameters * observations.ln() -} - #[cfg(test)] #[path = "craft_tests.rs"] mod tests; diff --git a/crates/processing/src/craft/diagnostics.rs b/crates/processing/src/craft/diagnostics.rs index ab652e0..4ccc5b5 100644 --- a/crates/processing/src/craft/diagnostics.rs +++ b/crates/processing/src/craft/diagnostics.rs @@ -17,21 +17,62 @@ pub struct CraftDiagnostics { pub normalized_residual: f64, /// `None` means at least one fitted design was rank deficient or unbounded. pub maximum_condition_number: Option, - /// One entry per internal fit window. A user region can require several - /// windows, but those windows never become user-visible region identities. - pub fit_windows: Vec, + /// One entry per protocol-owned modeling window. These are independent of + /// user-visible signal-region identities. + pub modeling_windows: Vec, pub warnings: Vec, + pub stability: CraftStabilityDiagnostics, +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftStabilityDiagnostics { + pub delta_ppm: f64, + pub regions: Vec, + pub ratio: Option, + pub passed: bool, + pub skipped: Vec, +} + +impl Default for CraftStabilityDiagnostics { + fn default() -> Self { + Self { + delta_ppm: 0.0, + regions: Vec::new(), + ratio: None, + passed: false, + skipped: Vec::new(), + } + } +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftStabilityRegion { + pub region: CraftRegionId, + pub metric: CraftStabilityMetric, + pub component_count_min: usize, + pub component_count_max: usize, + pub model_order_min: usize, + pub model_order_max: usize, +} + +#[derive(Clone, Copy, Debug, Default, PartialEq, Serialize, Deserialize)] +pub struct CraftStabilityMetric { + pub median: f64, + pub minimum: f64, + pub maximum: f64, + pub relative_dispersion: f64, } #[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)] #[serde(rename_all = "snake_case")] pub enum CraftWarningKind { - FitWindowFailure, + ModelingWindowFailure, ModelOrderLimit, EmptyRegion, LinewidthAtBound, UnboundedUncertainty, IllConditionedFit, + StabilityFailure, LowAmplitudeThreshold, SsfpQuantitation, InputAssessment, @@ -41,7 +82,7 @@ pub enum CraftWarningKind { pub struct CraftWarning { pub kind: CraftWarningKind, pub region: Option, - pub fit_window: Option, + pub modeling_window: Option, pub message: String, } @@ -49,27 +90,30 @@ impl CraftWarning { pub fn blocks_quantitation(&self) -> bool { matches!( self.kind, - CraftWarningKind::FitWindowFailure + CraftWarningKind::ModelingWindowFailure | CraftWarningKind::ModelOrderLimit | CraftWarningKind::EmptyRegion | CraftWarningKind::LinewidthAtBound | CraftWarningKind::UnboundedUncertainty | CraftWarningKind::IllConditionedFit + | CraftWarningKind::StabilityFailure ) } } #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] -pub struct CraftFitWindowDiagnostic { - pub region: CraftRegionId, - pub core_hz: (f64, f64), - pub padded_hz: (f64, f64), - pub actual_decimation: usize, - pub retained_samples: usize, +pub struct CraftModelingWindowDiagnostic { + pub retention_band_hz: (f64, f64), + pub modeling_band_hz: (f64, f64), + pub decimation_factor: usize, + pub modeled_sample_count: usize, pub evaluated_model_orders: usize, pub selected_model_order: usize, - pub bic: Option, + pub training_bic: Option, pub condition_number: Option, + pub modeled_duration_s: f64, + pub training_normalized_residual: f64, + pub validation_normalized_residual: f64, } #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] diff --git a/crates/processing/src/craft/fitting.rs b/crates/processing/src/craft/fitting.rs new file mode 100644 index 0000000..ae2669a --- /dev/null +++ b/crates/processing/src/craft/fitting.rs @@ -0,0 +1,496 @@ +use num_complex::Complex64; +use plotx_analysis::craft::{ + CraftFitBounds, CraftFitError, CraftFitOptions, DampedSinusoid, backward_linear_predict, + fit_damped_sinusoids_initialized_cancellable, matrix_pencil_estimates, +}; +use std::f64::consts::{PI, TAU}; + +use super::regions::ModelingWindow; +use super::{CraftError, CraftModelingPolicy, CraftParams, CraftProfile, CraftWarningKind}; + +pub(super) struct ModelingWindowResult { + pub(super) components: Vec, + pub(super) center_hz: f64, + pub(super) training_bic: Option, + pub(super) condition_number: f64, + pub(super) decimation: usize, + pub(super) modeled_sample_count: usize, + pub(super) evaluated_model_orders: usize, + pub(super) modeled_duration_s: f64, + pub(super) training_normalized_residual: f64, + pub(super) validation_normalized_residual: f64, + pub(super) warning: Option<(CraftWarningKind, String)>, +} + +pub(super) struct CraftModelingContext<'a> { + pub(super) input: &'a [Complex64], + pub(super) skipped_points: usize, + pub(super) group_delay_points: f64, + pub(super) spectral_width_hz: f64, + pub(super) params: &'a CraftParams, + pub(super) policy: CraftModelingPolicy, +} + +struct ValidatedCandidate { + order: usize, + components: Vec, + training_bic: f64, + condition_number: f64, + training_normalized_residual: f64, + validation_normalized_residual: f64, +} + +pub(super) fn fit_modeling_window( + context: &CraftModelingContext<'_>, + window: ModelingWindow, + cancelled: &impl Fn() -> bool, +) -> Result { + let CraftModelingContext { + input, + skipped_points, + group_delay_points, + spectral_width_hz: sw, + params, + policy, + } = context; + let center_hz = (window.modeling_band_hz.0 + window.modeling_band_hz.1) * 0.5; + let modeled_bandwidth_hz = window.modeling_band_hz.1 - window.modeling_band_hz.0; + // Include one filter length of guard samples so the modeled interval has + // complete centered-FIR support at both ends. + let modeled_points = (policy.modeling_duration_s * *sw).ceil().max(1.0) as usize; + let filter_input_len = input + .len() + .min(modeled_points.saturating_add(params.fir_filter_taps)); + let mixed: Vec = input[..filter_input_len] + .iter() + .enumerate() + .map(|(index, &value)| { + // Raw point numbers preserve the fractional group-delay time origin + // after the digital-filter transient has been skipped. + let time = (*skipped_points as f64 + index as f64 - *group_delay_points) / *sw; + value * Complex64::from_polar(1.0, -TAU * center_hz * time) + }) + .collect(); + let filtered = low_pass_fir( + &mixed, + *sw, + modeled_bandwidth_hz * 0.5, + params.fir_filter_taps, + cancelled, + )?; + let mut decimation = (*sw / (2.0 * modeled_bandwidth_hz).max(f64::MIN_POSITIVE)) + .floor() + .max(1.0) as usize; + decimation = decimation.max(filtered.len().div_ceil(params.maximum_modeled_sample_count)); + let filter_half = effective_filter_taps(params.fir_filter_taps, mixed.len()) / 2; + let valid_end = filtered.len().saturating_sub(filter_half); + // Match the established CRAFT digital-filter workflow: retain the early + // record through a phase-preserving FIR precharge, then replace the five + // boundary-dependent downsampled points by backward linear prediction. + let modeling_start = 0_usize; + let useful_end = modeling_start.saturating_add(modeled_points).min(valid_end); + let mut samples: Vec = filtered[modeling_start..useful_end] + .iter() + .step_by(decimation) + .copied() + .collect(); + let predicted_count = samples.len().min(5); + let training_count = samples.len().saturating_sub(predicted_count).min(256); + let configured_order = if samples.len() > 261 { 32 } else { 16 }; + let prediction_order = configured_order.min(training_count.saturating_sub(1) / 2); + if params.profile == CraftProfile::Conventional && predicted_count > 0 && prediction_order > 0 { + match backward_linear_predict( + &mut samples, + predicted_count, + training_count, + prediction_order, + ) { + Ok(()) => {} + // A rankless no-signal record has nothing to predict. Keep the + // phase-preserving FIR precharge so exploratory runs can report + // the empty window instead of failing the complete invocation. + Err(CraftFitError::Singular) => {} + Err(CraftFitError::Cancelled) => return Err(CraftError::Cancelled), + Err(error) => return Err(CraftError::Fit(error)), + } + } + let times: Vec = (0..samples.len()) + .map(|index| { + (*skipped_points as f64 + modeling_start as f64 + (index * decimation) as f64 + - *group_delay_points) + / *sw + }) + .collect(); + if samples.len() < 16 { + return Ok(ModelingWindowResult { + components: Vec::new(), + center_hz, + training_bic: None, + condition_number: 1.0, + decimation, + modeled_sample_count: samples.len(), + evaluated_model_orders: 0, + modeled_duration_s: samples.len() as f64 * decimation as f64 / *sw, + training_normalized_residual: 1.0, + validation_normalized_residual: 1.0, + warning: Some(( + CraftWarningKind::ModelingWindowFailure, + "too few samples remained after filtering".to_owned(), + )), + }); + } + + let relative_frequency_bounds = ( + window.modeling_band_hz.0 - center_hz, + window.modeling_band_hz.1 - center_hz, + ); + let validation_count = if samples.len() >= 32 { + ((samples.len() as f64 * policy.validation_tail_fraction).round() as usize).max(8) + } else { + 0 + }; + let validation_start = samples.len() - validation_count; + let training_samples = samples.as_slice(); + let training_times = times.as_slice(); + let validation_samples = &samples[validation_start..]; + let validation_times = ×[validation_start..]; + let training_energy = training_samples + .iter() + .map(Complex64::norm_sqr) + .sum::() + .max(f64::MIN_POSITIVE); + let validation_energy = validation_samples + .iter() + .map(Complex64::norm_sqr) + .sum::() + .max(f64::MIN_POSITIVE); + let observation_count = (training_samples.len() * 2) as f64; + let mut candidates = Vec::new(); + let mut warning = None; + + let fit_bounds = CraftFitBounds { + frequency_hz: relative_frequency_bounds, + linewidth_hz: policy.component_linewidth_bounds_hz, + }; + let dwell_s = decimation as f64 / *sw; + let minimum_separation_hz = *sw / input.len() as f64; + let maximum_order = params + .maximum_model_order + .min(training_samples.len() / 2 - 1); + let mut evaluated_model_orders = 0; + // Candidate generation is deliberately bounded because the Hankel SVD is + // cubic. Final amplitudes, evidence, covariance, and residuals use the + // complete sub-FID, matching the established CRAFT workflow. + let pencil_samples = &training_samples[..training_samples.len().min(256)]; + for order in 1..=maximum_order { + evaluated_model_orders += 1; + let Ok(candidate) = matrix_pencil_estimates(pencil_samples, dwell_s, order, fit_bounds) + else { + continue; + }; + if candidate.components.len() != order { + continue; + } + let fit = match fit_damped_sinusoids_initialized_cancellable( + training_samples, + training_times, + &candidate.components, + fit_bounds, + CraftFitOptions::default(), + cancelled, + ) { + Ok(fit) => Some(fit), + Err(CraftFitError::Cancelled) => return Err(CraftError::Cancelled), + Err(error) => { + warning = Some((CraftWarningKind::ModelingWindowFailure, error.to_string())); + None + } + }; + if let Some(mut fit) = fit { + fit.components + .sort_by(|left, right| left.frequency_hz.total_cmp(&right.frequency_hz)); + let candidate_bic = bic( + fit.rss, + observation_count, + (fit.components.len() * 4 + 1) as f64, + ); + let sufficiently_separated = fit.components.windows(2).all(|pair| { + (pair[1].frequency_hz - pair[0].frequency_hz).abs() >= minimum_separation_hz + }); + let linewidths_are_interior = fit.components.iter().all(|component| { + component.linewidth_hz > fit_bounds.linewidth_hz.0 + 1e-6 + && component.linewidth_hz < fit_bounds.linewidth_hz.1 - 1e-6 + }); + let uncertainties_are_bounded = fit.components.iter().all(|component| { + component.amplitude_std.is_some() + && component.frequency_std_hz.is_some() + && component.linewidth_std_hz.is_some() + && component.phase_std_rad.is_some() + }); + let model_amplitude_to_noise = coherent_amplitude(&fit.components) + / fit + .components + .iter() + .filter_map(|component| component.amplitude_std) + .map(|value| value * value) + .sum::() + .sqrt() + .max(f64::MIN_POSITIVE); + if sufficiently_separated + && linewidths_are_interior + && uncertainties_are_bounded + && model_amplitude_to_noise >= params.minimum_amplitude_to_noise + && fit.condition_number <= 1e8 + { + let validation_rss = if validation_samples.is_empty() { + fit.rss + } else { + validation_samples + .iter() + .zip(validation_times) + .map(|(&sample, &time)| { + (sample - damped_model_at(&fit.components, time)).norm_sqr() + }) + .sum() + }; + candidates.push(ValidatedCandidate { + order, + components: fit.components, + training_bic: candidate_bic, + condition_number: fit.condition_number, + training_normalized_residual: (fit.rss / training_energy).sqrt(), + validation_normalized_residual: if validation_samples.is_empty() { + (fit.rss / training_energy).sqrt() + } else { + (validation_rss / validation_energy).sqrt() + }, + }); + } + } + } + + let minimum_training_bic = candidates + .iter() + .map(|candidate| candidate.training_bic) + .min_by(f64::total_cmp) + .unwrap_or_else(|| bic(training_energy, observation_count, 1.0)); + // Bretthorst CRAFT selects model order from the evidence in the modeled + // record. BIC is the deterministic evidence approximation used here; a + // two-unit band retains the simplest statistically comparable model. + let comparable_bic_limit = minimum_training_bic + 2.0; + candidates.sort_by_key(|candidate| candidate.order); + let selected = candidates + .into_iter() + .find(|candidate| candidate.training_bic <= comparable_bic_limit); + let (mut components, training_bic, condition_number, training_residual, validation_residual) = + selected.map_or_else( + || { + ( + Vec::new(), + bic(training_energy, observation_count, 1.0), + 1.0, + 1.0, + 1.0, + ) + }, + |candidate| { + ( + candidate.components, + candidate.training_bic, + candidate.condition_number, + candidate.training_normalized_residual, + candidate.validation_normalized_residual, + ) + }, + ); + // Very small poles at a window edge are commonly transition-band leakage + // or a split of the dominant line, not an independently quantifiable + // resonance. Apply the threshold to the selected multiplet so weak lines + // are retained relative to their local partner rather than compared with + // a global raw-FID noise estimate. + if let Some(maximum_amplitude) = components + .iter() + .map(|component| component.amplitude) + .max_by(f64::total_cmp) + { + let minimum_amplitude = maximum_amplitude * 0.05; + components.retain(|component| component.amplitude >= minimum_amplitude); + } + if !components.is_empty() + && warning + .as_ref() + .is_some_and(|(kind, _)| *kind == CraftWarningKind::ModelingWindowFailure) + { + warning = None; + } + let actual_taps = effective_filter_taps(params.fir_filter_taps, mixed.len()); + for component in &mut components { + let gain = fir_response( + *sw, + modeled_bandwidth_hz * 0.5, + actual_taps, + component.frequency_hz, + component.decay_rate_s_inv, + ); + let gain_norm = gain.norm(); + if gain_norm > 0.1 { + component.amplitude /= gain_norm; + component.amplitude_std = component.amplitude_std.map(|value| value / gain_norm); + component.phase_rad -= gain.arg(); + } + } + if components.len() == maximum_order { + warning = Some(( + CraftWarningKind::ModelOrderLimit, + "model order reached the modeling-window limit; inspect the residual before quantitation" + .to_owned(), + )); + } + if !components.is_empty() && training_residual > 0.25 { + warning = Some(( + CraftWarningKind::ModelingWindowFailure, + format!( + "modeled-record residual {training_residual:.3} exceeded the quantitative limit" + ), + )); + } + Ok(ModelingWindowResult { + components, + center_hz, + training_bic: Some(training_bic), + condition_number, + decimation, + modeled_sample_count: samples.len(), + evaluated_model_orders, + modeled_duration_s: samples.len() as f64 * decimation as f64 / *sw, + training_normalized_residual: training_residual, + validation_normalized_residual: validation_residual, + warning, + }) +} + +fn damped_model_at(components: &[DampedSinusoid], time_s: f64) -> Complex64 { + components + .iter() + .fold(Complex64::new(0.0, 0.0), |sum, component| { + sum + Complex64::from_polar( + component.amplitude * (-component.decay_rate_s_inv * time_s).exp(), + component.phase_rad + TAU * component.frequency_hz * time_s, + ) + }) +} + +fn coherent_amplitude(components: &[DampedSinusoid]) -> f64 { + components + .iter() + .fold(Complex64::new(0.0, 0.0), |sum, component| { + sum + Complex64::from_polar(component.amplitude, component.phase_rad) + }) + .norm() +} + +fn low_pass_fir( + input: &[Complex64], + sample_rate_hz: f64, + cutoff_hz: f64, + requested_taps: usize, + cancelled: &impl Fn() -> bool, +) -> Result, CraftError> { + if cutoff_hz * 2.0 >= sample_rate_hz * 0.999 { + return Ok(input.to_vec()); + } + let taps = effective_filter_taps(requested_taps, input.len()); + if taps < 3 { + return Ok(input.to_vec()); + } + let half = taps / 2; + let kernel = fir_kernel(sample_rate_hz, cutoff_hz, taps); + let mut output = vec![Complex64::new(0.0, 0.0); input.len()]; + let first = input[0]; + let reflection = if first.norm_sqr() > f64::MIN_POSITIVE { + first / first.conj() + } else { + Complex64::new(1.0, 0.0) + }; + for (center, filtered) in output + .iter_mut() + .enumerate() + .take(input.len().saturating_sub(half)) + { + if center % 64 == 0 && cancelled() { + return Err(CraftError::Cancelled); + } + *filtered = kernel.iter().enumerate().fold( + Complex64::new(0.0, 0.0), + |sum, (index, &coefficient)| { + let source = center as isize + index as isize - half as isize; + let sample = if source >= 0 { + input[source as usize] + } else { + reflection * input[source.unsigned_abs()].conj() + }; + sum + sample * coefficient + }, + ); + } + Ok(output) +} + +fn fir_response( + sample_rate_hz: f64, + cutoff_hz: f64, + taps: usize, + frequency_hz: f64, + decay_rate_s_inv: f64, +) -> Complex64 { + if cutoff_hz * 2.0 >= sample_rate_hz * 0.999 || taps < 3 { + return Complex64::new(1.0, 0.0); + } + let half = taps / 2; + fir_kernel(sample_rate_hz, cutoff_hz, taps) + .iter() + .enumerate() + .fold(Complex64::new(0.0, 0.0), |sum, (index, coefficient)| { + let offset_s = (index as f64 - half as f64) / sample_rate_hz; + sum + Complex64::from_polar( + coefficient * (-decay_rate_s_inv * offset_s).exp(), + TAU * frequency_hz * offset_s, + ) + }) +} + +fn fir_kernel(sample_rate_hz: f64, cutoff_hz: f64, taps: usize) -> Vec { + let half = taps / 2; + let normalized = cutoff_hz / sample_rate_hz; + let mut kernel = (0..taps) + .map(|index| { + let x = index as isize - half as isize; + let sinc = if x == 0 { + 2.0 * normalized + } else { + (TAU * normalized * x as f64).sin() / (PI * x as f64) + }; + let window = 0.42 - 0.5 * (TAU * index as f64 / (taps - 1) as f64).cos() + + 0.08 * (2.0 * TAU * index as f64 / (taps - 1) as f64).cos(); + sinc * window + }) + .collect::>(); + let sum = kernel.iter().sum::(); + for coefficient in &mut kernel { + *coefficient /= sum; + } + kernel +} + +fn effective_filter_taps(requested_taps: usize, input_len: usize) -> usize { + let taps = requested_taps.min(input_len.saturating_sub(1)); + if taps.is_multiple_of(2) { + taps.saturating_sub(1) + } else { + taps + } +} + +fn bic(rss: f64, observations: f64, parameters: f64) -> f64 { + observations * (rss.max(f64::MIN_POSITIVE) / observations).ln() + parameters * observations.ln() +} diff --git a/crates/processing/src/craft/preflight.rs b/crates/processing/src/craft/preflight.rs index f2e034f..a2909ef 100644 --- a/crates/processing/src/craft/preflight.rs +++ b/crates/processing/src/craft/preflight.rs @@ -2,6 +2,7 @@ use plotx_analysis::peaks::{DetectParams, detect_peaks, estimate_noise}; use plotx_io::{Domain, NmrData}; use rustfft::FftPlanner; use serde::{Deserialize, Serialize}; +use std::f64::consts::PI; use super::{CraftDerivedPlan, CraftParams, CraftReference, CraftRegionId}; @@ -36,7 +37,7 @@ pub enum CraftIssueAction { CheckImport, CheckAcquisitionMetadata, CorrectReference, - ResetFitSettings, + ResetModelingSettings, AdjustRegions, ReduceSkippedPoints, ReviewAcquisition, @@ -65,7 +66,7 @@ pub struct CraftInputAssessment { pub point_count: usize, pub effective_point_count: usize, pub acquisition_duration_s: Option, - pub fit_window_count: usize, + pub modeling_window_count: usize, pub clear_signals: Vec, pub issues: Vec, } @@ -136,8 +137,8 @@ impl CraftInputAssessment { if params.validate().is_err() { error( CraftIssueCode::InvalidParameters, - "One or more explicit fit settings are invalid.", - CraftIssueAction::ResetFitSettings, + "One or more explicit component or acquisition settings are invalid.", + CraftIssueAction::ResetModelingSettings, ); } if regions_outside_bandwidth_or_overlap(data, reference, params) { @@ -155,14 +156,14 @@ impl CraftInputAssessment { ); } if plan - .fit_windows + .modeling_windows .iter() - .any(|window| window.planned_retained_samples < 16) + .any(|window| window.planned_modeled_sample_count < 16) { error( CraftIssueCode::TooFewEffectivePoints, - "Fewer than 16 samples remain in one or more fit windows after FIR filtering.", - CraftIssueAction::ResetFitSettings, + "Fewer than 16 samples remain in one or more modeling windows after FIR filtering.", + CraftIssueAction::ResetModelingSettings, ); } } @@ -193,9 +194,10 @@ impl CraftInputAssessment { if count == 0 && !clear_signals.is_empty() { issues.push(warning(CraftIssueCode::RegionWithoutClearSignal, Some(region.id), "A selected region contains no clear signal; adjust the region or confirm it with independent evidence.", CraftIssueAction::AdjustRegions)); } - if count >= params.max_components_per_fit_window { - issues.push(warning(CraftIssueCode::DenseSignalWindow, Some(region.id), "Detected peak density reaches the model-order limit; consider narrowing the region or increasing the limit.", CraftIssueAction::IncreaseModelLimit)); - } + // Peak-picking is only a preflight hint. The FFT can contain many + // transition-band extrema for one physical multiplet, so raw peak + // count must not be used as a model-capacity warning. Capacity is + // assessed after the bounded time-domain fit has selected a model. } Self { point_count: data.points.len(), @@ -203,7 +205,7 @@ impl CraftInputAssessment { acquisition_duration_s: (data.spectral_width_hz.is_finite() && data.spectral_width_hz > 0.0) .then(|| data.points.len() as f64 / data.spectral_width_hz), - fit_window_count: plan.fit_windows.len(), + modeling_window_count: plan.modeling_windows.len(), clear_signals, issues, } @@ -278,7 +280,7 @@ fn regions_outside_bandwidth_or_overlap( .any(|pair| pair[0].end_ppm > pair[1].start_ppm) } -fn detect_clear_signals( +pub(super) fn detect_clear_signals( data: &NmrData, reference: CraftReference, skip: usize, @@ -289,7 +291,12 @@ fn detect_clear_signals( } let fft_len = input.len().next_power_of_two(); let mut spectrum = vec![num_complex::Complex64::new(0.0, 0.0); fft_len]; - spectrum[..input.len()].copy_from_slice(input); + let duration_s = input.len() as f64 / data.spectral_width_hz; + let matched_line_broadening_hz = 1.0 / duration_s.max(f64::MIN_POSITIVE); + for (index, (&sample, output)) in input.iter().zip(&mut spectrum).enumerate() { + let time_s = index as f64 / data.spectral_width_hz; + *output = sample * (-PI * matched_line_broadening_hz * time_s).exp(); + } FftPlanner::::new() .plan_fft_forward(fft_len) .process(&mut spectrum); @@ -308,7 +315,14 @@ fn detect_clear_signals( &DetectParams { min_height: Some(6.0 * sigma), min_prominence: 5.0 * sigma, - min_spacing: None, + // Merge FFT extrema closer than one acquired spectral + // resolution element (1/acquisition time). Matched exponential + // apodization broadens a line and otherwise creates several + // equally significant extrema for a single resonance. + // `xs` is expressed in zero-padded FFT-bin indices. One acquired + // spectral-resolution element spans `fft_len / input.len()` of + // those bins when the FFT is zero-padded. + min_spacing: Some(fft_len as f64 / input.len() as f64), max_count: Some(64), }, ); diff --git a/crates/processing/src/craft/regions.rs b/crates/processing/src/craft/regions.rs index ed8affd..c5d3001 100644 --- a/crates/processing/src/craft/regions.rs +++ b/crates/processing/src/craft/regions.rs @@ -2,14 +2,13 @@ use plotx_io::NmrData; use super::{ CraftComponent, CraftError, CraftParams, CraftReference, CraftRegion, CraftRegionId, - CraftRegionRatio, CraftRegionSummary, + CraftRegionRatio, CraftRegionSummary, CraftSignalSuggestion, }; -#[derive(Clone, Copy)] -pub(super) struct HzRegion { - pub(super) selection: CraftRegion, - pub(super) core: (f64, f64), - pub(super) padded: (f64, f64), +#[derive(Clone, Copy, Debug)] +pub(super) struct ModelingWindow { + pub(super) retention_band_hz: (f64, f64), + pub(super) modeling_band_hz: (f64, f64), } pub(super) fn selections_are_valid(regions: &[CraftRegion]) -> bool { @@ -37,11 +36,12 @@ pub(super) fn selections_are_valid(regions: &[CraftRegion]) -> bool { .all(|pair| pair[0].end_ppm <= pair[1].start_ppm) } -pub(super) fn build_regions( +pub(super) fn build_modeling_windows( data: &NmrData, params: &CraftParams, reference: CraftReference, -) -> Result, CraftError> { + clear_signals: &[CraftSignalSuggestion], +) -> Result, CraftError> { let half_sw = data.spectral_width_hz * 0.5; let effective_carrier_ppm = reference.effective_carrier_ppm(); let requested: Vec<(CraftRegion, f64, f64)> = if params.regions.is_empty() { @@ -88,31 +88,101 @@ pub(super) fn build_regions( return Err(CraftError::InvalidParameters); } - let mut regions = Vec::new(); - for (selection, start, end) in requested_cores { - let pieces = ((end - start) / params.max_fit_window_width_hz) - .ceil() - .max(1.0) as usize; - let width = (end - start) / pieces as f64; - for index in 0..pieces { - let core = ( - start + index as f64 * width, - start + (index + 1) as f64 * width, - ); - let padding = (core.1 - core.0) * params.padding_fraction * 0.5; - regions.push(HzRegion { - selection, - core, - padded: ( - (core.0 - padding).max(-half_sw), - (core.1 + padding).min(half_sw), - ), + // Modeling windows are a profile-owned protocol. They tile the acquired + // bandwidth with a fixed physical width and therefore do not change when a + // user nudges a reporting region boundary. + let width = params + .profile + .modeling_bandwidth_hz() + .min(data.spectral_width_hz); + let signal_hz = clear_signals + .iter() + .filter_map(|signal| { + let frequency = + (signal.chemical_shift_ppm - effective_carrier_ppm) * data.observe_freq_mhz; + let weight = signal.prominence_sigma.max(f64::MIN_POSITIVE); + (frequency.is_finite() + && weight.is_finite() + && frequency >= -half_sw + && frequency <= half_sw) + .then_some((frequency, weight)) + }) + .collect::>(); + let mut centers = signal_cluster_centers(&signal_hz, width); + if centers.is_empty() { + let pieces = (data.spectral_width_hz / width).ceil().max(1.0) as usize; + centers.extend((0..pieces).map(|index| { + let start = -half_sw + index as f64 * width; + (start + (start + width).min(half_sw)) * 0.5 + })); + } + let mut regions = Vec::with_capacity(centers.len()); + for (index, center) in centers.iter().copied().enumerate() { + let nominal_start = (center - width * 0.5).max(-half_sw); + let nominal_end = (center + width * 0.5).min(half_sw); + let retention_start = centers + .get(index.wrapping_sub(1)) + .map_or(nominal_start, |previous| { + nominal_start.max((previous + center) * 0.5) }); - } + let retention_end = centers + .get(index + 1) + .map_or(nominal_end, |next| nominal_end.min((center + next) * 0.5)); + regions.push(ModelingWindow { + retention_band_hz: (retention_start, retention_end), + modeling_band_hz: (nominal_start, nominal_end), + }); } Ok(regions) } +fn signal_cluster_centers(signals: &[(f64, f64)], window_width_hz: f64) -> Vec { + let mut ranked = signals.to_vec(); + ranked.sort_by(|left, right| { + right + .1 + .total_cmp(&left.1) + .then_with(|| left.0.total_cmp(&right.0)) + }); + let minimum_center_spacing = window_width_hz * 0.5; + let mut seeds = Vec::new(); + for &(frequency, _) in &ranked { + if seeds + .iter() + .all(|seed: &f64| (frequency - *seed).abs() > minimum_center_spacing) + { + seeds.push(frequency); + } + } + + let mut clusters = vec![Vec::new(); seeds.len()]; + for signal in ranked { + if let Some((index, _)) = seeds.iter().enumerate().min_by(|(_, left), (_, right)| { + (signal.0 - **left) + .abs() + .total_cmp(&(signal.0 - **right).abs()) + }) { + clusters[index].push(signal); + } + } + let mut centers = clusters + .into_iter() + .filter_map(weighted_median_frequency) + .collect::>(); + centers.sort_by(f64::total_cmp); + centers +} + +fn weighted_median_frequency(mut signals: Vec<(f64, f64)>) -> Option { + signals.sort_by(|left, right| left.0.total_cmp(&right.0)); + let half_weight = signals.iter().map(|signal| signal.1).sum::() * 0.5; + let mut accumulated = 0.0; + signals.into_iter().find_map(|(frequency, weight)| { + accumulated += weight; + (accumulated >= half_weight).then_some(frequency) + }) +} + pub(super) fn summarize_regions( components: &[CraftComponent], selections: &[CraftRegion], diff --git a/crates/processing/src/craft/report.rs b/crates/processing/src/craft/report.rs new file mode 100644 index 0000000..33ad7e9 --- /dev/null +++ b/crates/processing/src/craft/report.rs @@ -0,0 +1,222 @@ +use super::{CraftComponent, CraftComponentId, CraftRegionId}; +use num_complex::Complex64; +use serde::{Deserialize, Serialize}; + +/// User-controlled definition of a derived CRAFT amplitude report. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftReportDefinition { + pub threshold_an: f64, + pub segment_width_hz: f64, + /// Empty selects every region in the source run. + #[serde(default)] + pub regions: Vec, +} + +impl Default for CraftReportDefinition { + fn default() -> Self { + Self { + threshold_an: 3.3, + segment_width_hz: 1.0, + regions: Vec::new(), + } + } +} + +impl CraftReportDefinition { + pub fn validate(&self) -> Result<(), CraftReportError> { + if !self.threshold_an.is_finite() || self.threshold_an <= 0.0 { + return Err(CraftReportError::InvalidThreshold); + } + if !self.segment_width_hz.is_finite() || self.segment_width_hz <= 0.0 { + return Err(CraftReportError::InvalidWidth); + } + Ok(()) + } +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftReportSegment { + pub center_hz: f64, + pub start_hz: f64, + pub end_hz: f64, + pub component_ids: Vec, + pub component_count: usize, + pub scalar_amplitude_sum_t0: f64, + pub coherent_amplitude_t0: f64, +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftAmplitudeReport { + pub schema_version: u32, + pub definition: CraftReportDefinition, + pub segments: Vec, +} + +impl CraftAmplitudeReport { + pub fn validate_against(&self, components: &[CraftComponent]) -> Result<(), CraftReportError> { + self.definition.validate()?; + let mut seen = std::collections::HashSet::new(); + for segment in &self.segments { + for id in &segment.component_ids { + if !seen.insert(*id) || !components.iter().any(|component| component.id == *id) { + return Err(CraftReportError::UnknownComponent); + } + } + } + Ok(()) + } +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq, thiserror::Error)] +pub enum CraftReportError { + #[error("CRAFT report threshold must be finite and positive")] + InvalidThreshold, + #[error("CRAFT report segment width must be finite and positive")] + InvalidWidth, + #[error("CRAFT report references an unknown component")] + UnknownComponent, +} + +/// Build a report from the complete retained component list. Components are +/// sorted by frequency, and overlapping windows are merged without duplicate +/// membership, making the result independent of input ordering. +pub fn calculate_craft_report( + components: &[CraftComponent], + definition: CraftReportDefinition, +) -> Result { + definition.validate()?; + let mut selected: Vec<&CraftComponent> = components + .iter() + .filter(|component| { + (definition.regions.is_empty() || definition.regions.contains(&component.region)) + && component.amplitude_to_noise >= definition.threshold_an + }) + .collect(); + selected.sort_by(|a, b| { + a.frequency_hz + .total_cmp(&b.frequency_hz) + .then(a.id.0.cmp(&b.id.0)) + }); + let half = definition.segment_width_hz * 0.5; + let mut segments = Vec::new(); + for component in selected { + let start = component.frequency_hz - half; + let end = component.frequency_hz + half; + let append = segments + .last() + .is_none_or(|segment: &CraftReportSegment| start > segment.end_hz); + if append { + segments.push(CraftReportSegment { + center_hz: component.frequency_hz, + start_hz: start, + end_hz: end, + component_ids: vec![component.id], + component_count: 1, + scalar_amplitude_sum_t0: component.amplitude_t0, + coherent_amplitude_t0: Complex64::from_polar( + component.amplitude_t0, + component.phase_rad, + ) + .norm(), + }); + } else { + let segment = segments.last_mut().expect("segment exists"); + segment.end_hz = segment.end_hz.max(end); + segment.center_hz = (segment.start_hz + segment.end_hz) * 0.5; + if !segment.component_ids.contains(&component.id) { + segment.component_ids.push(component.id); + segment.component_count += 1; + segment.scalar_amplitude_sum_t0 += component.amplitude_t0; + let phase_sum = segment + .component_ids + .iter() + .filter_map(|id| components.iter().find(|candidate| candidate.id == *id)) + .fold(Complex64::new(0.0, 0.0), |sum, item| { + sum + Complex64::from_polar(item.amplitude_t0, item.phase_rad) + }); + segment.coherent_amplitude_t0 = phase_sum.norm(); + } + } + } + Ok(CraftAmplitudeReport { + schema_version: 1, + definition, + segments, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + fn component( + id: u64, + frequency_hz: f64, + amplitude_to_noise: f64, + phase_rad: f64, + ) -> CraftComponent { + CraftComponent { + id: CraftComponentId(id), + region: CraftRegionId(1), + frequency_hz, + chemical_shift_ppm: 0.0, + amplitude_t0: 1.0, + phase_rad, + decay_rate_s_inv: 1.0, + linewidth_hz: 1.0, + amplitude_to_noise, + amplitude_std: None, + frequency_std_hz: None, + linewidth_std_hz: None, + phase_std_rad: None, + } + } + #[test] + fn filters_at_threshold_and_merges_without_duplicates() { + let components = vec![ + component(2, 1.4, 3.3, std::f64::consts::PI), + component(1, 1.0, 3.2, 0.0), + component(3, 1.8, 4.0, 0.0), + ]; + let report = calculate_craft_report( + &components, + CraftReportDefinition { + threshold_an: 3.3, + segment_width_hz: 1.0, + regions: vec![], + }, + ) + .unwrap(); + assert_eq!(report.segments.len(), 1); + assert_eq!( + report.segments[0].component_ids, + vec![CraftComponentId(2), CraftComponentId(3)] + ); + assert_eq!(report.segments[0].scalar_amplitude_sum_t0, 2.0); + assert!((report.segments[0].coherent_amplitude_t0 - 0.0).abs() < 1e-12); + } + #[test] + fn rejects_invalid_definition() { + assert_eq!( + calculate_craft_report( + &[], + CraftReportDefinition { + threshold_an: f64::NAN, + ..Default::default() + } + ) + .unwrap_err(), + CraftReportError::InvalidThreshold + ); + assert_eq!( + calculate_craft_report( + &[], + CraftReportDefinition { + segment_width_hz: 0.0, + ..Default::default() + } + ) + .unwrap_err(), + CraftReportError::InvalidWidth + ); + } +} diff --git a/crates/processing/src/craft/resolution.rs b/crates/processing/src/craft/resolution.rs index 5c81a66..b6e1d4d 100644 --- a/crates/processing/src/craft/resolution.rs +++ b/crates/processing/src/craft/resolution.rs @@ -1,6 +1,8 @@ use plotx_io::NmrData; use serde::{Deserialize, Serialize}; +use super::preflight::detect_clear_signals; +use super::regions::build_modeling_windows; use super::{ CraftInputAssessment, CraftInvocation, CraftParams, CraftProfile, CraftReference, CraftRegion, CraftRegionId, @@ -19,13 +21,11 @@ pub enum CraftParamSource { pub struct CraftParamOverrides { pub profile: Option, pub regions: Option>, - pub max_components_per_fit_window: Option, - pub min_amplitude_to_noise: Option, - pub linewidth_hz: Option<(f64, f64)>, - pub filter_taps: Option, - pub padding_fraction: Option, - pub max_fit_window_width_hz: Option, - pub max_downsampled_points: Option, + pub maximum_model_order: Option, + pub minimum_amplitude_to_noise: Option, + pub component_linewidth_bounds_hz: Option<(f64, f64)>, + pub fir_filter_taps: Option, + pub maximum_modeled_sample_count: Option, pub skip_duration_s: Option, pub reconstruction_duration_s: Option>, } @@ -36,13 +36,11 @@ impl CraftParamOverrides { Self { profile: Some(params.profile), regions: Some(params.regions), - max_components_per_fit_window: Some(params.max_components_per_fit_window), - min_amplitude_to_noise: Some(params.min_amplitude_to_noise), - linewidth_hz: Some(params.linewidth_hz), - filter_taps: Some(params.filter_taps), - padding_fraction: Some(params.padding_fraction), - max_fit_window_width_hz: Some(params.max_fit_window_width_hz), - max_downsampled_points: Some(params.max_downsampled_points), + maximum_model_order: Some(params.maximum_model_order), + minimum_amplitude_to_noise: Some(params.minimum_amplitude_to_noise), + component_linewidth_bounds_hz: Some(params.component_linewidth_bounds_hz), + fir_filter_taps: Some(params.fir_filter_taps), + maximum_modeled_sample_count: Some(params.maximum_modeled_sample_count), skip_duration_s: Some(params.skip_duration_s), reconstruction_duration_s: Some(params.reconstruction_duration_s), } @@ -64,13 +62,11 @@ impl CraftParamOverrides { pub struct CraftParameterSources { pub profile: CraftParamSource, pub regions: CraftParamSource, - pub max_components_per_fit_window: CraftParamSource, - pub min_amplitude_to_noise: CraftParamSource, - pub linewidth_hz: CraftParamSource, - pub filter_taps: CraftParamSource, - pub padding_fraction: CraftParamSource, - pub max_fit_window_width_hz: CraftParamSource, - pub max_downsampled_points: CraftParamSource, + pub maximum_model_order: CraftParamSource, + pub minimum_amplitude_to_noise: CraftParamSource, + pub component_linewidth_bounds_hz: CraftParamSource, + pub fir_filter_taps: CraftParamSource, + pub maximum_modeled_sample_count: CraftParamSource, pub skip_duration_s: CraftParamSource, pub reconstruction_duration_s: CraftParamSource, } @@ -80,13 +76,11 @@ impl CraftParameterSources { [ self.profile, self.regions, - self.max_components_per_fit_window, - self.min_amplitude_to_noise, - self.linewidth_hz, - self.filter_taps, - self.padding_fraction, - self.max_fit_window_width_hz, - self.max_downsampled_points, + self.maximum_model_order, + self.minimum_amplitude_to_noise, + self.component_linewidth_bounds_hz, + self.fir_filter_taps, + self.maximum_modeled_sample_count, self.skip_duration_s, self.reconstruction_duration_s, ] @@ -95,12 +89,44 @@ impl CraftParameterSources { } #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] -pub struct CraftDerivedWindow { - pub region: CraftRegionId, - pub core_hz: (f64, f64), - pub padded_hz: (f64, f64), - pub planned_decimation: usize, - pub planned_retained_samples: usize, +pub struct CraftDerivedModelingWindow { + pub retention_band_hz: (f64, f64), + pub modeling_band_hz: (f64, f64), + pub planned_decimation_factor: usize, + pub planned_modeled_sample_count: usize, + pub planned_modeled_duration_s: f64, +} + +#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] +pub struct CraftModelingPolicy { + pub modeling_bandwidth_hz: f64, + pub modeling_duration_s: f64, + pub validation_tail_fraction: f64, + pub boundary_stability_relative_tolerance: f64, + pub component_linewidth_bounds_hz: (f64, f64), +} + +impl Default for CraftModelingPolicy { + fn default() -> Self { + Self { + modeling_bandwidth_hz: 250.0, + modeling_duration_s: 1.0, + validation_tail_fraction: 0.2, + boundary_stability_relative_tolerance: 0.01, + component_linewidth_bounds_hz: (0.05, 20.0), + } + } +} + +impl CraftModelingPolicy { + pub(super) fn for_params(params: &CraftParams) -> Self { + Self { + modeling_bandwidth_hz: params.profile.modeling_bandwidth_hz(), + modeling_duration_s: params.profile.modeling_duration_s(), + component_linewidth_bounds_hz: params.component_linewidth_bounds_hz, + ..Self::default() + } + } } #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] @@ -108,10 +134,10 @@ pub struct CraftDerivedPlan { pub effective_skip_points: usize, pub effective_skip_source: CraftParamSource, pub available_points: usize, - pub actual_filter_taps: usize, + pub effective_fir_filter_taps: usize, pub reconstruction_points: usize, pub resolved_regions: Vec, - pub fit_windows: Vec, + pub modeling_windows: Vec, } pub fn resolve_craft_invocation( @@ -160,43 +186,40 @@ pub fn resolve_craft_invocation( regions = full_bandwidth_region(data, reference).into_iter().collect(); regions_source = CraftParamSource::InputDerived; } - let (max_components_per_fit_window, max_components_source) = - resolved!(max_components_per_fit_window); - let (min_amplitude_to_noise, min_amplitude_source) = resolved!(min_amplitude_to_noise); - let (linewidth_hz, linewidth_source) = resolved!(linewidth_hz); - let (filter_taps, filter_taps_source) = resolved!(filter_taps); - let (padding_fraction, padding_source) = resolved!(padding_fraction); - let (max_fit_window_width_hz, fit_window_source) = resolved!(max_fit_window_width_hz); - let (max_downsampled_points, downsample_source) = resolved!(max_downsampled_points); + let (maximum_model_order, maximum_model_order_source) = resolved!(maximum_model_order); + let (minimum_amplitude_to_noise, minimum_amplitude_source) = + resolved!(minimum_amplitude_to_noise); + let (component_linewidth_bounds_hz, component_linewidth_bounds_source) = + resolved!(component_linewidth_bounds_hz); + let (fir_filter_taps, fir_filter_taps_source) = resolved!(fir_filter_taps); + let (maximum_modeled_sample_count, maximum_modeled_sample_count_source) = + resolved!(maximum_modeled_sample_count); let (skip_duration_s, skip_source) = resolved!(skip_duration_s); let (reconstruction_duration_s, reconstruction_source) = resolved!(reconstruction_duration_s); let params = CraftParams { profile, regions, - max_components_per_fit_window, - min_amplitude_to_noise, - linewidth_hz, - filter_taps, - padding_fraction, - max_fit_window_width_hz, - max_downsampled_points, + maximum_model_order, + minimum_amplitude_to_noise, + component_linewidth_bounds_hz, + fir_filter_taps, + maximum_modeled_sample_count, skip_duration_s, reconstruction_duration_s, }; let sources = CraftParameterSources { profile: profile_source, regions: regions_source, - max_components_per_fit_window: max_components_source, - min_amplitude_to_noise: min_amplitude_source, - linewidth_hz: linewidth_source, - filter_taps: filter_taps_source, - padding_fraction: padding_source, - max_fit_window_width_hz: fit_window_source, - max_downsampled_points: downsample_source, + maximum_model_order: maximum_model_order_source, + minimum_amplitude_to_noise: minimum_amplitude_source, + component_linewidth_bounds_hz: component_linewidth_bounds_source, + fir_filter_taps: fir_filter_taps_source, + maximum_modeled_sample_count: maximum_modeled_sample_count_source, skip_duration_s: skip_source, reconstruction_duration_s: reconstruction_source, }; - let derived_plan = derive_plan(data, reference, ¶ms, &sources); + let modeling_policy = CraftModelingPolicy::for_params(¶ms); + let derived_plan = derive_plan(data, reference, ¶ms, &sources, modeling_policy); let assessment = CraftInputAssessment::assess(data, reference, ¶ms, &derived_plan); CraftInvocation { params, @@ -204,6 +227,7 @@ pub fn resolve_craft_invocation( reference, derived_plan, assessment, + modeling_policy, } } @@ -224,6 +248,7 @@ fn derive_plan( reference: CraftReference, params: &CraftParams, sources: &CraftParameterSources, + modeling_policy: CraftModelingPolicy, ) -> CraftDerivedPlan { let requested_skip = if data.spectral_width_hz.is_finite() && params.skip_duration_s.is_finite() { @@ -245,9 +270,16 @@ fn derive_plan( sources.skip_duration_s }; let available_points = data.points.len().saturating_sub(effective_skip_points); - let actual_filter_taps = effective_taps( - params.filter_taps, - available_points.min(6_000_usize.saturating_add(params.filter_taps)), + let fit_points = if data.spectral_width_hz.is_finite() && data.spectral_width_hz > 0.0 { + (modeling_policy.modeling_duration_s * data.spectral_width_hz) + .ceil() + .max(1.0) as usize + } else { + 0 + }; + let effective_fir_filter_taps = effective_taps( + params.fir_filter_taps, + available_points.min(fit_points.saturating_add(params.fir_filter_taps)), ); let acquired_duration = if data.spectral_width_hz.is_finite() && data.spectral_width_hz > 0.0 { data.points.len() as f64 / data.spectral_width_hz @@ -269,73 +301,50 @@ fn derive_plan( 0 }; let resolved_regions = params.regions.clone(); - let mut fit_windows = Vec::new(); + let mut modeling_windows = Vec::new(); if data.spectral_width_hz.is_finite() && data.spectral_width_hz > 0.0 && data.observe_freq_mhz.is_finite() && data.observe_freq_mhz > 0.0 && reference.effective_carrier_ppm().is_finite() - && params.max_fit_window_width_hz.is_finite() - && params.max_fit_window_width_hz > 0.0 + && params.profile.modeling_bandwidth_hz().is_finite() { - let half_sw = data.spectral_width_hz * 0.5; - let carrier = reference.effective_carrier_ppm(); - for selection in &resolved_regions { - let normalized = selection.normalized(); - let start = - ((normalized.start_ppm - carrier) * data.observe_freq_mhz).clamp(-half_sw, half_sw); - let end = - ((normalized.end_ppm - carrier) * data.observe_freq_mhz).clamp(-half_sw, half_sw); - if !start.is_finite() || !end.is_finite() || end <= start { - continue; - } - let pieces = ((end - start) / params.max_fit_window_width_hz) - .ceil() - .max(1.0) as usize; - let width = (end - start) / pieces as f64; - for index in 0..pieces { - let core = ( - start + index as f64 * width, - start + (index + 1) as f64 * width, - ); - let padding = width * params.padding_fraction.max(0.0) * 0.5; - let padded = ( - (core.0 - padding).max(-half_sw), - (core.1 + padding).min(half_sw), - ); - let padded_width = padded.1 - padded.0; - let filter_input = - available_points.min(6_000_usize.saturating_add(params.filter_taps)); - let mut decimation = (data.spectral_width_hz - / (2.0 * padded_width).max(f64::MIN_POSITIVE)) - .floor() - .max(1.0) as usize; - if params.max_downsampled_points > 0 { - decimation = - decimation.max(filter_input.div_ceil(params.max_downsampled_points)); - } - let retained = filter_input - .saturating_sub(actual_filter_taps) - .min(6_000) - .div_ceil(decimation); - fit_windows.push(CraftDerivedWindow { - region: selection.id, - core_hz: core, - padded_hz: padded, - planned_decimation: decimation, - planned_retained_samples: retained, - }); + let filter_input = available_points.min(fit_points.saturating_add(params.fir_filter_taps)); + let clear_signals = detect_clear_signals(data, reference, effective_skip_points); + for window in + build_modeling_windows(data, params, reference, &clear_signals).unwrap_or_default() + { + let modeled_bandwidth_hz = window.modeling_band_hz.1 - window.modeling_band_hz.0; + let mut decimation = (data.spectral_width_hz + / (2.0 * modeled_bandwidth_hz).max(f64::MIN_POSITIVE)) + .floor() + .max(1.0) as usize; + if params.maximum_modeled_sample_count > 0 { + decimation = + decimation.max(filter_input.div_ceil(params.maximum_modeled_sample_count)); } + let retained = filter_input + .saturating_sub(effective_fir_filter_taps) + .min(fit_points) + .div_ceil(decimation); + modeling_windows.push(CraftDerivedModelingWindow { + retention_band_hz: window.retention_band_hz, + modeling_band_hz: window.modeling_band_hz, + planned_decimation_factor: decimation, + planned_modeled_sample_count: retained, + planned_modeled_duration_s: retained as f64 * decimation as f64 + / data.spectral_width_hz, + }); } } CraftDerivedPlan { effective_skip_points, effective_skip_source, available_points, - actual_filter_taps, + effective_fir_filter_taps, reconstruction_points, resolved_regions, - fit_windows, + modeling_windows, } } diff --git a/crates/processing/src/craft/stability.rs b/crates/processing/src/craft/stability.rs new file mode 100644 index 0000000..f006e36 --- /dev/null +++ b/crates/processing/src/craft/stability.rs @@ -0,0 +1,193 @@ +use plotx_io::NmrData; + +use super::diagnostics::{ + CraftModelingWindowDiagnostic, CraftStabilityDiagnostics, CraftStabilityMetric, + CraftStabilityRegion, +}; +use super::regions::{region_ratio, selections_are_valid, summarize_regions}; +use super::{CraftComponent, CraftComponentId, CraftModelingPolicy, CraftReference, CraftRegion}; + +pub(super) fn stability_diagnostics( + all_components: &[CraftComponent], + selections: &[CraftRegion], + windows: &[CraftModelingWindowDiagnostic], + policy: CraftModelingPolicy, + reference: CraftReference, + data: &NmrData, +) -> CraftStabilityDiagnostics { + let delta_ppm = (0.01_f64).max(8.0 / data.observe_freq_mhz.max(f64::MIN_POSITIVE)); + let mut perturbations = vec![("original".to_owned(), selections.to_vec())]; + for (name, start_delta, end_delta) in [ + ("shift left", -delta_ppm, -delta_ppm), + ("shift right", delta_ppm, delta_ppm), + ("expand", -delta_ppm, delta_ppm), + ("contract", delta_ppm, -delta_ppm), + ] { + perturbations.push(( + name.to_owned(), + selections + .iter() + .map(|region| { + CraftRegion::new( + region.id, + region.start_ppm + start_delta, + region.end_ppm + end_delta, + ) + }) + .collect(), + )); + } + for (index, region) in selections.iter().enumerate() { + for (side, start_delta, end_delta) in [ + ("left edge left", -delta_ppm, 0.0), + ("left edge right", delta_ppm, 0.0), + ("right edge left", 0.0, -delta_ppm), + ("right edge right", 0.0, delta_ppm), + ] { + let mut moved = selections.to_vec(); + moved[index] = CraftRegion::new( + region.id, + region.start_ppm + start_delta, + region.end_ppm + end_delta, + ); + perturbations.push((format!("region {} {side}", region.id.0), moved)); + } + } + + let carrier = reference.effective_carrier_ppm(); + let half_ppm = data.spectral_width_hz / (2.0 * data.observe_freq_mhz); + let lower = carrier - half_ppm; + let upper = carrier + half_ppm; + let mut skipped = Vec::new(); + let mut observations = Vec::new(); + for (name, regions) in perturbations { + if !selections_are_valid(®ions) + || regions.iter().any(|region| { + let region = region.normalized(); + region.start_ppm < lower || region.end_ppm > upper + }) + { + skipped.push(format!("{name}: invalid or overlapping regions")); + continue; + } + let assigned = components_for_regions(all_components, ®ions); + let summaries = summarize_regions(&assigned, ®ions); + let ratio = region_ratio(&summaries).map(|ratio| ratio.value); + observations.push((summaries, ratio)); + } + + let total_model_order = windows + .iter() + .map(|window| window.selected_model_order) + .sum(); + let regions = selections + .iter() + .map(|selection| { + let summaries = observations + .iter() + .filter_map(|(summaries, _)| { + summaries + .iter() + .find(|summary| summary.region == selection.id) + }) + .collect::>(); + let values = summaries + .iter() + .map(|summary| summary.coherent_amplitude_t0) + .collect::>(); + CraftStabilityRegion { + region: selection.id, + metric: stability_metric(&values), + component_count_min: summaries + .iter() + .map(|summary| summary.component_count) + .min() + .unwrap_or(0), + component_count_max: summaries + .iter() + .map(|summary| summary.component_count) + .max() + .unwrap_or(0), + model_order_min: total_model_order, + model_order_max: total_model_order, + } + }) + .collect::>(); + let ratio_values = observations + .iter() + .filter_map(|(_, ratio)| *ratio) + .collect::>(); + let ratio = (selections.len() == 2).then(|| stability_metric(&ratio_values)); + let passed = !all_components.is_empty() + && !observations.is_empty() + && regions.iter().all(|region| { + region.metric.relative_dispersion <= policy.boundary_stability_relative_tolerance + && region.component_count_min == region.component_count_max + }) + && ratio.as_ref().is_none_or(|metric| { + ratio_values.len() == observations.len() + && metric.relative_dispersion <= policy.boundary_stability_relative_tolerance + }); + CraftStabilityDiagnostics { + delta_ppm, + regions, + ratio, + passed, + skipped, + } +} + +pub(super) fn components_for_regions( + components: &[CraftComponent], + regions: &[CraftRegion], +) -> Vec { + components + .iter() + .filter_map(|component| { + regions + .iter() + .find(|region| { + let region = region.normalized(); + component.chemical_shift_ppm >= region.start_ppm + && component.chemical_shift_ppm <= region.end_ppm + }) + .map(|region| { + let mut selected = component.clone(); + selected.region = region.id; + selected + }) + }) + .enumerate() + .map(|(id, mut component)| { + component.id = CraftComponentId(id as u64); + component + }) + .collect() +} + +fn stability_metric(values: &[f64]) -> CraftStabilityMetric { + if values.is_empty() { + return CraftStabilityMetric { + median: 0.0, + minimum: 0.0, + maximum: 0.0, + relative_dispersion: f64::MAX, + }; + } + let mut sorted = values.to_vec(); + sorted.sort_by(f64::total_cmp); + let middle = sorted.len() / 2; + let median = if sorted.len().is_multiple_of(2) { + (sorted[middle - 1] + sorted[middle]) * 0.5 + } else { + sorted[middle] + }; + let minimum = sorted[0]; + let maximum = sorted[sorted.len() - 1]; + CraftStabilityMetric { + median, + minimum, + maximum, + relative_dispersion: (maximum - minimum) / median.abs().max(f64::MIN_POSITIVE), + } +} diff --git a/crates/processing/src/craft_tests.rs b/crates/processing/src/craft_tests.rs index 03eff81..c735e03 100644 --- a/crates/processing/src/craft_tests.rs +++ b/crates/processing/src/craft_tests.rs @@ -1,4 +1,5 @@ use super::*; +use std::f64::consts::{PI, TAU}; fn data(components: &[(f64, f64, f64, f64)], count: usize, sw: f64) -> NmrData { let points = (0..count) @@ -35,8 +36,7 @@ fn complete_reduction_recovers_table_and_residual() { 2000.0, ); let params = CraftParams { - max_fit_window_width_hz: 2_000.0, - filter_taps: 127, + fir_filter_taps: 127, ..CraftParams::default() }; let invocation = CraftInvocation::acquisition(&input, params); @@ -44,18 +44,25 @@ fn complete_reduction_recovers_table_and_residual() { assert_eq!(result.components.len(), 2, "{:?}", result.components); assert!((result.components[0].frequency_hz + 75.0).abs() < 0.05); assert!((result.components[1].frequency_hz - 120.0).abs() < 0.05); - assert!(result.diagnostics.normalized_residual < 1e-4); + assert!((result.components[0].amplitude_t0 - 8.0).abs() / 8.0 < 0.01); + assert!((result.components[1].amplitude_t0 - 4.0).abs() / 4.0 < 0.01); + assert!( + result.diagnostics.normalized_residual < 0.01, + "components={:?} diagnostics={:?}", + result.components, + result.diagnostics + ); } #[test] -fn full_band_fit_treats_empty_internal_windows_as_valid_no_signal_results() { +fn full_band_modeling_treats_empty_windows_as_valid_no_signal_results() { let input = data( &[(-300.0, 4.0, 0.1, 2.0), (300.0, 2.0, -0.2, 2.0)], 4096, 2_000.0, ); let params = CraftParams { - filter_taps: 63, + fir_filter_taps: 63, ..CraftParams::default() }; @@ -69,7 +76,7 @@ fn full_band_fit_treats_empty_internal_windows_as_valid_no_signal_results() { assert_eq!(result.region_summaries.len(), 1); assert_eq!(result.region_summaries[0].region, CraftRegionId(0)); assert!(result.region_summaries[0].component_count >= 2); - assert_eq!(result.diagnostics.fit_windows.len(), 4); + assert_eq!(result.diagnostics.modeling_windows.len(), 2); } #[test] @@ -80,14 +87,18 @@ fn ssfp_skip_extrapolates_amplitude_to_time_zero() { profile: CraftProfile::Ssfp, skip_duration_s: 10.0 / input.spectral_width_hz, reconstruction_duration_s: Some(0.1), - max_fit_window_width_hz: 20_000.0, - filter_taps: 63, + fir_filter_taps: 63, ..CraftParams::default() }; let invocation = CraftInvocation::acquisition(&input, params); let result = process_craft_cancellable(&input, &invocation, &|| false).unwrap(); - assert_eq!(result.components.len(), 1); - assert!((result.components[0].amplitude_t0 - 5.0).abs() < 0.05); + assert_eq!(result.components.len(), 1, "{:?}", result.components); + assert!( + (result.components[0].amplitude_t0 - 5.0).abs() < 0.05, + "components={:?} diagnostics={:?}", + result.components, + result.diagnostics + ); assert_eq!(result.synthetic_fid.len(), 2000); assert!( result @@ -117,7 +128,7 @@ fn group_delay_uses_the_physical_fid_time_origin() { .collect(); let params = CraftParams { regions: vec![CraftRegion::new(CraftRegionId(1), 0.18, 0.22)], - filter_taps: 127, + fir_filter_taps: 127, ..CraftParams::default() }; @@ -149,10 +160,10 @@ fn default_model_limit_resolves_non_lorentzian_multiplet_quantitation() { let input = data(&components, 4096, 2_000.0); let params = CraftParams { regions: vec![ - CraftRegion::new(CraftRegionId(0), -0.23, -0.17), - CraftRegion::new(CraftRegionId(1), 0.14, 0.23), + CraftRegion::new(CraftRegionId(0), -0.25, -0.15), + CraftRegion::new(CraftRegionId(1), 0.12, 0.25), ], - filter_taps: 127, + fir_filter_taps: 127, ..CraftParams::default() }; @@ -163,9 +174,29 @@ fn default_model_limit_resolves_non_lorentzian_multiplet_quantitation() { ) .unwrap(); - assert!(result.diagnostics.fit_windows[1].selected_model_order > 7); + assert_eq!(result.region_summaries[1].component_count, 8); + assert!( + !result + .diagnostics + .warnings + .iter() + .any(|warning| warning.kind == CraftWarningKind::InputAssessment + && warning.message.contains("peak density")) + ); let ratio = result.region_ratio.unwrap().value; assert!((ratio - 1.5).abs() < 0.03, "ratio was {ratio}"); + assert!( + result.diagnostics.stability.passed, + "{:?}", + result.diagnostics.stability + ); + assert!( + result + .diagnostics + .stability + .ratio + .is_some_and(|metric| metric.relative_dispersion < 0.01) + ); } #[test] @@ -190,12 +221,11 @@ fn overlapping_requested_regions_are_rejected_as_ambiguous() { CraftRegion::new(CraftRegionId(10), -1.0, 1.0), CraftRegion::new(CraftRegionId(20), 0.5, 2.0), ], - max_fit_window_width_hz: 500.0, ..CraftParams::default() }; assert!(matches!( - build_regions(&input, ¶ms, CraftReference::acquisition(&input)), + build_modeling_windows(&input, ¶ms, CraftReference::acquisition(&input), &[],), Err(CraftError::InvalidParameters) )); } @@ -206,14 +236,24 @@ fn reference_maps_displayed_regions_and_reported_shifts_without_changing_frequen let reference = CraftReference::new(input.carrier_ppm, 0.15); let params = CraftParams { regions: vec![CraftRegion::new(CraftRegionId(7), 0.38, 0.40)], - max_fit_window_width_hz: 500.0, - filter_taps: 127, + fir_filter_taps: 127, ..CraftParams::default() }; - let regions = build_regions(&input, ¶ms, reference).unwrap(); - assert!((regions[0].core.0 - 115.0).abs() < 1e-9); - assert!((regions[0].core.1 - 125.0).abs() < 1e-9); + let clear_signals = preflight::detect_clear_signals(&input, reference, 0); + let regions = build_modeling_windows(&input, ¶ms, reference, &clear_signals).unwrap(); + assert_eq!(regions.len(), 1); + assert!( + ((regions[0].retention_band_hz.0 + regions[0].retention_band_hz.1) * 0.5 - 120.0).abs() + < 0.25, + "signals={clear_signals:?} window={:?}", + regions[0] + ); + assert!( + regions + .iter() + .all(|window| window.retention_band_hz.1 - window.retention_band_hz.0 <= 500.0) + ); let invocation = resolve_craft_invocation( &input, @@ -230,11 +270,10 @@ fn reference_maps_displayed_regions_and_reported_shifts_without_changing_frequen } #[test] -fn narrow_window_fit_is_independent_of_signal_phase() { +fn modeling_window_is_independent_of_signal_phase() { let params = CraftParams { regions: vec![CraftRegion::new(CraftRegionId(1), 0.18, 0.22)], - max_fit_window_width_hz: 500.0, - filter_taps: 127, + fir_filter_taps: 127, ..CraftParams::default() }; let fitted_frequency = |phase| { @@ -252,7 +291,7 @@ fn narrow_window_fit_is_independent_of_signal_phase() { } #[test] -fn fit_windows_preserve_user_region_identity_and_one_region_ratio() { +fn modeling_windows_are_independent_while_components_preserve_region_identity() { let input = data( &[(-300.0, 4.0, 0.1, 2.0), (300.0, 2.0, 0.1, 2.0)], 4096, @@ -263,24 +302,19 @@ fn fit_windows_preserve_user_region_identity_and_one_region_ratio() { CraftRegion::new(CraftRegionId(22), 0.1, 1.0), CraftRegion::new(CraftRegionId(11), -1.0, -0.1), ], - max_fit_window_width_hz: 150.0, - max_components_per_fit_window: 3, - filter_taps: 63, + maximum_model_order: 3, + fir_filter_taps: 63, ..CraftParams::default() }; let reference = CraftReference::acquisition(&input); - let windows = build_regions(&input, ¶ms, reference).unwrap(); - assert_eq!(windows.len(), 6); + let clear_signals = preflight::detect_clear_signals(&input, reference, 0); + let windows = build_modeling_windows(&input, ¶ms, reference, &clear_signals).unwrap(); + assert_eq!(windows.len(), 2); assert!( - windows[..3] + windows .iter() - .all(|window| window.selection.id == CraftRegionId(11)) - ); - assert!( - windows[3..] - .iter() - .all(|window| window.selection.id == CraftRegionId(22)) + .all(|window| window.retention_band_hz.1 - window.retention_band_hz.0 <= 500.0) ); let invocation = resolve_craft_invocation( @@ -302,6 +336,55 @@ fn fit_windows_preserve_user_region_identity_and_one_region_ratio() { assert!((result.region_ratio.unwrap().value - 0.5).abs() < 0.05); } +#[test] +fn overlapping_modeling_bands_retain_each_component_once() { + let input = data( + &[ + (-100.0, 10.0, 0.1, 2.0), + (0.0, 2.0, 0.1, 2.0), + (110.0, 8.0, 0.1, 2.0), + ], + 4096, + 2_000.0, + ); + let params = CraftParams { + maximum_model_order: 4, + fir_filter_taps: 127, + ..CraftParams::default() + }; + + let result = process_craft_cancellable( + &input, + &CraftInvocation::acquisition(&input, params), + &|| false, + ) + .unwrap(); + + assert_eq!(result.diagnostics.modeling_windows.len(), 2); + assert_eq!( + result + .diagnostics + .modeling_windows + .iter() + .map(|window| window.selected_model_order) + .sum::(), + 4, + "{:?}", + result.diagnostics.modeling_windows + ); + assert_eq!(result.components.len(), 3, "{:?}", result.components); + assert_eq!( + result + .components + .iter() + .filter(|component| component.frequency_hz.abs() < 1.0) + .count(), + 1, + "{:?}", + result.components + ); +} + #[test] fn rejects_non_finite_reference() { let input = data(&[], 128, 1_000.0); @@ -386,11 +469,11 @@ fn zero_overrides_resolve_complete_bandwidth_and_stable_sources() { fn resolver_applies_per_field_explicit_provenance_default_priority() { let input = data(&[(80.0, 10.0, 0.0, 2.0)], 1024, 1_000.0); let mut prior_params = CraftParams::ssfp(); - prior_params.min_amplitude_to_noise = 8.0; - prior_params.filter_taps = 127; + prior_params.minimum_amplitude_to_noise = 8.0; + prior_params.fir_filter_taps = 127; let prior = CraftInvocation::acquisition(&input, prior_params); let overrides = CraftParamOverrides { - min_amplitude_to_noise: Some(5.0), + minimum_amplitude_to_noise: Some(5.0), ..CraftParamOverrides::default() }; @@ -401,14 +484,14 @@ fn resolver_applies_per_field_explicit_provenance_default_priority() { Some(&prior), ); - assert_eq!(resolved.params.min_amplitude_to_noise, 5.0); + assert_eq!(resolved.params.minimum_amplitude_to_noise, 5.0); assert_eq!( - resolved.sources.min_amplitude_to_noise, + resolved.sources.minimum_amplitude_to_noise, CraftParamSource::ExplicitInput ); - assert_eq!(resolved.params.filter_taps, 127); + assert_eq!(resolved.params.fir_filter_taps, 127); assert_eq!( - resolved.sources.filter_taps, + resolved.sources.fir_filter_taps, CraftParamSource::ResultProvenance ); assert_eq!(resolved.params.profile, CraftProfile::Ssfp); @@ -419,7 +502,7 @@ fn selecting_profile_clears_profile_owned_overrides_but_keeps_regions() { let region = CraftRegion::new(CraftRegionId(9), -0.2, 0.2); let mut overrides = CraftParamOverrides { regions: Some(vec![region]), - min_amplitude_to_noise: Some(9.0), + minimum_amplitude_to_noise: Some(9.0), ..CraftParamOverrides::default() }; @@ -427,23 +510,18 @@ fn selecting_profile_clears_profile_owned_overrides_but_keeps_regions() { assert_eq!(overrides.profile, Some(CraftProfile::Ssfp)); assert_eq!(overrides.regions, Some(vec![region])); - assert_eq!(overrides.min_amplitude_to_noise, None); + assert_eq!(overrides.minimum_amplitude_to_noise, None); let input = data(&[(80.0, 10.0, 0.0, 2.0)], 1024, 1_000.0); - let mut conventional = CraftParams::conventional(); - conventional.max_fit_window_width_hz = 125.0; - let previous = CraftInvocation::acquisition(&input, conventional); + let previous = CraftInvocation::acquisition(&input, CraftParams::conventional()); let resolved = resolve_craft_invocation( &input, CraftReference::acquisition(&input), &overrides, Some(&previous), ); - assert_eq!(resolved.params.max_fit_window_width_hz, 2_000.0); - assert_eq!( - resolved.sources.max_fit_window_width_hz, - CraftParamSource::StableDefault - ); + assert_eq!(resolved.params.profile, CraftProfile::Ssfp); + assert_eq!(resolved.modeling_policy.modeling_bandwidth_hz, 2_000.0); } #[test] @@ -483,7 +561,7 @@ fn short_and_invalid_inputs_are_classified_before_execution() { } #[test] -fn derived_plan_matches_actual_fit_window_diagnostics() { +fn derived_plan_matches_actual_modeling_window_diagnostics() { let input = data(&[(100.0, 10.0, 0.1, 2.0)], 2048, 1_000.0); let invocation = resolve_craft_invocation( &input, @@ -494,18 +572,174 @@ fn derived_plan_matches_actual_fit_window_diagnostics() { let result = process_craft_cancellable(&input, &invocation, &|| false).unwrap(); assert_eq!( - result.diagnostics.fit_windows.len(), - invocation.derived_plan.fit_windows.len() + result.diagnostics.modeling_windows.len(), + invocation.derived_plan.modeling_windows.len() ); for (actual, planned) in result .diagnostics - .fit_windows + .modeling_windows + .iter() + .zip(&invocation.derived_plan.modeling_windows) + { + assert_eq!(actual.retention_band_hz, planned.retention_band_hz); + assert_eq!(actual.modeling_band_hz, planned.modeling_band_hz); + assert_eq!(actual.decimation_factor, planned.planned_decimation_factor); + } +} + +#[test] +fn user_boundaries_do_not_change_fixed_modeling_protocol() { + let input = data( + &[(-100.0, 6.0, 0.2, 2.0), (100.0, 4.0, 0.2, 2.0)], + 4096, + 2_000.0, + ); + let invocation = |regions| { + let params = CraftParams { + regions, + fir_filter_taps: 127, + ..CraftParams::default() + }; + CraftInvocation::acquisition(&input, params) + }; + let narrow = invocation(vec![CraftRegion::new(CraftRegionId(1), -0.24, -0.16)]); + let wide = invocation(vec![CraftRegion::new(CraftRegionId(1), -0.30, -0.10)]); + + assert_eq!(narrow.modeling_policy, wide.modeling_policy); + assert_eq!(narrow.modeling_policy.modeling_bandwidth_hz, 250.0); + assert_eq!(narrow.modeling_policy.modeling_duration_s, 1.0); + assert_eq!( + narrow.derived_plan.modeling_windows.len(), + wide.derived_plan.modeling_windows.len() + ); + for (left, right) in narrow + .derived_plan + .modeling_windows .iter() - .zip(&invocation.derived_plan.fit_windows) + .zip(&wide.derived_plan.modeling_windows) { - assert_eq!(actual.region, planned.region); - assert_eq!(actual.core_hz, planned.core_hz); - assert_eq!(actual.padded_hz, planned.padded_hz); - assert_eq!(actual.actual_decimation, planned.planned_decimation); + assert_eq!(left.retention_band_hz, right.retention_band_hz); + assert_eq!(left.modeling_band_hz, right.modeling_band_hz); + assert_eq!( + left.planned_decimation_factor, + right.planned_decimation_factor + ); + assert_eq!( + left.planned_modeled_sample_count, + right.planned_modeled_sample_count + ); + assert_eq!( + left.planned_modeled_duration_s, + right.planned_modeled_duration_s + ); + } +} + +#[test] +fn boundary_instability_marks_run_partial_but_keeps_components() { + let input = data(&[(100.0, 5.0, 0.2, 2.0)], 4096, 2_000.0); + let params = CraftParams { + regions: vec![CraftRegion::new(CraftRegionId(7), 0.195, 0.30)], + fir_filter_taps: 127, + ..CraftParams::default() + }; + + let result = process_craft_cancellable( + &input, + &CraftInvocation::acquisition(&input, params), + &|| false, + ) + .unwrap(); + + assert_eq!(result.components.len(), 1); + assert_eq!(result.diagnostics.status, CraftRunStatus::Partial); + assert!(!result.diagnostics.stability.passed); + assert!( + result + .diagnostics + .warnings + .iter() + .any(|warning| { warning.kind == CraftWarningKind::StabilityFailure }) + ); +} + +#[test] +fn global_zero_order_phase_does_not_change_coherent_amplitude() { + let input = data( + &[(-20.0, 3.0, 0.1, 1.5), (20.0, 2.0, 0.4, 1.8)], + 4096, + 2_000.0, + ); + let params = CraftParams { + regions: vec![CraftRegion::new(CraftRegionId(3), -0.10, 0.10)], + fir_filter_taps: 127, + ..CraftParams::default() + }; + let fit = |input: &NmrData| { + process_craft_cancellable( + input, + &CraftInvocation::acquisition(input, params.clone()), + &|| false, + ) + .unwrap() + .region_summaries[0] + .coherent_amplitude_t0 + }; + let expected = fit(&input); + let rotation = Complex64::from_polar(1.0, 1.1); + let mut rotated = input.clone(); + for point in &mut rotated.points { + *point *= rotation; } + + let actual = fit(&rotated); + assert!( + (actual - expected).abs() / expected < 1e-6, + "expected={expected} actual={actual}" + ); +} + +#[test] +fn no_clear_signal_allows_exploration_but_requires_review() { + let input = data(&[], 4096, 2_000.0); + let invocation = CraftInvocation::acquisition( + &input, + CraftParams { + fir_filter_taps: 63, + ..CraftParams::default() + }, + ); + + assert!(invocation.assessment.can_run()); + assert!(invocation.assessment.issues.iter().any(|issue| { + issue.code == CraftIssueCode::NoClearSignal && issue.severity == CraftIssueSeverity::Warning + })); + let result = process_craft_cancellable(&input, &invocation, &|| false).unwrap(); + + assert_eq!(result.diagnostics.status, CraftRunStatus::Partial); + assert!(result.components.is_empty()); + assert!(!result.diagnostics.stability.passed); +} + +#[test] +fn validation_selection_keeps_model_order_when_only_the_unused_tail_changes() { + let components = [(-30.0, 5.0, 0.2, 1.5), (42.0, 3.0, -0.3, 2.0)]; + let selected_orders = [4096, 4112].map(|count| { + let input = data(&components, count, 2_000.0); + let invocation = CraftInvocation::acquisition( + &input, + CraftParams { + fir_filter_taps: 127, + ..CraftParams::default() + }, + ); + process_craft_cancellable(&input, &invocation, &|| false) + .unwrap() + .diagnostics + .modeling_windows[0] + .selected_model_order + }); + + assert_eq!(selected_orders[0], 2); + assert_eq!(selected_orders[1], selected_orders[0]); } diff --git a/docs/src/content/docs/guides/craft.md b/docs/src/content/docs/guides/craft.md index 795369b..865b44b 100644 --- a/docs/src/content/docs/guides/craft.md +++ b/docs/src/content/docs/guides/craft.md @@ -50,27 +50,38 @@ fields under **Signal groups** are available when you need exact bounds. A signal group can contain several fitted components. A component is a fitted resonance contribution, not a compound identification or a guaranteed visible -multiplet. A wide group may be split into several calculation windows, but its -components are still reported under the group you selected. +multiplet. Fixed modeling windows determine how the FID is solved; the signal +group boundary only determines which completed components belong to the group. -## Advanced fit settings +## Advanced component settings -Leave **Advanced fit settings** closed for routine work. The conventional +Leave **Advanced component settings** closed for routine work. The conventional profile uses these defaults: - **Minimum A/N**: 3.3. Lower values retain weaker candidates but increase the chance of fitting noise; values below 3.3 are flagged for review. -- **Max components / fit window**: 15 (allowed range 1–64). Reaching the limit - is reported as a diagnostic warning. -- **Linewidth range (Hz)**: 0.05–10 Hz. -- **Fit window width (Hz)**: 500 Hz. This controls how wide groups are divided - for calculation; it does not create extra signal groups. +- **Maximum model order**: 15 (allowed range 1–64) for each modeling window. + Reaching the limit is reported as a diagnostic warning. +- **Component linewidth range (Hz)**: 0.05–20 Hz. The bound applies to each + component, not the frequency range modeled at once. The 20 Hz default is a + typical starting point rather than a universal constant. Change it only when + the acquisition and expected line shape justify a different bound. + +The fixed modeling bandwidth is 250 Hz for Conventional and 2000 Hz for SSFP. +This is the actual frequency width of a modeling window, not a linewidth and +not a quantitative tuning control. A modeling window can contain many component +lines, each still constrained by the separate component-linewidth range. Use **Reset** beside an edited value to restore the value inherited from the selected run or the profile default. Changing profiles keeps the selected groups and loads the new profile's settings. Conventional FID is always the default; PlotX does not infer SSFP from the waveform. +CRAFT keeps Bruker acquisition `GRPDLY` separate from its own FIR filtering. +The importer/FFT path uses `GRPDLY` to define the acquisition time origin; the +499-tap CRAFT FIR has an independent edge transient handled by phase-conjugate +precharge. Neither delay is silently folded into the other. + The **SSFP / interrupted FID** profile starts with **Skip initial** at 0.5 ms and **Extend reconstructed FID** enabled for 1.2 s. Skipping early points can remove fast-decaying background; reconstruction extends the modeled FID for @@ -106,11 +117,29 @@ or enabled **Reference** step changes; rerun it before interpreting the result. Runs with warnings or a partial fit are marked **Needs review** even when the calculation completes. +CRAFT performs deterministic boundary-perturbation checks around each selected +group. Small shifts, expansions, contractions, and one-sided moves must keep +amplitudes and ratios within 1%. A run that fails this stability gate keeps its +complete component table and residual for inspection, but cannot create or +export a quantitative amplitude report. + Choose **Export components…** to open the standard CSV, TSV, XLSX, or clipboard export dialog. Under **Signals**, **Create data table** creates a sortable PlotX table without leaving CRAFT. Choose **View data table** to inspect, chart, or export it, and **Add to board** only when you want the table on a board sheet. +## CRAFT amplitude reports + +The model's **Minimum A/N** is a trust criterion for retaining fitted components. +The **Reports** tab is a separate reporting layer: its **Report threshold** +selects which retained components are included, without refitting or changing +the complete component table. **Segment width** is the total fixed frequency +window (Hz) around each selected peak; overlapping windows are merged. Reports +show both the scalar sum of component amplitudes and the phase-aware coherent +amplitude. These segment amplitudes are summaries of fitted components, not +integrals of frequency-domain bins. A report whose source run changes is marked +for review rather than silently recalculated. + ## Interpretation CRAFT components describe the selected FID; they do not identify compounds or diff --git a/docs/src/content/docs/zh-cn/guides/craft.md b/docs/src/content/docs/zh-cn/guides/craft.md index 3e4be97..a0f6f85 100644 --- a/docs/src/content/docs/zh-cn/guides/craft.md +++ b/docs/src/content/docs/zh-cn/guides/craft.md @@ -37,25 +37,33 @@ CRAFT 直接拟合一维 NMR 采集中的原始复数 FID，并报告共振分 时每次移动十个）。需要精确边界时，可直接编辑 **Signal groups** 下的数值字段。 一个信号组可以包含多个拟合分量。分量是拟合得到的共振贡献，不是化合物鉴定，也不 -保证对应一个肉眼可见的多重峰。较宽的信号组可能被分成多个计算窗口，但分量仍会 -归在你选择的信号组下。 +保证对应一个肉眼可见的多重峰。固定建模窗口决定如何求解 FID；信号组边界只决定哪些 +已完成拟合的分量归入该信号组。 -## 高级拟合设置 +## 高级分量设置 -常规分析无需展开 **Advanced fit settings**。Conventional FID 的默认值为： +常规分析无需展开 **Advanced component settings**。Conventional FID 的默认值为： - **Minimum A/N**：3.3。较低值会保留更弱的候选分量，但也更容易拟合噪声；低于 3.3 时会标记为需要复核。 -- **Max components / fit window**：15（允许范围 1–64）。达到上限会在诊断中给出 - 警告。 -- **Linewidth range (Hz)**：0.05–10 Hz。 -- **Fit window width (Hz)**：500 Hz。它只决定宽信号组如何分段计算，不会增加 - 信号组。 +- **Maximum model order**：每个建模窗口默认为 15（允许范围 1–64）。达到上限会在 + 诊断中给出警告。 +- **Component linewidth range (Hz)**：0.05–20 Hz。该范围限制每个分量，而不是 + 一次建模所覆盖的频率范围。20 Hz 默认值是典型起点，并非对所有样品都固定不变； + 只有在采集条件和预期线形有独立依据时才应修改该范围。 + +固定建模带宽在 Conventional 中为 250 Hz，在 SSFP 中为 2000 Hz。这是一个建模窗口 +实际覆盖的频率宽度，不是线宽，也不是用于追逐定量结果的调参旋钮。一个建模窗口可以 +包含多个分量谱线，每个分量仍受独立的分量线宽范围约束。 编辑后的值可用旁边的 **Reset** 恢复为所选运行或配置的默认值。切换配置会保留已选 信号组，并载入新配置的设置。Conventional FID 始终是默认配置；PlotX 不会根据波形 猜测 SSFP。 +CRAFT 将 Bruker 采集的 `GRPDLY` 与自身 FIR 滤波严格分开。导入器/FFT 流程用 +`GRPDLY` 定义采集时间原点；499-tap CRAFT FIR 具有独立的边界瞬态，并通过相位共轭 +预充电处理。两种延迟不会互相重复补偿。 + **SSFP / interrupted FID** 配置的 **Skip initial** 默认值为 0.5 ms，并默认启用 时长 1.2 s 的 **Extend reconstructed FID**。跳过开头的数据点可以排除快速衰减的 背景；重建功能会为该配置延长模型 FID。除非实验已经完成定量验证，否则 SSFP 结果 @@ -83,11 +91,24 @@ CRAFT 直接拟合一维 NMR 采集中的原始复数 FID，并报告共振分 运行会标记为 **Stale**；请重新运行后再解读。即使计算完成，只要有警告或拟合不完整， 运行也会标记为 **Needs review**。 +CRAFT 会对每个所选信号组执行确定性的边界扰动检查，包括整体平移、两侧扩展或收缩以及 +单侧移动。振幅和比值的相对离散度必须保持在 1% 以内。未通过稳定性门槛的运行仍会保存 +完整分量表和残差供检查，但不能创建或导出可靠的定量振幅报告。 + 选择 **Export components…** 可打开标准 CSV、TSV、XLSX 或剪贴板导出对话框。在 **Signals** 中选择 **Create data table**，可在 CRAFT 卡片内创建可排序的 PlotX 表格。 创建后选择 **View data table** 查看、绘图或导出；只有需要把表格放到 board 的 sheet 上时，才选择 **Add to board**。 +## CRAFT 振幅报告 + +建模设置中的 **Minimum A/N** 是判断分量是否可信的标准。**Reports** 标签是独立的 +报告层；其中的 **Report threshold** 只决定哪些已保留分量进入报告，不会重新拟合，也 +不会改变完整分量表。**Segment width** 是围绕每个入选峰的固定总频宽（Hz）；相互重叠 +的窗口会合并。报告同时显示分量振幅的直接相加值和考虑相位的相干振幅。这些值是拟合 +分量的汇总，不是频域 bin 积分。来源运行发生变化时，报告会标记为需要复核，而不会 +静默地重新计算。 + ## 解释结果 CRAFT 分量描述的是所选 FID，不会鉴定化合物，也不能替代浓度校准。进行定量指纹分析