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/__init__.py b/ctis/instruments/__init__.py index bcc5400..6f82559 100644 --- a/ctis/instruments/__init__.py +++ b/ctis/instruments/__init__.py @@ -6,10 +6,12 @@ AbstractInstrument, AbstractLinearInstrument, IdealInstrument, + OptikaInstrument, ) __all__ = [ "AbstractInstrument", "AbstractLinearInstrument", "IdealInstrument", + "OptikaInstrument", ] diff --git a/ctis/instruments/_instruments.py b/ctis/instruments/_instruments.py index f65ae17..7cb2747 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", + "OptikaInstrument", ] @@ -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 @@ -68,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 @@ -86,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 @@ -106,14 +125,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): @@ -211,17 +222,11 @@ def _volume_scene(self) -> na.AbstractScalar: """ The volume of each voxel in :attr:`coordinates_scene`. """ - coords = self.coordinates_scene + # `.spectral_positional` accepts a Doppler scene (as used by + # `IdealInstrument`), which does not implement `volume_cell` directly. + coords = self.coordinates_scene.spectral_positional - 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 - - return dV + return coords.volume_cell((self.axis_wavelength, *self.axis_scene_xy)) @property def _energy_per_photon(self) -> u.Quantity | na.AbstractScalar: @@ -238,12 +243,49 @@ 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.min(axis), + coordinates.wavelength.max(axis), + ], + axis=axis, + ) + ) + + return outputs, coordinates + 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): @@ -270,22 +312,6 @@ 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, - ) - ) - return na.FunctionArray( inputs=coordinates, outputs=values_output, @@ -381,7 +407,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 +453,27 @@ 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}`. + + For this idealized instrument the quantum yield is a constant, so it can be + applied after integrating over wavelength. + """ - return _shot_noise + 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 + 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 +580,86 @@ 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, noise=noise, ) + # convert the measured photons into electrons + electrons = result.outputs * self.quantum_yield + coordinates = result.inputs + + if uncertainty: + # 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, noise=False).outputs * self.quantum_yield + ) + else: + expected = electrons + width = self._shot_noise(expected) + electrons = na.NormalUncertainScalarArray(nominal=electrons, width=width) + + 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, + ) + + 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 + 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, @@ -562,4 +667,153 @@ 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 +class OptikaInstrument( + 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) + + # 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, + integrate=integrate, + ) + + 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): + 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 + + # 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, + ) + + # the energy/photon conversion is handled by optika's backproject. + return self.system.backproject_from_weights( + self.weights_transpose, + image, + coordinates=self.coordinates_scene, + 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 9f1ad33..d706c37 100644 --- a/ctis/instruments/_instruments_test.py +++ b/ctis/instruments/_instruments_test.py @@ -1,10 +1,202 @@ import pytest import abc +import dataclasses import numpy as np import astropy.units as u +import astropy.constants 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_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, + ): + 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, + ): + # 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) + + 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, + 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) + + # 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, + 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, +): + pass + + velocity = na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s wavelength_rest = 171 * u.AA @@ -32,8 +224,6 @@ position=position_sensor, ) -gaussians = ctis.scenes.gaussians(coordinates_scene) - AA = dict( unit=u.AA, equivalencies=u.doppler_optical(wavelength_rest), @@ -65,83 +255,72 @@ ) -class AbstractTestAbstractInstrument( - abc.ABC, +@pytest.mark.parametrize( + argnames="a", + argvalues=[instrument_ideal], +) +class TestIdealInstrument( + AbstractTestAbstractLinearInstrument, ): + def _with_read_noise(self, a, read_noise): + return dataclasses.replace(a, read_noise=read_noise) - @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=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(16, 16) * u.pix, + ), + ), + 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="image", - argvalues=[ - instrument_ideal.image(gaussians.outputs).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"), ) - def test_uncertainty( - self, - a: ctis.instruments.AbstractInstrument, - image: na.ScalarArray, - ): - result = a.uncertainty(image) - - assert np.all(result >= 0 * u.photon) - - -class AbstractTestAbstractLinearInstrument( - AbstractTestAbstractInstrument, -): - pass @pytest.mark.parametrize( argnames="a", - argvalues=[instrument_ideal], + argvalues=[_instrument_optika()], ) -class TestIdealInstrument( +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), + ) 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/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), } diff --git a/pyproject.toml b/pyproject.toml index dfa2d13..6c1b3e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,8 @@ classifiers = [ ] dependencies = [ "astropy", - "named-arrays~=2.3", + "named-arrays~=2.4", + "optika~=2.1", ] dynamic = ["version"]