Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 84 additions & 0 deletions crates/analysis/src/craft.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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::<Vec<_>>();
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::<f64>::zeros(equation_count * 2, order * 2);
let mut observed = DVector::<f64>::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::<Vec<_>>();
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(
Expand Down
50 changes: 50 additions & 0 deletions crates/analysis/src/craft_tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
3 changes: 1 addition & 2 deletions crates/app/src/shot/craft_shot.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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}"))?;
Expand Down
142 changes: 142 additions & 0 deletions crates/app/src/ui/canvas/craft_results.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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::<Vec<_>>();
let regions = stored
.region_summaries
.iter()
.map(|region| (region.start_ppm, region.end_ppm))
.collect::<Vec<_>>();
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::<plotx_processing::craft::CraftAmplitudeReport>(
record.snapshot.clone(),
)
.ok()
})
.map(|report| {
report
.segments
.into_iter()
.map(|segment| {
(
carrier + segment.start_hz / observe,
carrier + segment.end_hz / observe,
)
})
.collect::<Vec<_>>()
})
.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(
&regions,
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);
}
}
}
2 changes: 1 addition & 1 deletion crates/app/src/ui/commands_craft_tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
Loading
Loading