From e541bcc5fbcb21fc327f5fd5c3385b46cb7e188b Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 14:43:12 -0600 Subject: [PATCH 01/13] Add electron-based instrument interface and `LinearOptikaInstrument` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Make the CTIS instrument interface electron-based: `image()` returns the electrons measured by the sensor and `backproject()` accepts electrons. `IdealInstrument` gains a `quantum_yield` parameter and converts photons<->electrons around its regridding model. `LinearOptikaInstrument` wraps an `optika.systems.LinearSystem` and delegates to its `image_from_weights`/`backproject_from_weights`, reusing the cached regridding weights so the forward/transpose model (including the sensor's `expose`/`photons_absorbed`) is not rebuilt each call. `image()` gains an `uncertainty` flag: when set, it attaches the standard deviation of the measurement noise as a `NormalUncertainScalarArray`. The variance is computed per wavelength and summed in quadrature during the wavelength integration, so it is exact even for the integrated image. The MART inverter uses this to weight its chi-squared, replacing the old per-instrument `uncertainty` callable. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ Co-Authored-By: Claude Opus 4.8 --- ctis/instruments/__init__.py | 2 + ctis/instruments/_instruments.py | 284 ++++++++++++++++++++--- ctis/instruments/_instruments_test.py | 17 +- ctis/inverters/_iterative/_iterative.py | 6 +- ctis/inverters/_iterative/_mart/_mart.py | 5 +- ctis/inverters/_results_test.py | 0 demo.ipynb | 38 +++ na.ipynb | 205 ++++++++++++++++ na.py | 0 pyproject.toml | 1 + 10 files changed, 518 insertions(+), 40 deletions(-) create mode 100644 ctis/inverters/_results_test.py create mode 100644 demo.ipynb create mode 100644 na.ipynb create mode 100644 na.py diff --git a/ctis/instruments/__init__.py b/ctis/instruments/__init__.py index bcc5400..95b43f5 100644 --- a/ctis/instruments/__init__.py +++ b/ctis/instruments/__init__.py @@ -6,10 +6,12 @@ AbstractInstrument, AbstractLinearInstrument, IdealInstrument, + LinearOptikaInstrument, ) __all__ = [ "AbstractInstrument", "AbstractLinearInstrument", "IdealInstrument", + "LinearOptikaInstrument", ] diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index f65ae17..3113bd4 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -6,11 +6,13 @@ import astropy.units as u import astropy.constants import named_arrays as na +import optika __all__ = [ "AbstractInstrument", "AbstractLinearInstrument", "IdealInstrument", + "LinearOptikaInstrument", ] @@ -42,10 +44,11 @@ def image( scene: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, noise: bool = True, + uncertainty: bool = False, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: r""" The forward model of this CTIS instrument, which maps spectral radiance - on the skyplane to photons measured by the instrument's sensor. + on the skyplane to the electrons measured by the instrument's sensor. Parameters ---------- @@ -61,6 +64,13 @@ def image( for demonstration purposes. noise Whether to include the effect of noise in the final image. + uncertainty + Whether to attach the standard deviation of the measurement noise + to the result, as a + :class:`~named_arrays.NormalUncertainScalarArray`. + The variance is computed for each wavelength *before* the + integration along the wavelength axis and summed in quadrature, so + it is exact even for the integrated image. """ @abc.abstractmethod @@ -106,14 +116,6 @@ def coordinates_sensor(self) -> na.AbstractSpectralPositionalVectorArray: A grid of wavelength and position coordinates on the detector plane. """ - @property - @abc.abstractmethod - def uncertainty(self) -> Callable[[na.ScalarArray], na.ScalarArray]: - """ - A function that returns the standard deviation of the uncertainty - for a given number of photons. - """ - @property @abc.abstractmethod def channel(self): @@ -238,6 +240,41 @@ def _energy_per_photon(self) -> u.Quantity | na.AbstractScalar: return energy_per_photon + def _integrate_wavelength( + self, + outputs: na.AbstractScalar, + coordinates: na.AbstractSpectralPositionalVectorArray, + ) -> tuple[na.AbstractScalar, na.AbstractSpectralPositionalVectorArray]: + """ + Integrate an image along the wavelength axis and collapse the + wavelength coordinates to the band edges. + + If `outputs` carries an uncertainty (a + :class:`~named_arrays.NormalUncertainScalarArray`), the per-wavelength + variances are summed in quadrature, which is exact because the noise in + each wavelength bin is independent. + """ + axis = self.axis_wavelength + + if isinstance(outputs, na.NormalUncertainScalarArray): + nominal = outputs.nominal.sum(axis) + width = np.sqrt(np.square(outputs.width).sum(axis)) + outputs = na.NormalUncertainScalarArray(nominal, width) + else: + outputs = outputs.sum(axis) + + coordinates = coordinates.replace( + wavelength=na.stack( + arrays=[ + coordinates.wavelength[{axis: +0}], + coordinates.wavelength[{axis: ~0}], + ], + axis=axis, + ) + ) + + return outputs, coordinates + def image( self, scene: na.AbstractScalar | na.AbstractFunctionArray, @@ -271,19 +308,9 @@ def image( coordinates = self.coordinates_sensor if integrate: - - axis = self.axis_wavelength - - values_output = values_output.sum(axis) - - coordinates = coordinates.replace( - wavelength=na.stack( - arrays=[ - coordinates.wavelength[{axis: +0}], - coordinates.wavelength[{axis: ~0}], - ], - axis=axis, - ) + values_output, coordinates = self._integrate_wavelength( + values_output, + coordinates, ) return na.FunctionArray( @@ -381,7 +408,7 @@ class IdealInstrument( position_ref: u.Quantity | na.AbstractScalar | na.Cartesian2dVectorArray """ The position on the sensor where center of the FOV lands at the reference - wavelength. + wavelength. """ coordinates_scene: na.AbstractSpectralPositionalVectorArray = dataclasses.MISSING @@ -427,12 +454,21 @@ class IdealInstrument( changing position coordinate. """ - @property - def uncertainty(self) -> Callable[[na.ScalarArray], na.ScalarArray]: - def _shot_noise(image: na.ScalarArray) -> na.ScalarArray: - return np.sqrt(image.to_value(u.ph)) * u.ph + quantum_yield: u.Quantity | na.AbstractScalar = 1 * u.electron / u.photon + r""" + The number of electrons generated in the sensor per incident photon, in + units equivalent to :math:`\text{electron} \, \text{photon}^{-1}`. - return _shot_noise + For this idealized instrument the quantum yield is a constant, so it can be + applied after integrating over wavelength. + """ + + def _shot_noise(self, image: na.ScalarArray) -> na.ScalarArray: + # photon shot noise, converted back into electrons to match the + # electron-valued image + photons = image / self.quantum_yield + uncertainty = np.sqrt(photons.to_value(u.ph)) * u.ph + return uncertainty * self.quantum_yield def distortion(self, coordinates: na.SpectralPositionalVectorArray): """ @@ -539,22 +575,47 @@ def image( scene: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, noise: bool = True, + uncertainty: bool = False, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: scene = scene * self.area_effective * self.timedelta_exposure - return super().image( + # keep the wavelength axis so the uncertainty can be computed for each + # wavelength before integrating. + result = super().image( scene=scene, - integrate=integrate, + integrate=False, noise=noise, ) + # convert the measured photons into electrons + electrons = result.outputs * self.quantum_yield + coordinates = result.inputs + + if uncertainty: + width = self._shot_noise(electrons) + electrons = na.NormalUncertainScalarArray(nominal=electrons, width=width) + + if integrate: + electrons, coordinates = self._integrate_wavelength(electrons, coordinates) + + return na.FunctionArray( + inputs=coordinates, + outputs=electrons, + ) + def backproject( self, image: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + # convert the measured electrons back into photons + if isinstance(image, na.AbstractFunctionArray): + image = image.replace(outputs=image.outputs / self.quantum_yield) + else: + image = image / self.quantum_yield + result = super().backproject( image=image, integrate=integrate, @@ -563,3 +624,166 @@ def backproject( result = result / (self.area_effective * self.timedelta_exposure) return result + + +@dataclasses.dataclass +class LinearOptikaInstrument( + AbstractLinearInstrument, +): + """ + A CTIS instrument whose forward model is an :mod:`optika` + :class:`~optika.systems.AbstractLinearSystem`. + + The optika system supplies the distortion, effective area, and vignetting; + this class adapts its regridding forward model to the + :class:`AbstractLinearInstrument` interface and adds the transpose + (:meth:`backproject`) used during inversion. The system may be + *channel-aware*: its component models can vary along :attr:`axis_channel` + to represent the different CTIS projections. + """ + + system: optika.systems.AbstractLinearSystem = dataclasses.MISSING + """A :mod:`optika` representation of a linear optical system.""" + + coordinates_scene: na.AbstractSpectralPositionalVectorArray = dataclasses.MISSING + """ + A grid of wavelength and position coordinates on the skyplane + which will be used to construct the inverted scene. + + Normally the pitch of this grid is chosen to be the average + plate scale of the instrument. + """ + + channel: str | na.AbstractScalar = dataclasses.MISSING + """ + Human-readable name of each independent CTIS channel. + """ + + axis_channel: str | tuple[str, ...] = dataclasses.MISSING + """ + The logical axis or axes of :attr:`system` corresponding to the different + CTIS channels. + """ + + axis_wavelength: str = dataclasses.MISSING + """ + The logical axis of :attr:`coordinates_scene` corresponding to changing + wavelength coordinate. + """ + + axis_scene_xy: tuple[str, str] = dataclasses.MISSING + """ + The logical axes of :attr:`coordinates_scene` corresponding to changing + position coordinate. + """ + + @property + def axis_sensor_xy(self) -> tuple[str, str]: + axis_pixel = self.system.sensor.axis_pixel + return (axis_pixel.x, axis_pixel.y) + + @property + def coordinates_sensor(self) -> na.AbstractSpectralPositionalVectorArray: + return na.SpectralPositionalVectorArray( + wavelength=self.coordinates_scene.wavelength, + position=self.system.coordinates_sensor, + ) + + @functools.cached_property + def weights(self) -> tuple[na.AbstractScalar, dict[str, int], dict[str, int]]: + return self.system.weights( + coordinates=self.coordinates_scene, + axis_wavelength=self.axis_wavelength, + axis_field=self.axis_scene_xy, + ) + + @functools.cached_property + def weights_transpose( + self, + ) -> tuple[na.AbstractScalar, dict[str, int], dict[str, int]]: + return self.system.weights_transposed( + weights=self.weights, + coordinates=self.coordinates_scene, + axis_wavelength=self.axis_wavelength, + axis_field=self.axis_scene_xy, + ) + + def image( + self, + scene: na.AbstractScalar | na.AbstractFunctionArray, + integrate: bool = True, + noise: bool = True, + uncertainty: bool = False, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + + if isinstance(scene, na.AbstractFunctionArray): + if not np.all(scene.inputs == self.coordinates_scene): + raise ValueError( + "`scene.inputs` and `self.coordinates_scene` are not equal." + ) + else: + scene = na.FunctionArray(inputs=self.coordinates_scene, outputs=scene) + + # apply the optika forward model (effective area, vignetting, and the + # sensor response) using the cached regridding weights, giving the + # electrons measured in each pixel (optionally with their uncertainty). + result = self.system.image_from_weights( + self.weights, + scene, + axis_wavelength=self.axis_wavelength, + axis_field=self.axis_scene_xy, + noise=noise, + uncertainty=uncertainty, + ) + + coordinates = result.inputs + values_output = result.outputs + + if integrate: + values_output, coordinates = self._integrate_wavelength( + values_output, + coordinates, + ) + + return na.FunctionArray( + inputs=coordinates, + outputs=values_output, + ) + + def backproject( + self, + image: na.AbstractScalar | na.AbstractFunctionArray, + integrate: bool = True, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + + if isinstance(image, na.AbstractFunctionArray): + if not np.all(image.inputs.position == self.coordinates_sensor.position): + raise ValueError( + "`image.inputs` and `self.coordinates_sensor` are not equal." + ) + values = image.outputs + else: + values = image + + axis_wavelength = self.axis_wavelength + num_wavelength = self.coordinates_scene.wavelength.shape[axis_wavelength] - 1 + + if integrate: + values = values / num_wavelength + + # rebuild the detector image over the full sensor wavelength grid so + # optika can invert the sensor response for each wavelength. + image = na.FunctionArray( + inputs=self.coordinates_sensor, + outputs=values, + ) + + # apply the transposed optika model with the cached transpose weights, + # inverting the sensor response and recovering the spectral radiance. + return self.system.backproject_from_weights( + self.weights_transpose, + image, + coordinates=self.coordinates_scene, + axis_wavelength=axis_wavelength, + axis_field=self.axis_scene_xy, + ) diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index 9f1ad33..0f76dbf 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -116,19 +116,24 @@ def test_num_channel( assert isinstance(result, int) @pytest.mark.parametrize( - argnames="image", + argnames="scene", argvalues=[ - instrument_ideal.image(gaussians.outputs).outputs, + gaussians.outputs, ], ) - def test_uncertainty( + @pytest.mark.parametrize("integrate", [False, True]) + def test_image_uncertainty( self, a: ctis.instruments.AbstractInstrument, - image: na.ScalarArray, + scene: na.AbstractScalar | na.AbstractFunctionArray, + integrate: bool, ): - result = a.uncertainty(image) + result = a.image(scene, integrate=integrate, noise=False, uncertainty=True) - assert np.all(result >= 0 * u.photon) + # the measurement noise is attached as a normal uncertain array + assert isinstance(result.outputs, na.NormalUncertainScalarArray) + assert result.outputs.width.unit.is_equivalent(u.electron) + assert np.all(result.outputs.width >= 0 * u.electron) class AbstractTestAbstractLinearInstrument( diff --git a/ctis/inverters/_iterative/_iterative.py b/ctis/inverters/_iterative/_iterative.py index 9d3a0f5..a1b1f76 100644 --- a/ctis/inverters/_iterative/_iterative.py +++ b/ctis/inverters/_iterative/_iterative.py @@ -45,6 +45,7 @@ def mean_chi_squared( self, images_observed: na.ScalarArray, images_predicted: na.ScalarArray, + uncertainty: na.ScalarArray, ) -> na.ScalarArray: r""" Evaluate :math:`\langle \chi^2 \rangle` for each observed/predicted @@ -56,10 +57,11 @@ def mean_chi_squared( The actual measured images. images_predicted The images predicted by the inversion. + uncertainty + The standard deviation of the measurement noise in the predicted + images. """ - uncertainty = self.instrument.uncertainty(images_predicted) - return ctis.inverters.merit.mean_chi_squared( observed=images_observed, expected=images_predicted, diff --git a/ctis/inverters/_iterative/_mart/_mart.py b/ctis/inverters/_iterative/_mart/_mart.py index e45c41b..77729f0 100644 --- a/ctis/inverters/_iterative/_mart/_mart.py +++ b/ctis/inverters/_iterative/_mart/_mart.py @@ -118,9 +118,10 @@ def __call__( if verbose: # pragma: nocover print(f"{i=}") - images_new = instrument.image(scene, noise=False).outputs + predicted = instrument.image(scene, noise=False, uncertainty=True).outputs + images_new = predicted.nominal - chi2_ij = self.mean_chi_squared(images, images_new) + chi2_ij = self.mean_chi_squared(images, images_new, predicted.width) r_ij = self.correlation_residual(images, images_new) chi2.append(chi2_ij) diff --git a/ctis/inverters/_results_test.py b/ctis/inverters/_results_test.py new file mode 100644 index 0000000..e69de29 diff --git a/demo.ipynb b/demo.ipynb new file mode 100644 index 0000000..f587e0b --- /dev/null +++ b/demo.ipynb @@ -0,0 +1,38 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "initial_id", + "metadata": { + "ExecuteTime": { + "end_time": "2025-04-28T14:54:05.839586Z", + "start_time": "2025-04-28T14:54:05.836086Z" + } + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/na.ipynb b/na.ipynb new file mode 100644 index 0000000..59b0c01 --- /dev/null +++ b/na.ipynb @@ -0,0 +1,205 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 13, + "id": "initial_id", + "metadata": { + "ExecuteTime": { + "start_time": "2025-05-21T21:46:24.866716Z" + }, + "jupyter": { + "is_executing": true + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import named_arrays as na" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "dc0632491b7ed3be", + "metadata": {}, + "outputs": [], + "source": [ + "a = na.linspace(0, 1, axis=\"x\", num=11)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9e8564ab", + "metadata": {}, + "outputs": [], + "source": [ + "b = na.linspace(0, 1, axis=\"y\", num=10)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "baf51200", + "metadata": {}, + "outputs": [], + "source": [ + "c = a + b" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "84026040", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'x': 11, 'y': 10}" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "1f503a04", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "named_arrays._scalars.scalars.ScalarArray" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "type(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "d9e3377a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "ScalarArray(\n", + " ndarray=[[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],\n", + " axes=('x', 'y'),\n", + ")" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "na.ScalarArray.zeros(c.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "5dfed7bc", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1, 2],\n", + " [3, 4]])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "d = np.array([[1, 2], [3, 4]])\n", + "d" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "646d8b33", + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "The number of axis names, ('a',), must match the number of dimensions, 2.", + "output_type": "error", + "traceback": [ + "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", + "\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)", + "Cell \u001B[1;32mIn[20], line 1\u001B[0m\n\u001B[1;32m----> 1\u001B[0m e \u001B[38;5;241m=\u001B[39m \u001B[43mna\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mScalarArray\u001B[49m\u001B[43m(\u001B[49m\u001B[43md\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43maxes\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43m(\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43ma\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 2\u001B[0m e\n", + "File \u001B[1;32m:5\u001B[0m, in \u001B[0;36m__init__\u001B[1;34m(self, ndarray, axes)\u001B[0m\n", + "File \u001B[1;32m~\\Kankelborg-Group\\named_arrays\\named_arrays\\_scalars\\scalars.py:798\u001B[0m, in \u001B[0;36mScalarArray.__post_init__\u001B[1;34m(self)\u001B[0m\n\u001B[0;32m 796\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes \u001B[38;5;241m=\u001B[39m (\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes, )\n\u001B[0;32m 797\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28mgetattr\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mndarray, \u001B[38;5;124m'\u001B[39m\u001B[38;5;124mndim\u001B[39m\u001B[38;5;124m'\u001B[39m, \u001B[38;5;241m0\u001B[39m) \u001B[38;5;241m!=\u001B[39m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes): \u001B[38;5;66;03m# pragma: nocover\u001B[39;00m\n\u001B[1;32m--> 798\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\n\u001B[0;32m 799\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mThe number of axis names, \u001B[39m\u001B[38;5;132;01m{\u001B[39;00m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m, \u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m 800\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mmust match the number of dimensions, \u001B[39m\u001B[38;5;132;01m{\u001B[39;00mnp\u001B[38;5;241m.\u001B[39mndim(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mndarray)\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m 801\u001B[0m )\n\u001B[0;32m 802\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes) \u001B[38;5;241m!=\u001B[39m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mset\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes)):\n\u001B[0;32m 803\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mEach axis name must be unique, got \u001B[39m\u001B[38;5;132;01m{\u001B[39;00m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;124m'\u001B[39m)\n", + "\u001B[1;31mValueError\u001B[0m: The number of axis names, ('a',), must match the number of dimensions, 2." + ] + } + ], + "source": [ + "e = na.ScalarArray(d, axes=(\"a\"))\n", + "e" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be015f1f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/na.py b/na.py new file mode 100644 index 0000000..e69de29 diff --git a/pyproject.toml b/pyproject.toml index dfa2d13..8bd551b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,6 +18,7 @@ classifiers = [ dependencies = [ "astropy", "named-arrays~=2.3", + "optika~=2.1", ] dynamic = ["version"] From a1b2bf30f7f1082057d3701e79d663bcf6424732 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 18:05:28 -0600 Subject: [PATCH 02/13] Rename to `OptikaInstrument` and generalize the instrument tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rename `LinearOptikaInstrument` to `OptikaInstrument` ("Linear" is already implied by `AbstractLinearInstrument`; the distinguishing trait is that it is backed by an `optika` system). Generalize `AbstractTestAbstractInstrument` so the shared tests derive the scene and coordinates from the instrument under test, and add `TestOptikaInstrument` with a channel-aware `LinearSystem` fixture whose scene grid spans the field of view. A new `test_backproject_conserves_flux` asserts that re-imaging a backprojection preserves the total measured electrons (the forward model's adjoint property) for every instrument. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ Co-Authored-By: Claude Opus 4.8 --- ctis/instruments/__init__.py | 4 +- ctis/instruments/_instruments.py | 4 +- ctis/instruments/_instruments_test.py | 226 +++++++++++++++++--------- 3 files changed, 156 insertions(+), 78 deletions(-) diff --git a/ctis/instruments/__init__.py b/ctis/instruments/__init__.py index 95b43f5..6f82559 100644 --- a/ctis/instruments/__init__.py +++ b/ctis/instruments/__init__.py @@ -6,12 +6,12 @@ AbstractInstrument, AbstractLinearInstrument, IdealInstrument, - LinearOptikaInstrument, + OptikaInstrument, ) __all__ = [ "AbstractInstrument", "AbstractLinearInstrument", "IdealInstrument", - "LinearOptikaInstrument", + "OptikaInstrument", ] diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 3113bd4..0284fc6 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -12,7 +12,7 @@ "AbstractInstrument", "AbstractLinearInstrument", "IdealInstrument", - "LinearOptikaInstrument", + "OptikaInstrument", ] @@ -627,7 +627,7 @@ def backproject( @dataclasses.dataclass -class LinearOptikaInstrument( +class OptikaInstrument( AbstractLinearInstrument, ): """ diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index 0f76dbf..be87825 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -3,8 +3,110 @@ import numpy as np import astropy.units as u import named_arrays as na +import optika import ctis + +def _scene( + a: ctis.instruments.AbstractInstrument, +) -> na.FunctionArray[na.AbstractSpectralPositionalVectorArray, na.AbstractScalar]: + """ + A compact spectral radiance defined on the scene grid of `a`. + + The signal occupies a single central field cell, so its dispersed image is + contained within the sensor and the backprojection is contained within the + scene grid. This lets the forward model conserve flux exactly through a + round trip. + """ + coordinates = a.coordinates_scene + axis_x, axis_y = a.axis_scene_xy + shape = na.shape(coordinates.position) + + values = np.zeros((shape[axis_x] - 1, shape[axis_y] - 1)) + values[values.shape[0] // 2, values.shape[1] // 2] = 1.0 + + radiance = 1e-11 * u.erg / u.s / u.cm**2 / u.arcsec**2 / u.nm + outputs = na.ScalarArray(values, axes=(axis_x, axis_y)) * radiance + return na.FunctionArray(inputs=coordinates, outputs=outputs) + + +class AbstractTestAbstractInstrument( + abc.ABC, +): + + def test_image( + self, + a: ctis.instruments.AbstractInstrument, + ): + scene = _scene(a) + result = a.image(scene.outputs, noise=False) + assert isinstance(result, na.FunctionArray) + assert np.all(result.inputs.position == a.coordinates_sensor.position) + assert result.outputs.sum() > 0 + + def test_backproject( + self, + a: ctis.instruments.AbstractInstrument, + ): + scene = _scene(a) + image = a.image(scene.outputs, noise=False) + result = a.backproject(image) + + assert isinstance(result, na.FunctionArray) + assert np.all(result.inputs == a.coordinates_scene) + assert np.all(np.isfinite(na.as_named_array(result.outputs).value)) + assert result.outputs.sum() > 0 + + def test_backproject_conserves_flux( + self, + a: ctis.instruments.AbstractInstrument, + ): + # backproject is the adjoint of the forward model, so for a scene + # contained within the field of view, re-imaging a backprojection + # preserves the total number of measured electrons. + scene = _scene(a) + image = a.image(scene.outputs, noise=False) + result = a.backproject(image) + image_check = a.image(result, noise=False) + assert np.allclose( + image.outputs.sum().ndarray.value, + image_check.outputs.sum().ndarray.value, + ) + + def test_num_channel( + self, + a: ctis.instruments.AbstractInstrument, + ): + result = a.num_channel + + assert isinstance(result, int) + + @pytest.mark.parametrize("integrate", [False, True]) + def test_image_uncertainty( + self, + a: ctis.instruments.AbstractInstrument, + integrate: bool, + ): + scene = _scene(a) + result = a.image( + scene.outputs, + integrate=integrate, + noise=False, + uncertainty=True, + ) + + # the measurement noise is attached as a normal uncertain array + assert isinstance(result.outputs, na.NormalUncertainScalarArray) + assert result.outputs.width.unit.is_equivalent(u.electron) + assert np.all(result.outputs.width >= 0 * u.electron) + + +class AbstractTestAbstractLinearInstrument( + AbstractTestAbstractInstrument, +): + pass + + velocity = na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s wavelength_rest = 171 * u.AA @@ -32,8 +134,6 @@ position=position_sensor, ) -gaussians = ctis.scenes.gaussians(coordinates_scene) - AA = dict( unit=u.AA, equivalencies=u.doppler_optical(wavelength_rest), @@ -65,88 +165,66 @@ ) -class AbstractTestAbstractInstrument( - abc.ABC, +@pytest.mark.parametrize( + argnames="a", + argvalues=[instrument_ideal], +) +class TestIdealInstrument( + AbstractTestAbstractLinearInstrument, ): + pass - @pytest.mark.parametrize( - argnames="scene", - argvalues=[ - gaussians.outputs, - ], - ) - def test_image( - self, - a: ctis.instruments.AbstractInstrument, - scene: na.AbstractScalar | na.AbstractFunctionArray, - ): - result = a.image(scene) - assert np.all(result.inputs.position == coordinates_sensor.position) - assert result.outputs.sum() > 0 - @pytest.mark.parametrize( - argnames="image", - argvalues=[ - instrument_ideal.image(gaussians, noise=False), - ], +def _instrument_optika() -> ctis.instruments.OptikaInstrument: + channel = na.linspace(0, 360, axis="channel", num=3, endpoint=False) * u.deg + system = optika.systems.LinearSystem( + area_effective=optika.radiometry.InterpolatedEffectiveAreaModel( + wavelength=na.linspace(400, 700, axis="wavelength", num=10) * u.nm, + area=na.linspace(1, 2, axis="wavelength", num=10) * u.cm**2, + axis_wavelength="wavelength", + ), + distortion=optika.distortion.SimpleDistortionModel( + plate_scale=50 * u.arcsec / u.mm, + dispersion=250 * u.nm / u.mm, + angle=channel, + reference=na.SpectralPositionalVectorArray( + wavelength=550 * u.nm, + position=na.Cartesian2dVectorArray(0, 0) * u.mm, + ), + ), + sensor=optika.sensors.ImagingSensor( + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("sensor_x", "sensor_y"), + timedelta_exposure=1 * u.s, + num_pixel=na.Cartesian2dVectorArray(32, 32), + ), ) - def test_backproject( - self, - a: ctis.instruments.AbstractInstrument, - image: na.AbstractScalar | na.AbstractFunctionArray, - ): - result = a.backproject(image) - - assert np.all(result.inputs == coordinates_scene) - assert result.outputs.sum() > 0 - - if isinstance(image, na.AbstractFunctionArray): - image = image.outputs - - image_check = a.image(result, noise=False).outputs - - assert np.allclose(image.sum(), image_check.sum()) - - def test_num_channel( - self, - a: ctis.instruments.AbstractInstrument, - ): - result = a.num_channel - - assert isinstance(result, int) - - @pytest.mark.parametrize( - argnames="scene", - argvalues=[ - gaussians.outputs, - ], + # the scene grid spans the field of view so the backprojection of a + # contained source is not clipped + coordinates_scene = na.SpectralPositionalVectorArray( + wavelength=na.linspace(530, 570, axis="wavelength", num=4) * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-20 * u.arcsec, + stop=+20 * u.arcsec, + axis=na.Cartesian2dVectorArray("scene_x", "scene_y"), + num=21, + ), + ) + return ctis.instruments.OptikaInstrument( + system=system, + coordinates_scene=coordinates_scene, + channel=channel, + axis_channel="channel", + axis_wavelength="wavelength", + axis_scene_xy=("scene_x", "scene_y"), ) - @pytest.mark.parametrize("integrate", [False, True]) - def test_image_uncertainty( - self, - a: ctis.instruments.AbstractInstrument, - scene: na.AbstractScalar | na.AbstractFunctionArray, - integrate: bool, - ): - result = a.image(scene, integrate=integrate, noise=False, uncertainty=True) - - # the measurement noise is attached as a normal uncertain array - assert isinstance(result.outputs, na.NormalUncertainScalarArray) - assert result.outputs.width.unit.is_equivalent(u.electron) - assert np.all(result.outputs.width >= 0 * u.electron) - - -class AbstractTestAbstractLinearInstrument( - AbstractTestAbstractInstrument, -): - pass @pytest.mark.parametrize( argnames="a", - argvalues=[instrument_ideal], + argvalues=[_instrument_optika()], ) -class TestIdealInstrument( +class TestOptikaInstrument( AbstractTestAbstractLinearInstrument, ): pass From fde855db403710d7e837bb3eff3d598cd6022c68 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 20:52:35 -0600 Subject: [PATCH 03/13] Apply read noise once per readout in the instruments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Delegate the wavelength integration to optika: `OptikaInstrument` now forwards `integrate` to `system.image_from_weights`/`backproject_from_weights` (which sum over wavelength and apply the sensor's read noise once), instead of integrating in ctis. This fixes the sqrt(N_wavelength) over-count that resulted from the sensor applying read noise per wavelength. `IdealInstrument` gains a `read_noise` field, applied once (Gaussian on the nominal, in quadrature on the uncertainty) after its own wavelength integration. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ Co-Authored-By: Claude Opus 4.8 --- ctis/instruments/_instruments.py | 55 +++++++++++++++----------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 0284fc6..3f3b00f 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -463,6 +463,12 @@ class IdealInstrument( applied after integrating over wavelength. """ + read_noise: u.Quantity | na.AbstractScalar = 0 * u.electron + """ + The standard deviation of the Gaussian read noise added to each pixel + once per readout (after integrating over wavelength), in electrons. + """ + def _shot_noise(self, image: na.ScalarArray) -> na.ScalarArray: # photon shot noise, converted back into electrons to match the # electron-valued image @@ -599,6 +605,16 @@ def image( if integrate: electrons, coordinates = self._integrate_wavelength(electrons, coordinates) + # add read noise once, at the integrated readout + if isinstance(electrons, na.NormalUncertainScalarArray): + nominal = electrons.nominal + width = np.sqrt(np.square(electrons.width) + np.square(self.read_noise)) + if noise: + nominal = na.random.normal(loc=nominal, scale=self.read_noise) + electrons = na.NormalUncertainScalarArray(nominal=nominal, width=width) + elif noise: + electrons = na.random.normal(loc=electrons, scale=self.read_noise) + return na.FunctionArray( inputs=coordinates, outputs=electrons, @@ -724,30 +740,17 @@ def image( else: scene = na.FunctionArray(inputs=self.coordinates_scene, outputs=scene) - # apply the optika forward model (effective area, vignetting, and the - # sensor response) using the cached regridding weights, giving the - # electrons measured in each pixel (optionally with their uncertainty). - result = self.system.image_from_weights( + # optika applies the effective area, vignetting, and sensor response + # (including the integration over wavelength and read noise) using the + # cached regridding weights. + return self.system.image_from_weights( self.weights, scene, axis_wavelength=self.axis_wavelength, axis_field=self.axis_scene_xy, noise=noise, uncertainty=uncertainty, - ) - - coordinates = result.inputs - values_output = result.outputs - - if integrate: - values_output, coordinates = self._integrate_wavelength( - values_output, - coordinates, - ) - - return na.FunctionArray( - inputs=coordinates, - outputs=values_output, + integrate=integrate, ) def backproject( @@ -765,25 +768,19 @@ def backproject( else: values = image - axis_wavelength = self.axis_wavelength - num_wavelength = self.coordinates_scene.wavelength.shape[axis_wavelength] - 1 - - if integrate: - values = values / num_wavelength - - # rebuild the detector image over the full sensor wavelength grid so - # optika can invert the sensor response for each wavelength. + # rebuild the detector image over the full sensor wavelength grid; optika + # spreads the integrated readout back over wavelength and inverts the + # sensor response with the cached transpose weights. image = na.FunctionArray( inputs=self.coordinates_sensor, outputs=values, ) - # apply the transposed optika model with the cached transpose weights, - # inverting the sensor response and recovering the spectral radiance. return self.system.backproject_from_weights( self.weights_transpose, image, coordinates=self.coordinates_scene, - axis_wavelength=axis_wavelength, + axis_wavelength=self.axis_wavelength, axis_field=self.axis_scene_xy, + integrate=integrate, ) From 4b54058a0a61f0e21118a7a0b43469d2159432ac Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 20:57:47 -0600 Subject: [PATCH 04/13] Test that read noise is applied once per readout MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add a shared `test_read_noise` (with a per-instrument `_with_read_noise` hook) verifying that a nonzero read noise raises the integrated uncertainty by exactly `read_noise` in quadrature, guarding against the sqrt(N_wavelength) over-count. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ Co-Authored-By: Claude Opus 4.8 --- ctis/instruments/_instruments_test.py | 33 +++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index be87825..0899155 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -1,5 +1,6 @@ import pytest import abc +import dataclasses import numpy as np import astropy.units as u import named_arrays as na @@ -100,6 +101,28 @@ def test_image_uncertainty( assert result.outputs.width.unit.is_equivalent(u.electron) assert np.all(result.outputs.width >= 0 * u.electron) + @abc.abstractmethod + def _with_read_noise( + self, + a: ctis.instruments.AbstractInstrument, + read_noise: u.Quantity, + ) -> ctis.instruments.AbstractInstrument: + """Return a copy of `a` with the given per-readout read noise.""" + + def test_read_noise( + self, + a: ctis.instruments.AbstractInstrument, + ): + # read noise is added once per readout, so it raises the integrated + # uncertainty by exactly `read_noise` in quadrature (not by + # sqrt(num_wavelength) * read_noise) + scene = _scene(a) + a_rn = self._with_read_noise(a, 10 * u.electron) + width_0 = a.image(scene.outputs, noise=False, uncertainty=True).outputs.width + width_1 = a_rn.image(scene.outputs, noise=False, uncertainty=True).outputs.width + contribution = np.sqrt(np.square(width_1) - np.square(width_0)) + assert np.allclose(contribution.to_value(u.electron), 10) + class AbstractTestAbstractLinearInstrument( AbstractTestAbstractInstrument, @@ -172,7 +195,8 @@ class AbstractTestAbstractLinearInstrument( class TestIdealInstrument( AbstractTestAbstractLinearInstrument, ): - pass + def _with_read_noise(self, a, read_noise): + return dataclasses.replace(a, read_noise=read_noise) def _instrument_optika() -> ctis.instruments.OptikaInstrument: @@ -227,4 +251,9 @@ def _instrument_optika() -> ctis.instruments.OptikaInstrument: class TestOptikaInstrument( AbstractTestAbstractLinearInstrument, ): - pass + def _with_read_noise(self, a, read_noise): + sensor = dataclasses.replace(a.system.sensor, read_noise=read_noise) + return dataclasses.replace( + a, + system=dataclasses.replace(a.system, sensor=sensor), + ) From 22f716c7cf125f3c98d37ebd3ad988c92dc4bada Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 21:23:50 -0600 Subject: [PATCH 05/13] Use pixel-coordinate distortion in the OptikaInstrument fixture MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `LinearSystem.coordinates_sensor` is now in pixel units, so the test fixture's distortion is expressed in arcsec/pix and nm/pix to match. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ Co-Authored-By: Claude Opus 4.8 --- ctis/instruments/_instruments_test.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index 0899155..b8ba1a7 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -208,12 +208,12 @@ def _instrument_optika() -> ctis.instruments.OptikaInstrument: axis_wavelength="wavelength", ), distortion=optika.distortion.SimpleDistortionModel( - plate_scale=50 * u.arcsec / u.mm, - dispersion=250 * u.nm / u.mm, + plate_scale=0.75 * u.arcsec / u.pix, + dispersion=3.75 * u.nm / u.pix, angle=channel, reference=na.SpectralPositionalVectorArray( wavelength=550 * u.nm, - position=na.Cartesian2dVectorArray(0, 0) * u.mm, + position=na.Cartesian2dVectorArray(16, 16) * u.pix, ), ), sensor=optika.sensors.ImagingSensor( From f3210c9a23eccd70ffed22b12b54c888471c9b48 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 27 Jul 2026 14:50:28 -0600 Subject: [PATCH 06/13] Add a `unit` parameter to `backproject` to match the scene units The forward model accepts a scene in photon or energy units, so let `backproject` return the spectral radiance in whichever the caller requests, so an inversion can be expressed in the same units as its input scene: instrument.backproject(image, unit=scene.outputs.unit) `OptikaInstrument` forwards `unit` to optika's `backproject`, where the wavelength and sensor physics live. `IdealInstrument`, a standalone analytic model with no optika system, converts locally: its backprojection is in energy units and is divided by the energy per photon when a photon unit is requested. The default (`unit=None`) preserves the previous behavior of each instrument. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/instruments/_instruments.py | 38 ++++++++++++++++++++++++++- ctis/instruments/_instruments_test.py | 31 ++++++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 3f3b00f..fdb7586 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -78,6 +78,7 @@ def backproject( self, image: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, + unit: None | u.UnitBase = None, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ The backward model of this CTIS instrument, which maps photons measured @@ -96,6 +97,14 @@ def backproject( in units of photons. integrate Complement of the `integrate` keyword of :meth:`image`. + unit + The unit of the backprojected spectral radiance. + The forward model, :meth:`image`, accepts a scene in either photon + or energy units, so the backprojection is expressed in whichever the + caller requests, converting between photon and energy units using + the energy per photon. + If :obj:`None` (the default), the radiance is left in the natural + units of the backprojection and is not converted. """ @property @@ -620,10 +629,32 @@ def image( outputs=electrons, ) + def _to_unit( + self, + radiance: na.AbstractScalar, + unit: None | u.UnitBase, + ) -> na.AbstractScalar: + """ + Express the (energy) backprojected radiance in `unit`, dividing by the + energy per photon (:attr:`_energy_per_photon`) when `unit` is a photon + rather than an energy unit. A `unit` compatible with neither raises a + :class:`~astropy.units.UnitConversionError`. + """ + if unit is None: + return radiance + + if unit.is_equivalent(na.unit_normalized(radiance)): + return radiance.to(unit) + + # the natural (energy) radiance is divided by the energy per photon to + # express it in photon units. + return (radiance / self._energy_per_photon).to(unit) + def backproject( self, image: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, + unit: None | u.UnitBase = None, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: # convert the measured electrons back into photons @@ -639,7 +670,9 @@ def backproject( result = result / (self.area_effective * self.timedelta_exposure) - return result + # the super() backprojection is in energy units; express it in the + # requested unit, converting to photons if necessary. + return result.replace(outputs=self._to_unit(result.outputs, unit)) @dataclasses.dataclass @@ -757,6 +790,7 @@ def backproject( self, image: na.AbstractScalar | na.AbstractFunctionArray, integrate: bool = True, + unit: None | u.UnitBase = None, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: if isinstance(image, na.AbstractFunctionArray): @@ -776,6 +810,7 @@ def backproject( outputs=values, ) + # the energy/photon conversion is handled by optika's backproject. return self.system.backproject_from_weights( self.weights_transpose, image, @@ -783,4 +818,5 @@ def backproject( axis_wavelength=self.axis_wavelength, axis_field=self.axis_scene_xy, integrate=integrate, + unit=unit, ) diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index b8ba1a7..a3656aa 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -3,6 +3,7 @@ import dataclasses import numpy as np import astropy.units as u +import astropy.constants import named_arrays as na import optika import ctis @@ -58,6 +59,36 @@ def test_backproject( assert np.all(np.isfinite(na.as_named_array(result.outputs).value)) assert result.outputs.sum() > 0 + def test_backproject_unit( + self, + a: ctis.instruments.AbstractInstrument, + ): + scene = _scene(a) + image = a.image(scene.outputs, noise=False) + + # the backprojection can be expressed in the (energy) units of the + # input scene, or in photon units, to match the forward model. + unit_energy = na.unit_normalized(scene.outputs) + unit_photon = u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm + + result_energy = a.backproject(image, unit=unit_energy) + result_photon = a.backproject(image, unit=unit_photon) + + assert na.unit_normalized(result_energy.outputs).is_equivalent(unit_energy) + assert na.unit_normalized(result_photon.outputs).is_equivalent(unit_photon) + + # the two expressions differ only by the energy per photon + wavelength = a.coordinates_scene.wavelength.cell_centers(a.axis_wavelength) + energy_per_photon = ( + astropy.constants.h * astropy.constants.c / wavelength / u.photon + ) + expected_energy = (result_photon.outputs * energy_per_photon).to(unit_energy) + assert np.allclose(result_energy.outputs, expected_energy) + + # a unit compatible with neither photon nor energy radiance is rejected + with pytest.raises(u.UnitConversionError): + a.backproject(image, unit=u.s) + def test_backproject_conserves_flux( self, a: ctis.instruments.AbstractInstrument, From ffb21b315ff9b1d9453ed7c3e159a0309dfa37fa Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 28 Jul 2026 15:51:18 -0600 Subject: [PATCH 07/13] Use `SpectralPositionalVectorArray.volume_cell` in `_volume_scene` Replace the hand-written `wavelength.volume_cell(...) * position.volume_cell(...)` (with the manual cell-center alignment) in `AbstractLinearInstrument._volume_scene` with `coordinates.volume_cell(...)` from named-arrays 2.3. The scene coordinates are coerced to a spectral-positional vector first so that `IdealInstrument`'s Doppler scene also gets the method. Bump the named-arrays pin to ~=2.3. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/instruments/_instruments.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index fdb7586..654d73c 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -224,15 +224,14 @@ def _volume_scene(self) -> na.AbstractScalar: """ coords = self.coordinates_scene - dw = coords.wavelength.volume_cell(self.axis_wavelength) - - dA = coords.position.volume_cell(self.axis_scene_xy) - dA = na.as_named_array(dA) - dA = dA.cell_centers(self.axis_wavelength) - - dV = dw * dA + # coerce to a spectral-positional vector so a Doppler scene (as used by + # `IdealInstrument`) also gets `volume_cell`. + coords = na.SpectralPositionalVectorArray( + wavelength=coords.wavelength, + position=coords.position, + ) - return dV + return coords.volume_cell((self.axis_wavelength, *self.axis_scene_xy)) @property def _energy_per_photon(self) -> u.Quantity | na.AbstractScalar: From 801315a7b9f2e9fdce4fa4ef7b995515d925e50e Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 30 Jul 2026 12:07:24 -0600 Subject: [PATCH 08/13] Remove accidentally-committed scratch files `na.py`, `demo.ipynb`, and `na.ipynb` are root-level scratch files that were committed by mistake; they are not part of the `ctis` package. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- demo.ipynb | 38 ---------- na.ipynb | 205 ----------------------------------------------------- na.py | 0 3 files changed, 243 deletions(-) delete mode 100644 demo.ipynb delete mode 100644 na.ipynb delete mode 100644 na.py diff --git a/demo.ipynb b/demo.ipynb deleted file mode 100644 index f587e0b..0000000 --- a/demo.ipynb +++ /dev/null @@ -1,38 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "initial_id", - "metadata": { - "ExecuteTime": { - "end_time": "2025-04-28T14:54:05.839586Z", - "start_time": "2025-04-28T14:54:05.836086Z" - } - }, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "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.11.1" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/na.ipynb b/na.ipynb deleted file mode 100644 index 59b0c01..0000000 --- a/na.ipynb +++ /dev/null @@ -1,205 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 13, - "id": "initial_id", - "metadata": { - "ExecuteTime": { - "start_time": "2025-05-21T21:46:24.866716Z" - }, - "jupyter": { - "is_executing": true - } - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "import named_arrays as na" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "dc0632491b7ed3be", - "metadata": {}, - "outputs": [], - "source": [ - "a = na.linspace(0, 1, axis=\"x\", num=11)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "9e8564ab", - "metadata": {}, - "outputs": [], - "source": [ - "b = na.linspace(0, 1, axis=\"y\", num=10)" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "baf51200", - "metadata": {}, - "outputs": [], - "source": [ - "c = a + b" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "id": "84026040", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{'x': 11, 'y': 10}" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "c.shape" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "1f503a04", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "named_arrays._scalars.scalars.ScalarArray" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "type(a)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "d9e3377a", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "ScalarArray(\n", - " ndarray=[[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", - " [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],\n", - " axes=('x', 'y'),\n", - ")" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "na.ScalarArray.zeros(c.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "id": "5dfed7bc", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[1, 2],\n", - " [3, 4]])" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "d = np.array([[1, 2], [3, 4]])\n", - "d" - ] - }, - { - "cell_type": "code", - "execution_count": 20, - "id": "646d8b33", - "metadata": {}, - "outputs": [ - { - "ename": "ValueError", - "evalue": "The number of axis names, ('a',), must match the number of dimensions, 2.", - "output_type": "error", - "traceback": [ - "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", - "\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)", - "Cell \u001B[1;32mIn[20], line 1\u001B[0m\n\u001B[1;32m----> 1\u001B[0m e \u001B[38;5;241m=\u001B[39m \u001B[43mna\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mScalarArray\u001B[49m\u001B[43m(\u001B[49m\u001B[43md\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43maxes\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43m(\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43ma\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 2\u001B[0m e\n", - "File \u001B[1;32m:5\u001B[0m, in \u001B[0;36m__init__\u001B[1;34m(self, ndarray, axes)\u001B[0m\n", - "File \u001B[1;32m~\\Kankelborg-Group\\named_arrays\\named_arrays\\_scalars\\scalars.py:798\u001B[0m, in \u001B[0;36mScalarArray.__post_init__\u001B[1;34m(self)\u001B[0m\n\u001B[0;32m 796\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes \u001B[38;5;241m=\u001B[39m (\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes, )\n\u001B[0;32m 797\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28mgetattr\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mndarray, \u001B[38;5;124m'\u001B[39m\u001B[38;5;124mndim\u001B[39m\u001B[38;5;124m'\u001B[39m, \u001B[38;5;241m0\u001B[39m) \u001B[38;5;241m!=\u001B[39m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes): \u001B[38;5;66;03m# pragma: nocover\u001B[39;00m\n\u001B[1;32m--> 798\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\n\u001B[0;32m 799\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mThe number of axis names, \u001B[39m\u001B[38;5;132;01m{\u001B[39;00m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m, \u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m 800\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mmust match the number of dimensions, \u001B[39m\u001B[38;5;132;01m{\u001B[39;00mnp\u001B[38;5;241m.\u001B[39mndim(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mndarray)\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m 801\u001B[0m )\n\u001B[0;32m 802\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes) \u001B[38;5;241m!=\u001B[39m \u001B[38;5;28mlen\u001B[39m(\u001B[38;5;28mset\u001B[39m(\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes)):\n\u001B[0;32m 803\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mEach axis name must be unique, got \u001B[39m\u001B[38;5;132;01m{\u001B[39;00m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39maxes\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;124m'\u001B[39m)\n", - "\u001B[1;31mValueError\u001B[0m: The number of axis names, ('a',), must match the number of dimensions, 2." - ] - } - ], - "source": [ - "e = na.ScalarArray(d, axes=(\"a\"))\n", - "e" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "be015f1f", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "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.11.1" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/na.py b/na.py deleted file mode 100644 index e69de29..0000000 From 7562d16b38a6db42682971ed087d49e9bb850284 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 11:18:35 -0600 Subject: [PATCH 09/13] Use `.spectral_positional` to accept Doppler grids Replace the hand-written `na.SpectralPositionalVectorArray(wavelength=..., position=...)` coercions in `ctis.regrid` and `AbstractLinearInstrument. _volume_scene` with the `.spectral_positional` property (named-arrays 2.4.0), and bump the pin. Add `test_regrid_doppler`, covering a Doppler-grid regrid. Behavior is unchanged. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/_regrid.py | 5 +++ ctis/_regrid_test.py | 54 ++++++++++++++++++++++++++++++++ ctis/instruments/_instruments.py | 11 ++----- pyproject.toml | 2 +- 4 files changed, 63 insertions(+), 9 deletions(-) diff --git a/ctis/_regrid.py b/ctis/_regrid.py index a7d7d68..bf96856 100644 --- a/ctis/_regrid.py +++ b/ctis/_regrid.py @@ -151,6 +151,11 @@ def regrid( axis = (axis_wavelength, *axis_position) + # normalize to plain spectral-positional vectors so grids defined on a + # Doppler vector (as produced by the `IdealInstrument`) are accepted. + coordinates_input = coordinates_input.spectral_positional + coordinates_output = coordinates_output.spectral_positional + # weight each input voxel by its volume so the conservative regridding # preserves the integral of the field rather than the per-voxel sum. values = values_input * coordinates_input.volume_cell(axis) diff --git a/ctis/_regrid_test.py b/ctis/_regrid_test.py index 55b5711..856c503 100644 --- a/ctis/_regrid_test.py +++ b/ctis/_regrid_test.py @@ -17,6 +17,22 @@ def _grid(num_wavelength: int, num_position: int) -> na.SpectralPositionalVector ) +def _grid_doppler( + num_wavelength: int, + num_position: int, +) -> na.DopplerPositionalVectorArray: + return na.DopplerPositionalVectorArray( + wavelength=na.linspace(500, 600, axis="wavelength", num=num_wavelength) * u.nm, + wavelength_rest=550 * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-10 * u.arcsec, + stop=+10 * u.arcsec, + axis=na.Cartesian2dVectorArray("x", "y"), + num=num_position, + ), + ) + + @pytest.mark.parametrize( argnames="coordinates_input,coordinates_output", argvalues=[ @@ -69,3 +85,41 @@ def test_regrid( integral_output = (result * coordinates_output.volume_cell(axis)).sum() ratio = float((integral_output / integral_input).ndarray) assert np.isclose(ratio, 1, rtol=0.05) + + +def test_regrid_doppler(): + # grids defined on a Doppler vector (as produced by the `IdealInstrument`) + # are accepted, normalized to spectral-positional for the volume weighting. + axis_wavelength = "wavelength" + axis_position = ("x", "y") + coordinates_input = _grid_doppler(9, 13) + coordinates_output = _grid_doppler(5, 7) + + num_input = coordinates_input.shape + values_input = na.random.uniform( + low=0, + high=1, + shape_random={ + axis_wavelength: num_input[axis_wavelength] - 1, + axis_position[0]: num_input[axis_position[0]] - 1, + axis_position[1]: num_input[axis_position[1]] - 1, + }, + ) + + result = ctis.regrid( + coordinates_input=coordinates_input, + coordinates_output=coordinates_output, + values_input=values_input, + axis_wavelength=axis_wavelength, + axis_position=axis_position, + ) + + assert isinstance(result, na.AbstractScalarArray) + assert np.all(np.isfinite(result)) + + num_output = coordinates_output.shape + assert na.shape(result) == { + axis_wavelength: num_output[axis_wavelength] - 1, + axis_position[0]: num_output[axis_position[0]] - 1, + axis_position[1]: num_output[axis_position[1]] - 1, + } diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 654d73c..723f735 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -222,14 +222,9 @@ def _volume_scene(self) -> na.AbstractScalar: """ The volume of each voxel in :attr:`coordinates_scene`. """ - coords = self.coordinates_scene - - # coerce to a spectral-positional vector so a Doppler scene (as used by - # `IdealInstrument`) also gets `volume_cell`. - coords = na.SpectralPositionalVectorArray( - wavelength=coords.wavelength, - position=coords.position, - ) + # `.spectral_positional` accepts a Doppler scene (as used by + # `IdealInstrument`), which does not implement `volume_cell` directly. + coords = self.coordinates_scene.spectral_positional return coords.volume_cell((self.axis_wavelength, *self.axis_scene_xy)) diff --git a/pyproject.toml b/pyproject.toml index 8bd551b..6c1b3e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,7 @@ classifiers = [ ] dependencies = [ "astropy", - "named-arrays~=2.3", + "named-arrays~=2.4", "optika~=2.1", ] dynamic = ["version"] From 4b0491ca2ae5c49a6c08f00e9d8a2234f9e7e4a7 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 16:03:20 -0600 Subject: [PATCH 10/13] Address review findings: remove stray test stub, robust band edges - Remove the empty `ctis/inverters/_results_test.py` that was committed by accident. - `_integrate_wavelength`: collapse the wavelength coordinates to their min/max rather than the first/last edge, so the reported band edges are correct for a non-ascending wavelength grid. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/instruments/_instruments.py | 4 ++-- ctis/inverters/_results_test.py | 0 2 files changed, 2 insertions(+), 2 deletions(-) delete mode 100644 ctis/inverters/_results_test.py diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 723f735..83f01e5 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -269,8 +269,8 @@ def _integrate_wavelength( coordinates = coordinates.replace( wavelength=na.stack( arrays=[ - coordinates.wavelength[{axis: +0}], - coordinates.wavelength[{axis: ~0}], + coordinates.wavelength.min(axis), + coordinates.wavelength.max(axis), ], axis=axis, ) diff --git a/ctis/inverters/_results_test.py b/ctis/inverters/_results_test.py deleted file mode 100644 index e69de29..0000000 From 9e2aeba2fa1f8e4ce437ccae4016e00f8e549f8a Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 16:41:16 -0600 Subject: [PATCH 11/13] Compute IdealInstrument shot-noise width from the expected signal When `noise=True`, the parent forward model returns a Poisson realization, so deriving the shot-noise width from it made the attached uncertainty depend on the noise draw (e.g. a pixel that realized zero photons got a width of zero). Recompute the noiseless image for the width when a realization was drawn, so the uncertainty is deterministic and matches the `noise=False` case. Add a test asserting the width is unchanged by noise. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/instruments/_instruments.py | 11 ++++++++++- ctis/instruments/_instruments_test.py | 13 +++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 83f01e5..5663d34 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -602,7 +602,16 @@ def image( coordinates = result.inputs if uncertainty: - width = self._shot_noise(electrons) + # the shot-noise width uses the *expected* electrons; recompute the + # noiseless image when `noise` replaced them with a realization. + if noise: + expected = ( + super().image(scene=scene, integrate=False, noise=False).outputs + * self.quantum_yield + ) + else: + expected = electrons + width = self._shot_noise(expected) electrons = na.NormalUncertainScalarArray(nominal=electrons, width=width) if integrate: diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index a3656aa..f25a3c2 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -132,6 +132,19 @@ def test_image_uncertainty( assert result.outputs.width.unit.is_equivalent(u.electron) assert np.all(result.outputs.width >= 0 * u.electron) + # the width is deterministic: it is computed from the expected signal, + # so drawing a noise realization must not change it. + result_noisy = a.image( + scene.outputs, + integrate=integrate, + noise=True, + uncertainty=True, + ) + assert np.allclose( + result_noisy.outputs.width.to_value(u.electron), + result.outputs.width.to_value(u.electron), + ) + @abc.abstractmethod def _with_read_noise( self, From 9391a8ad4af2c4ae11d008892036748de10539b8 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 2 Aug 2026 10:05:55 -0600 Subject: [PATCH 12/13] optika intersphinx --- docs/conf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 741929b..85629cd 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -102,5 +102,6 @@ 'numpy': ('https://numpy.org/doc/stable/', None), 'matplotlib': ('https://matplotlib.org/stable', None), 'astropy': ('https://docs.astropy.org/en/stable/', None), - 'named_arrays': ('https://named-arrays.readthedocs.io/en/stable/', None) + 'named_arrays': ('https://named-arrays.readthedocs.io/en/stable/', None), + 'optika': ('https://optika.readthedocs.io/en/stable/', None), } From 8fd5f2747977d225a8b4472a17b6f545de81f0a3 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sun, 2 Aug 2026 11:58:42 -0600 Subject: [PATCH 13/13] Cover the last uncovered lines in the instruments module - Drop the dead `integrate` handling from `AbstractLinearInstrument.image`: both concrete subclasses override `image` and integrate themselves, and `IdealInstrument` always calls `super().image(integrate=False)`, so the parent branch was unreachable. The parent is now a per-wavelength photon forward model with integration left to the subclass. - Add `test_axis_sensor_xy` and `test_backproject_scalar` to exercise `OptikaInstrument.axis_sensor_xy` and the bare-scalar `backproject` path. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- ctis/instruments/_instruments.py | 14 ++++---------- ctis/instruments/_instruments_test.py | 23 +++++++++++++++++++++++ 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index 5663d34..7cb2747 100644 --- a/ctis/instruments/_instruments.py +++ b/ctis/instruments/_instruments.py @@ -281,9 +281,11 @@ def _integrate_wavelength( def image( self, scene: na.AbstractScalar | na.AbstractFunctionArray, - integrate: bool = True, noise: bool = True, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + # this low-level forward model always returns the per-wavelength photons; + # integration over wavelength (and any conversion to electrons) is done + # by the concrete subclass's `image`. if isinstance(scene, na.AbstractFunctionArray): if not np.all(scene.inputs == self.coordinates_scene): @@ -310,12 +312,6 @@ def image( coordinates = self.coordinates_sensor - if integrate: - values_output, coordinates = self._integrate_wavelength( - values_output, - coordinates, - ) - return na.FunctionArray( inputs=coordinates, outputs=values_output, @@ -593,7 +589,6 @@ def image( # wavelength before integrating. result = super().image( scene=scene, - integrate=False, noise=noise, ) @@ -606,8 +601,7 @@ def image( # noiseless image when `noise` replaced them with a realization. if noise: expected = ( - super().image(scene=scene, integrate=False, noise=False).outputs - * self.quantum_yield + super().image(scene=scene, noise=False).outputs * self.quantum_yield ) else: expected = electrons diff --git a/ctis/instruments/_instruments_test.py b/ctis/instruments/_instruments_test.py index f25a3c2..d706c37 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -59,6 +59,19 @@ def test_backproject( assert np.all(np.isfinite(na.as_named_array(result.outputs).value)) assert result.outputs.sum() > 0 + def test_backproject_scalar( + self, + a: ctis.instruments.AbstractInstrument, + ): + # backproject also accepts the bare image outputs, not just a + # `FunctionArray` + scene = _scene(a) + image = a.image(scene.outputs, noise=False) + result = a.backproject(image.outputs) + + assert isinstance(result, na.FunctionArray) + assert np.all(np.isfinite(na.as_named_array(result.outputs).value)) + def test_backproject_unit( self, a: ctis.instruments.AbstractInstrument, @@ -113,6 +126,16 @@ def test_num_channel( assert isinstance(result, int) + def test_axis_sensor_xy( + self, + a: ctis.instruments.AbstractInstrument, + ): + result = a.axis_sensor_xy + + assert isinstance(result, tuple) + assert len(result) == 2 + assert all(isinstance(ax, str) for ax in result) + @pytest.mark.parametrize("integrate", [False, True]) def test_image_uncertainty( self,