Commit 2b0fd28e authored by Kenneth Rainer Corley's avatar Kenneth Rainer Corley Committed by Kenneth Corley
Browse files

Add cosmology feature to ligo-skymap-stats

parent f45db62d
......@@ -21,6 +21,9 @@ Changelog
``ligo-skymap-plot-pp-samples``, a tool for making P-P plots to compare a sky
map with posterior samples.
- Add ``--cosmology`` feature to ``ligo-skymap-stats`` to calculate comoving
volumes.
0.0.17 (2018-10-24)
===================
......
# -*- coding: utf-8 -*-
#
# Copyright (C) 2013-2018 Leo Singer, Rainer Corley
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
import numpy as np
import astropy.cosmology
import astropy.units as u
cosmo = astropy.cosmology.default_cosmology.get_cosmology_from_string(
'WMAP9')
def dVC_dVL_for_z(z):
"""Ratio between the comoving volume element dVC and a
DL**2 prior at redshift z."""
Ok0 = cosmo.Ok0
DH = cosmo.hubble_distance
DM_by_DH = (cosmo.comoving_transverse_distance(z) / DH).value
DC_by_DH = (cosmo.comoving_distance(z) / DH).value
zplus1 = z + 1.0
if Ok0 == 0.0:
ret = 1.0
elif Ok0 > 0.0:
ret = np.cosh(np.sqrt(Ok0) * DC_by_DH)
else: # Ok0 < 0.0 or Ok0 is nan
ret = np.cos(np.sqrt(-Ok0) * DC_by_DH)
ret *= zplus1
ret += DM_by_DH * cosmo.efunc(z)
ret *= np.square(zplus1)
return 1.0 / ret
@np.vectorize
def z_for_DL(DL):
return astropy.cosmology.z_at_value(
cosmo.luminosity_distance, DL * u.Mpc)
def dVC_dVL_for_DL(DL):
return dVC_dVL_for_z(z_for_DL(DL))
......@@ -28,6 +28,8 @@ from scipy.interpolate import interp1d
from .. import distance
from .. import moc
from .cosmology import dVC_dVL_for_DL
__all__ = ('find_injection_moc', 'FoundInjection')
......@@ -102,7 +104,8 @@ FoundInjection = namedtuple(
def find_injection_moc(sky_map, true_ra=None, true_dec=None, true_dist=None,
contours=(), areas=(), modes=False, nest=False):
contours=(), areas=(), modes=False, nest=False,
cosmology=False):
"""
Given a sky map and the true right ascension and declination (in radians),
find the smallest area in deg^2 that would have to be searched to find the
......@@ -240,6 +243,9 @@ def find_injection_moc(sky_map, true_ra=None, true_dec=None, true_dist=None,
dP[np.isnan(dP)] = 0 # Suppress invalid values
# Calculate probability density per unit volume.
if cosmology:
dV *= dVC_dVL_for_DL(r)
dP_dV = dP / dV
i = np.flipud(np.argsort(dP_dV.ravel()))
......
......@@ -126,11 +126,14 @@ def parser():
parser.add_argument(
'fitsfilenames', metavar='GLOB.fits[.gz]', nargs='+', action='glob',
help='Input FITS filenames and/or globs')
parser.add_argument(
'--cosmology', action='store_true',
help='Report volume localizations as comoving volumes.')
return parser
def startup(dbfilename, opts_contour, opts_modes, opts_area):
global db, contours, modes, areas
def startup(dbfilename, opts_contour, opts_modes, opts_area, opts_cosmology):
global db, contours, modes, areas, cosmology
if dbfilename is None:
db = None
else:
......@@ -138,6 +141,7 @@ def startup(dbfilename, opts_contour, opts_modes, opts_area):
contours = opts_contour
modes = opts_modes
areas = opts_area
cosmology = opts_cosmology
def process(fitsfilename):
......@@ -180,7 +184,7 @@ def process(fitsfilename):
searched_vol, searched_prob_vol, contour_vols
) = find_injection_moc(
sky_map, true_ra, true_dec, true_dist,
contours=contour_pvalues, areas=areas, modes=modes
contours=contour_pvalues, areas=areas, modes=modes, cosmology=cosmology
)
if snr is None:
......@@ -217,7 +221,7 @@ def main(args=None):
dbfilename = None
else:
dbfilename = sqlite.get_filename(opts.database)
args = (dbfilename, opts.contour, opts.modes, opts.area)
args = (dbfilename, opts.contour, opts.modes, opts.area, opts.cosmology)
if opts.jobs == 1:
pool_map = map
else:
......
......@@ -25,6 +25,8 @@ import astropy.units as u
import numpy as np
import scipy.misc
from ligo.skymap.postprocess.cosmology import dVC_dVL_for_z, z_for_DL, dVC_dVL_for_DL
parser = argparse.ArgumentParser()
parser.add_argument(
'cosmology', choices=astropy.cosmology.parameters.available,
......@@ -34,37 +36,6 @@ args = parser.parse_args()
cosmo = astropy.cosmology.default_cosmology.get_cosmology_from_string(
args.cosmology)
def dVC_dVL_for_z(z):
"""Ratio between the comoving volume element dVC and a
DL**2 prior at redshift z."""
Ok0 = cosmo.Ok0
DH = cosmo.hubble_distance
DM_by_DH = (cosmo.comoving_transverse_distance(z) / DH).value
DC_by_DH = (cosmo.comoving_distance(z) / DH).value
zplus1 = z + 1.0
if Ok0 == 0.0:
ret = 1.0
elif Ok0 > 0.0:
ret = np.cosh(np.sqrt(Ok0) * DC_by_DH)
else: # Ok0 < 0.0 or Ok0 is nan
ret = np.cos(np.sqrt(-Ok0) * DC_by_DH)
ret *= zplus1
ret += DM_by_DH * cosmo.efunc(z)
ret *= np.square(zplus1)
return 1.0 / ret
@np.vectorize
def z_for_DL(DL):
return astropy.cosmology.z_at_value(
cosmo.luminosity_distance, DL * u.Mpc)
def dVC_dVL_for_DL(DL):
return dVC_dVL_for_z(z_for_DL(DL))
DL = np.logspace(0, 6, 32)
log_DL = np.log(DL)
dVC_dVL = dVC_dVL_for_DL(DL)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment