Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #321 +/- ##
==========================================
+ Coverage 92.52% 92.56% +0.03%
==========================================
Files 18 18
Lines 2341 2367 +26
==========================================
+ Hits 2166 2191 +25
- Misses 175 176 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
python 3.11 still in CI when 3.12 is minimum required. |
… estimate noise in a 2D spectrum. This is for now used only in `specreduce.line_matching.find_arc_lines` as a fallback noise estimate if the spectrum doesn't contain uncertainties, but may be useful elsewhere as well.
fa7adef to
c9bc1b9
Compare
|
Shouldn't |
tepickering
left a comment
There was a problem hiding this comment.
The points I raised in specreduce/utils/utils.py are pretty significant and should be dealt with, though it may be an upstream issue as well. Definitely needs a deeper look. Once that's sorted, measure_noise should operate on the masked data. The docs issue is minor, but worth resolving.
As a heads-up: this PR and #320 both rewrite the "no uncertainty" fallback block in find_arc_lines, so whichever merges second will need a rebase. They complement each other: the second-difference noise estimate here doesn't depend on the baseline, which suits the baseline subtraction in #320. Do you have a preferred merge order?
| _, _, clipped_std = sigma_clipped_stats( | ||
| diff2, mask=invalid, sigma=sigma, maxiters=maxiters, stdfunc=mad_std, axis=-1 | ||
| ) |
There was a problem hiding this comment.
There's an issue here with NaN or masked input. With axis=-1 and stdfunc=mad_std, sigma_clipped_stats (astropy 8.0.1) returns the plain unclipped np.nanstd as soon as anything is masked or non-finite. A synthetic arc with 4 Gaussian lines plus σ = 2 noise:
| input | measure_noise |
|---|---|
| clean | 2.12 |
| one pixel set to NaN | 25.9 |
| three pixels masked (masked array) | 26.1 |
The same data with axis=None gives 2.10, so it's specific to the axis path. test_nan_pixels_are_ignored doesn't catch it, I think because it's line-free noise, where clipping makes no difference. In find_arc_lines, a single NaN pixel in a spectrum without uncertainties would raise the detection threshold about 12×.
A per-row loop over the finite values avoids it, something like (untested):
rows = diff2.reshape(-1, diff2.shape[-1])
noise = np.full(rows.shape[0], np.nan)
for i, row in enumerate(rows):
row = row[np.isfinite(row)]
if row.size > 0:
noise[i] = sigma_clipped_stats(
row, sigma=sigma, maxiters=maxiters, stdfunc=mad_std
)[2]
noise = noise.reshape(diff2.shape[:-1]) / np.sqrt(6.0)It would also be good for the NaN test to include a few strong lines so it exercises the clipping. This might be worth reporting upstream to astropy as well.
| _, _, clipped_std = sigma_clipped_stats( | ||
| diff2, mask=invalid, sigma=sigma, maxiters=maxiters, stdfunc=mad_std, axis=-1 | ||
| ) | ||
| noise = clipped_std / np.sqrt(6.0) |
There was a problem hiding this comment.
The MAD collapses to 0 whenever more than half of the second differences are identical, and then every positive pixel passes the detection threshold. Some cases I tried through find_arc_lines (no input uncertainty):
- Low-count integer Poisson arc (true σ ≈ 0.45): 51% of the differences are exactly 0, so the noise is 0 and the result is 73 "lines", some with centroids outside the frame (e.g. −24.3).
- Spectrum zero-padded over 60% of its length: noise 0, 47 junk lines.
- Noiseless synthetic arc on a 100-count pedestal: noise 0, 1 line found out of 4.
- All-NaN spectrum: noise
nan, 0 lines, and no error.
The raw-count case is a regression compared with the old sqrt(|flux|) fallback. Maybe fall back to the clipped standard deviation when the MAD is 0 (e.g. rerun with the default stdfunc="std" for those rows), and have find_arc_lines raise a clear error if the estimate is still non-finite or ≤ 0?
| if spectrum.uncertainty is None: | ||
| spectrum = deepcopy(spectrum) | ||
| spectrum.uncertainty = StdDevUncertainty(np.sqrt(np.abs(spectrum.flux.value))) | ||
| noise = measure_noise(spectrum.flux.value) |
There was a problem hiding this comment.
This ignores spectrum.mask, so masked pixels (cosmic rays, bad columns) still go into the noise estimate. measure_noise already handles masked arrays, so this could be:
flux = spectrum.flux.value
if spectrum.mask is not None:
flux = np.ma.masked_array(flux, mask=spectrum.mask)
noise = measure_noise(flux)This only helps once the masked-input issue in measure_noise (comment on utils.py) is fixed; at the moment, adding a mask makes the estimate worse.
| (`~astropy.nddata.StdDevUncertainty`, `~astropy.nddata.VarianceUncertainty`, or | ||
| `~astropy.nddata.InverseVariance`); it is converted to a standard deviation | ||
| before the line finding. If the spectrum has no uncertainty, a constant per-pixel | ||
| noise is estimated from the data with `~specreduce.utils.utils.measure_noise`. |
There was a problem hiding this comment.
specreduce.utils.utils isn't in docs/api.rst (only specreduce.utils.synth_data is), so this cross-reference won't resolve. With nitpicky = True that gives a docs warning, and measure_noise won't appear in the API docs at all. Could you add an automodapi entry for it?
.. automodapi:: specreduce.utils.utils
:no-inheritance-diagram:The module's __all__ also exposes the private _align_along_trace, so you might want :include: measure_noise, measure_cross_dispersion_profile (or similar) to keep that out of the docs.
This PR fixes two issues in how the
find_arc_lineshandles uncertainties:the
find_arc_linesfunction now allows uncertainties to be any of Astropy's uncertainty types (it allowed previously onlyStdDevUncertainty)If the spectrum doesn't include uncertainties, the uncertainties are estimated using a robust approach based on the DER_SNR algorithm by Stoerh et al. (2008). (with a small difference that we're using Astropy's sigma-clipped median instead of plain median).
I've also added the noise estimation function as
specreduce.utils.utils.measure_noise.AI/LLM Disclaimer: The same as in #319. I've used Claude Code with Fable 5.1 to co-develop the PR, but I understand what the code does and how it works, and can explain it when needed.