windhorst et al 2023 F200W mag Priors

.github/workflows/ci_tests.yml

@@ -18,8 +18,8 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python: [3.9]
-        deps: [current, numpy123, astropydev, numpydev, astropydev-numpydev]
+        python: [3.11]
+        deps: [current, numpy210, astropydev, numpydev, astropydev-numpydev]
 
     steps:
     - name: Check out repository
@@ -31,10 +31,10 @@ jobs:
     - name: Install base dependencies
       run: |
         python -m pip install --upgrade pip
-    - name: Test with numpy = 1.23
-      if: "contains(matrix.deps, 'numpy123')"
+    - name: Test with numpy = 2.1.0
+      if: "contains(matrix.deps, 'numpy210')"
       run: |
-        python -m pip install numpy==1.23
+        python -m pip install numpy==2.1.0
     - name: Test with dev version of numpy
       if: "contains(matrix.deps, 'numpydev')"
       run: |
@@ -65,7 +65,7 @@ jobs:
     - name: Python codestyle check
       uses: actions/setup-python@v2
       with:
-        python-version: 3.9
+        python-version: 3.11
     - name: Install base dependencies
       run: |
         python -m pip install --upgrade pip

astropath/chance.py

@@ -3,6 +3,9 @@
 
 from scipy import interpolate
 
+from importlib.resources import files
+import os
+
 from IPython import embed
 
 
@@ -12,7 +15,8 @@
 driver_spl = interpolate.UnivariateSpline._from_tck(driver_tck)
 
 
-def driver_sigma(rmag):
+
+def driver_sigma(mag):
     """
     Estimated incidence of galaxies per sq arcsec with r > rmag
     using Driver et al. 2016 number counts.
@@ -26,7 +30,32 @@ def driver_sigma(rmag):
         float or np.ndarray:  Galaxy number density
 
     """
-    return 10**driver_spl(rmag)
+    return 10**driver_spl(mag)
+
+def windhorst_sigma(mag):
+    """
+    Estimated incidence of galaxies per sq arcsec with F200W > mag
+    using Windhorst et al. 2024 number counts.
+
+    Spline parameters (globals) are for F200W vs Num counts
+
+    Args:
+        mag (float or np.ndarray): F200W band magnitude of galaxy
+
+    Returns:
+        float or np.ndarray:  Galaxy number density
+
+    """
+    data_path = os.path.join(files('astropath'),'data','galaxy_num_counts','windhorst2023_F200W.npz')
+    data = np.load(data_path)
+    mag_f200w = data['mag']
+    Num = data['Num(N/arcsec/0.5mag))']
+
+    winhorst_spline = interpolate.interp1d(mag_f200w, Num, kind='cubic', fill_value='extrapolate')
+
+    num_counts = winhorst_spline(mag)
+
+    return num_counts
 
 
 def bloom_sigma(rmag):

astropath/path.py

@@ -11,6 +11,8 @@
 from astropath import localization
 from astropath import priors
 
+from IPython import embed
+
 
 class PATH(object):
     """Convenience class to run PATH analysis
@@ -122,7 +124,7 @@ def init_localization(self, ltype:str, **kwargs):
         assert localization.vet_localization(self.localiz), 'Bad candidate prior input'
         logging.info("Localization is ready!")
 
-    def calc_priors(self):
+    def calc_priors(self, ifilter:str='r'):
         """Calculate and normalize the P(O) values for the candidates
 
         Raises:
@@ -138,7 +140,8 @@ def calc_priors(self):
         logging.info("Calculating priors")
         self.raw_prior_Oi = priors.raw_prior_Oi(
             self.cand_prior['P_O_method'], self.candidates['ang_size'], 
-            mag=self.candidates['mag'] if 'mag' in self.candidates.keys() else None)
+            mag=self.candidates['mag'] if 'mag' in self.candidates.keys() else None,
+            filter=ifilter)
 
         # Normalize
         logging.info("Normalizing priors")

astropath/priors.py

@@ -31,39 +31,50 @@
 }
 
 
-
 def raw_prior_Oi(method, ang_size, mag=None, filter='r'):
     """
-    Calculate raw prior for galaxy candidates
-    If Sigma(m) is required, we adopt the Driver et al. 2016 
-    evaluation 
+    Calculate the raw prior probability for galaxy candidates based on angular size and magnitude.
+
+    Parameters
+    ----------
+    method : str
+        Method to use for prior calculation. Options:
+        - 'inverse'       : 1 / Sigma_m
+        - 'inverse_ang'   : 1 / (Sigma_m * ang_size)
+        - 'inverse_ang2'  : 1 / (Sigma_m * ang_size^2)
+        - 'identical'     : All priors set equal (1)
+        - 'user'          : Use user-defined function (not implemented here)
+
+    ang_size : float or np.ndarray
+        Angular size(s) in arcseconds. Required for all methods except 'identical'.
+
+    mag : float or np.ndarray, optional
+        Apparent magnitudes (assumed extinction-corrected r-band). Required for all methods except 'identical'.
+
+    filter : str
+        Photometric filter to use for number counts. Supported: 'r', 'F200W'
+
+    Returns
+    -------
+    float or np.ndarray
+        Raw prior value(s)
+    """
+    if method == 'identical':
+        return 1.0 if np.isscalar(ang_size) else np.ones_like(ang_size)
 
-    Args:
-        method (str):
-            inverse :: Assign inverse to Sigma_m
-            inverse_ang :: Assign inverse to Sigma_m * half_light
-            inverse_ang2 :: Assign inverse to Sigma_m * half_light**2
-            identical :: All the same
-            user :: user-defined function
-        ang_size (float or np.ndarray):
-            Angular size of the galaxy in arcsec
-            Only required for several methods
-        mag (float or np.ndarray, optional):
-            Magnitudes of the sources;  assumes r-band and corrected
-            for Galactic extinction
+    if mag is None:
+        raise ValueError("Magnitude `mag` must be provided for non-identical methods.")
 
-    Returns:
-        float or np.ndarray:
+    supported_filters = {
+        'r': chance.driver_sigma,
+        'F200W': chance.windhorst_sigma
+    }
 
-    """
-    # allows a user to set this externally
-    global USR_raw_prior_Oi
-
-    # Convenience
-    if method not in ['identical']:
-        if filter != 'r':
-            raise IOError("Not ready for this.  Best to go with what you have that is closest to r-band")
-        Sigma_m = chance.driver_sigma(mag)
+    if filter not in supported_filters:
+        raise ValueError(f"Unsupported filter '{filter}'. Supported: {list(supported_filters.keys())}")
+
+    # Compute surface density
+    Sigma_m = supported_filters[filter](mag)
 
     # Do it
     if method == 'inverse':

astropath/tests/test_priors.py

@@ -1,8 +1,6 @@
 """ Test methods in bayesian.py """
 
 import os
-from pkg_resources import resource_filename
-
 import numpy as np
 import pandas
 import healpy as hp
@@ -12,6 +10,8 @@
 from astropy.io import fits
 from astropy.wcs import WCS
 
+from unittest.mock import patch
+
 from astropath import bayesian, priors
 from astropath import localization
 
@@ -20,53 +20,103 @@
 remote_data = pytest.mark.skipif(os.getenv('PATH_DATA') is None,
                                  reason='test requires dev suite')
 
-def test_raw_prior():
-    # Inverse
-    mag=np.array([21.,22.])
+
+# ================================
+# DRIVER TESTS (default filter='r')
+# ================================
+
+@patch('astropath.chance.driver_sigma')
+def test_raw_prior_inverse(mock_driver_sigma):
+    mag = np.array([21., 22.])
     ang_size = np.ones_like(mag)
-    raw_PO = priors.raw_prior_Oi('inverse', ang_size, mag=np.array([21.,22.]))
-    assert np.all(np.isclose(raw_PO, np.array([2736.80588898, 1074.04504448])))
-
-def user_func(a,b,c):
-    """
-    Arbitrary user-defined function
-    """
-    return a+b+c
-
-def test_user_prior():
-    # Inverse
-    mag=np.array([21.,22.])
+    mock_driver_sigma.return_value = np.array([0.0003655, 0.0009314])
+    result = priors.raw_prior_Oi('inverse', ang_size, mag=mag)
+    expected = 1. / mock_driver_sigma.return_value
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+@patch('astropath.chance.driver_sigma')
+def test_raw_prior_inverse_ang(mock_driver_sigma):
+    mag = np.array([21., 22.])
+    ang_size = np.array([2.0, 4.0])
+
+    mock_driver_sigma.return_value = np.array([0.0003655, 0.0009314])
+    result = priors.raw_prior_Oi('inverse_ang', ang_size, mag=mag)
+    expected = 1. / mock_driver_sigma.return_value / ang_size
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+@patch('astropath.chance.driver_sigma')
+def test_raw_prior_inverse_ang2(mock_driver_sigma):
+    mag = np.array([21., 22.])
+    ang_size = np.array([2.0, 4.0])
+
+    mock_driver_sigma.return_value = np.array([0.0003655, 0.0009314])
+    result = priors.raw_prior_Oi('inverse_ang2', ang_size, mag=mag)
+    expected = 1. / mock_driver_sigma.return_value / ang_size**2
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+
+# ================================
+# WINDHORST TESTS (filter='F200W')
+# ================================
+
+@patch('astropath.chance.windhorst_sigma')
+def test_raw_prior_windhorst(mock_windhorst_sigma):
+    mag = np.array([20., 22.])
     ang_size = np.ones_like(mag)
-    priors.USR_raw_prior_Oi = user_func
-    raw_PO = priors.raw_prior_Oi('user', ang_size, mag=np.array([21.,22.]))
-    assert np.all(np.isclose(raw_PO, np.array([22.00036539 23.00093106])))
+
+    mock_windhorst_sigma.return_value = np.array([0.001, 0.002])
+    result = priors.raw_prior_Oi('inverse', ang_size, mag=mag, filter='F200W')
+    expected = 1. / mock_windhorst_sigma.return_value
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+@patch('astropath.chance.windhorst_sigma')
+def test_raw_prior_windhorst_ang(mock_windhorst_sigma):
+    mag = np.array([20., 22.])
+    ang_size = np.array([2.0, 4.0])
+
+    mock_windhorst_sigma.return_value = np.array([0.001, 0.002])
+    result = priors.raw_prior_Oi('inverse_ang', ang_size, mag=mag, filter='F200W')
+    expected = 1. / mock_windhorst_sigma.return_value / ang_size
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+@patch('astropath.chance.windhorst_sigma')
+def test_raw_prior_windhorst_ang2(mock_windhorst_sigma):
+    mag = np.array([20., 22.])
+    ang_size = np.array([2.0, 4.0])
+
+    mock_windhorst_sigma.return_value = np.array([0.001, 0.002])
+    result = priors.raw_prior_Oi('inverse_ang2', ang_size, mag=mag, filter='F200W')
+    expected = 1. / mock_windhorst_sigma.return_value / ang_size**2
+    np.testing.assert_allclose(result, expected, rtol=1e-5)
+
+
+# ================================
+# GENERIC TESTS
+# ================================
+
+def test_raw_prior_identical():
+    ang_size = np.array([1.0, 2.0, 3.0])
+    result = priors.raw_prior_Oi('identical', ang_size)
+    np.testing.assert_array_equal(result, np.ones_like(ang_size))
 
 def test_renorm_priors():
-    # U=0
     raw_Oi = np.array([0.1, 0.2, 0.5])
     renorm_Oi = priors.renorm_priors(raw_Oi, 0.)
+    assert np.allclose(renorm_Oi, np.array([0.125, 0.25, 0.625]))
 
-    assert np.all(renorm_Oi == np.array([0.125, 0.25 , 0.625]))
-
-    # U != 0
     renorm_Oi = priors.renorm_priors(raw_Oi, 0.1)
-    assert np.all(renorm_Oi == np.array([0.1125, 0.225, 0.5625]))
+    assert np.allclose(renorm_Oi, np.array([0.1125, 0.225, 0.5625]))
 
 def test_norm_theta_priors():
-    # Grid me
     ngrid = 2000
     x = np.linspace(-30., 30., ngrid)
     grid_spacing_arcsec = x[1]-x[0]
 
-    xcoord, ycoord = np.meshgrid(x,x)
+    xcoord, ycoord = np.meshgrid(x, x)
     theta = np.sqrt(xcoord**2 + ycoord**2)
 
-    # Uniform
     for pdf in ['core', 'uniform', 'exp']:
         theta_prior = dict(max=6, PDF=pdf, scale=1.)
-        p_wOi = bayesian.pw_Oi(theta,
-                      2.0, # phi
-                      theta_prior)
-        assert np.isclose(np.sum(p_wOi)*grid_spacing_arcsec**2, 
-                          1., atol=1e-4), f'Failed normalization on {pdf}'
-
+        p_wOi = bayesian.pw_Oi(theta, 2.0, theta_prior)
+        assert np.isclose(np.sum(p_wOi) * grid_spacing_arcsec**2, 1., atol=1e-4), f'Failed normalization on {pdf}'
+