Commit 8a4c686b authored by Christos Karathanasis's avatar Christos Karathanasis
Browse files

Merge branch 'pending_review' into 'Speed_up_gwcosmo'

# Conflicts:
#   bin/gwcosmo_single_posterior
#   gwcosmo/gwcosmo.py
parents 5ac45bfc fd8c2e98
......@@ -54,6 +54,14 @@ parser = OptionParser(
help="Hard redshift cut to apply to the galaxy catalogue (default=None)"),
Option("--mth", default=None,
help="Override the apparent magnitude threshold of the catalogue, if provided (default=None)"),
Option("--schech_alpha", default=None,
help="Override the default value for slope of schechter function for given band, if provided (default=None)"),
Option("--schech_Mstar", default=None,
help="Override the default value for Mstar of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmin", default=None,
help="Override the default value for Mmin of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmax", default=None,
help="Override the default value for Mmax of schechter function for given band, if provided (default=None)"),
Option("--nside", default='32', type=int,
help="skymap nside choice for reading in galaxies from the overlap of catalogue and skymap (default=32)"),
Option("--sky_area", default='0.999', type=float,
......@@ -91,6 +99,26 @@ if opts.mth is None:
else:
mth_str = f"--mth {opts.mth}"
if opts.schech_alpha is None:
schech_alpha_str = ""
else:
schech_alpha_str = f"--schech_alpha {opts.schech_alpha}"
if opts.schech_Mstar is None:
schech_Mstar_str = ""
else:
schech_Mstar_str = f"--schech_Mstar {opts.schech_Mstar}"
if opts.schech_Mmin is None:
schech_Mmin_str = ""
else:
schech_Mmin_str = f"--schech_Mmin {opts.schech_Mmin}"
if opts.schech_Mmax is None:
schech_Mmax_str = ""
else:
schech_Mmax_str = f"--schech_Mmax {opts.schech_Mmax}"
if opts.zcut is None:
zcut_str = ""
else:
......@@ -151,6 +179,10 @@ common_args = f"\
--zmax {opts.zmax} \
{zcut_str} \
{mth_str} \
{schech_alpha_str} \
{schech_Mstar_str} \
{schech_Mmin_str} \
{schech_Mmax_str} \
--nside {opts.nside} \
--sky_area {opts.sky_area} \
--min_pixels {opts.min_pixels} \
......
#!/usr/bin/env python3
"""
This script computes posterior on H0 using a single gravitational wave event
and an electromagnetic counterpart.
......@@ -106,6 +107,14 @@ parser = OptionParser(
help="Hard redshift cut to apply to the galaxy catalogue (default=%default)"),
Option("--mth", default=None,
help="Override the apparent magnitude threshold of the catalogue, if provided (default=None)"),
Option("--schech_alpha", default=None,
help="Override the default value for slope of schechter function for given band, if provided (default=None)"),
Option("--schech_Mstar", default=None,
help="Override the default value for Mstar of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmin", default=None,
help="Override the default value for Mmin of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmax", default=None,
help="Override the default value for Mmax of schechter function for given band, if provided (default=None)"),
Option("--nside", default=32, type=int,
help="skymap nside choice for reading in galaxies from the overlap of catalogue and skymap (default=32)"),
Option("--sky_area", default=0.999, type=float,
......@@ -151,7 +160,9 @@ outputfile = str(opts.outputfile)
combine_pixels = str2bool(opts.combine_pixels)
if combine_pixels:
pixels = np.genfromtxt(outputfile+'_indices.txt',dtype=int)
low_res_nside = np.genfromtxt(outputfile+'_indices.txt',dtype=int,max_rows=1)[1]
opts.low_res_nside = low_res_nside
pixels = np.genfromtxt(outputfile+'_indices.txt',dtype=int, skip_header=1)
H0 = np.load(outputfile+'_pixel_{}.npz'.format(pixels[0]),allow_pickle=True)['arr_0'][0]
dH0 = H0[1] - H0[0]
min_H0 = np.amin(H0)
......@@ -304,8 +315,25 @@ else:
zcut = float(opts.zcut)
else:
zcut = None
if opts.schech_alpha is not None:
schech_alpha=float(opts.schech_alpha)
else:
schech_alpha = None
if opts.schech_Mstar is not None:
schech_Mstar=float(opts.schech_Mstar)
else:
schech_Mstar = None
if opts.schech_Mmin is not None:
schech_Mmin=float(opts.schech_Mmin)
else:
schech_Mmin = None
if opts.schech_Mmax is not None:
schech_Mmax=float(opts.schech_Mmax)
else:
schech_Mmax = None
sp = SchechterParams(band)
sp = SchechterParams(band, schech_alpha, schech_Mstar, schech_Mmin, schech_Mmax)
print("Schechter function with parameters: alpha={}, Mstar={}, Mmin={}, Mmax={}, ".format(sp.alpha, sp.Mstar, sp.Mmin, sp.Mmax))
p_M = SchechterMagFunction(Mstar_obs=sp.Mstar,alpha=sp.alpha)
if galaxy_weighting:
ps_M = gwcosmo.gwcosmo.LuminosityWeighting()
......@@ -314,20 +342,35 @@ 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, 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, numerical=numerical )
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, numerical=numerical )
if opts.method == 'pixel':
min_pixels = int(opts.min_pixels)
pixelated_samples = gwcosmo.likelihood.posterior_samples.make_pixel_px_function(samples, skymap, npixels=min_pixels, thresh=sky_area)
nside_low_res = pixelated_samples.nside
if nside_low_res > nside:
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:
np.savetxt(outputfile+'_indices.txt',pixelated_samples.indices,fmt='%i')
# 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()
nside_low_res = pixelated_samples.nside
# Select the coarse pixel from the skymap that we will use
catalog = catalog.select_pixel(nside_low_res, pixel_index, nested=skymap.nested)
......@@ -335,7 +378,7 @@ else:
print("Setting up p(x|z,H0)")
px_zH0 = pixelated_samples.make_los_px_function(samp_ind,H0,hyper_params_dict,name=mass_distribution,reweight_samples=reweight_samples)
me = gwcosmo.gwcosmo.SinglePixelGalaxyCatalogLikelihood(pixel_index, catalog, skymap, band, cosmo,\
me = gwcosmo.gwcosmo.SinglePixelGalaxyCatalogLikelihood(pixel_index, catalog, skymap, band, sp, cosmo,\
px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, \
ps_M, outputfile, Kcorr=Kcorr, mth=mth, zcut=zcut, \
zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, \
......
......@@ -190,6 +190,8 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
----------
skymap : gwcosmo.likelihood.skymap.skymap object
provides p(x|Omega) and skymap properties
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : gwcosmo.utilities.standard_cosmology.fast_cosmology object
Cosmological model
Kcorr : bool, optional
......@@ -199,8 +201,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.
......@@ -209,7 +209,7 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
"""
def __init__(self, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, zmax=10.):
def __init__(self, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, zmax=10.):
"""
Parameters
----------
......@@ -246,7 +246,6 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
self.Kcorr = Kcorr
self.band = observation_band
sp = SchechterParams(self.band)
self.Mmin_obs = sp.Mmin
self.Mmax_obs = sp.Mmax
......@@ -271,7 +270,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:
......@@ -742,6 +740,8 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The galaxy catalogue
skymap : gwcosmo.likelihood.skymap.skymap object
provides p(x|Omega) and skymap properties
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : gwcosmo.utilities.standard_cosmology.fast_cosmology object
Cosmological model
Kcorr : bool, optional
......@@ -761,7 +761,9 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
"""
def __init__(self, pixel_index, galaxy_catalog, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, outputfile, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, nside=32, nside_low_res = None, numerical=True):
def __init__(self, pixel_index, galaxy_catalog, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, outputfile, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, nside=32, nside_low_res = None, numerical=True):
"""
Parameters
----------
......@@ -804,7 +806,7 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The high-resolution value of nside to subdivide the current pixel into
"""
super().__init__(skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
super().__init__(skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
self.nside_low_res = nside_low_res
self.zcut = zcut
self.complete_catalog = complete_catalog
......@@ -812,11 +814,26 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
self.path = outputfile+'_'+str(pixel_index)+'_checkpoint.p'
self.numerical = numerical
# 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
......@@ -994,7 +1011,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))
......@@ -1037,7 +1057,8 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
catalogue method.
"""
def __init__(self, galaxy_catalog, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, sky_thresh = 0.999, nside=32, numerical=True):
def __init__(self, galaxy_catalog, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, sky_thresh = 0.999, nside=32, numerical=True):
"""
Parameters
----------
......@@ -1047,6 +1068,8 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The GW skymap
observation_band : str
Observation band (eg. 'B', 'K', 'u', 'g')
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : object
Fast cosmology
px_zH0 : object
......@@ -1079,7 +1102,7 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
TODO: UPDATE this for new catalog classes
"""
super().__init__(skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
super().__init__(skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
self.mth = mth
self.zcut = zcut
......@@ -1089,7 +1112,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:
......@@ -1188,7 +1225,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')
......@@ -16,7 +16,7 @@ class SchechterParams():
(note paper quotes M*=-23.55 but means M*=-23.55 + 5log10(h))
"""
def __init__(self, band):
def __init__(self, band, schech_alpha=None, schech_Mstar=None, schech_Mmin=None, schech_Mmax=None):
"""
Parameters
----------
......@@ -28,23 +28,32 @@ class SchechterParams():
self.Mmin = None
self.Mmax = None
self.alpha, self.Mstar, self.Mmin, self.Mmax = self.values(band)
self.alpha, self.Mstar, self.Mmin, self.Mmax = self.default_values(band) #initialize to default values
def values(self, band):
if schech_alpha is not None: #change to input value if provided
self.alpha=schech_alpha
if schech_Mstar is not None:
self.Mstar=schech_Mstar
if schech_Mmin is not None:
self.Mmin=schech_Mmin
if schech_Mmax is not None:
self.Mmax=schech_Mmax
def default_values(self, band):
if band == 'B':
return -1.21, -19.70, -22.96, -12.96
elif band == 'K':
return -1.02, -23.55, -27.0, -12.96
elif band == 'u':
return -0.92, -17.93, -21.93, -15.54 #TODO check Mmin and Mmax
elif band == 'u': #These values are actually u', g', r', i', zi, a.k.a. redshifted to the median redshift (0.1) of SDSS.
return -0.92, -17.93, -21.93, -15.54
elif band == 'g':