Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bug fix and added full run test #186

Merged
merged 20 commits into from
Sep 20, 2024
Merged
10 changes: 5 additions & 5 deletions .github/workflows/ci_tests.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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: |
Expand Down
16 changes: 16 additions & 0 deletions bin/frb_pzdm_mag
Original file line number Diff line number Diff line change
@@ -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)
5 changes: 3 additions & 2 deletions docs/installing.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:

Expand Down
Binary file added frb/data/DM/PDM_z.npz
Binary file not shown.
4 changes: 2 additions & 2 deletions frb/galaxies/eazy.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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!
Expand Down
7 changes: 6 additions & 1 deletion frb/galaxies/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from pkg_resources import resource_filename
import numpy as np
from scipy.interpolate import interp1d
import warnings

import pandas

Expand All @@ -20,7 +21,7 @@
from astropy import units

import pandas as pd
import extinction

from linetools.spectra import xspectrum1d

from frb import frb
Expand All @@ -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
Expand Down
73 changes: 60 additions & 13 deletions frb/scripts/limiting_mag.py
Original file line number Diff line number Diff line change
@@ -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

Expand All @@ -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()
Expand All @@ -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
Expand All @@ -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.))]
Expand All @@ -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

Expand Down
11 changes: 6 additions & 5 deletions frb/scripts/pz_dm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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.))]
Expand All @@ -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("")
Expand All @@ -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

Loading
Loading