From 5007513b01a114d5992519dd4a37af8676e1cabc Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 21 Apr 2025 21:10:14 -0600 Subject: [PATCH 1/8] Added `ctis.scenes.random_gaussians()`, a test pattern developed by Amy Winebarger. --- ctis/__init__.py | 6 + ctis/scenes/__init__.py | 9 ++ ctis/scenes/_gaussians.py | 198 +++++++++++++++++++++++++++++++++ ctis/scenes/_gaussians_test.py | 40 +++++++ 4 files changed, 253 insertions(+) create mode 100644 ctis/scenes/__init__.py create mode 100644 ctis/scenes/_gaussians.py create mode 100644 ctis/scenes/_gaussians_test.py diff --git a/ctis/__init__.py b/ctis/__init__.py index 41d44a5..8107466 100644 --- a/ctis/__init__.py +++ b/ctis/__init__.py @@ -2,3 +2,9 @@ A package for inverting imagery captured by a computed tomography imaging spectrograph. """ + +from . import scenes + +__all__ = [ + "scenes", +] diff --git a/ctis/scenes/__init__.py b/ctis/scenes/__init__.py new file mode 100644 index 0000000..dd5626c --- /dev/null +++ b/ctis/scenes/__init__.py @@ -0,0 +1,9 @@ +""" +Test patterns used to evaluate the quality of CTIS inversions. +""" + +from ._gaussians import random_gaussians + +__all__ = [ + "random_gaussians", +] diff --git a/ctis/scenes/_gaussians.py b/ctis/scenes/_gaussians.py new file mode 100644 index 0000000..b65f0d5 --- /dev/null +++ b/ctis/scenes/_gaussians.py @@ -0,0 +1,198 @@ +from typing import TypeVar +import numpy as np +import astropy.units as u +import named_arrays as na + +__all__ = [ + "random_gaussians", +] + +SpectralPositionalVectorT = TypeVar( + name="SpectralPositionalVectorT", + bound=na.AbstractSpectralPositionalVectorArray, +) + + +def _gaussian( + inputs: na.AbstractSpectralPositionalVectorArray, + amplitude: u.Quantity | na.AbstractScalar, + center: na.AbstractSpectralPositionalVectorArray, + width: na.AbstractSpectralPositionalVectorArray, +) -> na.AbstractScalar: + """ + Compute a Gaussian kernel. + + Parameters + ---------- + inputs + The grid on which to evaluate the Gaussian. + amplitude + The height of the Gausian distribution. + center + The mean of the Gaussian distribution. + width + The standard deviation of the Gaussian distribution. + """ + + arg = -np.square(((inputs - center) / width).length) / 2 + return amplitude * np.exp(arg) + + +def random_gaussians( + inputs: SpectralPositionalVectorT, + width: na.AbstractSpectralPositionalVectorArray, +) -> na.FunctionArray[SpectralPositionalVectorT, na.ScalarArray]: + r""" + A scene with a randomly-positioned Gaussian distributions originally + prepared by Amy Winebarger. + + Parameters + ---------- + axis + The names of the logical axes representing each + :math:`\lambda,x,y` coordinate, + num + The number of pixels along each coordinate: + + Examples + -------- + + Plot this scene for the input grid and Gaussian widths originally + prepared by Amy Winebarger. + + .. jupyter-execute:: + + import matplotlib.pyplot as plt + import astropy.units as u + import astropy.visualization + import named_arrays as na + import ctis + + # Define the spatial plate scale of the instrument + platescale = 0.59 * u.arcsec / u.pix + + # Define the grid of positions and velocities on which to evaluate the + # test pattern + inputs = na.SpectralPositionalVectorArray( + wavelength=na.linspace(-500, 500, axis="wavelength", num=101) * u.km / u.s, + position=na.Cartesian2dVectorLinearSpace( + start=-20 * platescale * u.pix, + stop=20 * platescale * u.pix, + axis=na.Cartesian2dVectorArray("x", "y"), + num=41, + ), + ) + + # Define the standard deviations of the Gaussians in space and velocity + width = na.SpectralPositionalVectorArray( + wavelength=27 * u.km / u.s, + position=2.4 / 2.35 * u.arcsec, + ) + + # Compute the scene of random Gaussians for the + # given input grid and standard deviations + scene = ctis.scenes.random_gaussians( + inputs=inputs, + width=width, + ) + + # Plot the result + with astropy.visualization.quantity_support(): + fig, axs = plt.subplots( + ncols=2, + gridspec_kw=dict(width_ratios=[.9,.1]), + constrained_layout=True, + ) + colorbar = na.plt.rgbmesh( + C=scene, + axis_wavelength="wavelength", + ax=axs[0], + vmin=0, + vmax=scene.outputs.max(), + ) + na.plt.pcolormesh( + C=colorbar, + axis_rgb="wavelength", + ax=axs[1], + ) + axs[1].yaxis.tick_right() + axs[1].yaxis.set_label_position("right") + """ + + amplitude = np.array( + [ + 3.38, + 6, + 4, + 3.38, + 5, + 7, + 3.7, + 5.5, + ] + ) + + center_x = np.array( + [ + -2.5, + 2.5, + 0, + 0, + 0, + -4, + 0, + 4, + ] + ) + + center_y = np.array( + [ + 0, + 0, + -2, + 2, + 4, + 4, + -4, + -4, + ] + ) + + center_v = np.array( + [ + +200, + -200, + -150, + +150, + -100, + +100, + -100, + -150, + ] + ) + + axis = "event" + + amplitude = na.ScalarArray(amplitude, axis) + center_x = na.ScalarArray(center_x, axis) << u.arcsec + center_y = na.ScalarArray(center_y, axis) << u.arcsec + center_v = na.ScalarArray(center_v, axis) << (u.km / u.s) + + center = na.SpectralPositionalVectorArray( + wavelength=center_v, + position=na.Cartesian2dVectorArray(center_x, center_y), + ) + + outputs = _gaussian( + inputs=inputs, + amplitude=amplitude, + center=center, + width=width, + ) + + outputs = outputs.sum(axis) + + return na.FunctionArray( + inputs=inputs, + outputs=outputs, + ) diff --git a/ctis/scenes/_gaussians_test.py b/ctis/scenes/_gaussians_test.py new file mode 100644 index 0000000..cd33dff --- /dev/null +++ b/ctis/scenes/_gaussians_test.py @@ -0,0 +1,40 @@ +import pytest +import astropy.units as u +import named_arrays as na +import ctis + + +@pytest.mark.parametrize( + argnames="inputs", + argvalues=[ + na.SpectralPositionalVectorArray( + wavelength=na.linspace(-500, 500, axis="wavelength", num=101) * u.km / u.s, + position=na.Cartesian2dVectorLinearSpace( + start=-10 * u.arcsec, + stop=10 * u.arcsec, + axis=na.Cartesian2dVectorArray("x", "y"), + num=41, + ), + ) + ], +) +@pytest.mark.parametrize( + argnames="width", + argvalues=[ + na.SpectralPositionalVectorArray( + wavelength=30 * u.km / u.s, + position=1 * u.arcsec, + ) + ], +) +def test_random_gaussians( + inputs: na.AbstractSpectralPositionalVectorArray, + width: na.AbstractSpectralPositionalVectorArray, +): + result = ctis.scenes.random_gaussians( + inputs=inputs, + width=width, + ) + + assert result.inputs is inputs + assert result.shape == inputs.shape From 7a1922eaabff50724a1371d0714470578bf5b549 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 21 Apr 2025 21:14:40 -0600 Subject: [PATCH 2/8] deps --- pyproject.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 1574239..bc4253a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,8 +17,7 @@ classifiers = [ ] dependencies = [ "astropy", - "named-arrays", - "optika", + "named-arrays==0.21.0", ] dynamic = ["version"] From 0eebb5b09ddd4010c487318bfe5d793d0f423119 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 21 Apr 2025 21:18:20 -0600 Subject: [PATCH 3/8] doc tweaks --- ctis/scenes/_gaussians.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctis/scenes/_gaussians.py b/ctis/scenes/_gaussians.py index b65f0d5..7a1ff1c 100644 --- a/ctis/scenes/_gaussians.py +++ b/ctis/scenes/_gaussians.py @@ -74,7 +74,7 @@ def random_gaussians( # Define the grid of positions and velocities on which to evaluate the # test pattern inputs = na.SpectralPositionalVectorArray( - wavelength=na.linspace(-500, 500, axis="wavelength", num=101) * u.km / u.s, + wavelength=na.linspace(-500, 500, axis="wavelength", num=31) * u.km / u.s, position=na.Cartesian2dVectorLinearSpace( start=-20 * platescale * u.pix, stop=20 * platescale * u.pix, From fc8b9aeec977fa07ca225b7f3571b84b98462ea9 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 21 Apr 2025 21:58:08 -0600 Subject: [PATCH 4/8] fixes --- ctis/scenes/__init__.py | 4 +- ctis/scenes/_gaussians.py | 131 +++++++++++++++++++------------------- 2 files changed, 69 insertions(+), 66 deletions(-) diff --git a/ctis/scenes/__init__.py b/ctis/scenes/__init__.py index dd5626c..13859de 100644 --- a/ctis/scenes/__init__.py +++ b/ctis/scenes/__init__.py @@ -2,8 +2,8 @@ Test patterns used to evaluate the quality of CTIS inversions. """ -from ._gaussians import random_gaussians +from ._gaussians import gaussians __all__ = [ - "random_gaussians", + "gaussians", ] diff --git a/ctis/scenes/_gaussians.py b/ctis/scenes/_gaussians.py index 7a1ff1c..e3b437e 100644 --- a/ctis/scenes/_gaussians.py +++ b/ctis/scenes/_gaussians.py @@ -4,7 +4,7 @@ import named_arrays as na __all__ = [ - "random_gaussians", + "gaussians", ] SpectralPositionalVectorT = TypeVar( @@ -34,31 +34,37 @@ def _gaussian( The standard deviation of the Gaussian distribution. """ + inputs = inputs.explicit + + if inputs.ndim != 3: # pragma: nocover + raise ValueError(f"The number of logical axes must be 3, got {inputs.ndim=}") + + inputs = inputs.cell_centers() + arg = -np.square(((inputs - center) / width).length) / 2 return amplitude * np.exp(arg) -def random_gaussians( +def gaussians( inputs: SpectralPositionalVectorT, width: na.AbstractSpectralPositionalVectorArray, ) -> na.FunctionArray[SpectralPositionalVectorT, na.ScalarArray]: r""" - A scene with a randomly-positioned Gaussian distributions originally - prepared by Amy Winebarger. + A scene with eight randomly-positioned Gaussian kernels originally + prepared by Amy R. Winebarger. Parameters ---------- - axis - The names of the logical axes representing each - :math:`\lambda,x,y` coordinate, - num - The number of pixels along each coordinate: + inputs + The grid of wavelengths and positions on which to evaluate the scene. + width + The standard deviation of the Gaussian kernels. Examples -------- - Plot this scene for the input grid and Gaussian widths originally - prepared by Amy Winebarger. + Plot this scene for an input grid and standard deviations similar to the + one originally tested by Amy R. Winebarger. .. jupyter-execute:: @@ -74,7 +80,7 @@ def random_gaussians( # Define the grid of positions and velocities on which to evaluate the # test pattern inputs = na.SpectralPositionalVectorArray( - wavelength=na.linspace(-500, 500, axis="wavelength", num=31) * u.km / u.s, + wavelength=na.linspace(-500, 500, axis="wavelength", num=21) * u.km / u.s, position=na.Cartesian2dVectorLinearSpace( start=-20 * platescale * u.pix, stop=20 * platescale * u.pix, @@ -91,7 +97,7 @@ def random_gaussians( # Compute the scene of random Gaussians for the # given input grid and standard deviations - scene = ctis.scenes.random_gaussians( + scene = ctis.scenes.gaussians( inputs=inputs, width=width, ) @@ -119,57 +125,54 @@ def random_gaussians( axs[1].yaxis.set_label_position("right") """ - amplitude = np.array( - [ - 3.38, - 6, - 4, - 3.38, - 5, - 7, - 3.7, - 5.5, - ] - ) - - center_x = np.array( - [ - -2.5, - 2.5, - 0, - 0, - 0, - -4, - 0, - 4, - ] - ) - - center_y = np.array( - [ - 0, - 0, - -2, - 2, - 4, - 4, - -4, - -4, - ] - ) - - center_v = np.array( - [ - +200, - -200, - -150, - +150, - -100, - +100, - -100, - -150, - ] - ) + amplitude = [ + 3.38, + 6, + 4, + 3.38, + 5, + 7, + 3.7, + 5.5, + ] + + center_x = [ + -2.5, + 2.5, + 0, + 0, + 0, + -4, + 0, + 4, + ] + + center_y = [ + 0, + 0, + -2, + 2, + 4, + 4, + -4, + -4, + ] + + center_v = [ + +200, + -200, + -150, + +150, + -100, + +100, + -100, + -150, + ] + + amplitude = np.array(amplitude) + center_x = np.array(center_x) + center_y = np.array(center_y) + center_v = np.array(center_v) axis = "event" From 2838a11a9395d1a724f6f53057dc4df6c24ba381 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 22 Apr 2025 11:11:59 -0600 Subject: [PATCH 5/8] fixes --- ctis/scenes/_gaussians_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctis/scenes/_gaussians_test.py b/ctis/scenes/_gaussians_test.py index cd33dff..2bbb9a3 100644 --- a/ctis/scenes/_gaussians_test.py +++ b/ctis/scenes/_gaussians_test.py @@ -31,7 +31,7 @@ def test_random_gaussians( inputs: na.AbstractSpectralPositionalVectorArray, width: na.AbstractSpectralPositionalVectorArray, ): - result = ctis.scenes.random_gaussians( + result = ctis.scenes.gaussians( inputs=inputs, width=width, ) From ebff774bc08e7c2b3f2540982c80e1ed33b20c03 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 22 Apr 2025 11:16:30 -0600 Subject: [PATCH 6/8] fix scaling --- ctis/scenes/_gaussians.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ctis/scenes/_gaussians.py b/ctis/scenes/_gaussians.py index e3b437e..f065a7a 100644 --- a/ctis/scenes/_gaussians.py +++ b/ctis/scenes/_gaussians.py @@ -174,6 +174,9 @@ def gaussians( center_y = np.array(center_y) center_v = np.array(center_v) + scale = 1000 * u.photon / (u.cm**2 * u.arcsec**2 * u.mAA * u.s) + amplitude = amplitude * scale + axis = "event" amplitude = na.ScalarArray(amplitude, axis) From 691b73b0b06154e4fb05063860777a05fde6f29f Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 22 Apr 2025 11:31:07 -0600 Subject: [PATCH 7/8] tests --- ctis/scenes/_gaussians_test.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ctis/scenes/_gaussians_test.py b/ctis/scenes/_gaussians_test.py index 2bbb9a3..2818055 100644 --- a/ctis/scenes/_gaussians_test.py +++ b/ctis/scenes/_gaussians_test.py @@ -1,4 +1,5 @@ import pytest +import numpy as np import astropy.units as u import named_arrays as na import ctis @@ -37,4 +38,5 @@ def test_random_gaussians( ) assert result.inputs is inputs - assert result.shape == inputs.shape + assert result.axes == inputs.axes + assert np.all(result >= 0 * u.photon / (u.cm**2 * u.sr * u.s * u.AA)) From a4ad911caf821e9e45721efff393e4a6f4acdfe2 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 22 Apr 2025 11:32:11 -0600 Subject: [PATCH 8/8] rename --- ctis/scenes/_gaussians_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctis/scenes/_gaussians_test.py b/ctis/scenes/_gaussians_test.py index 2818055..ae7397c 100644 --- a/ctis/scenes/_gaussians_test.py +++ b/ctis/scenes/_gaussians_test.py @@ -28,7 +28,7 @@ ) ], ) -def test_random_gaussians( +def test_gaussians( inputs: na.AbstractSpectralPositionalVectorArray, width: na.AbstractSpectralPositionalVectorArray, ):