Skip to content

research: noise recovery notebooks with physical ground truth#249

Merged
KedoKudo merged 12 commits intofeature/172-2d-imagingfrom
feature/172-noise-recovery-research
Feb 26, 2026
Merged

research: noise recovery notebooks with physical ground truth#249
KedoKudo merged 12 commits intofeature/172-2d-imagingfrom
feature/172-noise-recovery-research

Conversation

@KedoKudo
Copy link
Collaborator

@KedoKudo KedoKudo commented Feb 26, 2026

Summary

  • Rebuilt noise recovery research notebooks (NB0 + NB1) from scratch using physical ground truth from sample metadata instead of the flawed NNLS-on-clean reference
  • NB0 establishes ground truth, quantifies the NNLS normalization artifact (~63–218K ppm crosstalk), and characterizes the sparse data challenge (clipping, binning trade-offs)
  • NB1 compares 8 recovery methods against physical truth: per-pixel NNLS, T-domain fitting, spatial denoising (NLM/Wavelet/Gaussian), NMF neighborhoods, and hierarchical subdivision
  • All conclusions verified against executed output with deterministic seed

Key findings for March 4th presentation:

  • Wavelet+NNLS is best at moderate noise (n=10): 181K ppm avg MAE
  • Hierarchical bin=2 is the most robust method overall: 185K at n=10, 296K at n=2 (only method to beat raw NNLS at extreme noise)
  • Blind NMF does not outperform informed NNLS — motivates semi-supervised approaches
  • NNLS model bias floor (~100–335K ppm) limits all preprocessing methods

Test plan

  • pixi run jupyter nbconvert --execute --to notebook examples/Notebooks/pleiades_noise_recovery_nb0_ground_truth.ipynb runs without error
  • pixi run jupyter nbconvert --execute --to notebook examples/Notebooks/pleiades_noise_recovery_nb1_methods.ipynb runs without error
  • Verify truth values: U-235=0.9974, Pu-241=0.4863
  • Check _debug_images/NB0_*.png and NB1_*.png for reasonable plots

🤖 Generated with Claude Code

KedoKudo and others added 11 commits February 22, 2026 17:55
Add denoise_nmf parameter to PhysicsRecovery.recover_image() and
analyze_imaging() that applies NMF low-rank denoising before per-pixel
NNLS recovery. This exploits spatial redundancy across all pixels to
filter spectral noise while preserving resonance structure, combining
NMF's noise robustness with NNLS's physical accuracy from SAMMY
reference spectra.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
…ill timed-out slots

Two bugs in _collect_results_with_timeout:

1. The stall_rounds bailout cancelled all pending futures after 2 timeout
   windows without progress, even when no futures had started running yet.
   During slow ProcessPoolExecutor startup, this killed pixels that were
   never attempted. Fix: only increment stall_rounds when at least one
   future is in running_start_times.

2. Timed-out futures were removed from pending but submit_fn() was never
   called to replace the freed slot. With bounded submission, each timeout
   permanently shrank the in-flight window, causing remaining pixels to
   never be submitted. Fix: call submit_fn() after each timed-out future,
   mirroring the existing refill logic for completed futures.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
…lock queue

When all running futures time out, their worker processes keep running
(ProcessPoolExecutor cannot kill in-flight tasks). Replacement futures
stay queued forever because no worker slot is available. The previous
fix (only advance stall_rounds when futures are currently running)
prevented this case from ever bailing out, causing fit_pixels() to hang.

Fix: track whether any worker has ever been observed running. Advance the
stall counter when workers_have_started is True even if no pending future
is currently running — this correctly identifies zombie-occupied pools
while still protecting against spurious cancellation during genuine slow
pool startup (where workers_have_started remains False).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
NB0 (pleiades_noise_recovery_research.ipynb): Documents progressive failure
of NNLS, SG denoising, blind NMF, and Semi-NMF at low photon counts.
Establishes binning sweep, SAMMY strided fitting, and representative pixel
benchmarks as the testing framework.

NB1 (pleiades_noise_recovery_nb1_transmission_domain.ipynb): Tests fitting
in transmission domain to avoid log-clipping catastrophe, plus TV spatial
regularization. Key findings:
- T-domain+TV(w=0.1) achieves MAE 0.108 at n=2 (30% better than raw NNLS)
- T-domain has systematic bias vs NNLS on clean data (MAE 0.036)
- Iterative alternating provides no benefit (converges in 1 iteration)
- At n=2, per-pixel information is fundamentally inadequate; methods that
  borrow across similar spatial regions are needed

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…PM visualization

Add configurable vary flags (normalization, background, TZERO, thickness)
to ImagingConfig and InpDatasetMetadata, allowing users to fix nuisance
parameters for clean synthetic data while keeping them free for real data.
Defaults are True (vary all) for backwards compatibility.

Add units parameter ("fraction", "percent", "ppm") to all
AbundanceMapVisualizer plot methods, enabling instrument scientists to
view composition maps in parts-per-million as requested by VENUS team.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…onclusions

Switch all composition metrics from fractional (0-1) to ppm scale. Normalize
SAMMY validation tables. Add verified conclusions based on actual notebook
output: T-domain raw is the clear win for visualization (19-43% at n=10-50),
TV regularization is questionable due to blurring, and SAMMY is robust to
noise — preprocessing hurts characterization by erasing resonance fine structure.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
NB2 tests Gaussian, Non-Local Means, and Wavelet BayesShrink spatial
denoising on the raw transmission cube before T-domain fitting. All
three methods reduce composition map MAE by 27-29% at n=2. NLM wins
composition maps (59K ppm at n=10, +27%), Wavelet wins SAMMY fidelity
(mean |δ|=116K ppm). T-domain reconstruction universally hurts SAMMY
regardless of prior denoising — confirmed across all experiments.

Adds pywavelets dependency required by skimage wavelet denoising.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Replace 3 research notebooks (NB0-baseline, NB1-tdomain, NB2-spatial)
with 2 consolidated notebooks using correct physical ground truth from
sample metadata instead of NNLS-on-clean reference.

NB0 (ground truth): Derives truth from calculate_number_density(),
discovers 4 spectral classes, quantifies NNLS normalization artifact
(~250K ppm crosstalk), characterizes clipping catastrophe, and maps
the accuracy-resolution trade-off via binning sweep.

NB1 (methods): Compares T-domain fitting, spatial denoising (NLM,
wavelet, Gaussian), NMF neighborhood prototype, and hierarchical
subdivision against physical truth. Includes honest SAMMY assessment.
Key finding: blind NMF does not beat informed NNLS; spatial denoising
helps ~10-15% at moderate noise; hierarchical bin=2 is optimal.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Tightened all observation and conclusion cells against executed output:
- NB0 §3: precise crosstalk range (63-218K ppm)
- NB0 §5: per-class gate check numbers, model bias range (100-335K ppm)
- NB0 §6: binning only helps at n=5-10, not n>=50
- NB1 §2: exact spatial denoising improvement percentages
- NB1 §7: re-ranked methods table by n=10 avg MAE (Wavelet 1st, Hier 2nd)
- NB1 recap: updated to match verified NB0 findings

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
@KedoKudo
Copy link
Collaborator Author

Note: Poisson noise should be added to counts, not to transmission

The noise simulation in these notebooks adds Poisson noise directly to the transmission value T:

sigma = np.sqrt(T / n_incident)
T_noisy = T + noise  # or similar

This is physically incorrect. Transmission is not the Poisson random variable — raw counts are.

The correct simulation is:

counts = rng.poisson(n_incident * T_true)   # counts ~ Poisson(I₀ · T)
T_noisy = counts / n_incident               # normalize back to [0, 1]
sigma   = np.sqrt(np.maximum(counts, 1)) / n_incident  # propagated uncertainty

Why it matters:

  • Adding noise to T can produce values outside [0, 1] at low I₀ (the "clipping catastrophe" you document in NB0). This is a symptom of the wrong approach, not an inherent limitation.
  • Generating noise in count space then normalizing keeps T_noisy physical by construction, even at I₀ = 1.
  • At deep resonance dips where T → 0, the correct formulation gives counts ≈ 0 (a few Poisson-distributed zeros), whereas the T-domain approach blows up the relative uncertainty.

The distinction also matters for fitting: the correct sigma from the count-domain simulation feeds directly into a Gaussian chi-squared fitter (valid when counts >> 1) or a Poisson NLL fitter (valid at all count levels).

For reference, NEREIDS uses the count-domain simulation for its noise sensitivity analysis in examples/notebooks/applications/01_spatial_mapping_demo.ipynb, and is adding a TRINIDI-style Poisson NLL fitter (fitter='poisson' in spatial_map()) for the sparse-counts regime.

Address PR #249 comment: the noise simulation already uses the physically
correct count-domain approach (DataDegrader generates Poisson counts then
normalizes to T), but the notebook text was misleading.

Changes:
- NB0 §4: explicitly document count-domain Poisson model in code comments
- NB0 §4: reframe "clipping catastrophe" as "boundary values" — T=0 and
  T=1 are real Poisson outcomes, not simulation artifacts
- NB0: fix column headers and variable names (unphysical → boundary)
- NB1: update compute_poisson_sigma docstring to show equivalence:
  σ_T = √(counts)/I₀ = √(T_obs/I₀)
- Both: add DataDegrader source reference in code comments

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
@KedoKudo KedoKudo self-assigned this Feb 26, 2026
@KedoKudo KedoKudo merged commit 417f9ab into feature/172-2d-imaging Feb 26, 2026
3 checks passed
@KedoKudo KedoKudo deleted the feature/172-noise-recovery-research branch February 26, 2026 20:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant