diff --git a/.github/ci_orig_tests.yml b/.github/ci_orig_tests.yml new file mode 100644 index 00000000..14c06781 --- /dev/null +++ b/.github/ci_orig_tests.yml @@ -0,0 +1,85 @@ +name: CI Tests + +on: + push: + branches: + - master + tags: + pull_request: + +env: + SETUP_XVFB: True # avoid issues if mpl tries to open a GUI window + +jobs: + ci-tests: + name: Python-${{ matrix.python }}, deps=${{ matrix.deps }} + runs-on: ${{ matrix.os }} + if: "!(contains(github.event.head_commit.message, '[skip ci]') || contains(github.event.head_commit.message, '[ci skip]'))" + + strategy: + matrix: + os: [ubuntu-latest] + python: ['3.11', '3.12'] + deps: [current, numpy211, astropydev, numpydev, astropydev-numpydev] + + steps: + - name: Check out repository + uses: actions/checkout@v2 + - name: Set up Python ${{ matrix.python }} + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python }} + - name: Install base dependencies + run: | + python -m pip install --upgrade pip + - name: Test with numpy = 2.1.1 + if: "contains(matrix.deps, 'numpy211')" + run: | + python -m pip install numpy==2.1.1 + - name: Test with dev version of numpy + if: "contains(matrix.deps, 'numpydev')" + run: | + python -m pip install git+https://github.com/numpy/numpy.git#egg=numpy + - name: Test with dev version astropy + if: "contains(matrix.deps, 'astropydev')" + run: | + python -m pip install git+https://github.com/astropy/astropy.git#egg=astropy + - name: Install frb requirements + run: | + python -m pip install wheel scipy IPython + python -m pip install git+https://github.com/FRBs/ne2001.git#egg=ne2001 + python -m pip install git+https://github.com/FRBs/astropath.git#egg=astropath + python -m pip install git+https://github.com/linetools/linetools#egg=linetools + python -m pip install -r frb/requirements.txt + - name: Install numba + run: | + python -m pip install numba + - name: Install astroquery + run: | + python -m pip install astroquery + - name: Print Python, pip, astropy, numpy, and setuptools versions + run: | + python -c "import sys; print(f'Python {sys.version}')" + python -c "import pip; print(f'pip {pip.__version__}')" + python -c "import astropy; print(f'astropy {astropy.__version__}')" + python -c "import numpy; print(f'numpy {numpy.__version__}')" + python -c "import setuptools; print(f'setuptools {setuptools.__version__}')" + - name: Run tests + run: tox -e ${{ matrix.python }}-${{ matrix.toxenv }} + + codestyle: + runs-on: ubuntu-latest + if: "!contains(github.event.head_commit.message, '[ci skip]')" + steps: + - uses: actions/checkout@v2 + - name: Python codestyle check + uses: actions/setup-python@v2 + with: + python-version: 3.11 + - name: Install base dependencies + run: | + python -m pip install --upgrade pip + python -m pip install pycodestyle + - name: Check for runtime errors using pycodestyle + run: | + pycodestyle frb --count --select=E9 diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index 631e6dc6..368907f6 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -4,82 +4,33 @@ on: push: branches: - master - tags: pull_request: env: - SETUP_XVFB: True # avoid issues if mpl tries to open a GUI window + SETUP_XVFB: True # avoid issues if something tries to open a GUI window jobs: ci-tests: - name: Python-${{ matrix.python }}, deps=${{ matrix.deps }} + name: Tox env ${{ matrix.python }}-${{ matrix.toxenv }} runs-on: ${{ matrix.os }} - if: "!(contains(github.event.head_commit.message, '[skip ci]') || contains(github.event.head_commit.message, '[ci skip]'))" - strategy: matrix: os: [ubuntu-latest] python: ['3.11', '3.12'] - deps: [current, numpy211, astropydev, numpydev, astropydev-numpydev] - + toxenv: [test, test-astropydev] steps: - name: Check out repository - uses: actions/checkout@v2 + uses: actions/checkout@v3 - name: Set up Python ${{ matrix.python }} - uses: actions/setup-python@v2 + uses: actions/setup-python@v4 with: python-version: ${{ matrix.python }} - name: Install base dependencies run: | - python -m pip install --upgrade pip - - name: Test with numpy = 2.1.1 - if: "contains(matrix.deps, 'numpy211')" - run: | - python -m pip install numpy==2.1.1 - - name: Test with dev version of numpy - if: "contains(matrix.deps, 'numpydev')" - run: | - python -m pip install git+https://github.com/numpy/numpy.git#egg=numpy - - name: Test with dev version astropy - if: "contains(matrix.deps, 'astropydev')" - run: | - python -m pip install git+https://github.com/astropy/astropy.git#egg=astropy - - name: Install frb requirements - run: | - python -m pip install wheel scipy IPython + python -m pip install -r frb/requirements.txt --upgrade pip tox python -m pip install git+https://github.com/FRBs/ne2001.git#egg=ne2001 python -m pip install git+https://github.com/FRBs/astropath.git#egg=astropath python -m pip install git+https://github.com/linetools/linetools#egg=linetools - python -m pip install -r frb/requirements.txt - - name: Install numba - run: | - python -m pip install numba - - name: Install astroquery - run: | - python -m pip install astroquery - - name: Print Python, pip, astropy, numpy, and setuptools versions - run: | - python -c "import sys; print(f'Python {sys.version}')" - python -c "import pip; print(f'pip {pip.__version__}')" - python -c "import astropy; print(f'astropy {astropy.__version__}')" - python -c "import numpy; print(f'numpy {numpy.__version__}')" - python -c "import setuptools; print(f'setuptools {setuptools.__version__}')" - - name: Run tests - run: python setup.py test - - codestyle: - runs-on: ubuntu-latest - if: "!contains(github.event.head_commit.message, '[ci skip]')" - steps: - - uses: actions/checkout@v2 - - name: Python codestyle check - uses: actions/setup-python@v2 - with: - python-version: 3.11 - - name: Install base dependencies - run: | - python -m pip install --upgrade pip - python -m pip install pycodestyle - - name: Check for runtime errors using pycodestyle + - name: Test with tox run: | - pycodestyle frb --count --select=E9 + tox -e ${{ matrix.python }}-${{ matrix.toxenv }} diff --git a/bin/frb_mag_limit b/bin/frb_mag_limit deleted file mode 100755 index ba9ce1ef..00000000 --- a/bin/frb_mag_limit +++ /dev/null @@ -1,16 +0,0 @@ -#!/usr/bin/env python -# -# See top-level LICENSE file for Copyright information -# -# -*- coding: utf-8 -*- - - -""" -This script performs calculations related to the limiting mag -""" - -from frb.scripts import limiting_mag - -if __name__ == '__main__': - args = limiting_mag.parser() - limiting_mag.main(args) diff --git a/bin/frb_pz_dm b/bin/frb_pz_dm deleted file mode 100755 index ebc8aca9..00000000 --- a/bin/frb_pz_dm +++ /dev/null @@ -1,16 +0,0 @@ -#!/usr/bin/env python -# -# See top-level LICENSE file for Copyright information -# -# -*- coding: utf-8 -*- - - -""" -This script performs estimates P(z|DM) -""" - -from frb.scripts import pz_dm - -if __name__ == '__main__': - args = pz_dm.parser() - pz_dm.main(args) diff --git a/docs/installing.rst b/docs/installing.rst index 793e43f7..b1fe416b 100644 --- a/docs/installing.rst +++ b/docs/installing.rst @@ -32,7 +32,7 @@ to install and/or update these packages. * `healpy `_ version 1.15 or later * `pandas `_ version 1.5 or later * `requests `_ version 2.18 or later -* `extinction `_ version 0.4.6 or greater +* `dust_extinction `_ * `matplotlib `_ version 3.7 or greater * `linetools `_ version 0.3 or later * `astropath `_ version 0.1 or later diff --git a/frb/builds/build_hosts.py b/frb/builds/build_hosts.py index eb71ee33..b4e927ad 100644 --- a/frb/builds/build_hosts.py +++ b/frb/builds/build_hosts.py @@ -26,17 +26,13 @@ except: print('WARNING: ppxf not installed') from frb.galaxies import nebular +from frb.galaxies import utils as galaxy_utils from frb.galaxies import hosts from frb.galaxies import defs as galaxy_defs from frb.surveys import survey_utils from frb import utils import pandas -try: - import extinction -except ImportError: - print("extinction package not loaded. Extinction corrections will fail") - try: from frb.galaxies import cigale except ModuleNotFoundError: @@ -386,14 +382,9 @@ def run(host_input:pandas.core.series.Series, # Correct for Galactic extinction ebv = float(nebular.get_ebv(Host.coord)['meanValue']) print(f'Correcting the spectrum for Galactic extinction with reddening E(B-V)={ebv}') - AV = ebv * 3.1 # RV - Al = extinction.ccm89(spectrum.wavelength.value, AV, 3.1) - # New spec - new_flux = spectrum.flux * 10**(Al/2.5) - new_sig = spectrum.sig * 10**(Al/2.5) - new_spec = XSpectrum1D.from_tuple((spectrum.wavelength, new_flux, new_sig)) - - # + new_spec = galaxy_utils.deredden_spec(spectrum, ebv) + + # Run it ppxf.run(new_spec, R, host_input.z, results_file=ppxf_results_file, spec_fit=spec_file, diff --git a/frb/data/DM/CHIME_pzdm_repeaters.npz b/frb/data/DM/CHIME_pzdm_repeaters.npz new file mode 100644 index 00000000..7c324ad6 Binary files /dev/null and b/frb/data/DM/CHIME_pzdm_repeaters.npz differ diff --git a/frb/data/DM/CRAFT_ICS_1300_pzdm.npz b/frb/data/DM/CRAFT_ICS_1300_pzdm.npz new file mode 100644 index 00000000..8f76a6d0 Binary files /dev/null and b/frb/data/DM/CRAFT_ICS_1300_pzdm.npz differ diff --git a/frb/data/DM/CRAFT_ICS_1632_pzdm.npz b/frb/data/DM/CRAFT_ICS_1632_pzdm.npz new file mode 100644 index 00000000..c6923b97 Binary files /dev/null and b/frb/data/DM/CRAFT_ICS_1632_pzdm.npz differ diff --git a/frb/data/DM/CRAFT_ICS_892_pzdm.npz b/frb/data/DM/CRAFT_ICS_892_pzdm.npz new file mode 100644 index 00000000..edbb8733 Binary files /dev/null and b/frb/data/DM/CRAFT_ICS_892_pzdm.npz differ diff --git a/frb/data/DM/CRAFT_class_I_and_II_pzdm.npz b/frb/data/DM/CRAFT_class_I_and_II_pzdm.npz new file mode 100644 index 00000000..d863d8b5 Binary files /dev/null and b/frb/data/DM/CRAFT_class_I_and_II_pzdm.npz differ diff --git a/frb/data/DM/DSA_pzdm.npz b/frb/data/DM/DSA_pzdm.npz new file mode 100644 index 00000000..c77ea20c Binary files /dev/null and b/frb/data/DM/DSA_pzdm.npz differ diff --git a/frb/data/DM/FAST_pzdm.npz b/frb/data/DM/FAST_pzdm.npz new file mode 100644 index 00000000..41386524 Binary files /dev/null and b/frb/data/DM/FAST_pzdm.npz differ diff --git a/frb/data/DM/convert_npy_to_npz.py b/frb/data/DM/convert_npy_to_npz.py new file mode 100644 index 00000000..e2b024d9 --- /dev/null +++ b/frb/data/DM/convert_npy_to_npz.py @@ -0,0 +1,36 @@ +""" Simple script to convert .npy files to .npz files. """ + +import numpy as np + +from frb.dm import prob_dmz + +def main(): + + npy_files = ['DSA_pzdm.npy', + 'parkes_mb_class_I_and_II_pzdm.npy', + 'CRAFT_class_I_and_II_pzdm.npy', + 'CRAFT_ICS_1300_pzdm.npy', + 'CRAFT_ICS_892_pzdm.npy', + 'CRAFT_ICS_1632_pzdm.npy', + 'FAST_pzdm.npy'] + + # Load of z, DM_EG from CHIME + chime_dict = prob_dmz.grab_repo_grid(prob_dmz.telescope_dict['CHIME']) + + # Parse + z = chime_dict['z'] + DM = chime_dict['DM'] + + # Loop on the files + for npy_file in npy_files: + # Load + pzdm = prob_dmz.grab_repo_grid(npy_file) + # Save + npz_file = npy_file.replace('.npy', '.npz') + np.savez(npz_file, pzdm=pzdm, z=z, DM=DM) + # Report + print(f"Wrote: {npz_file}") + +# Command line execution +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/frb/data/DM/parkes_mb_class_I_and_II_pzdm.npz b/frb/data/DM/parkes_mb_class_I_and_II_pzdm.npz new file mode 100644 index 00000000..bb15f934 Binary files /dev/null and b/frb/data/DM/parkes_mb_class_I_and_II_pzdm.npz differ diff --git a/frb/dm/prob_dmz.py b/frb/dm/prob_dmz.py index 7a2becb0..e1316f37 100644 --- a/frb/dm/prob_dmz.py +++ b/frb/dm/prob_dmz.py @@ -14,6 +14,20 @@ from IPython import embed +# Dict to resolve telescope/survey into the appropriate grid filename +telescope_dict = { + 'CHIME_repeaters': 'CHIME_pzdm_repeaters.npz', + 'CHIME': 'CHIME_pzdm.npz', + 'DSA': 'DSA_pzdm.npz', + 'Parkes': 'parkes_mb_class_I_and_II_pzdm.npz', + 'CRAFT': 'CRAFT_class_I_and_II_pzdm.npz', + 'CRAFT_ICS_1300': 'CRAFT_ICS_1300_pzdm.npz', + 'CRAFT_ICS_892': 'CRAFT_ICS_892_pzdm.npz', + 'CRAFT_ICS_1632': 'CRAFT_ICS_1632_pzdm.npz', + 'FAST': 'FAST_pzdm.npz', + 'perfect': 'PDM_z.npz' +} + class P_DM_z(object): pass diff --git a/frb/galaxies/nebular.py b/frb/galaxies/nebular.py index d3924647..fd258fe7 100644 --- a/frb/galaxies/nebular.py +++ b/frb/galaxies/nebular.py @@ -12,10 +12,7 @@ from astropy.table import Table from astropy import units -try: - import extinction -except ImportError: - warnings.warn("extinction package not loaded. Extinction corrections will fail") +import dust_extinction try: from linetools.lists import linelist @@ -33,7 +30,7 @@ def calc_dust_extinct(neb_lines, method): """ Estimate the Visual extinction A_V based on input nebular emission lines - Uses the Fitzpatrick & Massa (2007) extinction law + Uses the Gordon 2024 Args: neb_lines (dict): Line fluxes @@ -71,8 +68,11 @@ def calc_dust_extinct(neb_lines, method): raise IOError("Not ready for this mode") # Extinction - a1AV = extinction.fm07(np.atleast_1d(wave1), 1.0)[0] - a2AV = extinction.fm07(np.atleast_1d(wave2), 1.0)[0] + #a1AV = extinction.fm07(np.atleast_1d(wave1), 1.0)[0] + #a2AV = extinction.fm07(np.atleast_1d(wave2), 1.0)[0] + extmod = dust_extinction.parameter_averages.G23(Rv=3.1) + a1AV = extmod(np.atleast_1d(wave1*units.AA)) # *units.AA) + a2AV = extmod(np.atleast_1d(wave2*units.AA)) # *units.AA) # Observed ratio fratio_obs = F1_obs/F2_obs @@ -153,7 +153,10 @@ def calc_lum(neb_lines, line, z, cosmo, AV=None): # Dust correct? if AV is not None: - al = extinction.fm07(np.atleast_1d(wave.to('Angstrom').value), AV)[0] + #al = extinction.fm07(np.atleast_1d(wave.to('Angstrom').value), AV)[0] + extmod = dust_extinction.parameter_averages.G23(Rv=3.1) + AlAV = extmod(np.atleast_1d(wave)*units.AA)[0] + al = AlAV * AV else: al = 0. diff --git a/frb/galaxies/photom.py b/frb/galaxies/photom.py index 61452452..ea091681 100644 --- a/frb/galaxies/photom.py +++ b/frb/galaxies/photom.py @@ -2,6 +2,7 @@ import os import warnings +import dust_extinction.parameter_averages import numpy as np from pkg_resources import resource_filename @@ -22,10 +23,7 @@ from frb.galaxies import defs -try: - import extinction -except ImportError: - print("extinction package not loaded. Extinction corrections will fail") +import dust_extinction # Photometry globals table_format = 'ascii.fixed_width' @@ -118,7 +116,7 @@ def extinction_correction(filt, EBV, RV=3.1, max_wave=None, required=True): """ calculate MW extinction correction for given filter - Uses the Fitzpatrick & Massa (2007) extinction law + Uses the Gordon 2024 extinction model Args: filt (str): @@ -173,7 +171,14 @@ def extinction_correction(filt, EBV, RV=3.1, max_wave=None, required=True): #get MW extinction correction AV = EBV * RV #AlAV = nebular.load_extinction('MW') - Alambda = extinction.fm07(wave, AV) + #Alambda = extinction.fm07(wave, AV) + + # Gordon 2024 + extmod = dust_extinction.parameter_averages.G23(Rv=RV) + AlAV = extmod(wave*units.AA) + Alambda = AlAV * AV + + source_flux = 1. #calculate linear correction delta = np.trapz(throughput * source_flux * diff --git a/frb/galaxies/utils.py b/frb/galaxies/utils.py index 5a8b3597..1604e34c 100644 --- a/frb/galaxies/utils.py +++ b/frb/galaxies/utils.py @@ -22,6 +22,8 @@ import pandas as pd +import dust_extinction + from linetools.spectra import xspectrum1d from frb import frb @@ -38,12 +40,13 @@ 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) + extmod = dust_extinction.parameter_averages.G23(Rv=3.1) + AlAV = extmod(spectrum.wavelength)#*units.AA) + Al = AlAV * AV + #Al = extinction.ccm89(spectrum.wavelength.value, AV, 3.1) + # New spec new_flux = spectrum.flux * 10**(Al/2.5) new_sig = spectrum.sig * 10**(Al/2.5) diff --git a/frb/requirements.txt b/frb/requirements.txt index a2061883..c486a83f 100644 --- a/frb/requirements.txt +++ b/frb/requirements.txt @@ -9,10 +9,13 @@ pytest>=6.0.0 pandas>=1.2 photutils>=1.1 requests>=2.18 -extinction>=0.4.2 +dust_extinction>=1.4 matplotlib>=3.3 healpy>=1.12 linetools>=0.3.1 PyYAML>=5.1 numba>=0.50.0 +setuptools>=75.0.0 +tox>=4.15.0 +pyvo>=1.5.3 diff --git a/frb/scripts/galaxies.py b/frb/scripts/galaxies.py index b5a107af..1d867eca 100644 --- a/frb/scripts/galaxies.py +++ b/frb/scripts/galaxies.py @@ -30,6 +30,7 @@ def plot_spec(specDB, meta, frb, dust_correct=False): from astropy.coordinates import SkyCoord from linetools import utils as ltu from frb.galaxies import nebular + from frb.galaxies import utils as galaxy_utils """ Plot galaxy spectra """ # Grab spectra @@ -37,15 +38,16 @@ def plot_spec(specDB, meta, frb, dust_correct=False): # Dust? if dust_correct: - import extinction EBV = nebular.get_ebv(frb.coord)['meanValue'] # AV = EBV * 3.1 # RV for kk in range(spec.nspec): spec.select = kk - Al = extinction.fm07(spec.wavelength.value, AV)#, 3.1) - spec.flux = spec.flux * 10 ** (Al / 2.5) - spec.sig = spec.sig * 10 ** (Al / 2.5) + embed(header='plot_spec: THIS IS BROKEN') + new_spec = galaxy_utils.deredden_spec(spec.select, AV) + #Al = extinction.fm07(spec.wavelength.value, AV)#, 3.1) + #spec.flux = spec.flux * 10 ** (Al / 2.5) + #spec.sig = spec.sig * 10 ** (Al / 2.5) # Add labels lbls = ['None']*len(meta) diff --git a/frb/scripts/limiting_mag.py b/frb/scripts/limiting_mag.py deleted file mode 100644 index f33fa392..00000000 --- a/frb/scripts/limiting_mag.py +++ /dev/null @@ -1,148 +0,0 @@ -#!/usr/bin/env python -""" -Estimate the limiting luminosity given the magnitude limit and DMs and coord and telescope -""" -from IPython import embed - -def parser(options=None): - import argparse - # Parse - parser = argparse.ArgumentParser(description='Script to print a summary of an FRB to the screen [v1.0]') - parser.add_argument("coord", type=str, help="Coordinates, e.g. J081240.7+320809 or 122.223,-23.2322 or 07:45:00.47,34:17:31.1 or FRB name (FRB180924)") - 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_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() - else: - pargs = parser.parse_args(options) - return pargs - - -def main(pargs): - """ Run - """ - import numpy as np - - from linetools import utils as ltu - from linetools.scripts.utils import coord_arg_to_coord - - from frb.galaxies import nebular - from frb.galaxies import photom - from frb.galaxies import utils as frb_gal_u - from frb import mw - from frb.dm import prob_dmz - - # Deal with coord - icoord = ltu.radec_to_coord(coord_arg_to_coord(pargs.coord)) - - # EBV - EBV = nebular.get_ebv(icoord)['meanValue'] # - print(f"EBV = {EBV}") - - # NE 2001 - DM_ISM = mw.ismDM(icoord) - print(f"NE2001 = {DM_ISM}") - - # 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 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'] - - # 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 = 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.))] - - # Setup Luminosity - - # Extinction correct - dust_correct = photom.extinction_correction(pargs.filter, EBV) - mag_dust = 2.5 * np.log10(1. / dust_correct) - mag_corr = pargs.mag_limit + mag_dust - - # ##########################3 - # Convert to L - - # Load f_mL - f_mL = frb_gal_u.load_f_mL() - # m_r(L*) - m_r_Lstar_min = float(f_mL(z_min)) - m_r_Lstar_max = float(f_mL(z_max)) - - frac_Lstar_min = 10**(-0.4*(mag_corr-m_r_Lstar_min)) - frac_Lstar_max = 10**(-0.4*(mag_corr-m_r_Lstar_max)) - - # Finish - print("-----------------------------------------------------") - 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_mag_limit J151849.52+122235.8 200. 23. - \ No newline at end of file diff --git a/frb/scripts/pz_dm.py b/frb/scripts/pz_dm.py deleted file mode 100644 index f8c39a45..00000000 --- a/frb/scripts/pz_dm.py +++ /dev/null @@ -1,139 +0,0 @@ - -""" -Estimate p(z|DM) for an assumed location on the sky and DM_FRB -Defaults to using a perfect telescope model for the DM-z grid -""" -from IPython import embed - - -def parser(options=None): - import argparse - # Parse - parser = argparse.ArgumentParser(description='Script to print a summary of an FRB to the screen [v1.0]') - parser.add_argument("coord", type=str, help="Coordinates, e.g. J081240.7+320809 or 122.223,-23.2322 or 07:45:00.47,34:17:31.1 or FRB name (FRB180924)") - 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=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") - 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("--fig_title", type=str, help="title for the figure; e.g., FRBXXXXX") - - if options is None: - pargs = parser.parse_args() - else: - pargs = parser.parse_args(options) - return pargs - - -def main(pargs): - """ Run - """ - import numpy as np - - from linetools import utils as ltu - from linetools.scripts.utils import coord_arg_to_coord - - from frb import mw - from frb.dm import prob_dmz - from frb.galaxies import mag_dm - - - # Deal with coord - icoord = ltu.radec_to_coord(coord_arg_to_coord(pargs.coord)) - - # NE 2001 - DM_ISM = mw.ismDM(icoord) - print("") - print("-----------------------------------------------------") - print(f"NE2001 = {DM_ISM:.2f}") - - # 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 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'] - - # 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 = [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.))] - - z_50 = z[np.argmin(np.abs(cum_sum-50./100.))] - z_mode = z[np.argmax(PzDM)] - - # Finish - 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("") - print("") - print(f"The mean redshift value is: {z_50:.3f}") - print(f"The mode redshift value is: {z_mode:.3f}") - print("") - print(f"The redshift range for your confidence interval {pargs.cl} is:") - print(f"z = [{z_min:.3f}, {z_max:.3f}]") - 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 (Macqurt+2020)") - else: - print("This assumes the "+(str(pargs.telescope))+" telescope and a model of the scatter in DM_cosmic (Macqurt+2020)") - print("-----------------------------------------------------") - - - - # make the magnitude vs redshift plot with z-range if requested - if pargs.magdm_plot: - mag_dm.r_vs_dm_figure(z_min, z_max, z, PzDM, outfile='fig_r_vs_z.png', - flipy=True, known_hosts=False, title=pargs.fig_title, logz_scale=False) - - return z_min, z_max, z_50, z_mode - diff --git a/frb/scripts/pzdm_mag.py b/frb/scripts/pzdm_mag.py index c7190b7a..c7314e11 100644 --- a/frb/scripts/pzdm_mag.py +++ b/frb/scripts/pzdm_mag.py @@ -66,52 +66,26 @@ def main(pargs): DM_cosmic = DM_extragalactic - pargs.dm_host - # Redshift estimates - # 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' - } + telescope_dict = prob_dmz.telescope_dict # 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'] - - # 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']) + if not pargs.telescope or pargs.telescope == 'perfect': + sdict = prob_dmz.grab_repo_grid(telescope_dict['perfect']) + DM = sdict['DM'] + PDM_z = sdict['PDM_z'] + PDM_z = PDM_z.T # Perfect is opposite of the rest + iDM = np.argmin(np.abs(DM - DM_cosmic)) + else: # Grab a non-perfect telescope grid, which use DM_extragalactic + sdict = prob_dmz.grab_repo_grid(telescope_dict[pargs.telescope]) 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]) + # Grab the right entry + z = sdict['z'] + PzDM = PDM_z[:,iDM] / np.sum(PDM_z[:,iDM]) + cum_sum = np.cumsum(PzDM) limits = [float(item) for item in pargs.cl.split(',')] diff --git a/frb/tests/test_frbgalaxies.py b/frb/tests/test_frbgalaxies.py index ca93ddaa..013be8d1 100644 --- a/frb/tests/test_frbgalaxies.py +++ b/frb/tests/test_frbgalaxies.py @@ -18,6 +18,7 @@ from frb.frb import FRB from frb.galaxies import mag_dm +from frb.dm import prob_dmz def data_path(filename): data_dir = os.path.join(os.path.dirname(__file__), 'files') @@ -154,19 +155,19 @@ def test_mag_dm_figure(): def test_pzdm_telescopes(): - telescope_dict = { - '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', - } + telescope_dict = prob_dmz.telescope_dict + #telescope_dict = { + # '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', + #} # Load the CHIME grid - from frb.dm import prob_dmz - sdict = prob_dmz.grab_repo_grid('CHIME_pzdm.npz') + sdict = prob_dmz.grab_repo_grid(telescope_dict['CHIME'])#'CHIME_pzdm.npz') PDM_z = sdict['pzdm'] z = sdict['z'] DM = sdict['DM'] @@ -177,7 +178,7 @@ def test_pzdm_telescopes(): assert PDM_z.shape == (500, 1400) # check the perfect grid - sdict = prob_dmz.grab_repo_grid('PDM_z.npz') + sdict = prob_dmz.grab_repo_grid(telescope_dict['perfect']) PDM_z = sdict['PDM_z'] z = sdict['z'] DM = sdict['DM'] @@ -189,8 +190,10 @@ def test_pzdm_telescopes(): # check the other grids for key in telescope_dict.keys(): - PDM_z = prob_dmz.grab_repo_grid(telescope_dict[key]) - # Test + if key == 'perfect': + continue + grid = prob_dmz.grab_repo_grid(telescope_dict[key]) + PDM_z = grid['pzdm'] assert PDM_z.shape == (500, 1400) @@ -202,8 +205,9 @@ def test_pzdm_telescopes(): assert isinstance(z_max, float) assert isinstance(z_50, float) assert isinstance(z_mode, float) - assert isinstance(frac_Lstar_min, float) - assert isinstance(frac_Lstar_max, float) + + assert np.isclose(frac_Lstar_min, 0.00804, atol=1e-4) + assert np.isclose(frac_Lstar_max, 2.88776, atol=1e-4) assert z_min < z_max assert 0 <= z_50 <= 1 assert 0 <= z_mode <= 1 diff --git a/frb/tests/test_photom.py b/frb/tests/test_photom.py index 52cfe403..37407271 100644 --- a/frb/tests/test_photom.py +++ b/frb/tests/test_photom.py @@ -20,13 +20,10 @@ from frb.surveys.catalog_utils import convert_mags_to_flux -# TODO -- Turn this test back on when we replace extinction -extinction = pytest.mark.skipif(True, reason='Extinction module needs to be replaced.') -@extinction def test_dust_correct(): correct = photom.extinction_correction('GMOS_S_r', 0.138) - assert np.isclose(correct, 1.3818590723917497) + assert np.isclose(correct, 1.3869201954307397) def test_flux_conversion(): diff --git a/frb/tests/test_scripts.py b/frb/tests/test_scripts.py index 68fc7469..58ac4330 100644 --- a/frb/tests/test_scripts.py +++ b/frb/tests/test_scripts.py @@ -11,6 +11,8 @@ from frb.scripts import pzdm_mag from frb.scripts import tns +remote_data = pytest.mark.skipif(os.getenv('FRB_GDB') is None, + reason='test requires dev suite') def test_frb_summary(): pargs = frb_summary.parser(['180924']) @@ -25,8 +27,9 @@ def test_frb_pzdm_mag(): assert np.isclose(zmax, 0.16080402010050251) assert np.isclose(z_50, 0.10050251256281408) assert np.isclose(z_mode, 0.12060301507537688) - assert np.isclose(Lmax, 0.023780306538345886) + assert np.isclose(Lmax, 0.023809, atol=1e-4) +@remote_data def test_tns(): tns.parser.units = 'deg' tns.parser.radius = 0.5 diff --git a/tox.ini b/tox.ini new file mode 100644 index 00000000..4e46082f --- /dev/null +++ b/tox.ini @@ -0,0 +1,77 @@ +[tox] +envlist = + {3.10,3.11,3.12}-test{,-alldeps} + {3.10,3.11,3.12}-test-numpy{125,126} + {3.10,3.11,3.12}-test-{numpy,astropy}dev + codestyle +requires = + setuptools >= 50.3.0 + pip >= 19.3.1 +isolated_build = true + +[testenv] +# pytest +allowlist_externals = pytest +# Suppress display of matplotlib plots generated during docs build +setenv = + MPLBACKEND=agg + numpydev: PIP_EXTRA_INDEX_URL = https://pypi.anaconda.org/scipy-wheels-nightly/simple + PYTHONPATH = {toxinidir} + + +# Pass through the following environment variables which may be needed for the CI +passenv = HOME,WINDIR,LC_ALL,LC_CTYPE,CC,CI + +# Run the tests in a temporary directory to make sure that we don't import +# this package from the source tree +#changedir = .tmp/{envname} + +# tox environments are constructed with so-called 'factors' (or terms) +# separated by hyphens, e.g. test-devdeps-cov. Lines below starting with factor: +# will only take effect if that factor is included in the environment name. To +# see a list of example environments that can be run, along with a description, +# run: +# +# tox -l -v +# +description = + run tests + alldeps: with all optional dependencies + devdeps: with the latest developer version of key dependencies + oldestdeps: with the oldest supported version of key dependencies + cov: and test coverage + numpy123: with numpy 1.23.* + numpy124: with numpy 1.24.* + numpy125: with numpy 1.25.* + numpy126: with numpy 1.26.* + +# The following provides some specific pinnings for key packages +deps = + + cov: coverage + numpy123: numpy==1.23.* + numpy124: numpy==1.24.* + numpy125: numpy==1.25.* + numpy126: numpy==1.26.* + + numpydev: numpy>=0.0.dev0 + astropydev: git+https://github.com/astropy/astropy.git#egg=astropy + + linetoolsdev: git+https://github.com/linetools/linetools.git#egg=linetools + +# The following indicates which extras_require from setup.cfg will be installed +extras = + test + alldeps + +commands = + pip freeze + !cov: pytest --pyargs frb {posargs} + cov: pytest --pyargs frb --cov frb --cov-config={toxinidir}/setup.cfg {posargs} + cov: coverage xml -o {toxinidir}/coverage.xml + +[testenv:conda] +description = run tests in environment created via conda +requires = tox-conda +conda_env = {toxinidir}/environment.yml +commands = pytest --pyargs frb {posargs}