From 05d239c1291769dfe529846f77dd1bb9a4d9f5c7 Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Fri, 9 Jan 2026 20:32:58 +0000 Subject: [PATCH 1/4] Added uncertainty propagation in `BoxcarExtract`, enabled support for various uncertainty types (`VarianceUncertainty`, `StdDevUncertainty`, `InverseVariance`), and added associated unit tests. --- specreduce/extract.py | 54 +++++++++++++++++++++++++------- specreduce/tests/test_extract.py | 30 ++++++++++++++++++ 2 files changed, 73 insertions(+), 11 deletions(-) diff --git a/specreduce/extract.py b/specreduce/extract.py index 8553124f..376a4225 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u from astropy.modeling import Model, models, fitting -from astropy.nddata import NDData, VarianceUncertainty +from astropy.nddata import NDData, VarianceUncertainty, StdDevUncertainty, InverseVariance from numpy import ndarray from scipy.integrate import trapezoid from scipy.interpolate import RectBivariateSpline @@ -178,7 +178,6 @@ class BoxcarExtract(SpecreduceOperation): width: float = 5 disp_axis: int = 1 crossdisp_axis: int = 0 - # TODO: should disp_axis and crossdisp_axis be defined in the Trace object? mask_treatment: MaskingOption = "apply" _valid_mask_treatment_methods = ( "apply", @@ -246,19 +245,52 @@ def __call__( # the window multiplied by the window width. window_weights = _ap_weight_image(trace, width, disp_axis, cdisp_axis, self.image.shape) + # Extract variance for uncertainty propagation (if available) + if self.image.uncertainty is not None: + variance = self.image.uncertainty.represent_as(VarianceUncertainty).array + orig_uncty_type = type(self.image.uncertainty) + else: + variance = None + orig_uncty_type = None + if self.mask_treatment == "apply": - image_cleaned = np.where(~self.image.mask, self.image.data * window_weights, 0.0) + flux = np.where(~self.image.mask, self.image.data * window_weights, 0.0) weights = np.where(~self.image.mask, window_weights, 0.0) - spectrum = ( - image_cleaned.sum(axis=cdisp_axis) - / weights.sum(axis=cdisp_axis) - * window_weights.sum(axis=cdisp_axis) - ) + weights_sum = weights.sum(axis=cdisp_axis) + window_sum = window_weights.sum(axis=cdisp_axis) + extracted_flux = flux.sum(axis=cdisp_axis) / weights_sum * window_sum + + if variance is not None: + variance = np.where(~self.image.mask, variance * window_weights**2, 0.0) + extracted_variance = variance.sum(axis=cdisp_axis) / weights_sum**2 * window_sum**2 + else: + flux = np.where(window_weights, self.image.data * window_weights, 0.0) + extracted_flux = flux.sum(axis=cdisp_axis) + + if variance is not None: + variance = np.where(window_weights, variance * window_weights**2, 0.0) + extracted_variance = np.sum(variance, axis=cdisp_axis) + + # Create output uncertainty (preserving original type) + if variance is not None: + # Convert variance to the original uncertainty type with proper units + if orig_uncty_type == VarianceUncertainty: + output_uncertainty = VarianceUncertainty(extracted_variance * self.image.unit**2) + elif orig_uncty_type == StdDevUncertainty: + output_uncertainty = StdDevUncertainty(np.sqrt(extracted_variance) * self.image.unit) + elif orig_uncty_type == InverseVariance: + output_uncertainty = InverseVariance(1.0 / extracted_variance / self.image.unit**2) + else: + # Fallback to VarianceUncertainty for unknown types + output_uncertainty = VarianceUncertainty(extracted_variance * self.image.unit**2) else: - image_windowed = np.where(window_weights, self.image.data * window_weights, 0.0) - spectrum = np.sum(image_windowed, axis=cdisp_axis) + output_uncertainty = None - return Spectrum(spectrum * self.image.unit, spectral_axis=self.image.spectral_axis) + return Spectrum( + extracted_flux * self.image.unit, + spectral_axis=self.image.spectral_axis, + uncertainty=output_uncertainty + ) @dataclass diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index 069996e3..0d029b0b 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -119,6 +119,36 @@ def test_boxcar_array_trace(mk_test_img): assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 75.0)) +def test_boxcar_uncertainty_propagation(): + """Test that BoxcarExtract propagates uncertainties correctly.""" + # Create image with known uniform flux and uncertainty + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + variance = np.full((nrows, ncols), 4.0) + + img = Spectrum( + flux * u.DN, + uncertainty=VarianceUncertainty(variance), + spectral_axis=np.arange(ncols) * u.pix + ) + + trace = FlatTrace(img, nrows // 2) + width = 3 + extracted = BoxcarExtract(img, trace, width=width)() + + # Check uncertainty exists and is correct type + assert extracted.uncertainty is not None + assert isinstance(extracted.uncertainty, VarianceUncertainty) + + # For width=3, summing 3 pixels each with variance=4 gives variance=12 + # (weights are 1.0 for full pixels in aperture) + expected_variance = width * 4.0 + np.testing.assert_allclose(extracted.uncertainty.array, expected_variance, rtol=0.01) + + # Check units are correct (flux_unit^2) + assert extracted.uncertainty.unit == u.DN**2 + + def test_horne_image_validation(mk_test_img): # # Test HorneExtract scenarios specific to its use with an image of From 6ef366818d3074f9f0299d924072c9ade19214ef Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Fri, 9 Jan 2026 21:56:56 +0000 Subject: [PATCH 2/4] Slightly refactored the uncertainty propagation code in `BoxcarExtract` and added more unit tests. --- specreduce/extract.py | 15 +++- specreduce/tests/test_extract.py | 134 ++++++++++++++++++++++++++++++- 2 files changed, 141 insertions(+), 8 deletions(-) diff --git a/specreduce/extract.py b/specreduce/extract.py index 376a4225..1b7a12d5 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -261,8 +261,13 @@ def __call__( extracted_flux = flux.sum(axis=cdisp_axis) / weights_sum * window_sum if variance is not None: - variance = np.where(~self.image.mask, variance * window_weights**2, 0.0) - extracted_variance = variance.sum(axis=cdisp_axis) / weights_sum**2 * window_sum**2 + extracted_variance = ( + np.where(~self.image.mask, variance * window_weights**2, 0.0).sum( + axis=cdisp_axis + ) + / weights_sum**2 + * window_sum**2 + ) else: flux = np.where(window_weights, self.image.data * window_weights, 0.0) extracted_flux = flux.sum(axis=cdisp_axis) @@ -277,7 +282,9 @@ def __call__( if orig_uncty_type == VarianceUncertainty: output_uncertainty = VarianceUncertainty(extracted_variance * self.image.unit**2) elif orig_uncty_type == StdDevUncertainty: - output_uncertainty = StdDevUncertainty(np.sqrt(extracted_variance) * self.image.unit) + output_uncertainty = StdDevUncertainty( + np.sqrt(extracted_variance) * self.image.unit + ) elif orig_uncty_type == InverseVariance: output_uncertainty = InverseVariance(1.0 / extracted_variance / self.image.unit**2) else: @@ -289,7 +296,7 @@ def __call__( return Spectrum( extracted_flux * self.image.unit, spectral_axis=self.image.spectral_axis, - uncertainty=output_uncertainty + uncertainty=output_uncertainty, ) diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index 0d029b0b..fc79a120 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -2,7 +2,13 @@ import pytest from astropy import units as u from astropy.modeling import models -from astropy.nddata import NDData, VarianceUncertainty, UnknownUncertainty +from astropy.nddata import ( + NDData, + VarianceUncertainty, + StdDevUncertainty, + InverseVariance, + UnknownUncertainty, +) from astropy.tests.helper import assert_quantity_allclose from specreduce.background import Background @@ -121,7 +127,7 @@ def test_boxcar_array_trace(mk_test_img): def test_boxcar_uncertainty_propagation(): """Test that BoxcarExtract propagates uncertainties correctly.""" - # Create image with known uniform flux and uncertainty + # Create an image with known uniform flux and uncertainty nrows, ncols = 10, 20 flux = np.full((nrows, ncols), 100.0) variance = np.full((nrows, ncols), 4.0) @@ -129,14 +135,14 @@ def test_boxcar_uncertainty_propagation(): img = Spectrum( flux * u.DN, uncertainty=VarianceUncertainty(variance), - spectral_axis=np.arange(ncols) * u.pix + spectral_axis=np.arange(ncols) * u.pix, ) trace = FlatTrace(img, nrows // 2) width = 3 extracted = BoxcarExtract(img, trace, width=width)() - # Check uncertainty exists and is correct type + # Check uncertainty exists and is the correct type assert extracted.uncertainty is not None assert isinstance(extracted.uncertainty, VarianceUncertainty) @@ -149,6 +155,126 @@ def test_boxcar_uncertainty_propagation(): assert extracted.uncertainty.unit == u.DN**2 +def test_boxcar_uncertainty_stddev_type(): + """Test that StdDevUncertainty type is preserved.""" + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + stddev = np.full((nrows, ncols), 2.0) + + img = Spectrum( + flux * u.DN, uncertainty=StdDevUncertainty(stddev), spectral_axis=np.arange(ncols) * u.pix + ) + + trace = FlatTrace(img, nrows // 2) + extracted = BoxcarExtract(img, trace, width=3)() + + assert isinstance(extracted.uncertainty, StdDevUncertainty) + # For width=3, variance=12, so stddev=sqrt(12) + expected_stddev = np.sqrt(3 * 4.0) + np.testing.assert_allclose(extracted.uncertainty.array, expected_stddev, rtol=0.01) + assert extracted.uncertainty.unit == u.DN + + +def test_boxcar_uncertainty_inverse_variance_type(): + """Test that InverseVariance type is preserved.""" + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + ivar = np.full((nrows, ncols), 0.25) + + img = Spectrum( + flux * u.DN, uncertainty=InverseVariance(ivar), spectral_axis=np.arange(ncols) * u.pix + ) + + trace = FlatTrace(img, nrows // 2) + extracted = BoxcarExtract(img, trace, width=3)() + + assert isinstance(extracted.uncertainty, InverseVariance) + # For width=3, variance=12, so ivar=1/12 + expected_ivar = 1.0 / (3 * 4.0) + np.testing.assert_allclose(extracted.uncertainty.array, expected_ivar, rtol=0.01) + assert extracted.uncertainty.unit == 1 / u.DN**2 + + +def test_boxcar_uncertainty_with_nonfinite(): + """Test uncertainty propagation with nonfinite values in data.""" + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + variance = np.full((nrows, ncols), 4.0) + flux[5, 5] = np.nan + + img = Spectrum( + flux * u.DN, + uncertainty=VarianceUncertainty(variance), + spectral_axis=np.arange(ncols) * u.pix, + ) + + trace = FlatTrace(img, 5) + + # With mask_treatment="apply", NaN pixels excluded from both flux and variance + extracted = BoxcarExtract(img, trace, width=3, mask_treatment="apply")() + assert extracted.uncertainty is not None + # Column 5 has only 2 valid pixels instead of 3 + # Variance formula: (sum w^2 sigma^2) / (sum w)^2 * W^2 = (2*4) / 4 * 9 = 18 + np.testing.assert_allclose(extracted.uncertainty.array[5], 18.0, rtol=0.01) + # Other columns should have normal variance=12 + np.testing.assert_allclose(extracted.uncertainty.array[0], 12.0, rtol=0.01) + + # With mask_treatment="ignore", NaN in flux propagates to flux but variance + # is computed from the variance array which is finite + extracted_ignore = BoxcarExtract(img, trace, width=3, mask_treatment="ignore")() + assert np.isnan(extracted_ignore.flux.value[5]) + # Variance is still computed from a finite variance array + np.testing.assert_allclose(extracted_ignore.uncertainty.array[5], 12.0, rtol=0.01) + + # Test that NaN in a variance array does propagate + variance_with_nan = variance.copy() + variance_with_nan[5, 5] = np.nan + img_var_nan = Spectrum( + np.full((nrows, ncols), 100.0) * u.DN, + uncertainty=VarianceUncertainty(variance_with_nan), + spectral_axis=np.arange(ncols) * u.pix, + ) + extracted_var_nan = BoxcarExtract(img_var_nan, trace, width=3, mask_treatment="ignore")() + assert np.isnan(extracted_var_nan.uncertainty.array[5]) + + +def test_boxcar_uncertainty_with_mask(): + """Test uncertainty propagation with masked pixels.""" + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + variance = np.full((nrows, ncols), 4.0) + mask = np.zeros((nrows, ncols), dtype=bool) + + # Mask one pixel in the extraction aperture at column 3 + mask[5, 3] = True + + img = Spectrum( + flux * u.DN, + uncertainty=VarianceUncertainty(variance), + mask=mask, + spectral_axis=np.arange(ncols) * u.pix, + ) + + trace = FlatTrace(img, 5) + extracted = BoxcarExtract(img, trace, width=3, mask_treatment="apply")() + assert extracted.uncertainty is not None + np.testing.assert_allclose(extracted.uncertainty.array[3], 18.0, rtol=0.01) + np.testing.assert_allclose(extracted.uncertainty.array[0], 12.0, rtol=0.01) + + +def test_boxcar_no_input_uncertainty(): + """Test that no uncertainty is returned when input has none.""" + from astropy.nddata import CCDData + + nrows, ncols = 10, 20 + flux = np.full((nrows, ncols), 100.0) + + img = CCDData(flux, unit=u.DN) + trace = FlatTrace(img, nrows // 2) + extracted = BoxcarExtract(img, trace, width=3)() + assert extracted.uncertainty is None + + def test_horne_image_validation(mk_test_img): # # Test HorneExtract scenarios specific to its use with an image of From 3e67b76a4583efa01adbe9bc88de53c9b119a120 Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Mon, 12 Jan 2026 22:39:48 +0000 Subject: [PATCH 3/4] - Simplified uncertainty handling in `BoxcarExtract`. - Cleaned up the related unit tests a bit. --- specreduce/extract.py | 25 ++++++++----------------- specreduce/tests/test_extract.py | 2 -- 2 files changed, 8 insertions(+), 19 deletions(-) diff --git a/specreduce/extract.py b/specreduce/extract.py index 1b7a12d5..61ccd713 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u from astropy.modeling import Model, models, fitting -from astropy.nddata import NDData, VarianceUncertainty, StdDevUncertainty, InverseVariance +from astropy.nddata import NDData, VarianceUncertainty from numpy import ndarray from scipy.integrate import trapezoid from scipy.interpolate import RectBivariateSpline @@ -253,6 +253,7 @@ def __call__( variance = None orig_uncty_type = None + extracted_variance = None if self.mask_treatment == "apply": flux = np.where(~self.image.mask, self.image.data * window_weights, 0.0) weights = np.where(~self.image.mask, window_weights, 0.0) @@ -276,27 +277,17 @@ def __call__( variance = np.where(window_weights, variance * window_weights**2, 0.0) extracted_variance = np.sum(variance, axis=cdisp_axis) - # Create output uncertainty (preserving original type) - if variance is not None: - # Convert variance to the original uncertainty type with proper units - if orig_uncty_type == VarianceUncertainty: - output_uncertainty = VarianceUncertainty(extracted_variance * self.image.unit**2) - elif orig_uncty_type == StdDevUncertainty: - output_uncertainty = StdDevUncertainty( - np.sqrt(extracted_variance) * self.image.unit - ) - elif orig_uncty_type == InverseVariance: - output_uncertainty = InverseVariance(1.0 / extracted_variance / self.image.unit**2) - else: - # Fallback to VarianceUncertainty for unknown types - output_uncertainty = VarianceUncertainty(extracted_variance * self.image.unit**2) + if extracted_variance is not None: + spectrum_uncty = VarianceUncertainty( + extracted_variance * self.image.unit**2 + ).represent_as(orig_uncty_type) else: - output_uncertainty = None + spectrum_uncty = None return Spectrum( extracted_flux * self.image.unit, spectral_axis=self.image.spectral_axis, - uncertainty=output_uncertainty, + uncertainty=spectrum_uncty, ) diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index fc79a120..9df36c6c 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -169,7 +169,6 @@ def test_boxcar_uncertainty_stddev_type(): extracted = BoxcarExtract(img, trace, width=3)() assert isinstance(extracted.uncertainty, StdDevUncertainty) - # For width=3, variance=12, so stddev=sqrt(12) expected_stddev = np.sqrt(3 * 4.0) np.testing.assert_allclose(extracted.uncertainty.array, expected_stddev, rtol=0.01) assert extracted.uncertainty.unit == u.DN @@ -189,7 +188,6 @@ def test_boxcar_uncertainty_inverse_variance_type(): extracted = BoxcarExtract(img, trace, width=3)() assert isinstance(extracted.uncertainty, InverseVariance) - # For width=3, variance=12, so ivar=1/12 expected_ivar = 1.0 / (3 * 4.0) np.testing.assert_allclose(extracted.uncertainty.array, expected_ivar, rtol=0.01) assert extracted.uncertainty.unit == 1 / u.DN**2 From 2126702f2a1989a8bcc3e3b88af74fc2c479f1cb Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Tue, 13 Jan 2026 19:34:31 +0000 Subject: [PATCH 4/4] Updated CHANGES.rst to mention added uncertainty propagation in `BoxcarExtract`. --- CHANGES.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 533fc56f..1bbdaad1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,11 @@ -1.7.1 (2025-11-2x) +1.8.0 (2026-01-xx) ------------------ +New Features +^^^^^^^^^^^^ + +- Added uncertainty propagation to ``specreduce.extract.BoxcarExtract``. + Other changes ^^^^^^^^^^^^^