Skip to content

Add error propagation estimates for BoxcarExtract and HorneExtract #286

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add error propagation estimates for BoxcarExtract and HorneExtract #286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
326 changes: 326 additions & 0 deletions notebook_sandbox/horne_extract/horne_uncertainty_testing.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,326 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Horne Extraction - Toy Tests\n",
"\n",
"This notebook recreates some simple tests to run for `HorneExtract`.\n",
"Each section builds a synthetic 2‑D spectrum and checks how uncertainties are propagated under different conditions.\n",
"\n",
"**Key ideas:**\n",
"- The input is a 2‑D long‑slit spectrum (spatial x spectral).\n",
"- We supply per‑pixel variance via `NDData.uncertainty` in a few different flavors.\n",
"- `HorneExtract` uses a spatial profile and inverse‑variance weighting: \n",
" $F(\\lambda)=\\frac{\\sum P\\,D/\\sigma^2}{\\sum P^2/\\sigma^2}$ and $\\mathrm{Var}[F]=\\frac{1}{\\sum P^2/\\sigma^2}$.\n",
"\n",
"You should be able to run each cell top‑to‑bottom."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1e8f4f8",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import astropy.units as u\n",
"from astropy.nddata import NDData, VarianceUncertainty, StdDevUncertainty\n",
"\n",
"from specreduce.extract import HorneExtract\n",
"from specreduce.tracing import FlatTrace"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Global synthetic image settings\n",
"ny, nx = 20, 60 # spatial, spectral\n",
"true_flux = 100.0 # DN per pixel\n",
"sigma_per_pixel = 5.0 # DN stddev per pixel\n",
"rng = np.random.default_rng(42) # reproducible noise\n",
"\n",
"def make_base_data(ny=ny, nx=nx, mu=true_flux, sigma=sigma_per_pixel):\n",
" return mu + rng.normal(0, sigma, size=(ny, nx))\n",
"\n",
"def set_spectral_axis(image, ndisp=nx):\n",
" # Specutils wants a spectral axis that matches a flux axis.\n",
" image.spectral_axis = np.arange(ndisp) * u.pix\n",
" return image\n",
"\n",
"def run_case(title, image, trace_pos=None, profile=\"gaussian\"):\n",
" print(f\"=== {title} ===\")\n",
" # Provide a spectral axis on the NDData\n",
" set_spectral_axis(image, ndisp=image.data.shape[1])\n",
"\n",
" # Choose a flat, centered trace unless requested otherwise\n",
" if trace_pos is None:\n",
" trace_pos = image.data.shape[0] // 2\n",
" trace = FlatTrace(image, trace_pos=trace_pos)\n",
"\n",
" # Build extractor and run\n",
" extractor = HorneExtract(\n",
" image=image,\n",
" trace_object=trace,\n",
" spatial_profile=profile,\n",
" disp_axis=1,\n",
" crossdisp_axis=0,\n",
" )\n",
" try:\n",
" spec = extractor()\n",
" arr = np.asarray(spec.uncertainty.array) if getattr(spec, 'uncertainty', None) is not None else None\n",
" print(\"flux shape:\", spec.flux.shape)\n",
" if arr is None:\n",
" print(\"No uncertainty returned\")\n",
" else:\n",
" print(\"first 5 uncertainties:\", arr[:5])\n",
" print(\"finite, positive:\", np.isfinite(arr).all() and np.all(arr > 0))\n",
" except Exception as e:\n",
" print(\"ERROR:\", e)\n",
" print()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case A - VarianceUncertainty (correct units)\n",
"We pass per‑pixel variance as `VarianceUncertainty` with units of `DN^2`. \n",
"Expected: uncertainties are finite and roughly constant across wavelength for stationary noise."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataA = make_base_data()\n",
"var_unc = VarianceUncertainty(np.full((ny, nx), sigma_per_pixel**2) * u.DN**2)\n",
"imgA = NDData(data=dataA * u.DN, uncertainty=var_unc)\n",
"run_case(\"A: VarianceUncertainty (correct units)\", imgA)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case B - StdDevUncertainty (correct units)\n",
"We pass standard deviation with units (`DN`). Internally it will be squared to variance. \n",
"Expected: same result as Case A."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataB = make_base_data()\n",
"std_unc = StdDevUncertainty(np.full((ny, nx), sigma_per_pixel) * u.DN)\n",
"imgB = NDData(data=dataB * u.DN, uncertainty=std_unc)\n",
"run_case(\"B: StdDevUncertainty (correct units)\", imgB)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case C - StdDevUncertainty (no units)\n",
"We pass a unitless stddev array; the parser will assume it's in the same units as the data. \n",
"Expected: matches Case B."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataC = make_base_data()\n",
"std_unc_nounits = StdDevUncertainty(np.full((ny, nx), sigma_per_pixel))\n",
"imgC = NDData(data=dataC * u.DN, uncertainty=std_unc_nounits)\n",
"run_case(\"C: StdDevUncertainty (no units)\", imgC)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case D - No uncertainty\n",
"We omit uncertainties; `HorneExtract` requires them. This should raise an error."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataD = make_base_data()\n",
"imgD = NDData(data=dataD * u.DN)\n",
"run_case(\"D: No uncertainty\", imgD)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case E - Negative variance\n",
"Construct a variance array with a negative element, should be rejected."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataE = make_base_data()\n",
"v = np.full((ny, nx), sigma_per_pixel**2)\n",
"v[0,0] = -1.0\n",
"var_unc_bad = VarianceUncertainty(v * u.DN**2)\n",
"imgE = NDData(data=dataE * u.DN, uncertainty=var_unc_bad)\n",
"run_case(\"E: Negative variance\", imgE)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case F - Zero variance everywhere\n",
"All zeros are treated as an unweighted case inside the implementation (variances replaced so extraction can proceed)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataF = make_base_data()\n",
"var_unc_zero = VarianceUncertainty(np.zeros((ny, nx)) * u.DN**2)\n",
"imgF = NDData(data=dataF * u.DN, uncertainty=var_unc_zero)\n",
"run_case(\"F: Zero variance everywhere\", imgF)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case G - Some zero variances\n",
"Only the first spatial row has zero variance; those pixels are masked out for the variance calculation.\n",
"Expected: finite uncertainties; first elements may differ slightly."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataG = make_base_data()\n",
"v = np.full((ny, nx), sigma_per_pixel**2)\n",
"v[0, :] = 0.0\n",
"var_unc_mixed = VarianceUncertainty(v * u.DN**2)\n",
"imgG = NDData(data=dataG * u.DN, uncertainty=var_unc_mixed)\n",
"run_case(\"G: Some zero variances\", imgG)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case H/I - Trace near detector edges\n",
"We move the flat trace near the top/bottom rows to simulate partial apertures hitting the edge.\n",
"Expected: still finite uncertainties (the spatial profile normalization handles partial coverage)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataH = make_base_data()\n",
"var_unc_H = VarianceUncertainty(np.full((ny, nx), sigma_per_pixel**2) * u.DN**2)\n",
"imgH = NDData(data=dataH * u.DN, uncertainty=var_unc_H)\n",
"run_case(\"H: Trace near top edge\", imgH, trace_pos=1)\n",
"\n",
"dataI = make_base_data()\n",
"var_unc_I = VarianceUncertainty(np.full((ny, nx), sigma_per_pixel**2) * u.DN**2)\n",
"imgI = NDData(data=dataI * u.DN, uncertainty=var_unc_I)\n",
"run_case(\"I: Trace near bottom edge\", imgI, trace_pos=ny-2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case J - Fully masked input\n",
"We feed an image where all pixels are NaN; the input parser should reject it as fully masked."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataJ = np.full((ny, nx), np.nan)\n",
"var_unc_J = VarianceUncertainty(np.full((ny, nx), sigma_per_pixel**2) * u.DN**2)\n",
"imgJ = NDData(data=dataJ * u.DN, uncertainty=var_unc_J)\n",
"run_case(\"J: Fully masked input\", imgJ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Case K - Zero profile sumP in a column\n",
"We force one spectral column to be entirely NaN, so the spatial profile normalization `sumP` would be zero there.\n",
"Implementation guards against division by zero by replacing `sumP==0` with 1.\n",
"Expected: uncertainties remain finite."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataK = make_base_data()\n",
"dataK[:, nx//2] = np.nan # one dead spectral column\n",
"var_unc_K = VarianceUncertainty(np.full((ny, nx), sigma_per_pixel**2) * u.DN**2)\n",
"imgK = NDData(data=dataK * u.DN, uncertainty=var_unc_K)\n",
"run_case(\"K: Zero profile sumP\", imgK)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "testing",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading
Loading