Commit de9d108e authored by Riccardo Barbieri's avatar Riccardo Barbieri
Browse files

Merge branch 'pending_review' into 'schech_params_commandline'

# Conflicts:
#   bin/gwcosmo_single_posterior
parents 2df59052 0c204878
......@@ -340,9 +340,7 @@ else:
if opts.method == 'statistical':
#catalog = gwcosmo.prior.catalog.galaxyCatalog(catalog_file=catalog_path,skymap=skymap,band=band,thresh=sky_area,Kcorr=Kcorr,nside=nside)
me = gwcosmo.gwcosmo.WholeSkyGalaxyCatalogLikelihood(catalog, skymap, band, sp, cosmo, px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, ps_M, Kcorr=Kcorr, mth=mth, zcut=zcut, zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, sky_thresh = opts.sky_area, nside=opts.nside)
me = gwcosmo.gwcosmo.WholeSkyGalaxyCatalogLikelihood(catalog, skymap, band, cosmo, px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, ps_M, Kcorr=Kcorr, mth=mth, zcut=zcut, zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, sky_thresh = opts.sky_area, nside=opts.nside)
if opts.method == 'pixel':
......@@ -355,10 +353,15 @@ else:
raise ValueError(f'Low resolution nside {nside_low_res} is higher than high resolution nside {nside}. Try decreasing min_pixels or increasing nside.')
if return_skymap_indices:
# Make the catalog cache file
if not catalog.read_pixel_index_cache(nside_low_res):
print(f'No pixel index found for nside {nside_low_res}, generating now.')
catalog.build_pixel_index_file(nside_low_res)
with open(outputfile+f'_indices.txt','w') as file:
file.write(f'nside: {nside_low_res}\n')
for i in pixelated_samples.indices:
file.write(f'{i}\n')
exit()
# Select the coarse pixel from the skymap that we will use
......
......@@ -144,8 +144,6 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
mth : float, optional
Specify an apparent magnitude threshold for the galaxy catalogue
(default=None). If none, mth is estimated from the galaxy catalogue.
zcut : float, optional
An artificial redshift cut to the galaxy catalogue (default=None)
zmax : float, optional
The upper redshift limit for integrals (default=10.). Should be well
beyond the highest redshift reachable by GW data or selection effects.
......@@ -246,7 +244,6 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
numerator and denominator
"""
# TODO: Move into the catalog class
if self.Kcorr:
Kcorr = self.full_catalog.get_k_correction(self.band, sampz, color_names[self.band], sampcolor)
else:
......@@ -553,11 +550,26 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
self.full_catalog = galaxy_catalog
self.path = outputfile+'_'+str(pixel_index)+'_checkpoint.p'
# Set redshift and colour limits based on whether Kcorrections are applied
# Set redshift and colour limits based on whether Kcorrections are applied
if Kcorr == True:
if zcut is None:
self.zcut = 0.5
if observation_band == 'W1':
# Polynomial k corrections out to z=1
self.zcut = 1.0
else:
# color-based k corrections valid to z=0.5
self.zcut = 0.5
else:
if observation_band == 'W1' and zcut > 1.0:
print(f"Warning, your requested zcut {zcut} is greater than the valid range (1.0) for W1-band k corrections")
elif zcut > 0.5:
print(f"Warning, your requested zcut {zcut} is greater than the valid range (0.5) for k corrections")
else:
# zcut is < valid for k-corr, do nothing
pass
self.full_catalog = self.full_catalog.apply_color_limit(observation_band,
*color_limits[color_names[observation_band]])
*color_limits[color_names[observation_band]])
else:
if zcut is None:
self.zcut = self.zmax
......@@ -731,7 +743,10 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
dec = subcatalog['dec']
m = subcatalog.get_magnitudes(self.band)
sigmaz = subcatalog['sigmaz']
color = subcatalog.get_color(self.band) # TODO: Fix based on Kcorr
if self.Kcorr:
color = subcatalog.get_color(self.band)
else:
color = np.zeros(len(m))
if len(subcatalog)==0:
self.pxG[:,i] = np.zeros(len(H0))
......@@ -828,7 +843,21 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
# Set redshift and colour limits based on whether Kcorrections are applied
if Kcorr == True:
if zcut is None:
self.zcut = 0.5
if observation_band == 'W1':
# Polynomial k corrections out to z=1
self.zcut = 1.0
else:
# color-based k corrections valid to z=0.5
self.zcut = 0.5
else:
if observation_band == 'W1' and zcut > 1.0:
print(f"Warning, your requested zcut {zcut} is greater than the valid range (1.0) for W1-band k corrections")
elif zcut > 0.5:
print(f"Warning, your requested zcut {zcut} is greater than the valid range (0.5) for k corrections")
else:
# zcut is < valid for k-corr, do nothing
pass
self.full_catalog = self.full_catalog.apply_color_limit(observation_band,
*color_limits[color_names[observation_band]])
else:
......@@ -922,7 +951,11 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
galdec = subcatalog['dec']
galm = subcatalog.get_magnitudes(self.band)
galsigmaz = subcatalog['sigmaz']
galcolor = subcatalog.get_color(self.band) # TODO: Fix based on Kcorr
if self.Kcorr:
color = subcatalog.get_color(self.band)
else:
color = np.zeros(len(galm))
print('Computing the in-catalogue part')
self.pxG, self.pDG = self.pxD_GH0_multi(H0, galz, galsigmaz, galm, galra,
galdec, galcolor, nfine=self.nfine,
......
import os
import glob
import pickle
import math
import array
......@@ -14,12 +15,14 @@ from ..utilities.cache import get_cachedir
from ..utilities.calc_kcor import calc_kcor
DEG2RAD = math.pi/180.0
color_names = {'B': None, 'K': None, 'u': 'u - r', 'g': 'g - r', 'r': 'g - r',
'i': 'g - i', 'z': 'r - z', 'W1': None}
# Taken from limits of K corrections from fig 4 of
# https://ui.adsabs.harvard.edu/abs/2010MNRAS.405.1409C/abstract
color_limits = {'u - r': [-0.1, 2.9], 'g - r': [-0.1, 1.9], 'g - i': [0, 3],
'r - z': [0, 1.5], None : [-np.inf, np.inf]}
Kcorr_bands = {'B': 'B', 'K': 'K', 'u': 'r', 'g': 'r', 'r': 'g', 'i': 'g', 'z': 'r'}
Kcorr_signs = {'B': 1, 'K': 1, 'u': 1, 'g': 1, 'r': -1, 'i': -1, 'z': -1}
# istarmap.py for Python 3.8+
import multiprocessing.pool as mpp
......@@ -98,17 +101,6 @@ def load_catalog_from_opts(opts):
band = opts.catalog_band
return load_catalog(name, band)
# Load original data
# Load index files
# Filter by pixel
# Discard non-galaxies
# SDSS-standardise (see https://git.ligo.org/cbc-cosmo/O2-H0/-/wikis/K-corrections-and-galaxy-catalog-review)
# Discard missing z vals
# Filter by z
# Filter by color
# Apply K corrections
# Return catalog
from abc import ABC, abstractmethod
class GalaxyCatalog:
......@@ -118,7 +110,7 @@ class GalaxyCatalog:
colnames = {'z','sigmaz','ra','dec'}
def __init__(self, data = None, name = 'Unknown Catalog',
supported_bands=None, cachedir=None, Kcorr = False):
supported_bands=None, cachedir=None):
self.data = data
self.name = name
self.supported_bands = supported_bands
......@@ -126,7 +118,6 @@ class GalaxyCatalog:
self.colnames.union([f'm_{band}' for band in supported_bands])
# Cache for pixel index to array index lookup
self.pixmap = {}
self.Kcorr = Kcorr # Whether this catalog has k-corrections applied
def __getitem__(self, *args, **kwargs):
return self.data.__getitem__(*args, **kwargs)
......@@ -169,6 +160,20 @@ class GalaxyCatalog:
else:
return False
def clean_cache(self, mtime, cachedir=None):
"""
Remove any index cache files older than mtime
"""
if cachedir is None:
cachedir = get_cachedir()
for f in glob.glob(cachedir+'/'+self.name+'_*'):
try:
if os.path.getmtime(f) < mtime:
os.remove(f)
except FileNotFoundError:
# Here in case a parallel job has removed the file
pass
def magnitude_thresh(self, band, ra=None, dec=None):
"""
Return the magnitude threshold for a specific
......@@ -205,8 +210,7 @@ class GalaxyCatalog:
idx = pixmap[pixel_index]
new = GalaxyCatalog(data = self.data[idx], name = self.name+f'_nside{nside}_pixel{pixel_index}',
supported_bands = self.supported_bands,
Kcorr = self.Kcorr)
supported_bands = self.supported_bands)
return new
def idx2pixdict(self, nside, idx, nested=True):
......@@ -223,18 +227,17 @@ class GalaxyCatalog:
def get_color(self, band):
"""
Return color index (??)
TODO: CHECK THIS AND UPDATE
Return color index for K corrections
"""
Kcorr_bands = {'B': 'B', 'K': 'K', 'u': 'r', 'g': 'r', 'r': 'g', 'i': 'g', 'z': 'r', 'W1': 'W1'}
if self.Kcorr:
Kcorr_bands = {'B': 'B', 'K': 'K', 'u': 'r', 'g': 'r', 'r': 'g', 'i': 'g', 'z': 'r'}
Kcorr_signs = {'B': 1, 'K': 1, 'u': 1, 'g': 1, 'r': -1, 'i': -1, 'z': -1}
if not band in Kcorr_bands.keys():
print('Band {band} not supported for color-based K corrections')
color = np.zeros(len(self))
else:
m = self.get_magnitudes(band)
m_K = self.get_magnitudes(Kcorr_bands[band])
color = Kcorr_signs[band]*(m - m_K)
else:
color = np.zeros(len(self))
return color
def apply_redshift_cut(self, zcut):
......@@ -246,9 +249,13 @@ class GalaxyCatalog:
m = self.get_magnitudes(band)
idx = np.where(np.isfinite(m))
return GalaxyCatalog(data = self.data[idx], name = self.name+f'_valid_{band}_mags',
supported_bands = [band])
supported_bands = self.supported_bands)
def apply_color_limit(self, band, cmin, cmax):
"""
Apply the given color limits to the given band.
Will filter out any galaxies whose magnitude is missing in given band.
"""
if band == 'W1':
print('Not applying color limits for W1 band, as we use the z-dependent k-correction')
return self
......@@ -300,66 +307,22 @@ def pixelate(cat, nside, allowed_pixels=None, nested=True):
pixlists[k]=np.sort(v)
return pixlists
class DESI(GalaxyCatalog):
"""
DESI data from data release 8
"""
supports_kcorrections = False # K corrections not implemented yet for DESI
supported_bands = {'G', 'R', 'Z', 'W1', 'W2'}
def __init__(self,
catalog_file='LS_DR8_total_csp.fits',
band='G',
Kcorr = False
):
super().__init__(band = band, Kcorr=Kcorr, name = "DESI")
self.band = band
assert band in self.supported_bands
self.filename = catalog_file
self.OmegaG = 1.0
self.px_OmegaG = 1.0
#self.gal_indices_per_pixel = self.pixelate()
self.populate()
def populate(self):
from astropy.io.fits import open
with open(self.filename, 'readonly', memmap=True) as desifile:
names = desifile[1].data.names
self.data = np.rec.fromarrays([desifile[1].data[n] if (n!='ra' and n!='dec') else desifile[1].data[n]*DEG2RAD for n in names], names)
def __getstate__(self):
state = self.__dict__.copy()
del state['data']
return state
def __setstate__(self, state):
self.__dict__ = state
self.populate()
class OldStyleCatalog(GalaxyCatalog):
"""
Catalog in the old GWCosmo format. Must have been
preprocessed.
preprocessed into HDF5 files.
"""
filename = None
def __init__(self,
catalog_file=None,
Kcorr = False,
name = None):
self.filename = get_catalogfile(catalog_file)
self.OmegaG = 1.0
self.px_OmegaG = 1.0
super().__init__(supported_bands = self.supported_bands, Kcorr=Kcorr, name = name)
super().__init__(supported_bands = self.supported_bands, name = name)
self.populate()
self.clean_cache(os.path.getmtime(self.filename))
def populate(self):
"""
......@@ -385,31 +348,27 @@ class OldStyleCatalog(GalaxyCatalog):
class OldStyleDESI(OldStyleCatalog):
supported_bands = {'g', 'W1'}
supports_kcorrections = True
def __init__(self, catalog_file = 'DESI.hdf5', band='W1', Kcorr=True):
supported_bands = {'W1'}
def __init__(self, catalog_file = 'DESI.hdf5', band='W1'):
print('WARNING: DESI catalog is not fully supported yet')
self.colnames = set(self.colnames).union([f'm_{b}' for b in self.supported_bands])
super().__init__(catalog_file = catalog_file, Kcorr=Kcorr, name = 'DESI')
super().__init__(catalog_file = catalog_file, name = 'DESI')
class OldStyleGLADEPlus(OldStyleCatalog):
supported_bands = {'B', 'K', 'W1'}
supports_kcorrections = False
def __init__(self, catalog_file = 'glade+.hdf5', band='W1', Kcorr=False):
def __init__(self, catalog_file = 'glade+.hdf5', band='W1'):
self.colnames = set(self.colnames).union([f'm_{b}' for b in self.supported_bands])
super().__init__(catalog_file = catalog_file, Kcorr=Kcorr, name = 'GladePlus')
super().__init__(catalog_file = catalog_file, name = 'GladePlus')
class OldStyleGLADE(OldStyleCatalog):
supported_bands = {'B','K'}
supports_kcorrections = False
def __init__(self, catalog_file = 'glade.hdf5', band='B', Kcorr=False):
def __init__(self, catalog_file = 'glade.hdf5', band='B'):
self.colnames = set(self.colnames).union([f'm_{b}' for b in self.supported_bands])
super().__init__(catalog_file = catalog_file, Kcorr=Kcorr, name = 'Glade')
super().__init__(catalog_file = catalog_file, name = 'Glade')
class OldStyleDES(OldStyleCatalog):
supported_bands = {'G','I','R','Z'}
supports_kcorrections = True
def __init__(self, catalog_file = 'des.hdf5', band='G', Kcorr=True):
supported_bands = {'g','i','r','z'}
def __init__(self, catalog_file = 'des.hdf5', band='G'):
self.colnames = set(self.colnames).union([f'm_{b.lower()}' for b in self.supported_bands])
print(self.colnames)
super().__init__(catalog_file = catalog_file, Kcorr=Kcorr, name = 'DES')
super().__init__(catalog_file = catalog_file, name = 'DES')
Supports Markdown
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