Merge pull request #65 from CosmoStat/feat/size-conversions

feat: size module — Gaussian size-conversion web (T, sigma, r50, FWHM)

Author: cailmdaley

cs_util/size.py

@@ -0,0 +1,212 @@
+"""SIZE.
+
+:Name: size.py
+
+:Description: Conversions between the Gaussian object-size measures
+    used across the UNIONS / ShapePipe stack.
+
+    All conversions assume a circular Gaussian profile, for which the
+    measures are related through the scale parameter ``sigma``:
+
+    - ``T = 2 sigma^2`` — the ngmix / DES area parameter (arcsec^2),
+    - ``r50 = sqrt(2 ln 2) sigma = 1.17741 sigma`` — the half-light
+      radius (the primary size in the UNIONS shape-catalogue papers),
+    - ``FWHM = 2 sqrt(2 ln 2) sigma = 2.35482 sigma``.
+
+    These functions are the single source of truth for the size web;
+    producers (ShapePipe) and consumers (sp_validation) should import
+    from here rather than re-deriving the factors locally.
+
+:Author: Cail Daley <cailmdaley@gmail.com>
+
+"""
+
+import numpy as np
+
+# Half-light radius of a circular Gaussian per unit sigma:
+# r50 = sqrt(2 ln 2) * sigma
+SIGMA_TO_R50 = np.sqrt(2 * np.log(2))
+
+# Full width at half maximum of a Gaussian per unit sigma:
+# FWHM = 2 sqrt(2 ln 2) * sigma
+SIGMA_TO_FWHM = 2 * np.sqrt(2 * np.log(2))
+
+
+def T_to_sigma(T):
+    """T To Sigma.
+
+    Gaussian scale ``sigma = sqrt(T / 2)`` from the ngmix area
+    parameter ``T = 2 sigma^2``.
+
+    Parameters
+    ----------
+    T : float or numpy.ndarray
+        ngmix area parameter, ``T = 2 sigma^2``
+
+    Returns
+    -------
+    float or numpy.ndarray
+        Gaussian scale ``sigma``, in units of ``sqrt(T)``
+
+    """
+    return np.sqrt(T / 2)
+
+
+def sigma_to_T(sigma):
+    """Sigma To T.
+
+    ngmix area parameter ``T = 2 sigma^2`` from the Gaussian scale.
+
+    Parameters
+    ----------
+    sigma : float or numpy.ndarray
+        Gaussian scale
+
+    Returns
+    -------
+    float or numpy.ndarray
+        ngmix area parameter ``T``, in units of ``sigma^2``
+
+    """
+    return 2 * sigma ** 2
+
+
+def sigma_to_r50(sigma):
+    """Sigma To R50.
+
+    Half-light radius ``r50 = sqrt(2 ln 2) sigma`` of a circular
+    Gaussian.
+
+    Parameters
+    ----------
+    sigma : float or numpy.ndarray
+        Gaussian scale
+
+    Returns
+    -------
+    float or numpy.ndarray
+        half-light radius, in units of ``sigma``
+
+    """
+    return SIGMA_TO_R50 * sigma
+
+
+def r50_to_sigma(r50):
+    """R50 To Sigma.
+
+    Gaussian scale from the half-light radius.
+
+    Parameters
+    ----------
+    r50 : float or numpy.ndarray
+        half-light radius
+
+    Returns
+    -------
+    float or numpy.ndarray
+        Gaussian scale ``sigma``, in units of ``r50``
+
+    """
+    return r50 / SIGMA_TO_R50
+
+
+def sigma_to_fwhm(sigma):
+    """Sigma To Fwhm.
+
+    Full width at half maximum ``FWHM = 2 sqrt(2 ln 2) sigma`` of a
+    Gaussian.
+
+    Parameters
+    ----------
+    sigma : float or numpy.ndarray
+        Gaussian scale
+
+    Returns
+    -------
+    float or numpy.ndarray
+        FWHM, in units of ``sigma``
+
+    """
+    return SIGMA_TO_FWHM * sigma
+
+
+def fwhm_to_sigma(fwhm):
+    """Fwhm To Sigma.
+
+    Gaussian scale from the full width at half maximum.
+
+    Parameters
+    ----------
+    fwhm : float or numpy.ndarray
+        full width at half maximum
+
+    Returns
+    -------
+    float or numpy.ndarray
+        Gaussian scale ``sigma``, in units of ``fwhm``
+
+    """
+    return fwhm / SIGMA_TO_FWHM
+
+
+def T_to_r50(T):
+    """T To R50.
+
+    Half-light radius ``r50 = sqrt(2 ln 2) * sqrt(T / 2)
+    = sqrt(ln 2 * T)`` from the ngmix area parameter.
+
+    Parameters
+    ----------
+    T : float or numpy.ndarray
+        ngmix area parameter, ``T = 2 sigma^2``
+
+    Returns
+    -------
+    float or numpy.ndarray
+        half-light radius, in units of ``sqrt(T)``
+
+    """
+    return sigma_to_r50(T_to_sigma(T))
+
+
+def r50_to_T(r50):
+    """R50 To T.
+
+    ngmix area parameter ``T = 2 (r50 / sqrt(2 ln 2))^2 = r50^2 / ln 2``
+    from the half-light radius.
+
+    Parameters
+    ----------
+    r50 : float or numpy.ndarray
+        half-light radius
+
+    Returns
+    -------
+    float or numpy.ndarray
+        ngmix area parameter ``T``, in units of ``r50^2``
+
+    """
+    return sigma_to_T(r50_to_sigma(r50))
+
+
+def T_to_fwhm(T):
+    """T To Fwhm.
+
+    Full width at half maximum ``FWHM = 2 sqrt(2 ln 2) * sqrt(T / 2)``
+    from the ngmix area parameter.
+
+    Note that ``T`` is an *area*: the conversion to the length ``FWHM``
+    involves a square root, not a linear factor.
+
+    Parameters
+    ----------
+    T : float or numpy.ndarray
+        ngmix area parameter, ``T = 2 sigma^2``
+
+    Returns
+    -------
+    float or numpy.ndarray
+        FWHM, in units of ``sqrt(T)``
+
+    """
+    return sigma_to_fwhm(T_to_sigma(T))

cs_util/tests/test_size.py

@@ -0,0 +1,93 @@
+"""UNIT TESTS FOR SIZE SUBPACKAGE.
+
+This module contains unit tests for the size subpackage.
+
+"""
+
+import numpy as np
+from numpy import testing as npt
+
+from unittest import TestCase
+
+from cs_util import size
+
+
+class SizeTestCase(TestCase):
+    """Test case for the ``size`` module."""
+
+    def setUp(self):
+        """Set test parameter values."""
+        # Unit-sigma Gaussian: T = 2, r50 = 1.17741, FWHM = 2.35482
+        self._sigma = 1.0
+        self._T = 2.0
+        self._r50 = 1.1774100226
+        self._fwhm = 2.3548200450
+
+    def tearDown(self):
+        """Unset test parameter values."""
+        self._sigma = None
+        self._T = None
+        self._r50 = None
+        self._fwhm = None
+
+    def test_constants(self):
+        """Test the module constants against their closed forms."""
+        npt.assert_almost_equal(size.SIGMA_TO_R50, np.sqrt(2 * np.log(2)))
+        npt.assert_almost_equal(size.SIGMA_TO_FWHM, 2 * np.sqrt(2 * np.log(2)))
+        npt.assert_almost_equal(size.SIGMA_TO_R50, 1.17741, decimal=5)
+        npt.assert_almost_equal(size.SIGMA_TO_FWHM, 2.35482, decimal=5)
+
+    def test_unit_sigma_values(self):
+        """Test all conversions on the unit-sigma Gaussian."""
+        npt.assert_almost_equal(size.T_to_sigma(self._T), self._sigma)
+        npt.assert_almost_equal(size.sigma_to_T(self._sigma), self._T)
+        npt.assert_almost_equal(size.sigma_to_r50(self._sigma), self._r50)
+        npt.assert_almost_equal(size.sigma_to_fwhm(self._sigma), self._fwhm)
+        npt.assert_almost_equal(size.T_to_r50(self._T), self._r50)
+        npt.assert_almost_equal(size.T_to_fwhm(self._T), self._fwhm)
+
+    def test_T_to_r50_closed_form(self):
+        """Test ``T_to_r50(T) == sqrt(ln 2 * T)``."""
+        T = np.array([0.05, 0.18, 2.0, 10.0])
+        npt.assert_allclose(size.T_to_r50(T), np.sqrt(np.log(2) * T))
+
+    def test_inverse_direct_values(self):
+        """Test each inverse converter directly on the unit-sigma case."""
+        npt.assert_almost_equal(size.r50_to_sigma(self._r50), self._sigma)
+        npt.assert_almost_equal(size.fwhm_to_sigma(self._fwhm), self._sigma)
+        npt.assert_almost_equal(size.r50_to_T(self._r50), self._T)
+
+    def test_nonpositive_T(self):
+        """Test that T = 0 maps to size 0 and negative T to NaN."""
+        npt.assert_equal(size.T_to_sigma(0.0), 0.0)
+        npt.assert_equal(size.T_to_r50(0.0), 0.0)
+        npt.assert_equal(size.T_to_fwhm(0.0), 0.0)
+        with np.errstate(invalid="ignore"):
+            self.assertTrue(np.isnan(size.T_to_r50(-1.0)))
+
+    def test_round_trips(self):
+        """Test that inverse pairs compose to the identity."""
+        values = np.array([0.01, 0.1, 1.0, 7.5])
+        npt.assert_allclose(size.r50_to_T(size.T_to_r50(values)), values)
+        npt.assert_allclose(size.T_to_r50(size.r50_to_T(values)), values)
+        npt.assert_allclose(
+            size.sigma_to_T(size.T_to_sigma(values)), values
+        )
+        npt.assert_allclose(
+            size.r50_to_sigma(size.sigma_to_r50(values)), values
+        )
+        npt.assert_allclose(
+            size.fwhm_to_sigma(size.sigma_to_fwhm(values)), values
+        )
+
+    def test_fwhm_is_twice_r50(self):
+        """Test ``FWHM = 2 r50`` for a Gaussian, via both paths."""
+        T = np.array([0.05, 0.18, 2.0])
+        npt.assert_allclose(size.T_to_fwhm(T), 2 * size.T_to_r50(T))
+
+    def test_array_input(self):
+        """Test that array input returns arrays of the same shape."""
+        T = np.linspace(0.01, 5.0, 11)
+        r50 = size.T_to_r50(T)
+        self.assertEqual(r50.shape, T.shape)
+        self.assertTrue(np.all(np.diff(r50) > 0))