Added ctis.inverters.MartInverter, an implementation of the Multiplicative Algebraic Reconstruction Technique (MART).

ctis/__init__.py

@@ -5,8 +5,10 @@
 
 from . import scenes
 from . import instruments
+from . import inverters
 
 __all__ = [
     "scenes",
     "instruments",
+    "inverters",
 ]

ctis/instruments/_instruments.py

@@ -114,12 +114,19 @@ def uncertainty(self) -> Callable[[na.ScalarArray], na.ScalarArray]:
         for a given number of photons.
         """
 
+    @property
+    @abc.abstractmethod
+    def channel(self):
+        """
+        Human-readable name of each independent CTIS channel.
+        """
+
     @property
     @abc.abstractmethod
     def axis_channel(self) -> str | tuple[str, ...]:
         """
         The logical axis or axes of this instrument corresponding to
-        the different dispersion magnitudes and angles.
+        the different CTIS channels.
         """
 
     @property
@@ -391,6 +398,11 @@ class IdealInstrument(
     A grid of wavelength and position coordinates on the sensor plane.
     """
 
+    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 this instrument corresponding to

ctis/instruments/_instruments_test.py

@@ -40,16 +40,19 @@
 dispersion = 200 * u.km / u.s
 dispersion = (dispersion.to(**AA) - wavelength_rest) / u.pix
 
+angle = na.linspace(0, 360, axis="channel", num=3, endpoint=False)
+
 instrument_ideal = ctis.instruments.IdealInstrument(
     area_effective=1 * u.cm**2,
     timedelta_exposure=10 * u.s,
     plate_scale=2 * u.arcsec / u.pix,
     dispersion=dispersion,
-    angle=na.linspace(0, 360, axis="channel", num=3, endpoint=False),
+    angle=angle,
     wavelength_ref=wavelength_rest,
     position_ref=32 * u.pix,
     coordinates_scene=coordinates_scene,
     coordinates_sensor=coordinates_sensor,
+    channel=angle,
     axis_channel="channel",
     axis_wavelength="wavelength",
     axis_scene_xy=("scene_x", "scene_y"),

ctis/inverters/__init__.py

@@ -0,0 +1,19 @@
+"""Inversion algorithms which can reconstruct scenes from observed images."""
+
+from . import merit
+from ._results import InversionResult
+from ._inverters import AbstractInverter
+from ._iterative import (
+    AbstractIterativeInverter,
+    MartInverter,
+    IterativeInversionResult,
+)
+
+__all__ = [
+    "merit",
+    "AbstractInverter",
+    "AbstractIterativeInverter",
+    "MartInverter",
+    "InversionResult",
+    "IterativeInversionResult",
+]

ctis/inverters/_inverters.py

@@ -0,0 +1,48 @@
+import abc
+import dataclasses
+import named_arrays as na
+import ctis
+from ._results import InversionResult
+
+__all__ = [
+    "AbstractInverter",
+]
+
+
+@dataclasses.dataclass
+class AbstractInverter(
+    abc.ABC,
+):
+    """
+    An interface describing an algorithm which can invert CTIS observations
+    to yield a reconstruction of the observed scene.
+    """
+
+    @property
+    @abc.abstractmethod
+    def instrument(self) -> ctis.instruments.AbstractInstrument:
+        """
+        A model of a CTIS instrument which transforms the radiance of an observed
+        scene to photons measured by the sensors.
+        """
+
+    @abc.abstractmethod
+    def __call__(
+        self,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> InversionResult:
+        """
+        Reconstruct a scene using the observed images.
+
+        Parameters
+        ----------
+        images
+            The observed images used to calculate the reconstruction.
+            Must be evaluated on the same coordinates as
+            :attr:`~ctis.instruments.AbstractInstrument.coordinates_sensor`
+            attribute of :attr:`instrument`.
+        kwargs
+            Additional keyword arguments which can be used by subclass
+            implementations.
+        """

ctis/inverters/_inverters_test.py

@@ -0,0 +1,31 @@
+import abc
+import numpy as np
+import named_arrays as na
+import ctis
+
+
+class AbstractTestAbstractInverter(
+    abc.ABC,
+):
+
+    def test_instrument(self, a: ctis.inverters.AbstractInverter):
+        result = a.instrument
+        assert isinstance(result, ctis.instruments.AbstractInstrument)
+
+    def test__call__(
+        self,
+        a: ctis.inverters.AbstractInverter,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> ctis.inverters.InversionResult:
+        result = a(images, **kwargs)
+
+        assert isinstance(result, ctis.inverters.InversionResult)
+
+        assert result.solution.sum() > 0
+        assert isinstance(result.success, bool)
+        assert isinstance(result.message, str)
+        assert np.all(result.images == images)
+        assert result.inverter == a
+
+        return result

ctis/inverters/_iterative/__init__.py

@@ -0,0 +1,8 @@
+from ._iterative import AbstractIterativeInverter, IterativeInversionResult
+from ._mart import MartInverter
+
+__all__ = [
+    "AbstractIterativeInverter",
+    "IterativeInversionResult",
+    "MartInverter",
+]

ctis/inverters/_iterative/_iterative.py

@@ -0,0 +1,118 @@
+from typing import ClassVar
+import abc
+import dataclasses
+import numpy as np
+import astropy.units as u
+import named_arrays as na
+import ctis
+from .. import AbstractInverter, InversionResult
+
+__all__ = [
+    "AbstractIterativeInverter",
+    "IterativeInversionResult",
+]
+
+
+@dataclasses.dataclass
+class AbstractIterativeInverter(
+    AbstractInverter,
+):
+    """
+    An abstract inversion algorithm which reconstructs an observed scene
+    using iterative methods.
+
+    These methods will apply some operation repeatedly until a specified
+    convergence criteria is met.
+    """
+
+    axis_iteration: ClassVar[str] = "iteration"
+    """The logical axis associated with changing iteration index."""
+
+    @property
+    @abc.abstractmethod
+    def num_iteration(self) -> int:
+        """
+        The maximum number of iterations to perform.
+
+        If convergence is not reached before this number is exceeded,
+        a warning is raised and an unsuccessful result is returned.
+        """
+
+    def mean_chi_squared(
+        self,
+        images_observed: na.ScalarArray,
+        images_predicted: na.ScalarArray,
+    ) -> na.ScalarArray:
+        r"""
+        Evaluate :math:`\langle \chi^2 \rangle` for each observed/predicted
+        image pair.
+
+        Parameters
+        ----------
+        images_observed
+            The actual measured images.
+        images_predicted
+            The images predicted by the inversion.
+        """
+
+        uncertainty = self.instrument.uncertainty(images_predicted)
+
+        uncertainty = np.maximum(uncertainty, 1 * u.photon)
+
+        return ctis.inverters.merit.mean_chi_squared(
+            observed=images_observed,
+            expected=images_predicted,
+            uncertainty=uncertainty,
+            axis=self.instrument.axis_sensor_xy,
+        )
+
+    def correlation_residual(
+        self,
+        images_observed: na.ScalarArray,
+        images_predicted: na.ScalarArray,
+    ) -> na.ScalarArray:
+        """
+        Evaluate the correlation between the predicted images and the residual.
+
+        Parameters
+        ----------
+        images_observed
+            The actual measured images.
+        images_predicted
+            The images predicted by the inversion.
+        """
+        return ctis.inverters.merit.correlation_residual(
+            observed=images_observed,
+            expected=images_predicted,
+            axis=self.instrument.axis_sensor_xy,
+        )
+
+
+@dataclasses.dataclass
+class IterativeInversionResult(
+    InversionResult,
+):
+    """The results of an iterative inversion attempt."""
+
+    inverter: AbstractIterativeInverter
+
+    num_iteration: int
+    """The number of iterations performed by the inverter."""
+
+    mean_chi_squared: na.ScalarArray
+    """The mean chi squared statistic for each iteration."""
+
+    correlation_residual: na.ScalarArray
+    """
+    The correlation between the predicted images and the residuals 
+    for each iteration.
+    """
+
+    @property
+    def iteration(self) -> na.ScalarArray:
+        """The iteration value for each iteration."""
+        return na.arange(
+            start=0,
+            stop=self.num_iteration,
+            axis=self.inverter.axis_iteration,
+        )

ctis/inverters/_iterative/_iterative_test.py

@@ -0,0 +1,34 @@
+import ctis
+import named_arrays as na
+from .._inverters_test import AbstractTestAbstractInverter
+
+
+class AbstractTestAbstractIterativeInverter(
+    AbstractTestAbstractInverter,
+):
+
+    def test__call__(
+        self,
+        a: ctis.inverters.AbstractIterativeInverter,
+        images: na.FunctionArray[na.SpectralPositionalVectorArray, na.ScalarArray],
+        **kwargs,
+    ) -> ctis.inverters.IterativeInversionResult:
+
+        result = super().test__call__(
+            a=a,
+            images=images,
+            **kwargs,
+        )
+
+        axis_iteration = result.inverter.axis_iteration
+
+        assert result.iteration.size == result.num_iteration
+        assert result.mean_chi_squared.shape[axis_iteration] == result.num_iteration
+        assert result.correlation_residual.shape[axis_iteration] == result.num_iteration
+
+        return result
+
+    def test_num_iteration(self, a: ctis.inverters.AbstractIterativeInverter):
+        result = a.num_iteration
+        assert isinstance(result, int)
+        assert result > 0

ctis/inverters/_iterative/_mart/__init__.py

@@ -0,0 +1,5 @@
+from ._mart import MartInverter
+
+__all__ = [
+    "MartInverter",
+]