WavelengthSolution1D.resample bug fixes and improvements

[hpparvi]

This PR fixes several bugs in the WavelengthSolution1D.resample method and improves the method's overall functionality.

Bug fixes

• The flux was multiplied by the pixel wavelength width before dividing by the bin width, so the output was labeled as a flux density but was off by the local dispersion and did not conserve the integrated flux. The density and its uncertainty are now computed solely from the fractional pixel overlaps, with the variance using the squared weights. (I somehow remembered I had fixed this a long time ago, so I'm not sure what happened here...)
• Spectra whose spectral axis does not start at pixel zero were indexed from the wrong pixels.
• Solutions where the wavelength decreases with pixel number silently produced zero flux. The output is now always ascending in wavelength.

Changes

• The input mask is propagated (a bin is flagged if any contributing pixel is masked) and the metadata is copied.
• A spectrum without uncertainty returns no uncertainty instead of a fabricated zero-valued one.
• The output spectral axis carries the bin edges actually used, so spectral_axis.bin_edges is exact for non-uniform grids and integrating the flux over those edges recovers the input counts exactly.

AI and LLM use disclaimer: I'm systematically going through the specreduce codebase with Claude Code (Fable 5.1) to identify bugs, issues, and any room for (non-breaking) improvements in the existing functionality before the specreduce 2.10 release. The first, rather serious, bug I was aware of beforehand, but I remembered I had it already fixed several months ago. I also used Claude Code to implement most of the changes. I understand the code and the changes and could have implemented them myself, but it would have taken significantly longer than the 2-3 hours of continuous work it took to create this PR.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.54%. Comparing base (4e330bc) to head (e4708d5).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #319      +/-   ##
==========================================
+ Coverage   92.52%   92.54%   +0.02%     
==========================================
  Files          18       18              
  Lines        2341     2348       +7     
==========================================
+ Hits         2166     2173       +7     
  Misses        175      175              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[tepickering]

looks like python 3.12 is now required, but the CI not updated to remove 3.11. the windows failure is the old tcl bug. i think that can be fixed by forcing matplotlib to use the Agg head-less backend in the test fixtures.

[tepickering]

The most substantive changes needed are to address cpu-dependent handling of bins outside the solution range and masking of bins with partial/missing coverage. The checking of spectral_axis and bin_edges would be good to cover as well. The rest are recommended, but not necessary. Vectorizing of python loops is always a good idea if it can be managed straightforwardly.

[tepickering]

w2p is an interp1d with fill_value=np.nan that only spans bounds_pix ± 2 pixels, so any bin edge outside that range comes back as NaN. np.clip passes the NaN through, and then np.floor(nan).astype(int) is platform-dependent: on arm64 it gives 0, so the bin quietly becomes NaN with a RuntimeWarning: invalid value encountered in cast. On x86_64 it gives INT_MIN, which should end in an IndexError instead.

Reproducer: a linear solution starting at 5000 Å over (0, 100), then ws.resample(spectrum, wlbounds=(4990, 5100), nbins=10).

This predates the PR, but since this indexing is being rewritten anyway, could we handle it explicitly? Either raise a ValueError when the requested range extends beyond the solution, or treat NaN edges as "no coverage" (e.g. np.nan_to_num before the clip, then mask those bins; see the next comment).

[tepickering]

Related to the above: bins that extend past the input pixels are clipped, so their flux density is underestimated (or zero if they fall completely outside), and nothing marks them. The mask is only built when the input already has one, so a spectrum without a mask gets mask=None back even when some bins are NaN or only partially covered.

Since the PR now propagates masks, how about always returning one that also flags bins with incomplete pixel coverage? Something like:

x_raw = self.w2p(bin_edges_wav) + 0.5 - pixels[0]
uncovered = ~((np.minimum(x_raw[:-1], x_raw[1:]) >= 0) & (np.maximum(x_raw[:-1], x_raw[1:]) <= npix))
# NaN edges compare False, so they end up flagged as well
mask_wl = uncovered | propagated_mask

That would let users drop the edge bins downstream without having to recompute coverage themselves.

[tepickering]

The docstring now says the input spectral axis "must be in pixels with unit spacing", but nothing checks this. A spectrum whose spectral axis is already in Å runs without error and returns all-NaN flux, and a 2D flux array fails with "setting an array element with a sequence". A short guard would turn these into clear errors:

if spectrum.flux.ndim != 1:
    raise ValueError("resample only supports 1D spectra.")
if spectrum.spectral_axis.unit != u.pix:
    raise ValueError("The spectral axis of the input spectrum must be in pixels.")
if npix > 1 and not np.allclose(np.diff(pixels), 1.0):
    raise ValueError("The spectral axis of the input spectrum must have unit spacing.")

[tepickering]

np.sort quietly reorders whatever is passed in, which can hide a user error. Duplicate edges give zero-width bins that come out as NaN with "invalid value encountered in divide" warnings, and a single edge returns an empty spectrum without complaint. I'd suggest validating instead:

bin_edges_wav = np.asarray(bin_edges, dtype=float)
if bin_edges_wav.ndim != 1 or bin_edges_wav.size < 2:     
    raise ValueError("bin_edges must be a 1D array with at least two edges.")
if np.all(np.diff(bin_edges_wav) < 0):
    bin_edges_wav = bin_edges_wav[::-1]
if not np.all(np.diff(bin_edges_wav) > 0):
    raise ValueError("bin_edges must be strictly monotonic.")

The wlbounds path could get a matching check that l1 != l2.

[tepickering]

Small housekeeping item: the four new entries still have [#XXX] placeholders. These should be [#319].

[tepickering]

This claims a bit more precision than the code delivers. w2p is a linear interp1d of the polynomial on an integer pixel grid, so the bin edges in pixel space (and therefore the conserved flux) are only approximate for nonlinear solutions. The updated test_resample comment says the same thing and uses rtol=1e-4.

The paragraph also says the binning is exact and flux-conserving twice. Maybe something like:

Suggested change

	The bin edges are stored on the output spectral axis, so the integrated flux is
	The bin edges are stored on the output spectral axis, so the integrated flux can be
	recovered as ``(flux * np.diff(spectral_axis.bin_edges)).sum()`` for any bin grid.
	Flux is conserved up to the accuracy of the interpolated wavelength-to-pixel
	inverse, which is exact for linear solutions. The variance is propagated

[tepickering]

This text isn't changed by the PR, but since the docstring is being edited: with edges [e_0, ..., e_n] there are n bins, not n = nbins + 1, and the last bin should be [e_n-1, e_n]. Maybe also mention the ordering and validation behaviour here, whatever ends up being decided in the earlier bin_edges comment:

Suggested change

	bin_edges
	Explicit bin edges in the wavelength space, given as a 1D array-like
	``[e_0, e_1, ..., e_n]`` that defines the ``n`` bins ``[e_0, e_1], [e_1, e_2], ...,
	[e_n-1, e_n]``. If provided, ``nbins`` and ``wlbounds`` are ignored.

[tepickering]

Now that the output carries its bin edges, this could use them directly. That also exercises the new SpectralAxis(..., bin_specification="edges") path instead of assuming a uniform grid:

Suggested change

	f1 = (resampled.flux.value * np.diff(resampled.spectral_axis.value)[0]).sum()
	f1 = (resampled.flux.value * np.diff(resampled.spectral_axis.bin_edges.value)).sum()

It might also be worth adding a non-uniform case, e.g. resampling with bin_edges=np.geomspace(...) and checking that the integrated flux matches the input to the same tolerance.

[tepickering]

Optional, and fine to leave for a follow-up: the per-bin Python loop could be vectorized with cumulative sums, which would help for long spectra or fine output grids. Since the weights are 1 for interior pixels and fractional only at the two ends, everything reduces to cumulative-sum lookups. Here's a rough sketch (untested, and it assumes the clipped x_left/x_right/i_left/i_right from above):

fl = x_left - i_left            # fractional position of the left edge within its pixel
fr = x_right - i_right          # fractional position of the right edge within its pixel
same = i_left == i_right

# flux: the cumulative integral evaluated at the edges
cflux = np.concatenate([[0.0], np.cumsum(flux)])
flux_wl = np.interp(x_right, np.arange(npix + 1), cflux) - np.interp(
    x_left, np.arange(npix + 1), cflux
)

# variance: squared weights at the end pixels plus the plain sum over interior pixels
cvar = np.concatenate([[0.0], np.cumsum(ucty)])
interior = cvar[i_right] - cvar[np.minimum(i_left + 1, i_right)]
ucty_wl = np.where(
    same,
    (x_right - x_left) ** 2 * ucty[i_left],
    (1 - fl) ** 2 * ucty[i_left] + fr**2 * ucty[i_right] + interior,
)

# mask: any masked pixel with a non-zero weight
cmask = np.concatenate([[0], np.cumsum(mask)])
mask_wl = (cmask[i_right] - cmask[np.minimum(i_left + 1, i_right)] > 0) | (
    mask[i_left] & (fl < 1)
) | (mask[i_right] & (fr > 0))

The existing tests (offset axis, decreasing dispersion, variance, mask) should be enough to verify it matches the loop.

[tepickering]

Also, look into the failure on the dev dependencies test. That may take care of itself once #324 is merged.