Skip to content

Commit

Permalink
updated interface and allow for new feature in CosmoInterp() class
Browse files Browse the repository at this point in the history
  • Loading branch information
@sibirrer
sibirrer committed Sep 4, 2023
1 parent ac61393 commit 9642a12
Show file tree
Hide file tree
Showing 8 changed files with 129 additions and 90 deletions.
7 changes: 7 additions & 0 deletions HISTORY.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,3 +21,10 @@ History
------------------

* using CosmoInterp module from lenstronomy

1.1.2 (2023-09-04)
------------------

* double source plane likelihood
* Pantheon+ likelihood
* improved API
16 changes: 14 additions & 2 deletions hierarc/Likelihood/cosmo_likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_likelih
:param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space
(and thus flat prior in log)
:param custom_prior: None or a definition that takes the keywords from the CosmoParam conventions and returns a
log likelihood value (e.g. prior)
log likelihood value (e.g. prior)
:param interpolate_cosmo: bool, if True, uses interpolated comoving distance in the calculation for speed-up
:param num_redshift_interp: int, number of redshift interpolation steps
:param cosmo_fixed: astropy.cosmology instance to be used and held fixed throughout the sampling
Expand Down Expand Up @@ -122,10 +122,13 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_likelih
z_max = kwargs_lens['z_source_2']
self._z_max = z_max

def likelihood(self, args):
def likelihood(self, args, kwargs_cosmo_interp=None):
"""
:param args: list of sampled parameters
:param kwargs_cosmo_interp: interpolated angular diameter distances with
'ang_diameter_distances' and 'redshifts', and optionally 'ok' and 'K' in none-flat scenarios
:type kwargs_cosmo_interp: dict
:return: log likelihood of the combined lenses
"""
for i in range(0, len(args)):
Expand All @@ -149,6 +152,8 @@ def likelihood(self, args):
# make sure that Omega_DE is not negative...
if 1.0 - om - ok <= 0:
return -np.inf
if kwargs_cosmo_interp is not None:
kwargs_cosmo = {**kwargs_cosmo, **kwargs_cosmo_interp}
cosmo = self.cosmo_instance(kwargs_cosmo)
log_l = self._likelihoodLensSample.log_likelihood(cosmo=cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin,
kwargs_source=kwargs_source)
Expand All @@ -175,6 +180,13 @@ def cosmo_instance(self, kwargs_cosmo):
:param kwargs_cosmo: cosmology parameter keyword argument list
:return: astropy.cosmology (or equivalent interpolation scheme class)
"""
if hasattr(kwargs_cosmo, 'ang_diameter_distances') and hasattr(kwargs_cosmo, 'redshifts'):
# in that case we use directly the interpolation mode to approximate angular diameter distances
cosmo = CosmoInterp(ang_dist_list=kwargs_cosmo['ang_diameter_distances'],
z_list=kwargs_cosmo['redshifts'],
Ok0=kwargs_cosmo.get('ok', 0),
K=kwargs_cosmo.get('K', None))
return cosmo
if self._cosmo_fixed is None:
cosmo = self.param.cosmo(kwargs_cosmo)
if self._interpolate_cosmo is True:
Expand Down
167 changes: 86 additions & 81 deletions hierarc/Sampling/ParamManager/cosmo_param.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ def __init__(self, cosmology, ppn_sampling=False, kwargs_fixed=None):
kwargs_fixed = {}
self._kwargs_fixed = kwargs_fixed
self._ppn_sampling = ppn_sampling
self._supported_cosmologies = ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]
self._supported_cosmologies = ["FLCDM", "FwCDM", "w0waCDM", "oLCDM", "NONE"]
if cosmology not in self._supported_cosmologies:
raise ValueError(
'cosmology %s not supported!. Please chose among %s ' % (cosmology, self._supported_cosmologies))
Expand All @@ -29,40 +29,41 @@ def param_list(self, latex_style=False):
:return: list of the free parameters being sampled in the same order as the sampling
"""
list = []
if 'h0' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$H_0$')
else:
list.append('h0')
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$\Omega_{\rm m}$')
else:
list.append('om')
if self._cosmology in ["FwCDM"]:
if 'w' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w$')
else:
list.append('w')
if self._cosmology in ["w0waCDM"]:
if 'w0' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w_0$')
else:
list.append('w0')
if 'wa' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w_{\rm a}$')
else:
list.append('wa')
if self._cosmology in ["oLCDM"]:
if 'ok' not in self._kwargs_fixed:
if self._cosmology not in ["NONE"]:
if 'h0' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$\Omega_{\rm k}$')
list.append(r'$H_0$')
else:
list.append('ok')
list.append('h0')
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$\Omega_{\rm m}$')
else:
list.append('om')
if self._cosmology in ["FwCDM"]:
if 'w' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w$')
else:
list.append('w')
if self._cosmology in ["w0waCDM"]:
if 'w0' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w_0$')
else:
list.append('w0')
if 'wa' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$w_{\rm a}$')
else:
list.append('wa')
if self._cosmology in ["oLCDM"]:
if 'ok' not in self._kwargs_fixed:
if latex_style is True:
list.append(r'$\Omega_{\rm k}$')
else:
list.append('ok')
if self._ppn_sampling is True:
if 'gamma_ppn' not in self._kwargs_fixed:
if latex_style is True:
Expand All @@ -78,40 +79,41 @@ def args2kwargs(self, args, i=0):
:return: keyword argument list with parameter names
"""
kwargs = {}
if 'h0' in self._kwargs_fixed:
kwargs['h0'] = self._kwargs_fixed['h0']
else:
kwargs['h0'] = args[i]
i += 1
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' in self._kwargs_fixed:
kwargs['om'] = self._kwargs_fixed['om']
else:
kwargs['om'] = args[i]
i += 1
if self._cosmology in ["FwCDM"]:
if 'w' in self._kwargs_fixed:
kwargs['w'] = self._kwargs_fixed['w']
if self._cosmology not in ["NONE"]:
if 'h0' in self._kwargs_fixed:
kwargs['h0'] = self._kwargs_fixed['h0']
else:
kwargs['w'] = args[i]
i += 1
if self._cosmology in ["w0waCDM"]:
if 'w0' in self._kwargs_fixed:
kwargs['w0'] = self._kwargs_fixed['w0']
else:
kwargs['w0'] = args[i]
i += 1
if 'wa' in self._kwargs_fixed:
kwargs['wa'] = self._kwargs_fixed['wa']
else:
kwargs['wa'] = args[i]
i += 1
if self._cosmology in ["oLCDM"]:
if 'ok' in self._kwargs_fixed:
kwargs['ok'] = self._kwargs_fixed['ok']
else:
kwargs['ok'] = args[i]
kwargs['h0'] = args[i]
i += 1
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' in self._kwargs_fixed:
kwargs['om'] = self._kwargs_fixed['om']
else:
kwargs['om'] = args[i]
i += 1
if self._cosmology in ["FwCDM"]:
if 'w' in self._kwargs_fixed:
kwargs['w'] = self._kwargs_fixed['w']
else:
kwargs['w'] = args[i]
i += 1
if self._cosmology in ["w0waCDM"]:
if 'w0' in self._kwargs_fixed:
kwargs['w0'] = self._kwargs_fixed['w0']
else:
kwargs['w0'] = args[i]
i += 1
if 'wa' in self._kwargs_fixed:
kwargs['wa'] = self._kwargs_fixed['wa']
else:
kwargs['wa'] = args[i]
i += 1
if self._cosmology in ["oLCDM"]:
if 'ok' in self._kwargs_fixed:
kwargs['ok'] = self._kwargs_fixed['ok']
else:
kwargs['ok'] = args[i]
i += 1
if self._ppn_sampling is True:
if 'gamma_ppn' in self._kwargs_fixed:
kwargs['gamma_ppn'] = self._kwargs_fixed['gamma_ppn']
Expand All @@ -127,22 +129,23 @@ def kwargs2args(self, kwargs):
:return: sampling argument list in specified order
"""
args = []
if 'h0' not in self._kwargs_fixed:
args.append(kwargs['h0'])
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' not in self._kwargs_fixed:
args.append(kwargs['om'])
if self._cosmology in ["FwCDM"]:
if 'w' not in self._kwargs_fixed:
args.append(kwargs['w'])
if self._cosmology in ["w0waCDM"]:
if 'w0' not in self._kwargs_fixed:
args.append(kwargs['w0'])
if 'wa' not in self._kwargs_fixed:
args.append(kwargs['wa'])
if self._cosmology in ["oLCDM"]:
if 'ok' not in self._kwargs_fixed:
args.append(kwargs['ok'])
if self._cosmology not in ["NONE"]:
if 'h0' not in self._kwargs_fixed:
args.append(kwargs['h0'])
if self._cosmology in ["FLCDM", "FwCDM", "w0waCDM", "oLCDM"]:
if 'om' not in self._kwargs_fixed:
args.append(kwargs['om'])
if self._cosmology in ["FwCDM"]:
if 'w' not in self._kwargs_fixed:
args.append(kwargs['w'])
if self._cosmology in ["w0waCDM"]:
if 'w0' not in self._kwargs_fixed:
args.append(kwargs['w0'])
if 'wa' not in self._kwargs_fixed:
args.append(kwargs['wa'])
if self._cosmology in ["oLCDM"]:
if 'ok' not in self._kwargs_fixed:
args.append(kwargs['ok'])
if self._ppn_sampling is True:
if 'gamma_ppn' not in self._kwargs_fixed:
args.append(kwargs['gamma_ppn'])
Expand All @@ -162,6 +165,8 @@ def cosmo(self, kwargs):
cosmo = w0waCDM(H0=kwargs['h0'], Om0=kwargs['om'], Ode0=1.0 - kwargs['om'], w0=kwargs['w0'], wa=kwargs['wa'])
elif self._cosmology == "oLCDM":
cosmo = LambdaCDM(H0=kwargs['h0'], Om0=kwargs['om'], Ode0=1.0 - kwargs['om'] - kwargs['ok'])
elif self._cosmology == "NONE":
cosmo = None
else:
raise ValueError("Cosmology %s is not supported" % self._cosmology)
return cosmo
2 changes: 1 addition & 1 deletion hierarc/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,4 @@

__author__ = """Simon Birrer"""
__email__ = '[email protected]'
__version__ = '1.1.1'
__version__ = '1.1.2'
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ mpmath
h5py
matplotlib
pandas
lenstronomy>=1.9.2
lenstronomy>=1.11.3
# -e git://github.com/sibirrer/lenstronomy.git@main#egg=lenstronomy
scikit-learn>=0.22
pandas
4 changes: 2 additions & 2 deletions setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ def run_tests(self):
setup(
author="Simon Birrer",
author_email='[email protected]',
python_requires='>=3.6',
python_requires='>=3.8',
classifiers=[
'Development Status :: 5 - Production/Stable',
'Intended Audience :: Developers',
Expand All @@ -61,7 +61,7 @@ def run_tests(self):
test_suite='test',
tests_require=test_requirements,
url='https://github.com/sibirrer/hierarc',
version='1.1.1',
version='1.1.2',
zip_safe=False,
cmdclass={'test': PyTest}
)
13 changes: 13 additions & 0 deletions test/test_Likelihood/test_cosmo_likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,12 +112,25 @@ def test_cosmo_instance(self):
kwargs_cosmo_wrong = {'h0': 10, 'om': .3, 'ok': 0}
cosmo_fixed = cosmoL.cosmo_instance(kwargs_cosmo_wrong)

# and here we inject pre-computed angular diameter distances as a function of redshift
redshifts = np.linspace(start=0, stop=5, num=100)
ang_diameter_distances = cosmo_astropy.angular_diameter_distance(redshifts)
Ok0 = kwargs_cosmo['ok']
K = cosmo_interp.k # compute the curvature from the interpolation class (as easily available)
kwargs_cosmo_interp = {'ang_diameter_distances': ang_diameter_distances,
'redshifts': redshifts, 'ok': Ok0, 'K': K}

cosmo_interp_input = cosmoL.cosmo_instance(kwargs_cosmo_interp)

z = 1
dd_astropy = cosmo_astropy.angular_diameter_distance(z=z).value
dd_interp = cosmo_interp.angular_diameter_distance(z=z).value
dd_interp_input = cosmo_interp_input.angular_diameter_distance(z=z).value
dd_fixed = cosmo_fixed.angular_diameter_distance(z=z).value
dd_fixed_interp = cosmo_fixed_interp.angular_diameter_distance(z=z).value

npt.assert_almost_equal(dd_astropy, dd_interp, decimal=1)
npt.assert_almost_equal(dd_astropy, dd_interp_input, decimal=1)
npt.assert_almost_equal(dd_astropy, dd_fixed, decimal=1)
npt.assert_almost_equal(dd_astropy, dd_fixed_interp, decimal=1)

Expand Down
8 changes: 5 additions & 3 deletions test/test_Sampling/test_ParamManager/test_cosmo_param.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,16 +7,17 @@
class TestCosmoParamFLCDM(object):

def setup(self):
cosmology_list = ['FLCDM', "FwCDM", "w0waCDM", "oLCDM"]
cosmology_list = ['FLCDM', "FwCDM", "w0waCDM", "oLCDM", "NONE"]
kwargs_fixed = {'h0': 70, 'om': 0.3, 'ok': 0., 'w': -1, 'wa': -0, 'w0': -0, 'gamma_ppn': 1}
self._param_list = []
self._param_list_fixed = []
for cosmology in cosmology_list:
self._param_list.append(CosmoParam(cosmology, ppn_sampling=True, kwargs_fixed=None))
self._param_list_fixed.append(CosmoParam(cosmology, ppn_sampling=True, kwargs_fixed=kwargs_fixed))
self.cosmology_list = cosmology_list

def test_param_list(self):
num_param_list = [3, 4, 5, 4] # number of parameters for the cosmological models cosmology_list
num_param_list = [3, 4, 5, 4, 1] # number of parameters for the cosmological models cosmology_list
for i, param in enumerate(self._param_list):
param_list = param.param_list(latex_style=False)
assert len(param_list) == num_param_list[i]
Expand Down Expand Up @@ -45,7 +46,8 @@ def test_cosmo(self):
kwargs_cosmo = {'h0': 70, 'om': 0.3, 'ok': 0., 'w': -1, 'wa': -0, 'w0': -0, 'gamma_ppn': 1}
for i, param in enumerate(self._param_list):
cosmo = param.cosmo(kwargs_cosmo)
assert hasattr(cosmo, 'H0')
if self.cosmology_list[i] != "NONE":
assert hasattr(cosmo, 'H0')


class TestRaise(unittest.TestCase):
Expand Down

0 comments on commit 9642a12

Please sign in to comment.