diff --git a/HISTORY.rst b/HISTORY.rst index fc7c8a8..ba4f73f 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -21,3 +21,10 @@ History ------------------ * using CosmoInterp module from lenstronomy + +1.1.2 (2023-09-04) +------------------ + +* double source plane likelihood +* Pantheon+ likelihood +* improved API diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index 6827ef3..4d0d4a9 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -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 @@ -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)): @@ -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) @@ -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: diff --git a/hierarc/Sampling/ParamManager/cosmo_param.py b/hierarc/Sampling/ParamManager/cosmo_param.py index f88cdb1..d3a86d4 100644 --- a/hierarc/Sampling/ParamManager/cosmo_param.py +++ b/hierarc/Sampling/ParamManager/cosmo_param.py @@ -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)) @@ -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: @@ -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'] @@ -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']) @@ -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 diff --git a/hierarc/__init__.py b/hierarc/__init__.py index 78f5e05..7307ef0 100644 --- a/hierarc/__init__.py +++ b/hierarc/__init__.py @@ -2,4 +2,4 @@ __author__ = """Simon Birrer""" __email__ = 'sibirrer@gmail.com' -__version__ = '1.1.1' +__version__ = '1.1.2' diff --git a/requirements.txt b/requirements.txt index 5a07dad..00e04f0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -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 diff --git a/setup.py b/setup.py index 1ad043b..63e55c9 100644 --- a/setup.py +++ b/setup.py @@ -36,7 +36,7 @@ def run_tests(self): setup( author="Simon Birrer", author_email='sibirrer@gmail.com', - python_requires='>=3.6', + python_requires='>=3.8', classifiers=[ 'Development Status :: 5 - Production/Stable', 'Intended Audience :: Developers', @@ -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} ) diff --git a/test/test_Likelihood/test_cosmo_likelihood.py b/test/test_Likelihood/test_cosmo_likelihood.py index acc4663..09f7e65 100644 --- a/test/test_Likelihood/test_cosmo_likelihood.py +++ b/test/test_Likelihood/test_cosmo_likelihood.py @@ -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) diff --git a/test/test_Sampling/test_ParamManager/test_cosmo_param.py b/test/test_Sampling/test_ParamManager/test_cosmo_param.py index 2b5ac5e..fa4c25d 100644 --- a/test/test_Sampling/test_ParamManager/test_cosmo_param.py +++ b/test/test_Sampling/test_ParamManager/test_cosmo_param.py @@ -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] @@ -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):