diff --git a/cs_util/shape.py b/cs_util/shape.py new file mode 100644 index 0000000..d3c38bf --- /dev/null +++ b/cs_util/shape.py @@ -0,0 +1,82 @@ +"""SHAPE. + +:Name: shape.py + +:Description: Conversions between the two ellipticity conventions used + across the UNIONS / ShapePipe stack. + + A source's shape is a spin-2 quantity that can be expressed either + as the *distortion* (e-type) or the *reduced shear* (g-type). For + an ellipse of axis ratio ``q = b / a`` they have magnitudes + + - ``|e| = (1 - q^2) / (1 + q^2)`` — distortion, and + - ``|g| = (1 - q) / (1 + q)`` — reduced shear, + + related (at fixed position angle) by + + - ``g = e / (1 + sqrt(1 - e^2))`` and + - ``e = 2 g / (1 + g^2)``. + + Both components share the orientation, so the conversion scales the + complex ellipticity ``(comp1, comp2)`` by the magnitude ratio; the + two-component functions below apply that ratio component-wise. + + Convention note: ShapePipe's HSM moments (``*_E*_*_HSM`` columns) + store the distortion ``e``, while ngmix stores the reduced shear + ``g``; this module is the single bridge between them, so producers + (ShapePipe) and consumers (sp_validation) import from here rather + than re-deriving the factors locally. + +:Author: Cail Daley + +""" + +import numpy as np + + +def e_to_g(e1, e2): + """E To G. + + Reduced shear ``g`` (g-type) from the distortion ``e`` (e-type), + component-wise: ``g_i = e_i / (1 + sqrt(1 - |e|^2))`` with + ``|e|^2 = e1^2 + e2^2``. + + Parameters + ---------- + e1 : float or numpy.ndarray + First distortion component + e2 : float or numpy.ndarray + Second distortion component + + Returns + ------- + tuple of float or numpy.ndarray + Reduced-shear components ``(g1, g2)`` + + """ + factor = 1 / (1 + np.sqrt(1 - e1 ** 2 - e2 ** 2)) + return e1 * factor, e2 * factor + + +def g_to_e(g1, g2): + """G To E. + + Distortion ``e`` (e-type) from the reduced shear ``g`` (g-type), + component-wise: ``e_i = 2 g_i / (1 + |g|^2)`` with + ``|g|^2 = g1^2 + g2^2``. + + Parameters + ---------- + g1 : float or numpy.ndarray + First reduced-shear component + g2 : float or numpy.ndarray + Second reduced-shear component + + Returns + ------- + tuple of float or numpy.ndarray + Distortion components ``(e1, e2)`` + + """ + factor = 2 / (1 + g1 ** 2 + g2 ** 2) + return g1 * factor, g2 * factor diff --git a/cs_util/tests/test_shape.py b/cs_util/tests/test_shape.py new file mode 100644 index 0000000..ddc5906 --- /dev/null +++ b/cs_util/tests/test_shape.py @@ -0,0 +1,97 @@ +"""UNIT TESTS FOR SHAPE SUBPACKAGE. + +This module contains unit tests for the shape subpackage. + +""" + +import numpy as np +from numpy import testing as npt + +from unittest import TestCase + +from cs_util import shape + + +class ShapeTestCase(TestCase): + """Test case for the ``shape`` module.""" + + def setUp(self): + """Set test parameter values. + + Axis ratio q = 1/3 gives the clean pair + |e| = (1 - q^2) / (1 + q^2) = 0.8, + |g| = (1 - q) / (1 + q) = 0.5. + """ + self._e = 0.8 + self._g = 0.5 + + def tearDown(self): + """Unset test parameter values.""" + self._e = None + self._g = None + + def test_known_pair(self): + """Test the q = 1/3 distortion/shear pair, both directions.""" + g1, g2 = shape.e_to_g(self._e, 0.0) + npt.assert_almost_equal(g1, self._g) + npt.assert_almost_equal(g2, 0.0) + e1, e2 = shape.g_to_e(self._g, 0.0) + npt.assert_almost_equal(e1, self._e) + npt.assert_almost_equal(e2, 0.0) + + def test_axis_ratio_magnitudes(self): + """Test |e| and |g| against their axis-ratio closed forms.""" + q = np.array([0.1, 0.3, 0.5, 0.9]) + e_mag = (1 - q ** 2) / (1 + q ** 2) + g_mag = (1 - q) / (1 + q) + g1, g2 = shape.e_to_g(e_mag, np.zeros_like(e_mag)) + npt.assert_allclose(g1, g_mag) + e1, e2 = shape.g_to_e(g_mag, np.zeros_like(g_mag)) + npt.assert_allclose(e1, e_mag) + + def test_round_trips(self): + """Test that the converters compose to the identity, component-wise.""" + g1 = np.array([-0.4, -0.05, 0.0, 0.1, 0.45]) + g2 = np.array([0.2, 0.3, 0.0, -0.25, 0.1]) + e1, e2 = shape.g_to_e(g1, g2) + g1b, g2b = shape.e_to_g(e1, e2) + npt.assert_allclose(g1b, g1) + npt.assert_allclose(g2b, g2) + + e1 = np.array([-0.6, 0.0, 0.3, 0.7]) + e2 = np.array([0.1, 0.0, -0.4, 0.2]) + g1, g2 = shape.e_to_g(e1, e2) + e1b, e2b = shape.g_to_e(g1, g2) + npt.assert_allclose(e1b, e1) + npt.assert_allclose(e2b, e2) + + def test_orientation_preserved(self): + """Test that the conversion scales both components by one factor.""" + e1, e2 = 0.36, 0.48 # |e| = 0.6, position angle fixed + g1, g2 = shape.e_to_g(e1, e2) + # g1 / g2 == e1 / e2: orientation unchanged by the magnitude rescale. + npt.assert_almost_equal(g1 / g2, e1 / e2) + + def test_zero_maps_to_zero(self): + """Test that a round object stays round under both converters.""" + npt.assert_array_equal(shape.e_to_g(0.0, 0.0), (0.0, 0.0)) + npt.assert_array_equal(shape.g_to_e(0.0, 0.0), (0.0, 0.0)) + + def test_small_ellipticity_limit(self): + """Test ``e -> 2 g`` and ``g -> e / 2`` in the round-object limit.""" + g = np.array([1e-4, 5e-4, 1e-3]) + e1, _ = shape.g_to_e(g, np.zeros_like(g)) + npt.assert_allclose(e1, 2 * g, rtol=1e-5) + e = np.array([1e-4, 5e-4, 1e-3]) + g1, _ = shape.e_to_g(e, np.zeros_like(e)) + npt.assert_allclose(g1, e / 2, rtol=1e-5) + + def test_array_input(self): + """Test that array input returns arrays of the same shape.""" + e1 = np.linspace(-0.7, 0.7, 11) + e2 = np.linspace(0.0, 0.5, 11) + g1, g2 = shape.e_to_g(e1, e2) + self.assertEqual(g1.shape, e1.shape) + self.assertEqual(g2.shape, e2.shape) + # |g| < |e| < 1 for every input ellipticity. + self.assertTrue(np.all(np.hypot(g1, g2) < np.hypot(e1, e2) + 1e-12))