diff --git a/cs_util/size.py b/cs_util/size.py new file mode 100644 index 0000000..9c12760 --- /dev/null +++ b/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 + +""" + +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)) diff --git a/cs_util/tests/test_size.py b/cs_util/tests/test_size.py new file mode 100644 index 0000000..3ab01a6 --- /dev/null +++ b/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))