Merge pull request #186 from FRBs/mag_vs_dm

bug fix and added full run test

.github/workflows/ci_tests.yml

@@ -19,8 +19,8 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python: ['3.10', '3.11']
-        deps: [current, numpy124, astropydev, numpydev, astropydev-numpydev]
+        python: ['3.11', '3.12']
+        deps: [current, numpy211, astropydev, numpydev, astropydev-numpydev]
 
     steps:
     - name: Check out repository
@@ -32,10 +32,10 @@ jobs:
     - name: Install base dependencies
       run: |
         python -m pip install --upgrade pip
-    - name: Test with numpy = 1.24
-      if: "contains(matrix.deps, 'numpy124')"
+    - name: Test with numpy = 2.1.1
+      if: "contains(matrix.deps, 'numpy211')"
       run: |
-        python -m pip install numpy==1.24
+        python -m pip install numpy==2.1.1
     - name: Test with dev version of numpy
       if: "contains(matrix.deps, 'numpydev')"
       run: |

bin/frb_pzdm_mag

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+#
+# See top-level LICENSE file for Copyright information
+#
+# -*- coding: utf-8 -*-
+
+
+"""
+This script performs estimates P(z|DM) and limiting mag
+"""
+
+from frb.scripts import pzdm_mag
+
+if __name__ == '__main__':
+    args = pzdm_mag.parser()
+    pzdm_mag.main(args)

docs/installing.rst

@@ -26,13 +26,13 @@ We recommend that you use `Anaconda <https://www.continuum.io/downloads/>`_
 to install and/or update these packages.
 
 * `python <http://www.python.org/>`_ versions 3.8 or later
-* `numpy <http://www.numpy.org/>`_ version 1.22 or later
+* `numpy <http://www.numpy.org/>`_ version 2.1 or later
 * `astropy <http://www.astropy.org/>`_ version 5.1 or later
 * `scipy <http://www.scipy.org/>`_ version 1.11 or later
 * `healpy <https://healpy.readthedocs.io/en/latest/index.html>`_ version 1.15 or later
 * `pandas <https://pandas.pydata.org/>`_ version 1.5 or later
 * `requests <https://pillow.readthedocs.io/en/5.3.x/>`_  version 2.18 or later
-* `extinction <https://extinction.readthedocs.io/en/latest/>`_ version 0.4.2 or greater
+* `extinction <https://extinction.readthedocs.io/en/latest/>`_ version 0.4.6 or greater
 * `matplotlib <https://matplotlib.org/>`_ version 3.7 or greater
 * `linetools <https://github.com/linetools/linetools>`_  version 0.3 or later
 * `astropath <https://github.com/FRBs/astropath>`_  version 0.1 or later
@@ -58,6 +58,7 @@ The following package(s) is/are required to access FRB galaxy spectra:
 The following package is required to map a slit onto a finder chart (frb.figures.finder):
 
 * `photutils <https://photutils.readthedocs.io/en/stable/>`_  version 1.11.0 or later
+* `scikit-image <https://scikit-image.org/>`_  version 0.21.0 or later
 
 The following are required to use our KCWI datacube handling tools:
 

frb/galaxies/eazy.py

@@ -3,8 +3,8 @@
 import os
 import warnings
 from pkg_resources import resource_filename
-from distutils import spawn
 import subprocess
+import shutil
 
 import numpy as np
 import pandas
@@ -373,7 +373,7 @@ def run_eazy(input_dir, name, logfile):
     _, param_file, translate_file = eazy_filenames(input_dir, name)
 
     # Find the eazy executable
-    path_to_eazy = spawn.find_executable('eazy')
+    path_to_eazy = shutil.which('eazy')
     if path_to_eazy is None:
         raise ValueError("You must have eazy in your Unix path..")
     # Run it!

frb/galaxies/utils.py

@@ -6,6 +6,7 @@
 from pkg_resources import resource_filename
 import numpy as np
 from scipy.interpolate import interp1d
+import warnings
 
 import pandas
 
@@ -20,7 +21,7 @@
 from astropy import units
 
 import pandas as pd
-import extinction
+
 from linetools.spectra import xspectrum1d
 
 from frb import frb
@@ -37,6 +38,10 @@ def deredden_spec(spectrum:xspectrum1d.XSpectrum1D, ebv:float):
     """
 
     # Correct for Galactic extinction
+    # Hidnig this here to avoid a hard dependency
+    # TODO
+    #   Need to replace it 
+    import extinction
     AV = ebv * 3.1  # RV
     Al = extinction.ccm89(spectrum.wavelength.value, AV, 3.1)
     # New spec

frb/scripts/limiting_mag.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 """
-Estimate the limiting luminosity given the magnitude limit and DM and coord
+Estimate the limiting luminosity given the magnitude limit and DMs and coord and telescope
 """
 from IPython import embed
 
@@ -12,8 +12,13 @@ def parser(options=None):
     parser.add_argument("DM_FRB", type=float, help="FRB DM")
     parser.add_argument("mag_limit", type=float, help="Magnitude limit in filter *without* extinction correction")
     parser.add_argument("--filter", type=str, default='DECaL_r', help="Filter -- only used for extinction correction.  Must be a Repo approved choice")
-    parser.add_argument("--dm_hostmw", type=float, default=100., help="Assumed DM contribution from MW and Host")
+    parser.add_argument("--dm_host", type=float, default=50., help="Assumed DM contribution from the Host. Default = 50")
+    parser.add_argument("--dm_mwhalo", type=float, default=50., help="Assumed DM contribution from the MW halo. Default = 50")
     #parser.add_argument("-v", "--verbose", default=False, action="store_true", help="Overwhelm the screen?")
+    parser.add_argument("--telescope", type=str, default='perfect', help="telescope model for the DM-z grid: CHIME, DSA, Parkes, FAST, CRAFT, \
+                        CRAFT_ICS_892/1300/1632, perfect. Default = perfect")
+    parser.add_argument("--cl", type=tuple, default=(2.5,97.5), 
+                        help="Confidence limits for the z estimate [default is a 95 percent c.l., (2.5,97.5)]")
 
     if options is None:
         pargs = parser.parse_args()
@@ -25,10 +30,7 @@ def parser(options=None):
 def main(pargs):
     """ Run
     """
-    import json
-    import os
     import numpy as np
-    from pkg_resources import resource_filename
 
     from linetools import utils as ltu
     from linetools.scripts.utils import coord_arg_to_coord
@@ -50,23 +52,58 @@ def main(pargs):
     DM_ISM = mw.ismDM(icoord)
     print(f"NE2001 = {DM_ISM}")
 
-    # DM Cosmic
-    DM_cosmic = pargs.DM_FRB - DM_ISM.value - pargs.dm_hostmw
+    # DM cosmic and EG
+    DM_extragalactic = pargs.DM_FRB - DM_ISM.value - pargs.dm_mwhalo
+    DM_cosmic = DM_extragalactic - pargs.dm_host
 
     # Redshift estimates
 
-    # Load
-    sdict = prob_dmz.grab_repo_grid()
+    # Load the telescope specific grid
+    telescope_dict = {
+        'CHIME': 'CHIME_pzdm.npz',
+        'DSA': 'DSA_pzdm.npy',
+        'Parkes': 'parkes_mb_class_I_and_II_pzdm.npy',
+        'CRAFT': 'CRAFT_class_I_and_II_pzdm.npy',
+        'CRAFT_ICS_1300': 'CRAFT_ICS_1300_pzdm.npy',
+        'CRAFT_ICS_892': 'CRAFT_ICS_892_pzdm.npy',
+        'CRAFT_ICS_1632': 'CRAFT_ICS_1632_pzdm.npy',
+        'FAST': 'FAST_pzdm.npy',
+        'perfect': 'PDM_z.npz'
+    }
+
+    # Get the perfect telescope grid (default)
+    sdict = prob_dmz.grab_repo_grid(telescope_dict['perfect'])
     PDM_z = sdict['PDM_z']
     z = sdict['z']
     DM = sdict['DM']
 
-    # Do it
+    # Grab the right entry
     iDM = np.argmin(np.abs(DM - DM_cosmic))
     PzDM = PDM_z[iDM, :] / np.sum(PDM_z[iDM, :])
 
+
+    # Get the telescope specific PZDM grid
+    if pargs.telescope and pargs.telescope != 'CHIME' and pargs.telescope != 'perfect':
+        if pargs.telescope not in telescope_dict:
+            raise ValueError(f"Unknown telescope: {pargs.telescope}")
+        zdict = prob_dmz.grab_repo_grid(telescope_dict['CHIME'])
+        z = zdict['z']
+        DM = zdict['DM']
+        PDM_z = prob_dmz.grab_repo_grid(telescope_dict[pargs.telescope])
+        iDM = np.argmin(np.abs(DM - DM_extragalactic))
+        PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM])
+
+
+    if pargs.telescope and pargs.telescope == 'CHIME':
+        sdict = prob_dmz.grab_repo_grid(telescope_dict['CHIME'])
+        PDM_z = sdict['pzdm']
+        z = sdict['z']
+        DM = sdict['DM']
+        iDM = np.argmin(np.abs(DM - DM_extragalactic))
+        PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM])
+
     cum_sum = np.cumsum(PzDM)
-    limits = [10, 90]
+    limits = pargs.cl
 
     z_min = z[np.argmin(np.abs(cum_sum-limits[0]/100.))]
     z_max = z[np.argmin(np.abs(cum_sum-limits[1]/100.))]
@@ -92,8 +129,18 @@ def main(pargs):
 
     # Finish
     print("-----------------------------------------------------")
-    print(f"For z_{limits[0]}={z_min:.2f}, the limiting magnitude corresponds to L={frac_Lstar_min:.5f}L*")
-    print(f"For z_{limits[1]}={z_max:.2f}, the limiting magnitude corresponds to L={frac_Lstar_max:.5f}L*")
+    print("")
+    print(f"Allowing for the MW halo, DM_MW_halo = {int(pargs.dm_mwhalo)} pc/cm^3")
+    print(f"Allowing for the Host, DM_host = {int(pargs.dm_host)} pc/cm^3")
+    print("")
+    if not pargs.telescope or pargs.telescope == 'perfect':
+        print("WARNING: This all assumes a perfect telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
+    else:
+        print("This assumes the "+(str(pargs.telescope))+" telescope and a model of the scatter in DM_cosmic (Macquart+2020)")
+    print("")
+
+    print(f"For z_({limits[0]} %)={z_min:.2f}, the limiting magnitude corresponds to L={frac_Lstar_min:.5f}L*")
+    print(f"For z_({limits[1]} %)={z_max:.2f}, the limiting magnitude corresponds to L={frac_Lstar_max:.5f}L*")
 
     return frac_Lstar_min, frac_Lstar_max
 

frb/scripts/pz_dm.py

@@ -14,7 +14,7 @@ def parser(options=None):
     parser.add_argument("DM_FRB", type=float, help="FRB DM (pc/cm^3)")
     parser.add_argument("--dm_host", type=float, default=50., help="Assumed DM contribution from the Host. Default = 50")
     parser.add_argument("--dm_mwhalo", type=float, default=50., help="Assumed DM contribution from the MW halo. Default = 50")
-    parser.add_argument("--cl", type=tuple, default=(2.5,97.5), 
+    parser.add_argument("--cl", type=str, default="2.5,97.5", 
                         help="Confidence limits for the z estimate [default is a 95 percent c.l., (2.5,97.5)]")
     parser.add_argument("--magdm_plot", default=False, action='store_true', 
                         help="Plot the host redshift range given DM on the magnitude vs redshift evolution")
@@ -52,8 +52,8 @@ def main(pargs):
     print(f"NE2001 = {DM_ISM:.2f}")
 
     # DM cosmic and EG
-    DM_cosmic = pargs.DM_FRB - DM_ISM.value - pargs.dm_mwhalo
-    DM_extragalactic = DM_cosmic + pargs.dm_host  
+    DM_extragalactic = pargs.DM_FRB - DM_ISM.value - pargs.dm_mwhalo
+    DM_cosmic = DM_extragalactic - pargs.dm_host
 
     # Redshift estimates
 
@@ -102,7 +102,7 @@ def main(pargs):
         PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM])
 
     cum_sum = np.cumsum(PzDM)
-    limits = pargs.cl
+    limits = [float(item) for item in pargs.cl.split(',')]
 
     z_min = z[np.argmin(np.abs(cum_sum-limits[0]/100.))]
     z_max = z[np.argmin(np.abs(cum_sum-limits[1]/100.))]
@@ -112,7 +112,7 @@ def main(pargs):
 
     # Finish
     print("")
-    print(f"Allowing for the MW halo, DM_MW_halo = {int(pargs.dm_mwhalo))} pc/cm^3")
+    print(f"Allowing for the MW halo, DM_MW_halo = {int(pargs.dm_mwhalo)} pc/cm^3")
     print(f"Allowing for the Host, DM_host = {int(pargs.dm_host)} pc/cm^3")
     print("")
     print("")
@@ -136,3 +136,4 @@ def main(pargs):
                flipy=True, known_hosts=False, title=pargs.fig_title, logz_scale=False)
 
     return z_min, z_max, z_50, z_mode
+