From 6f59c3bf391282dad66c188500ef0dd918ab7775 Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Sun, 28 Apr 2024 18:21:02 -0400 Subject: [PATCH 1/4] new line of sight distribution handling for multiple distributions and separate list of dictionary, including generalized extreme value distribution --- hierarc/Diagnostics/goodness_of_fit.py | 18 ++- .../KDELikelihood/kde_likelihood.py | 4 +- .../LensLikelihood/double_source_plane.py | 3 +- hierarc/Likelihood/cosmo_likelihood.py | 20 ++-- hierarc/Likelihood/hierarchy_likelihood.py | 79 ++++++------ hierarc/Likelihood/lens_sample_likelihood.py | 11 +- hierarc/Likelihood/los_distributions.py | 80 +++++++++++++ hierarc/Sampling/ParamManager/lens_param.py | 34 ------ hierarc/Sampling/ParamManager/los_param.py | 113 ++++++++++++++++++ .../Sampling/ParamManager/param_manager.py | 37 ++++-- test/test_Diagnostics/test_goodness_of_fit.py | 9 +- .../test_hierarchy_likelihood.py | 35 +++--- .../test_lens_sample_likelihood.py | 2 +- test/test_Likelihood/test_los_distribution.py | 103 ++++++++++++++++ .../test_ParamManager/test_lens_param.py | 16 +-- .../test_ParamManager/test_los_param.py | 75 ++++++++++++ .../test_ParamManager/test_param_manager.py | 31 +++-- 17 files changed, 524 insertions(+), 146 deletions(-) create mode 100644 hierarc/Likelihood/los_distributions.py create mode 100644 hierarc/Sampling/ParamManager/los_param.py create mode 100644 test/test_Likelihood/test_los_distribution.py create mode 100644 test/test_Sampling/test_ParamManager/test_los_param.py diff --git a/hierarc/Diagnostics/goodness_of_fit.py b/hierarc/Diagnostics/goodness_of_fit.py index cdfd82cf..8f14cb14 100644 --- a/hierarc/Diagnostics/goodness_of_fit.py +++ b/hierarc/Diagnostics/goodness_of_fit.py @@ -28,6 +28,7 @@ def plot_ddt_fit( cosmo, kwargs_lens, kwargs_kin, + kwargs_los, color_measurement=None, color_prediction=None, redshift_trend=False, @@ -64,7 +65,7 @@ def plot_ddt_fit( ddt_model_sigma, dd_model_mean, dd_model_sigma, - ) = likelihood.ddt_dd_model_prediction(cosmo, kwargs_lens=kwargs_lens) + ) = likelihood.ddt_dd_model_prediction(cosmo, kwargs_lens=kwargs_lens, kwargs_los=kwargs_los) ddt_name_list.append(name) ddt_model_mean_list.append(ddt_model_mean) @@ -134,13 +135,14 @@ def plot_ddt_fit( ax.legend() return f, ax - def kin_fit(self, cosmo, kwargs_lens, kwargs_kin): + def kin_fit(self, cosmo, kwargs_lens, kwargs_kin, kwargs_los): """Plots the prediction and the uncorrelated error bars on the individual lenses currently works for likelihood classes 'TDKinGaussian', 'KinGaussian'. :param cosmo: astropy.cosmology instance :param kwargs_lens: lens model parameter keyword arguments :param kwargs_kin: kinematics model keyword arguments + :param kwargs_los: line of sight list of dictionaries :return: list of name, measurement, measurement errors, model prediction, model prediction error """ @@ -160,7 +162,7 @@ def kin_fit(self, cosmo, kwargs_lens, kwargs_kin): sigma_v_predict_mean, cov_error_predict, ) = likelihood.sigma_v_measured_vs_predict( - cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) if sigma_v_measurement is not None: @@ -188,6 +190,7 @@ def plot_kin_fit( cosmo, kwargs_lens, kwargs_kin, + kwargs_los, color_measurement=None, color_prediction=None, ): @@ -197,11 +200,12 @@ def plot_kin_fit( :param cosmo: astropy.cosmology instance :param kwargs_lens: lens model parameter keyword arguments :param kwargs_kin: kinematics model keyword arguments + :param kwargs_los: line of sight list of dictionaries :param color_measurement: color of measurement :param color_prediction: color of model prediction :return: fig, axes of matplotlib instance """ - logL = self._sample_likelihood.log_likelihood(cosmo, kwargs_lens, kwargs_kin) + logL = self._sample_likelihood.log_likelihood(cosmo, kwargs_lens, kwargs_kin, kwargs_los=kwargs_los) print(logL, "log likelihood") ( sigma_v_name_list, @@ -209,7 +213,7 @@ def plot_kin_fit( sigma_v_measurement_error_list, sigma_v_model_list, sigma_v_model_error_list, - ) = self.kin_fit(cosmo, kwargs_lens, kwargs_kin) + ) = self.kin_fit(cosmo, kwargs_lens, kwargs_kin, kwargs_los) f, ax = plt.subplots(1, 1, figsize=(int(len(sigma_v_name_list) / 2), 4)) ax.errorbar( @@ -256,6 +260,7 @@ def plot_ifu_fit( cosmo, kwargs_lens, kwargs_kin, + kwargs_los, lens_index, bin_edges, show_legend=True, @@ -268,6 +273,7 @@ def plot_ifu_fit( :param cosmo: astropy.cosmology instance :param kwargs_lens: lens model parameter keyword arguments :param kwargs_kin: kinematics model keyword arguments + :param kwargs_los: line of sight list of dictionaries :param lens_index: int, index in kwargs_lens to be plotted (needs to be of type 'IFUKinCov') :param bin_edges: radial bin edges in arc seconds. If number, then uniform @@ -293,7 +299,7 @@ def plot_ifu_fit( sigma_v_predict_mean, cov_error_predict, ) = likelihood.sigma_v_measured_vs_predict( - cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) if len(np.atleast_1d(bin_edges)) < 2: diff --git a/hierarc/Likelihood/KDELikelihood/kde_likelihood.py b/hierarc/Likelihood/KDELikelihood/kde_likelihood.py index 809bd63e..7357c1a9 100644 --- a/hierarc/Likelihood/KDELikelihood/kde_likelihood.py +++ b/hierarc/Likelihood/KDELikelihood/kde_likelihood.py @@ -64,9 +64,9 @@ def init_loglikelihood(self): ) self.loglikelihood = self.kdelikelihood() else: - raise ValueError( + raise NameError( "likelihood_type %s not supported! Supported are %s." - % (likelihood_type, LIKELIHOOD_TYPES) + % (self.loglikelihood_type, LIKELIHOOD_TYPES) ) def kdelikelihood(self): diff --git a/hierarc/Likelihood/LensLikelihood/double_source_plane.py b/hierarc/Likelihood/LensLikelihood/double_source_plane.py index d45dd009..b742183a 100644 --- a/hierarc/Likelihood/LensLikelihood/double_source_plane.py +++ b/hierarc/Likelihood/LensLikelihood/double_source_plane.py @@ -30,7 +30,7 @@ def __init__( self._normalized = normalized def lens_log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None + self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, kwargs_los=None ): """Log likelihood of the data given a model. @@ -38,6 +38,7 @@ def lens_log_likelihood( :param kwargs_lens: keyword arguments of lens :param kwargs_kin: keyword arguments of kinematics :param kwargs_source: keyword argument of source + :param kwargs_los: keyword argument list of line of sight :return: log likelihood of data given model """ beta_model = self._beta_model(cosmo) diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index 5be38d85..46c9c2ae 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -27,8 +27,8 @@ def __init__( gamma_in_distribution="NONE", log_m2l_sampling=False, log_m2l_distribution="NONE", - kappa_ext_sampling=False, - kappa_ext_distribution="NONE", + los_sampling=False, + los_distributions=None, alpha_lambda_sampling=False, beta_lambda_sampling=False, alpha_gamma_in_sampling=False, @@ -73,8 +73,10 @@ def __init__( according to the lens posterior kwargs 'lambda_scaling_property' :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly according to the lens posterior kwargs 'lambda_scaling_property_beta' - :param kappa_ext_sampling: bool, if True samples a global external convergence parameter - :param kappa_ext_distribution: string, distribution function of the kappa_ext parameter + :param los_sampling: if sampling of the parameters should be done + :type los_sampling: bool + :param los_distributions: what distribution to be sampled + :type los_distributions: list of str :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens kinematic prediction :param anisotropy_model: string, specifies the stellar anisotropy model @@ -106,7 +108,8 @@ def __init__( if sigma_v_systematics is True: normalized = True self._likelihoodLensSample = LensSampleLikelihood( - kwargs_likelihood_list, normalized=normalized + kwargs_likelihood_list, normalized=normalized, + los_distributions=los_distributions ) self.param = ParamManager( cosmology, @@ -127,8 +130,8 @@ def __init__( sne_distribution=sne_distribution, z_apparent_m_anchor=z_apparent_m_anchor, sigma_v_systematics=sigma_v_systematics, - kappa_ext_sampling=kappa_ext_sampling, - kappa_ext_distribution=kappa_ext_distribution, + los_sampling=los_sampling, + los_distributions=los_distributions, anisotropy_sampling=anisotropy_sampling, anisotropy_model=anisotropy_model, anisotropy_distribution=anisotropy_distribution, @@ -188,7 +191,7 @@ def likelihood(self, args, kwargs_cosmo_interp=None): if args[i] < self._lower_limit[i] or args[i] > self._upper_limit[i]: return -np.inf - kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source = self.param.args2kwargs( + kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los = self.param.args2kwargs( args ) if self._cosmology == "oLCDM": @@ -215,6 +218,7 @@ def likelihood(self, args, kwargs_cosmo_interp=None): kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_source=kwargs_source, + kwargs_los=kwargs_los, ) if self._sne_evaluate is True: diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py index af2127fd..3294023b 100644 --- a/hierarc/Likelihood/hierarchy_likelihood.py +++ b/hierarc/Likelihood/hierarchy_likelihood.py @@ -1,13 +1,13 @@ from hierarc.Likelihood.transformed_cosmography import TransformedCosmography from hierarc.Likelihood.LensLikelihood.base_lens_likelihood import LensLikelihoodBase from hierarc.Likelihood.parameter_scaling import ParameterScalingIFU -from hierarc.Util.distribution_util import PDFSampling +from hierarc.Likelihood.los_distributions import LOSDistribution import numpy as np import copy class LensLikelihood(TransformedCosmography, LensLikelihoodBase, ParameterScalingIFU): - """Master class containing the likelihood definitions of different analysis for s + """Master class containing the likelihood definitions of different analysis for a single lens.""" def __init__( @@ -24,7 +24,7 @@ def __init__( log_m2l_array=None, param_scaling_grid_list=None, num_distribution_draws=50, - kappa_ext_bias=False, + global_los_distribution=False, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, @@ -34,6 +34,7 @@ def __init__( kwargs_lens_properties=None, gamma_in_prior_mean=None, gamma_in_prior_std=None, + los_distributions=None, **kwargs_likelihood ): """ @@ -52,8 +53,9 @@ def __init__( anisotropy, gamma_in, and log_m2l. In that case, gamma_in_array and log_m2l_array need to be provided. :param num_distribution_draws: int, number of distribution draws from the likelihood that are being averaged over - :param kappa_ext_bias: bool, if True incorporates the global external selection function into the likelihood. - If False, the likelihood needs to incorporate the individual selection function with sufficient accuracy. + :param global_los_distribution: if integer, will draw from the global kappa distribution specified in that + integer. If False, will instead draw from the distribution specified in kappa_pdf. + :type global_los_distribution: bool or integer :param kappa_pdf: array of probability density function of the external convergence distribution binned according to kappa_bin_edges :param kappa_bin_edges: array of length (len(kappa_pdf)+1), bin edges of the kappa PDF @@ -69,7 +71,9 @@ def __init__( :param gamma_in_prior_mean: prior mean for inner power-law slope of the NFW profile, if available :param gamma_in_prior_std: standard deviation of the Gaussian prior for gamma_in :param kwargs_likelihood: keyword arguments specifying the likelihood function, - see individual classes for their use + see individual classes for their use + :param los_distributions: list of all line of sight distributions parameterized + :type los_distributions: list of str or None """ TransformedCosmography.__init__(self, z_lens=z_lens, z_source=z_source) if ani_scaling_array_list is None and ani_scaling_array is not None: @@ -122,15 +126,11 @@ def __init__( **kwargs_likelihood ) self._num_distribution_draws = int(num_distribution_draws) - self._kappa_ext_bias = kappa_ext_bias + + self._los = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=global_los_distribution, + los_distributions=los_distributions) self._mst_ifu = mst_ifu - if kappa_pdf is not None and kappa_bin_edges is not None: - self._kappa_dist = PDFSampling( - bin_edges=kappa_bin_edges, pdf_array=kappa_pdf - ) - self._draw_kappa = True - else: - self._draw_kappa = False self._lambda_scaling_property = lambda_scaling_property self._lambda_scaling_property_beta = lambda_scaling_property_beta self._gamma_in_array = gamma_in_array @@ -140,7 +140,7 @@ def __init__( self._gamma_in_prior_std = gamma_in_prior_std def lens_log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None + self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, kwargs_los=None, ): """Log likelihood of the data of a lens given a model (defined with hyper- parameters) and cosmology. @@ -149,6 +149,7 @@ def lens_log_likelihood( :param kwargs_lens: keywords of the hyper parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model hyper parameters :param kwargs_source: keyword argument of the source model (such as SNe) + :param kwargs_los: list of keyword arguments of global line of sight distributions :return: log likelihood of the data given the model """ @@ -168,6 +169,7 @@ def lens_log_likelihood( kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_source=kwargs_source, + kwargs_los=kwargs_los, cosmo=cosmo, ) return a @@ -180,6 +182,7 @@ def hyper_param_likelihood( kwargs_lens=None, kwargs_kin=None, kwargs_source=None, + kwargs_los=None, cosmo=None, ): """Log likelihood of the data of a lens given a model (defined with hyper- @@ -191,6 +194,7 @@ def hyper_param_likelihood( :param kwargs_lens: keywords of the hyper parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model hyper parameters :param kwargs_source: keyword argument of the source model (such as SNe) + :param kwargs_los: list of keyword arguments of global line of sight distributions :param cosmo: astropy.cosmology instance :return: log likelihood given the single lens analysis for the given hyper parameter @@ -202,7 +206,7 @@ def hyper_param_likelihood( sigma_v_sys_error = kwargs_kin_copy.pop("sigma_v_sys_error", None) if self.check_dist( - kwargs_lens, kwargs_kin, kwargs_source + kwargs_lens, kwargs_kin, kwargs_source, kwargs_los ): # sharp distributions return self.log_likelihood_single( ddt, @@ -211,6 +215,7 @@ def hyper_param_likelihood( kwargs_lens, kwargs_kin_copy, kwargs_source, + kwargs_los, sigma_v_sys_error=sigma_v_sys_error, ) else: @@ -223,6 +228,7 @@ def hyper_param_likelihood( kwargs_lens, kwargs_kin_copy, kwargs_source, + kwargs_los, sigma_v_sys_error=sigma_v_sys_error, ) exp_logl = np.exp(logl) @@ -240,6 +246,7 @@ def log_likelihood_single( kwargs_lens, kwargs_kin, kwargs_source, + kwargs_los=None, sigma_v_sys_error=None, ): """Log likelihood of the data of a lens given a specific model (as a draw from @@ -251,12 +258,14 @@ def log_likelihood_single( :param kwargs_lens: keywords of the hyper parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model hyper parameters :param kwargs_source: keyword arguments of source brightness + :param kwargs_los: line of sight list of dictionaries :param sigma_v_sys_error: unaccounted uncertainty in the velocity dispersion measurement :return: log likelihood given the single lens analysis for a single (random) realization of the hyper parameter distribution """ - lambda_mst, kappa_ext, gamma_ppn = self.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = self.draw_lens(**kwargs_lens) + kappa_ext = self._los.draw_los(kwargs_los) # draw intrinsic source magnitude mag_source = self.draw_source(lum_dist=delta_lum_dist, **kwargs_source) ddt_, dd_, mag_source_ = self.displace_prediction( @@ -420,37 +429,35 @@ def luminosity_distance_modulus(self, cosmo, z_apparent_m_anchor): delta_lum_dist = lum_dists - lum_dist_anchor return delta_lum_dist - def check_dist(self, kwargs_lens, kwargs_kin, kwargs_source): + def check_dist(self, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los): """Checks if the provided keyword arguments describe a distribution function of hyper parameters or are single values. :param kwargs_lens: lens model hyper-parameter keywords :param kwargs_kin: kinematic model hyper-parameter keywords :param kwargs_source: source brightness hyper-parameter keywords + :param kwargs_los: list of dictionaries for line of sight hyper-parameters :return: bool, True if delta function, else False """ lambda_mst_sigma = kwargs_lens.get("lambda_mst_sigma", 0) # scatter in MST - kappa_ext_sigma = kwargs_lens.get("kappa_ext_sigma", 0) + draw_kappa_bool = self._los.draw_bool(kwargs_los) a_ani_sigma = kwargs_kin.get("a_ani_sigma", 0) beta_inf_sigma = kwargs_kin.get("beta_inf_sigma", 0) sne_sigma = kwargs_source.get("sigma_sne", 0) if ( a_ani_sigma == 0 and lambda_mst_sigma == 0 - and kappa_ext_sigma == 0 and beta_inf_sigma == 0 and sne_sigma == 0 + and not draw_kappa_bool ): - if self._draw_kappa is False: - return True + return True return False def draw_lens( self, lambda_mst=1, lambda_mst_sigma=0, - kappa_ext=0, - kappa_ext_sigma=0, gamma_ppn=1, lambda_ifu=1, lambda_ifu_sigma=0, @@ -468,8 +475,6 @@ def draw_lens( :param lambda_mst: MST transform :param lambda_mst_sigma: spread in the distribution - :param kappa_ext: external convergence mean in distribution - :param kappa_ext_sigma: spread in the distribution :param gamma_ppn: Post-Newtonian parameter :param lambda_ifu: secondary lambda_mst parameter for subset of lenses specified for @@ -503,13 +508,7 @@ def draw_lens( + beta_lambda * self._lambda_scaling_property_beta ) lambda_mst_draw = np.random.normal(lambda_lens, lambda_mst_sigma) - if self._draw_kappa is True: - kappa_ext_draw = self._kappa_dist.draw_one - elif self._kappa_ext_bias is True: - kappa_ext_draw = np.random.normal(kappa_ext, kappa_ext_sigma) - else: - kappa_ext_draw = 0 - return lambda_mst_draw, kappa_ext_draw, gamma_ppn + return lambda_mst_draw, gamma_ppn @staticmethod def draw_source(mu_sne=1, sigma_sne=0, lum_dist=0, **kwargs): @@ -530,13 +529,14 @@ def draw_source(mu_sne=1, sigma_sne=0, lum_dist=0, **kwargs): # return linear amplitude with base log 10 return mag_source - def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None): + def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_los=None): """Mean and error covariance of velocity dispersion measurement mean and error covariance of velocity dispersion predictions. :param cosmo: astropy.cosmology instance :param kwargs_lens: keywords of the hyper parameters of the lens model - :param kwargs_kin: keyword arguments of the kinematic model hyper parameters + :param kwargs_kin: keyword arguments of the kinematic model hyper-parameters + :param kwargs_los: line of sight parapers :return: sigma_v_measurement, cov_error_measurement, sigma_v_predict_mean, cov_error_predict """ @@ -557,7 +557,8 @@ def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None): sigma_v_predict_mean = np.zeros_like(sigma_v_measurement) cov_error_predict = np.zeros_like(cov_error_measurement) for i in range(self._num_distribution_draws): - lambda_mst, kappa_ext, gamma_ppn = self.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = self.draw_lens(**kwargs_lens) + kappa_ext = self._los.draw_los(kwargs_los) ddt_, dd_, _ = self.displace_prediction( ddt, dd, gamma_ppn=gamma_ppn, lambda_mst=lambda_mst, kappa_ext=kappa_ext ) @@ -586,12 +587,13 @@ def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None): cov_error_predict, ) - def ddt_dd_model_prediction(self, cosmo, kwargs_lens=None): + def ddt_dd_model_prediction(self, cosmo, kwargs_lens=None, kwargs_los=None): """Predicts the model uncertainty corrected ddt prediction of the applied model (e.g. power-law) :param cosmo: astropy.cosmology instance :param kwargs_lens: keywords of the hyper parameters of the lens model + :param kwargs_los: line of slight list of dictionaries :return: ddt_model mean, ddt_model sigma, dd_model mean, dd_model sigma """ if kwargs_lens is None: @@ -600,7 +602,8 @@ def ddt_dd_model_prediction(self, cosmo, kwargs_lens=None): ddt_draws = [] dd_draws = [] for i in range(self._num_distribution_draws): - lambda_mst, kappa_ext, gamma_ppn = self.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = self.draw_lens(**kwargs_lens) + kappa_ext = self._los.draw_los(kwargs_los) ddt_, dd_, _ = self.displace_prediction( ddt, dd, gamma_ppn=gamma_ppn, lambda_mst=lambda_mst, kappa_ext=kappa_ext ) diff --git a/hierarc/Likelihood/lens_sample_likelihood.py b/hierarc/Likelihood/lens_sample_likelihood.py index 7385aa3c..3b0bc539 100644 --- a/hierarc/Likelihood/lens_sample_likelihood.py +++ b/hierarc/Likelihood/lens_sample_likelihood.py @@ -9,12 +9,14 @@ class LensSampleLikelihood(object): diameter posteriors Currently this class does not include possible covariances between the lens samples.""" - def __init__(self, kwargs_lens_list, normalized=False): + def __init__(self, kwargs_lens_list, normalized=False, los_distributions=None): """ :param kwargs_lens_list: keyword argument list specifying the arguments of the LensLikelihood class :param normalized: bool, if True, returns the normalized likelihood, if False, separates the constant prefactor (in case of a Gaussian 1/(sigma sqrt(2 pi)) ) to compute the reduced chi2 statistics + :param los_distributions: list of all line of sight distributions parameterized + :type los_distributions: list of str or None """ self._lens_list = [] for kwargs_lens in kwargs_lens_list: @@ -26,11 +28,12 @@ def __init__(self, kwargs_lens_list, normalized=False): ) else: self._lens_list.append( - LensLikelihood(normalized=normalized, **kwargs_lens) + LensLikelihood(normalized=normalized, los_distributions=los_distributions, **kwargs_lens) ) def log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None + self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, + kwargs_los=None, ): """ @@ -38,6 +41,7 @@ def log_likelihood( :param kwargs_lens: keywords of the parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model :param kwargs_source: keyword argument of the source model (such as SNe) + :param kwargs_los: line of sight keyword argument list :return: log likelihood of the combined lenses """ log_likelihood = 0 @@ -47,6 +51,7 @@ def log_likelihood( kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_source=kwargs_source, + kwargs_los=kwargs_los, ) return log_likelihood diff --git a/hierarc/Likelihood/los_distributions.py b/hierarc/Likelihood/los_distributions.py new file mode 100644 index 00000000..9a7bd8a5 --- /dev/null +++ b/hierarc/Likelihood/los_distributions.py @@ -0,0 +1,80 @@ +import numpy as np +from hierarc.Util.distribution_util import PDFSampling + + +class LOSDistribution(object): + """ + line of sight distribution drawing + """ + def __init__(self, kappa_pdf=None, kappa_bin_edges=None, global_los_distribution=False, + los_distributions=None): + """ + + :param global_los_distribution: if integer, will draw from the global kappa distribution specified in that + integer. If False, will instead draw from the distribution specified in kappa_pdf. + :type global_los_distribution: bool or int + :param kappa_pdf: array of probability density function of the external convergence distribution + binned according to kappa_bin_edges + :param kappa_bin_edges: array of length (len(kappa_pdf)+1), bin edges of the kappa PDF + :param los_distributions: list of all line of sight distributions parameterized + :type los_distributions: list of str or None + """ + + self._global_los_distribution = global_los_distribution + if isinstance(self._global_los_distribution, int) and self._global_los_distribution is not False: + self._draw_kappa_global = True + self._los_distribution = los_distributions[global_los_distribution] + else: + self._draw_kappa_global = False + if kappa_pdf is not None and kappa_bin_edges is not None and not self._draw_kappa_global: + print("test kappa pdf sampling") + self._kappa_dist = PDFSampling( + bin_edges=kappa_bin_edges, pdf_array=kappa_pdf + ) + self._draw_kappa_individual = True + else: + self._draw_kappa_individual = False + + def draw_los(self, kwargs_los, size=1): + """ + Draw from the distribution of line of sight convergence + + :param kwargs_los: line of sight parameters + :type kwargs_los: list of dict + :param size: how many samples to be drawn + :type size: int>0 + :return: external convergence draw + """ + if self._draw_kappa_individual is True: + kappa_ext_draw = self._kappa_dist.draw(n=size) + elif self._draw_kappa_global: + kwargs_los_i = kwargs_los[self._global_los_distribution] + if self._los_distribution == "GAUSSIAN": + los_mean = kwargs_los_i["mean"] + los_sigma = kwargs_los_i["sigma"] + kappa_ext_draw = np.random.normal(los_mean, los_sigma, size=size) + elif self._los_distribution == "GEV": + mean = kwargs_los_i["mean"] + sigma = kwargs_los_i["sigma"] + xi = kwargs_los_i["xi"] + from scipy.stats import genextreme + kappa_ext_draw = genextreme.rvs(c=xi, loc=mean, scale=sigma, size=size) + else: + raise ValueError("line of sight distribution %s not valid." % self._los_distribution) + else: + kappa_ext_draw = 0 + return kappa_ext_draw + + def draw_bool(self, kwargs_los): + """ + whether single-valued or extended distribution (need to draw from) + + :param kwargs_los: list of keyword arguments for line of sight distributions + :return: boolean, True with samples need to be drawn, else False + """ + if self._draw_kappa_individual is True: + return True + elif self._draw_kappa_global is True: + if kwargs_los[self._global_los_distribution]["sigma"] != 0: + return True + return False diff --git a/hierarc/Sampling/ParamManager/lens_param.py b/hierarc/Sampling/ParamManager/lens_param.py index 7cda2414..9b2e03af 100644 --- a/hierarc/Sampling/ParamManager/lens_param.py +++ b/hierarc/Sampling/ParamManager/lens_param.py @@ -14,8 +14,6 @@ def __init__( gamma_in_distribution="NONE", log_m2l_sampling=False, log_m2l_distribution="NONE", - kappa_ext_sampling=False, - kappa_ext_distribution="NONE", alpha_lambda_sampling=False, beta_lambda_sampling=False, alpha_gamma_in_sampling=False, @@ -58,8 +56,6 @@ def __init__( self._gamma_in_distribution = gamma_in_distribution self._log_m2l_sampling = log_m2l_sampling self._log_m2l_distribution = log_m2l_distribution - self._kappa_ext_sampling = kappa_ext_sampling - self._kappa_ext_distribution = kappa_ext_distribution self._alpha_lambda_sampling = alpha_lambda_sampling self._beta_lambda_sampling = beta_lambda_sampling self._alpha_gamma_in_sampling = alpha_gamma_in_sampling @@ -137,18 +133,6 @@ def param_list(self, latex_style=False): list.append(r"$\sigma(\Upsilon_{\rm stars})$") else: list.append("log_m2l_sigma") - if self._kappa_ext_sampling is True: - if "kappa_ext" not in self._kwargs_fixed: - if latex_style is True: - list.append(r"$\overline{\kappa}_{\rm ext}$") - else: - list.append("kappa_ext") - if self._kappa_ext_distribution == "GAUSSIAN": - if "kappa_ext_sigma" not in self._kwargs_fixed: - if latex_style is True: - list.append(r"$\sigma(\kappa_{\rm ext})$") - else: - list.append("kappa_ext_sigma") if self._alpha_lambda_sampling is True: if "alpha_lambda" not in self._kwargs_fixed: if latex_style is True: @@ -242,18 +226,6 @@ def args2kwargs(self, args, i=0): else: kwargs["log_m2l_sigma"] = args[i] i += 1 - if self._kappa_ext_sampling is True: - if "kappa_ext" in self._kwargs_fixed: - kwargs["kappa_ext"] = self._kwargs_fixed["kappa_ext"] - else: - kwargs["kappa_ext"] = args[i] - i += 1 - if self._kappa_ext_distribution == "GAUSSIAN": - if "kappa_ext_sigma" in self._kwargs_fixed: - kwargs["kappa_ext_sigma"] = self._kwargs_fixed["kappa_ext_sigma"] - else: - kwargs["kappa_ext_sigma"] = args[i] - i += 1 if self._alpha_lambda_sampling is True: if "alpha_lambda" in self._kwargs_fixed: kwargs["alpha_lambda"] = self._kwargs_fixed["alpha_lambda"] @@ -324,12 +296,6 @@ def kwargs2args(self, kwargs): args.append(np.log10(kwargs["log_m2l_sigma"])) else: args.append(kwargs["log_m2l_sigma"]) - if self._kappa_ext_sampling is True: - if "kappa_ext" not in self._kwargs_fixed: - args.append(kwargs["kappa_ext"]) - if self._kappa_ext_distribution == "GAUSSIAN": - if "kappa_ext_sigma" not in self._kwargs_fixed: - args.append(kwargs["kappa_ext_sigma"]) if self._alpha_lambda_sampling is True: if "alpha_lambda" not in self._kwargs_fixed: args.append(kwargs["alpha_lambda"]) diff --git a/hierarc/Sampling/ParamManager/los_param.py b/hierarc/Sampling/ParamManager/los_param.py new file mode 100644 index 00000000..3ed09df3 --- /dev/null +++ b/hierarc/Sampling/ParamManager/los_param.py @@ -0,0 +1,113 @@ +_LOS_DISTRIBUTIONS = ["GEV", "GAUSSIAN", "NONE"] + + +class LOSParam(object): + """Manager for the source property parameters (currently particularly source + magnitudes for SNe)""" + + def __init__( + self, + los_sampling=False, + los_distributions=None, + kwargs_fixed=None, + + ): + """ + + :param los_sampling: if sampling of the parameters should be done + :type los_sampling: bool + :param los_distributions: what distribution to be sampled + :type los_distributions: list of str + :param kwargs_fixed: fixed arguments in sampling + :type kwargs_fixed: list of dictionaries or None + """ + self._los_sampling = los_sampling + if los_distributions is None: + los_distributions = [] + for los_distribution in los_distributions: + if los_distribution not in _LOS_DISTRIBUTIONS: + raise ValueError( + "SNE distribution %s not supported. Please chose among %s." + % (los_distribution, _LOS_DISTRIBUTIONS) + ) + self._los_distributions = los_distributions + if kwargs_fixed is None: + kwargs_fixed = [{} for _ in range(len(los_distributions))] + self._kwargs_fixed = kwargs_fixed + + def param_list(self, latex_style=False): + """ + + :param latex_style: bool, if True returns strings in latex symbols, else in the convention of the sampler + :return: list of the free parameters being sampled in the same order as the sampling + """ + name_list = [] + if self._los_sampling is True: + for i, los_distribution in enumerate(self._los_distributions): + if los_distribution in ["GEV", "GAUSSIAN"]: + if "mean" not in self._kwargs_fixed[i]: + if latex_style is True: + name_list.append(r"$\mu_{\rm los %s}$" %i) + else: + name_list.append(str("mean_los_"+str(i))) + if "sigma" not in self._kwargs_fixed[i]: + if latex_style is True: + name_list.append(r"$\sigma_{\rm los %s}$" %i) + else: + name_list.append(str("sigma_los_"+str(i))) + if los_distribution in ["GEV"]: + if "xi" not in self._kwargs_fixed[i]: + if latex_style is True: + name_list.append(r"$\xi_{\rm los} %s$" %i) + else: + name_list.append(str("xi_los_" + str(i))) + #str(name + "_" + type + str(k)) + return name_list + + def args2kwargs(self, args, i=0): + """ + + :param args: sampling argument list + :param i: index of argument list to start reading out + :return: keyword argument list with parameter names + """ + kwargs = [{} for _ in range(len(self._los_distributions))] + if self._los_sampling is True: + for k, los_distribution in enumerate(self._los_distributions): + if los_distribution in ["GEV", "GAUSSIAN"]: + if "mean" in self._kwargs_fixed[k]: + kwargs[k]["mean"] = self._kwargs_fixed[k]["mean"] + else: + kwargs[k]["mean"] = args[i] + i += 1 + if "sigma" in self._kwargs_fixed[k]: + kwargs[k]["sigma"] = self._kwargs_fixed[k]["sigma"] + else: + kwargs[k]["sigma"] = args[i] + i += 1 + if los_distribution in ["GEV"]: + if "xi" in self._kwargs_fixed[k]: + kwargs[k]["xi"] = self._kwargs_fixed[k]["xi"] + else: + kwargs[k]["xi"] = args[i] + i += 1 + return kwargs, i + + def kwargs2args(self, kwargs): + """ + + :param kwargs: keyword argument list of parameters + :return: sampling argument list in specified order + """ + args = [] + if self._los_sampling is True: + for k, los_distribution in enumerate(self._los_distributions): + if los_distribution in ["GEV", "GAUSSIAN"]: + if "mean" not in self._kwargs_fixed[k]: + args.append(kwargs[k]["mean"]) + if "sigma" not in self._kwargs_fixed[k]: + args.append(kwargs[k]["sigma"]) + if los_distribution in ["GEV"]: + if "xi" not in self._kwargs_fixed[k]: + args.append(kwargs[k]["xi"]) + return args diff --git a/hierarc/Sampling/ParamManager/param_manager.py b/hierarc/Sampling/ParamManager/param_manager.py index aa56ec6b..34c00b51 100644 --- a/hierarc/Sampling/ParamManager/param_manager.py +++ b/hierarc/Sampling/ParamManager/param_manager.py @@ -2,6 +2,7 @@ from hierarc.Sampling.ParamManager.cosmo_param import CosmoParam from hierarc.Sampling.ParamManager.lens_param import LensParam from hierarc.Sampling.ParamManager.source_param import SourceParam +from hierarc.Sampling.ParamManager.los_param import LOSParam class ParamManager(object): @@ -20,8 +21,6 @@ def __init__( gamma_in_distribution="NONE", log_m2l_sampling=False, log_m2l_distribution="NONE", - kappa_ext_sampling=False, - kappa_ext_distribution="NONE", lambda_ifu_sampling=False, lambda_ifu_distribution="NONE", alpha_lambda_sampling=False, @@ -33,6 +32,8 @@ def __init__( sne_distribution="GAUSSIAN", z_apparent_m_anchor=0.1, log_scatter=False, + los_sampling=False, + los_distributions=None, kwargs_lower_cosmo=None, kwargs_upper_cosmo=None, kwargs_fixed_cosmo=None, @@ -45,6 +46,9 @@ def __init__( kwargs_lower_source=None, kwargs_upper_source=None, kwargs_fixed_source=None, + kwargs_lower_los=None, + kwargs_upper_los=None, + kwargs_fixed_los=None, ): """ @@ -78,6 +82,12 @@ def __init__( :param sigma_v_systematics: bool, if True samples paramaters relative to systematics in the velocity dispersion measurement :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space (and thus flat prior in log) + :param los_sampling: if sampling of the parameters should be done + :type los_sampling: bool + :param los_distributions: list of line of sight distributions to be sampled + :type los_distributions: list of str + :param kwargs_fixed_los: fixed arguments in sampling + :type kwargs_fixed_los: list of dictionaries for each los distribution """ self._kin_param = KinParam( anisotropy_sampling=anisotropy_sampling, @@ -101,8 +111,6 @@ def __init__( gamma_in_distribution=gamma_in_distribution, log_m2l_sampling=log_m2l_sampling, log_m2l_distribution=log_m2l_distribution, - kappa_ext_sampling=kappa_ext_sampling, - kappa_ext_distribution=kappa_ext_distribution, alpha_lambda_sampling=alpha_lambda_sampling, beta_lambda_sampling=beta_lambda_sampling, alpha_gamma_in_sampling=alpha_gamma_in_sampling, @@ -116,6 +124,9 @@ def __init__( z_apparent_m_anchor=z_apparent_m_anchor, kwargs_fixed=kwargs_fixed_source, ) + self._los_param = LOSParam(los_sampling=los_sampling, + los_distributions=los_distributions, + kwargs_fixed=kwargs_fixed_los,) self._kwargs_upper_cosmo, self._kwargs_lower_cosmo = ( kwargs_upper_cosmo, kwargs_lower_cosmo, @@ -132,6 +143,10 @@ def __init__( kwargs_upper_source, kwargs_lower_source, ) + self._kwargs_upper_los, self._kwargs_lower_los = ( + kwargs_upper_los, + kwargs_lower_los, + ) @property def num_param(self): @@ -152,6 +167,7 @@ def param_list(self, latex_style=False): list_param += self._lens_param.param_list(latex_style=latex_style) list_param += self._kin_param.param_list(latex_style=latex_style) list_param += self._source_param.param_list(latex_style=latex_style) + list_param += self._los_param.param_list(latex_style=latex_style) return list_param def args2kwargs(self, args): @@ -165,10 +181,12 @@ def args2kwargs(self, args): kwargs_lens, i = self._lens_param.args2kwargs(args, i=i) kwargs_kin, i = self._kin_param.args2kwargs(args, i=i) kwargs_source, i = self._source_param.args2kwargs(args, i=i) - return kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source + kwargs_los, i = self._los_param.args2kwargs(args, i=i) + return kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los def kwargs2args( - self, kwargs_cosmo=None, kwargs_lens=None, kwargs_kin=None, kwargs_source=None + self, kwargs_cosmo=None, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, + kwargs_los=None ): """ @@ -176,6 +194,7 @@ def kwargs2args( :param kwargs_lens: keyword argument list of parameters for lens model sampling :param kwargs_kin: keyword argument list of parameters for kinematic sampling :param kwargs_source: keyword arguments of parameters of source brightness + :param kwargs_los: keyword arguments of parameters of the line of sight :return: sampling argument list in specified order """ args = [] @@ -183,13 +202,15 @@ def kwargs2args( args += self._lens_param.kwargs2args(kwargs_lens) args += self._kin_param.kwargs2args(kwargs_kin) args += self._source_param.kwargs2args(kwargs_source) + args += self._los_param.kwargs2args(kwargs_los) return args def cosmo(self, kwargs_cosmo): """ :param kwargs_cosmo: keyword arguments of parameters (can include others not used for the cosmology) - :return: astropy.cosmology instance + :return: cosmology + :rtype: ~astropy.cosmology instance """ return self._cosmo_param.cosmo(kwargs_cosmo) @@ -204,11 +225,13 @@ def param_bounds(self): kwargs_lens=self._kwargs_lower_lens, kwargs_kin=self._kwargs_lower_kin, kwargs_source=self._kwargs_lower_source, + kwargs_los=self._kwargs_lower_los, ) upper_limit = self.kwargs2args( kwargs_cosmo=self._kwargs_upper_cosmo, kwargs_lens=self._kwargs_upper_lens, kwargs_kin=self._kwargs_upper_kin, kwargs_source=self._kwargs_upper_source, + kwargs_los=self._kwargs_upper_los, ) return lower_limit, upper_limit diff --git a/test/test_Diagnostics/test_goodness_of_fit.py b/test/test_Diagnostics/test_goodness_of_fit.py index b832c8ab..82ec6748 100644 --- a/test/test_Diagnostics/test_goodness_of_fit.py +++ b/test/test_Diagnostics/test_goodness_of_fit.py @@ -104,18 +104,18 @@ def test_plot_ddt_fit(self): kwargs_lens = {"lambda_mst": 1} kwargs_kin = {} f, ax = self.goodnessofFit.plot_ddt_fit( - self.cosmo, kwargs_lens, kwargs_kin, redshift_trend=False + self.cosmo, kwargs_lens, kwargs_kin, kwargs_los=None, redshift_trend=False ) plt.close() f, ax = self.goodnessofFit.plot_ddt_fit( - self.cosmo, kwargs_lens, kwargs_kin, redshift_trend=True + self.cosmo, kwargs_lens, kwargs_kin, kwargs_los=None, redshift_trend=True ) plt.close() def test_plot_kin_fit(self): kwargs_lens = {"lambda_mst": 1} kwargs_kin = {} - f, ax = self.goodnessofFit.plot_kin_fit(self.cosmo, kwargs_lens, kwargs_kin) + f, ax = self.goodnessofFit.plot_kin_fit(self.cosmo, kwargs_lens, kwargs_kin, kwargs_los=None) plt.close() def test_plot_ifu_fit(self): @@ -127,6 +127,7 @@ def test_plot_ifu_fit(self): self.cosmo, kwargs_lens, kwargs_kin, + kwargs_los=None, lens_index=self.ifu_index, bin_edges=1.0, show_legend=True, @@ -152,12 +153,14 @@ def test_raise(self): f, ax = plt.subplots(1, 1, figsize=(4, 4)) kwargs_lens = {"lambda_mst": 1} kwargs_kin = {} + kwargs_los = None cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.0) goodness_of_fit.plot_ifu_fit( ax, cosmo, kwargs_lens, kwargs_kin, + kwargs_los=kwargs_los, lens_index=0, bin_edges=1, show_legend=True, diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py index 86bbafd7..9dc3e676 100644 --- a/test/test_Likelihood/test_hierarchy_likelihood.py +++ b/test/test_Likelihood/test_hierarchy_likelihood.py @@ -1,7 +1,6 @@ from hierarc.Likelihood.hierarchy_likelihood import LensLikelihood from astropy.cosmology import FlatLambdaCDM import pytest -import unittest import numpy as np import numpy.testing as npt @@ -35,7 +34,8 @@ def setup_method(self): ani_scaling_array_list=None, ani_scaling_array=ani_scaling_array, num_distribution_draws=200, - kappa_ext_bias=True, + los_distributions=["GAUSSIAN"], + global_los_distribution=0, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=True, @@ -51,7 +51,8 @@ def setup_method(self): ani_scaling_array_list=None, ani_scaling_array=ani_scaling_array, num_distribution_draws=200, - kappa_ext_bias=False, + los_distributions=["GAUSSIAN"], + global_los_distribution=0, # testing previously set to =False kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, @@ -67,7 +68,8 @@ def setup_method(self): ani_scaling_array_list=None, ani_scaling_array=ani_scaling_array, num_distribution_draws=0, - kappa_ext_bias=True, + los_distributions=["GAUSSIAN"], + global_los_distribution=0, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=True, @@ -86,7 +88,8 @@ def setup_method(self): ani_scaling_array_list=None, ani_scaling_array=ani_scaling_array, num_distribution_draws=200, - kappa_ext_bias=True, + # los_distributions=["GAUSSIAN"], + global_los_distribution=False, kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, mst_ifu=False, @@ -111,7 +114,6 @@ def setup_method(self): gamma_in_array=gamma_in_array, log_m2l_array=log_m2l_array, num_distribution_draws=200, - kappa_ext_bias=False, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, @@ -130,7 +132,6 @@ def setup_method(self): gamma_in_array=gamma_in_array, log_m2l_array=log_m2l_array, num_distribution_draws=200, - kappa_ext_bias=False, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, @@ -149,7 +150,6 @@ def setup_method(self): gamma_in_array=gamma_in_array, log_m2l_array=log_m2l_array, num_distribution_draws=200, - kappa_ext_bias=False, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, @@ -162,43 +162,42 @@ def test_lens_log_likelihood(self): kwargs_lens = { "lambda_mst": 1, "lambda_mst_sigma": 0.01, - "kappa_ext": 0, - "kappa_ext_sigma": 0.03, "lambda_ifu": 1, "lambda_ifu_sigma": 0.01, } + kwargs_los = [{"mean": 0, "sigma": 0.03}] + kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0.1} ln_likelihood = self.likelihood.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) npt.assert_almost_equal(ln_likelihood, -0.5, decimal=1) ln_likelihood_zero = self.likelihood_zero_dist.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) assert ln_likelihood_zero == -np.inf ln_likelihood_kappa_ext = self.likelihood_kappa_ext.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) npt.assert_almost_equal(ln_likelihood, ln_likelihood_kappa_ext, decimal=1) kwargs_lens = { "lambda_mst": 1000000, "lambda_mst_sigma": 0, - "kappa_ext": 0, - "kappa_ext_sigma": 0, "lambda_ifu": 1, "lambda_ifu_sigma": 0, } + kwargs_los = [{"mean": 0, "sigma": 0}] kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0} ln_inf = self.likelihood_single.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) assert ln_inf < -10000000 ln_inf = self.likelihood_single.lens_log_likelihood( - self.cosmo, kwargs_lens=None, kwargs_kin=kwargs_kin + self.cosmo, kwargs_lens=None, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los ) npt.assert_almost_equal(ln_inf, 0.0, decimal=1) @@ -250,7 +249,7 @@ def test_lens_log_likelihood(self): ddt = (1.0 + z_lens) * dd * ds / dds ln_likelihood = self.likelihood_gamma_in_fail_case.log_likelihood_single( - ddt, dd, delta_lum_dist, kwargs_lens, kwargs_kin, kwargs_source + ddt, dd, delta_lum_dist, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los ) assert ln_likelihood < -10000000 diff --git a/test/test_Likelihood/test_lens_sample_likelihood.py b/test/test_Likelihood/test_lens_sample_likelihood.py index 194490ed..8bd3acc0 100644 --- a/test/test_Likelihood/test_lens_sample_likelihood.py +++ b/test/test_Likelihood/test_lens_sample_likelihood.py @@ -54,7 +54,7 @@ def setup_method(self): self.likelihood = LensSampleLikelihood(kwargs_lens_list=self.kwargs_lens_list) def test_log_likelihood(self): - kwargs_lens = {"kappa_ext": 0, "gamma_ppn": 1} + kwargs_lens = {"gamma_ppn": 1} kwargs_kin = {"a_ani": 1} logl = self.likelihood.log_likelihood( self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin diff --git a/test/test_Likelihood/test_los_distribution.py b/test/test_Likelihood/test_los_distribution.py new file mode 100644 index 00000000..c0df55b7 --- /dev/null +++ b/test/test_Likelihood/test_los_distribution.py @@ -0,0 +1,103 @@ +from hierarc.Likelihood.los_distributions import LOSDistribution +from scipy.stats import genextreme +import numpy as np +import numpy.testing as npt + + +class TestLOSDistribution(object): + + def setup_method(self): + pass + + def test_gev(self): + + xi = -0.1 + mean_gev = 0.02 + sigma_gev = np.exp(-5.46) + + mean_gauss = 0.1 + sigma_gauss = 0.2 + + kappa_ext_draw = genextreme.rvs(c=xi, loc=mean_gev, scale=sigma_gev, size=10000) + npt.assert_almost_equal(np.mean(kappa_ext_draw), mean_gev, decimal=2) + npt.assert_almost_equal(np.std(kappa_ext_draw), sigma_gev, decimal=2) + + kappa_pdf, kappa_bin_edges = np.histogram(kappa_ext_draw, bins=100) + kappa_pdf = np.array(kappa_pdf, dtype=float) / np.sum(kappa_pdf) + + los_distribution = ["GAUSSIAN", "GEV"] + + kwargs_los = [{"mean": mean_gauss, "sigma": sigma_gauss}, + {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}] + + # here we draw from the scipy function + dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=1, + los_distributions=los_distribution) + + kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) + npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gev, decimal=2) + npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gev, decimal=2) + + # here we draw from the distribution + dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=False, + los_distributions=los_distribution) + + kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) + npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gev, decimal=2) + npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gev, decimal=2) + + # draw from Gaussian + dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=0, + los_distributions=los_distribution) + + kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) + npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gauss, decimal=2) + npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gauss, decimal=2) + + + def test_draw_bool(self): + xi = -0.1 + mean_gev = 0.02 + sigma_gev = np.exp(-5.46) + + mean_gauss = 0.1 + sigma_gauss = 0 + + kappa_ext_draw = genextreme.rvs(c=xi, loc=mean_gev, scale=sigma_gev, size=10000) + npt.assert_almost_equal(np.mean(kappa_ext_draw), mean_gev, decimal=2) + npt.assert_almost_equal(np.std(kappa_ext_draw), sigma_gev, decimal=2) + + kappa_pdf, kappa_bin_edges = np.histogram(kappa_ext_draw, bins=100) + kappa_pdf = np.array(kappa_pdf, dtype=float) / np.sum(kappa_pdf) + + los_distribution = ["GAUSSIAN", "GEV"] + + kwargs_los = [{"mean": mean_gauss, "sigma": sigma_gauss}, + {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}] + + dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=1, + los_distributions=los_distribution) + bool_draw = dist.draw_bool(kwargs_los) + assert bool_draw is True + + dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=0, + los_distributions=los_distribution) + bool_draw = dist.draw_bool(kwargs_los) + assert bool_draw is False + + dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, + global_los_distribution=False, + los_distributions=los_distribution) + bool_draw = dist.draw_bool(kwargs_los) + assert bool_draw is True + + + + + + diff --git a/test/test_Sampling/test_ParamManager/test_lens_param.py b/test/test_Sampling/test_ParamManager/test_lens_param.py index dd4b96c8..9267d3e2 100644 --- a/test/test_Sampling/test_ParamManager/test_lens_param.py +++ b/test/test_Sampling/test_ParamManager/test_lens_param.py @@ -7,8 +7,6 @@ def setup_method(self): self._param = LensParam( lambda_mst_sampling=True, lambda_mst_distribution="GAUSSIAN", - kappa_ext_sampling=True, - kappa_ext_distribution="GAUSSIAN", lambda_ifu_sampling=True, lambda_ifu_distribution="GAUSSIAN", alpha_lambda_sampling=True, @@ -21,16 +19,12 @@ def setup_method(self): "lambda_mst_sigma": 0.1, "lambda_ifu": 1.1, "lambda_ifu_sigma": 0.2, - "kappa_ext": 0.01, - "kappa_ext_sigma": 0.03, "alpha_lambda": 0, "beta_lambda": 0, } self._param_fixed = LensParam( lambda_mst_sampling=True, lambda_mst_distribution="GAUSSIAN", - kappa_ext_sampling=True, - kappa_ext_distribution="GAUSSIAN", lambda_ifu_sampling=True, lambda_ifu_distribution="GAUSSIAN", alpha_lambda_sampling=True, @@ -40,8 +34,6 @@ def setup_method(self): self._param_log_scatter = LensParam( lambda_mst_sampling=True, lambda_mst_distribution="GAUSSIAN", - kappa_ext_sampling=True, - kappa_ext_distribution="GAUSSIAN", lambda_ifu_sampling=True, lambda_ifu_distribution="GAUSSIAN", alpha_lambda_sampling=True, @@ -52,14 +44,14 @@ def setup_method(self): def test_param_list(self): param_list = self._param.param_list(latex_style=False) - assert len(param_list) == 8 + assert len(param_list) == 6 param_list = self._param.param_list(latex_style=True) - assert len(param_list) == 8 + assert len(param_list) == 6 param_list = self._param_log_scatter.param_list(latex_style=False) - assert len(param_list) == 8 + assert len(param_list) == 6 param_list = self._param_log_scatter.param_list(latex_style=True) - assert len(param_list) == 8 + assert len(param_list) == 6 param_list = self._param_fixed.param_list(latex_style=False) assert len(param_list) == 0 diff --git a/test/test_Sampling/test_ParamManager/test_los_param.py b/test/test_Sampling/test_ParamManager/test_los_param.py new file mode 100644 index 00000000..e86ad012 --- /dev/null +++ b/test/test_Sampling/test_ParamManager/test_los_param.py @@ -0,0 +1,75 @@ +from hierarc.Sampling.ParamManager.los_param import LOSParam +import numpy.testing as npt +import pytest + + +class TestLOSParam(object): + def setup_method(self): + self._param = LOSParam( + los_sampling=True, + los_distributions=["GEV"], + kwargs_fixed=None, + ) + + self._param_gauss = LOSParam( + los_sampling=True, + los_distributions=["GAUSSIAN"], + kwargs_fixed=None, + ) + + kwargs_fixed = [{ + "mean": 0, + "sigma": 0.05, + "xi": 0.1, + }] + self._param_fixed = LOSParam( + los_sampling=True, + los_distributions=["GEV"], + kwargs_fixed=kwargs_fixed, + ) + + def test_param_list(self): + param_list = self._param.param_list(latex_style=False) + assert len(param_list) == 3 + param_list = self._param.param_list(latex_style=True) + assert len(param_list) == 3 + + param_list = self._param_gauss.param_list(latex_style=False) + assert len(param_list) == 2 + param_list = self._param_gauss.param_list(latex_style=True) + assert len(param_list) == 2 + + param_list = self._param_fixed.param_list(latex_style=False) + assert len(param_list) == 0 + param_list = self._param_fixed.param_list(latex_style=True) + assert len(param_list) == 0 + + def test_args2kwargs(self): + kwargs = [{ + "mean": 0.1, + "sigma": 0.1, + "xi": 0.2, + }] + + kwargs_gauss = [{ + "mean": 0.1, + "sigma": 0.1, + }] + args = self._param.kwargs2args(kwargs) + kwargs_new, i = self._param.args2kwargs(args, i=0) + args_new = self._param.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args) + + args = self._param_gauss.kwargs2args(kwargs_gauss) + kwargs_new, i = self._param_gauss.args2kwargs(args, i=0) + args_new = self._param_gauss.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args) + + args = self._param_fixed.kwargs2args(kwargs) + kwargs_new, i = self._param_fixed.args2kwargs(args, i=0) + args_new = self._param_fixed.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args) + + +if __name__ == "__main__": + pytest.main() diff --git a/test/test_Sampling/test_ParamManager/test_param_manager.py b/test/test_Sampling/test_ParamManager/test_param_manager.py index 85f0ec58..37f54568 100644 --- a/test/test_Sampling/test_ParamManager/test_param_manager.py +++ b/test/test_Sampling/test_ParamManager/test_param_manager.py @@ -19,9 +19,8 @@ def setup_method(self): kwargs_lower_lens = { "lambda_mst": 0, "lambda_mst_sigma": 0.1, - "kappa_ext": -0.2, - "kappa_ext_sigma": 0, } + kwargs_lower_los = [{"mean": -0.2, "sigma": 0.0}] kwargs_lower_kin = {"a_ani": 0.1, "a_ani_sigma": 0.1} kwargs_lower_source = {"mu_sne": 0, "sigma_sne": 0} @@ -37,9 +36,8 @@ def setup_method(self): kwargs_upper_lens = { "lambda_mst": 2, "lambda_mst_sigma": 0.1, - "kappa_ext": 0.2, - "kappa_ext_sigma": 1, } + kwargs_upper_los = [{"mean": 0.2, "sigma": 0.5}] kwargs_upper_kin = {"a_ani": 0.1, "a_ani_sigma": 0.1} kwargs_upper_source = {"mu_sne": 100, "sigma_sne": 10} @@ -55,9 +53,8 @@ def setup_method(self): kwargs_fixed_lens = { "lambda_mst": 1, "lambda_mst_sigma": 0.1, - "kappa_ext": 0, - "kappa_ext_sigma": 0, } + kwargs_fixed_los = [{"mean": 0, "sigma": 0.0}] kwargs_fixed_kin = {"a_ani": 0.1, "a_ani_sigma": 0.1} kwargs_fixed_source = {"mu_sne": 1, "sigma_sne": 0.1} @@ -73,8 +70,8 @@ def setup_method(self): anisotropy_sampling=True, anisotropy_model="OM", kwargs_lower_cosmo=kwargs_lower_cosmo, - kappa_ext_sampling=True, - kappa_ext_distribution="GAUSSIAN", + los_sampling=True, + los_distributions=["GAUSSIAN"], sne_apparent_m_sampling=True, sne_distribution="GAUSSIAN", kwargs_upper_cosmo=kwargs_upper_cosmo, @@ -88,6 +85,9 @@ def setup_method(self): kwargs_fixed_source=kwargs_fixed_source, kwargs_lower_source=kwargs_lower_source, kwargs_upper_source=kwargs_upper_source, + kwargs_fixed_los=kwargs_fixed_los, + kwargs_lower_los=kwargs_lower_los, + kwargs_upper_los=kwargs_upper_los ) ) @@ -100,8 +100,8 @@ def setup_method(self): anisotropy_distribution="GAUSSIAN", anisotropy_sampling=True, anisotropy_model="OM", - kappa_ext_sampling=True, - kappa_ext_distribution="GAUSSIAN", + los_sampling=True, + los_distributions=["GAUSSIAN"], sne_apparent_m_sampling=True, sne_distribution="GAUSSIAN", kwargs_lower_cosmo=kwargs_lower_cosmo, @@ -116,6 +116,9 @@ def setup_method(self): kwargs_lower_source=kwargs_lower_source, kwargs_upper_source=kwargs_upper_source, kwargs_fixed_source=None, + kwargs_lower_los=kwargs_lower_los, + kwargs_upper_los=kwargs_upper_los, + kwargs_fixed_los=None ) ) self.param_list = param_list @@ -142,9 +145,9 @@ def test_kwargs2args(self): kwargs_lens = { "lambda_mst": 1, "lambda_mst_sigma": 0, - "kappa_ext": 0, - "kappa_ext_sigma": 0, + } + kwargs_los = [{"mean": 0, "sigma": 0.05}] kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0.3} kwargs_source = {"mu_sne": 2, "sigma_sne": 0.2} for param in self.param_list: @@ -153,15 +156,17 @@ def test_kwargs2args(self): kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_source=kwargs_source, + kwargs_los=kwargs_los ) ( kwargs_cosmo_new, kwargs_lens_new, kwargs_kin_new, kwargs_source_new, + kwargs_los_new ) = param.args2kwargs(args) args_new = param.kwargs2args( - kwargs_cosmo_new, kwargs_lens_new, kwargs_kin_new, kwargs_source_new + kwargs_cosmo_new, kwargs_lens_new, kwargs_kin_new, kwargs_source_new, kwargs_los_new ) npt.assert_almost_equal(args_new, args) From f1cf577e411d69443fadf2503619e017037f30f5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 28 Apr 2024 22:22:01 +0000 Subject: [PATCH 2/4] [pre-commit.ci] auto fixes from pre-commit.com hooks --- hierarc/Diagnostics/goodness_of_fit.py | 13 +++- .../LensLikelihood/double_source_plane.py | 7 +- hierarc/Likelihood/cosmo_likelihood.py | 9 ++- hierarc/Likelihood/hierarchy_likelihood.py | 26 +++++-- hierarc/Likelihood/lens_sample_likelihood.py | 12 ++- hierarc/Likelihood/los_distributions.py | 36 ++++++--- hierarc/Sampling/ParamManager/los_param.py | 13 ++-- .../Sampling/ParamManager/param_manager.py | 16 ++-- test/test_Diagnostics/test_goodness_of_fit.py | 4 +- .../test_hierarchy_likelihood.py | 20 ++++- test/test_Likelihood/test_los_distribution.py | 75 +++++++++++-------- .../test_ParamManager/test_los_param.py | 34 +++++---- .../test_ParamManager/test_param_manager.py | 15 ++-- 13 files changed, 184 insertions(+), 96 deletions(-) diff --git a/hierarc/Diagnostics/goodness_of_fit.py b/hierarc/Diagnostics/goodness_of_fit.py index 8f14cb14..c65158d0 100644 --- a/hierarc/Diagnostics/goodness_of_fit.py +++ b/hierarc/Diagnostics/goodness_of_fit.py @@ -65,7 +65,9 @@ def plot_ddt_fit( ddt_model_sigma, dd_model_mean, dd_model_sigma, - ) = likelihood.ddt_dd_model_prediction(cosmo, kwargs_lens=kwargs_lens, kwargs_los=kwargs_los) + ) = likelihood.ddt_dd_model_prediction( + cosmo, kwargs_lens=kwargs_lens, kwargs_los=kwargs_los + ) ddt_name_list.append(name) ddt_model_mean_list.append(ddt_model_mean) @@ -162,7 +164,10 @@ def kin_fit(self, cosmo, kwargs_lens, kwargs_kin, kwargs_los): sigma_v_predict_mean, cov_error_predict, ) = likelihood.sigma_v_measured_vs_predict( - cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los + cosmo, + kwargs_lens=kwargs_lens, + kwargs_kin=kwargs_kin, + kwargs_los=kwargs_los, ) if sigma_v_measurement is not None: @@ -205,7 +210,9 @@ def plot_kin_fit( :param color_prediction: color of model prediction :return: fig, axes of matplotlib instance """ - logL = self._sample_likelihood.log_likelihood(cosmo, kwargs_lens, kwargs_kin, kwargs_los=kwargs_los) + logL = self._sample_likelihood.log_likelihood( + cosmo, kwargs_lens, kwargs_kin, kwargs_los=kwargs_los + ) print(logL, "log likelihood") ( sigma_v_name_list, diff --git a/hierarc/Likelihood/LensLikelihood/double_source_plane.py b/hierarc/Likelihood/LensLikelihood/double_source_plane.py index b742183a..c0b5fa12 100644 --- a/hierarc/Likelihood/LensLikelihood/double_source_plane.py +++ b/hierarc/Likelihood/LensLikelihood/double_source_plane.py @@ -30,7 +30,12 @@ def __init__( self._normalized = normalized def lens_log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, kwargs_los=None + self, + cosmo, + kwargs_lens=None, + kwargs_kin=None, + kwargs_source=None, + kwargs_los=None, ): """Log likelihood of the data given a model. diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index 46c9c2ae..c627dbe2 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -108,8 +108,9 @@ def __init__( if sigma_v_systematics is True: normalized = True self._likelihoodLensSample = LensSampleLikelihood( - kwargs_likelihood_list, normalized=normalized, - los_distributions=los_distributions + kwargs_likelihood_list, + normalized=normalized, + los_distributions=los_distributions, ) self.param = ParamManager( cosmology, @@ -191,8 +192,8 @@ def likelihood(self, args, kwargs_cosmo_interp=None): if args[i] < self._lower_limit[i] or args[i] > self._upper_limit[i]: return -np.inf - kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los = self.param.args2kwargs( - args + kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los = ( + self.param.args2kwargs(args) ) if self._cosmology == "oLCDM": # assert we are not in a crazy cosmological situation that prevents computing the angular distance integral diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py index 3294023b..5ff14904 100644 --- a/hierarc/Likelihood/hierarchy_likelihood.py +++ b/hierarc/Likelihood/hierarchy_likelihood.py @@ -127,9 +127,12 @@ def __init__( ) self._num_distribution_draws = int(num_distribution_draws) - self._los = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=global_los_distribution, - los_distributions=los_distributions) + self._los = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=global_los_distribution, + los_distributions=los_distributions, + ) self._mst_ifu = mst_ifu self._lambda_scaling_property = lambda_scaling_property self._lambda_scaling_property_beta = lambda_scaling_property_beta @@ -140,7 +143,12 @@ def __init__( self._gamma_in_prior_std = gamma_in_prior_std def lens_log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, kwargs_los=None, + self, + cosmo, + kwargs_lens=None, + kwargs_kin=None, + kwargs_source=None, + kwargs_los=None, ): """Log likelihood of the data of a lens given a model (defined with hyper- parameters) and cosmology. @@ -149,7 +157,8 @@ def lens_log_likelihood( :param kwargs_lens: keywords of the hyper parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model hyper parameters :param kwargs_source: keyword argument of the source model (such as SNe) - :param kwargs_los: list of keyword arguments of global line of sight distributions + :param kwargs_los: list of keyword arguments of global line of sight + distributions :return: log likelihood of the data given the model """ @@ -194,7 +203,8 @@ def hyper_param_likelihood( :param kwargs_lens: keywords of the hyper parameters of the lens model :param kwargs_kin: keyword arguments of the kinematic model hyper parameters :param kwargs_source: keyword argument of the source model (such as SNe) - :param kwargs_los: list of keyword arguments of global line of sight distributions + :param kwargs_los: list of keyword arguments of global line of sight + distributions :param cosmo: astropy.cosmology instance :return: log likelihood given the single lens analysis for the given hyper parameter @@ -529,7 +539,9 @@ def draw_source(mu_sne=1, sigma_sne=0, lum_dist=0, **kwargs): # return linear amplitude with base log 10 return mag_source - def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_los=None): + def sigma_v_measured_vs_predict( + self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_los=None + ): """Mean and error covariance of velocity dispersion measurement mean and error covariance of velocity dispersion predictions. diff --git a/hierarc/Likelihood/lens_sample_likelihood.py b/hierarc/Likelihood/lens_sample_likelihood.py index 3b0bc539..dedf24e9 100644 --- a/hierarc/Likelihood/lens_sample_likelihood.py +++ b/hierarc/Likelihood/lens_sample_likelihood.py @@ -28,11 +28,19 @@ def __init__(self, kwargs_lens_list, normalized=False, los_distributions=None): ) else: self._lens_list.append( - LensLikelihood(normalized=normalized, los_distributions=los_distributions, **kwargs_lens) + LensLikelihood( + normalized=normalized, + los_distributions=los_distributions, + **kwargs_lens + ) ) def log_likelihood( - self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, + self, + cosmo, + kwargs_lens=None, + kwargs_kin=None, + kwargs_source=None, kwargs_los=None, ): """ diff --git a/hierarc/Likelihood/los_distributions.py b/hierarc/Likelihood/los_distributions.py index 9a7bd8a5..2e02b3e8 100644 --- a/hierarc/Likelihood/los_distributions.py +++ b/hierarc/Likelihood/los_distributions.py @@ -3,11 +3,15 @@ class LOSDistribution(object): - """ - line of sight distribution drawing - """ - def __init__(self, kappa_pdf=None, kappa_bin_edges=None, global_los_distribution=False, - los_distributions=None): + """Line of sight distribution drawing.""" + + def __init__( + self, + kappa_pdf=None, + kappa_bin_edges=None, + global_los_distribution=False, + los_distributions=None, + ): """ :param global_los_distribution: if integer, will draw from the global kappa distribution specified in that @@ -21,12 +25,19 @@ def __init__(self, kappa_pdf=None, kappa_bin_edges=None, global_los_distribution """ self._global_los_distribution = global_los_distribution - if isinstance(self._global_los_distribution, int) and self._global_los_distribution is not False: + if ( + isinstance(self._global_los_distribution, int) + and self._global_los_distribution is not False + ): self._draw_kappa_global = True self._los_distribution = los_distributions[global_los_distribution] else: self._draw_kappa_global = False - if kappa_pdf is not None and kappa_bin_edges is not None and not self._draw_kappa_global: + if ( + kappa_pdf is not None + and kappa_bin_edges is not None + and not self._draw_kappa_global + ): print("test kappa pdf sampling") self._kappa_dist = PDFSampling( bin_edges=kappa_bin_edges, pdf_array=kappa_pdf @@ -36,8 +47,7 @@ def __init__(self, kappa_pdf=None, kappa_bin_edges=None, global_los_distribution self._draw_kappa_individual = False def draw_los(self, kwargs_los, size=1): - """ - Draw from the distribution of line of sight convergence + """Draw from the distribution of line of sight convergence. :param kwargs_los: line of sight parameters :type kwargs_los: list of dict @@ -58,16 +68,18 @@ def draw_los(self, kwargs_los, size=1): sigma = kwargs_los_i["sigma"] xi = kwargs_los_i["xi"] from scipy.stats import genextreme + kappa_ext_draw = genextreme.rvs(c=xi, loc=mean, scale=sigma, size=size) else: - raise ValueError("line of sight distribution %s not valid." % self._los_distribution) + raise ValueError( + "line of sight distribution %s not valid." % self._los_distribution + ) else: kappa_ext_draw = 0 return kappa_ext_draw def draw_bool(self, kwargs_los): - """ - whether single-valued or extended distribution (need to draw from) + """Whether single-valued or extended distribution (need to draw from) :param kwargs_los: list of keyword arguments for line of sight distributions :return: boolean, True with samples need to be drawn, else False diff --git a/hierarc/Sampling/ParamManager/los_param.py b/hierarc/Sampling/ParamManager/los_param.py index 3ed09df3..e2292215 100644 --- a/hierarc/Sampling/ParamManager/los_param.py +++ b/hierarc/Sampling/ParamManager/los_param.py @@ -10,7 +10,6 @@ def __init__( los_sampling=False, los_distributions=None, kwargs_fixed=None, - ): """ @@ -47,21 +46,21 @@ def param_list(self, latex_style=False): if los_distribution in ["GEV", "GAUSSIAN"]: if "mean" not in self._kwargs_fixed[i]: if latex_style is True: - name_list.append(r"$\mu_{\rm los %s}$" %i) + name_list.append(r"$\mu_{\rm los %s}$" % i) else: - name_list.append(str("mean_los_"+str(i))) + name_list.append(str("mean_los_" + str(i))) if "sigma" not in self._kwargs_fixed[i]: if latex_style is True: - name_list.append(r"$\sigma_{\rm los %s}$" %i) + name_list.append(r"$\sigma_{\rm los %s}$" % i) else: - name_list.append(str("sigma_los_"+str(i))) + name_list.append(str("sigma_los_" + str(i))) if los_distribution in ["GEV"]: if "xi" not in self._kwargs_fixed[i]: if latex_style is True: - name_list.append(r"$\xi_{\rm los} %s$" %i) + name_list.append(r"$\xi_{\rm los} %s$" % i) else: name_list.append(str("xi_los_" + str(i))) - #str(name + "_" + type + str(k)) + # str(name + "_" + type + str(k)) return name_list def args2kwargs(self, args, i=0): diff --git a/hierarc/Sampling/ParamManager/param_manager.py b/hierarc/Sampling/ParamManager/param_manager.py index 34c00b51..11413773 100644 --- a/hierarc/Sampling/ParamManager/param_manager.py +++ b/hierarc/Sampling/ParamManager/param_manager.py @@ -124,9 +124,11 @@ def __init__( z_apparent_m_anchor=z_apparent_m_anchor, kwargs_fixed=kwargs_fixed_source, ) - self._los_param = LOSParam(los_sampling=los_sampling, - los_distributions=los_distributions, - kwargs_fixed=kwargs_fixed_los,) + self._los_param = LOSParam( + los_sampling=los_sampling, + los_distributions=los_distributions, + kwargs_fixed=kwargs_fixed_los, + ) self._kwargs_upper_cosmo, self._kwargs_lower_cosmo = ( kwargs_upper_cosmo, kwargs_lower_cosmo, @@ -185,8 +187,12 @@ def args2kwargs(self, args): return kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los def kwargs2args( - self, kwargs_cosmo=None, kwargs_lens=None, kwargs_kin=None, kwargs_source=None, - kwargs_los=None + self, + kwargs_cosmo=None, + kwargs_lens=None, + kwargs_kin=None, + kwargs_source=None, + kwargs_los=None, ): """ diff --git a/test/test_Diagnostics/test_goodness_of_fit.py b/test/test_Diagnostics/test_goodness_of_fit.py index 82ec6748..0c9be7d8 100644 --- a/test/test_Diagnostics/test_goodness_of_fit.py +++ b/test/test_Diagnostics/test_goodness_of_fit.py @@ -115,7 +115,9 @@ def test_plot_ddt_fit(self): def test_plot_kin_fit(self): kwargs_lens = {"lambda_mst": 1} kwargs_kin = {} - f, ax = self.goodnessofFit.plot_kin_fit(self.cosmo, kwargs_lens, kwargs_kin, kwargs_los=None) + f, ax = self.goodnessofFit.plot_kin_fit( + self.cosmo, kwargs_lens, kwargs_kin, kwargs_los=None + ) plt.close() def test_plot_ifu_fit(self): diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py index 9dc3e676..a150d70c 100644 --- a/test/test_Likelihood/test_hierarchy_likelihood.py +++ b/test/test_Likelihood/test_hierarchy_likelihood.py @@ -169,17 +169,26 @@ def test_lens_log_likelihood(self): kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0.1} ln_likelihood = self.likelihood.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los + self.cosmo, + kwargs_lens=kwargs_lens, + kwargs_kin=kwargs_kin, + kwargs_los=kwargs_los, ) npt.assert_almost_equal(ln_likelihood, -0.5, decimal=1) ln_likelihood_zero = self.likelihood_zero_dist.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los + self.cosmo, + kwargs_lens=kwargs_lens, + kwargs_kin=kwargs_kin, + kwargs_los=kwargs_los, ) assert ln_likelihood_zero == -np.inf ln_likelihood_kappa_ext = self.likelihood_kappa_ext.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los + self.cosmo, + kwargs_lens=kwargs_lens, + kwargs_kin=kwargs_kin, + kwargs_los=kwargs_los, ) npt.assert_almost_equal(ln_likelihood, ln_likelihood_kappa_ext, decimal=1) @@ -192,7 +201,10 @@ def test_lens_log_likelihood(self): kwargs_los = [{"mean": 0, "sigma": 0}] kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0} ln_inf = self.likelihood_single.lens_log_likelihood( - self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los + self.cosmo, + kwargs_lens=kwargs_lens, + kwargs_kin=kwargs_kin, + kwargs_los=kwargs_los, ) assert ln_inf < -10000000 diff --git a/test/test_Likelihood/test_los_distribution.py b/test/test_Likelihood/test_los_distribution.py index c0df55b7..4fdef5be 100644 --- a/test/test_Likelihood/test_los_distribution.py +++ b/test/test_Likelihood/test_los_distribution.py @@ -27,37 +27,47 @@ def test_gev(self): los_distribution = ["GAUSSIAN", "GEV"] - kwargs_los = [{"mean": mean_gauss, "sigma": sigma_gauss}, - {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}] + kwargs_los = [ + {"mean": mean_gauss, "sigma": sigma_gauss}, + {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}, + ] # here we draw from the scipy function - dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=1, - los_distributions=los_distribution) + dist_gev = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=1, + los_distributions=los_distribution, + ) kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gev, decimal=2) npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gev, decimal=2) # here we draw from the distribution - dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=False, - los_distributions=los_distribution) + dist_gev = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=False, + los_distributions=los_distribution, + ) kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gev, decimal=2) npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gev, decimal=2) # draw from Gaussian - dist_gev = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=0, - los_distributions=los_distribution) + dist_gev = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=0, + los_distributions=los_distribution, + ) kappa_dist_drawn = dist_gev.draw_los(kwargs_los, size=10000) npt.assert_almost_equal(np.mean(kappa_dist_drawn), mean_gauss, decimal=2) npt.assert_almost_equal(np.std(kappa_dist_drawn), sigma_gauss, decimal=2) - def test_draw_bool(self): xi = -0.1 mean_gev = 0.02 @@ -75,29 +85,34 @@ def test_draw_bool(self): los_distribution = ["GAUSSIAN", "GEV"] - kwargs_los = [{"mean": mean_gauss, "sigma": sigma_gauss}, - {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}] - - dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=1, - los_distributions=los_distribution) + kwargs_los = [ + {"mean": mean_gauss, "sigma": sigma_gauss}, + {"mean": mean_gev, "sigma": sigma_gev, "xi": xi}, + ] + + dist = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=1, + los_distributions=los_distribution, + ) bool_draw = dist.draw_bool(kwargs_los) assert bool_draw is True - dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=0, - los_distributions=los_distribution) + dist = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=0, + los_distributions=los_distribution, + ) bool_draw = dist.draw_bool(kwargs_los) assert bool_draw is False - dist = LOSDistribution(kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, - global_los_distribution=False, - los_distributions=los_distribution) + dist = LOSDistribution( + kappa_pdf=kappa_pdf, + kappa_bin_edges=kappa_bin_edges, + global_los_distribution=False, + los_distributions=los_distribution, + ) bool_draw = dist.draw_bool(kwargs_los) assert bool_draw is True - - - - - - diff --git a/test/test_Sampling/test_ParamManager/test_los_param.py b/test/test_Sampling/test_ParamManager/test_los_param.py index e86ad012..85b95ab1 100644 --- a/test/test_Sampling/test_ParamManager/test_los_param.py +++ b/test/test_Sampling/test_ParamManager/test_los_param.py @@ -17,11 +17,13 @@ def setup_method(self): kwargs_fixed=None, ) - kwargs_fixed = [{ - "mean": 0, - "sigma": 0.05, - "xi": 0.1, - }] + kwargs_fixed = [ + { + "mean": 0, + "sigma": 0.05, + "xi": 0.1, + } + ] self._param_fixed = LOSParam( los_sampling=True, los_distributions=["GEV"], @@ -45,16 +47,20 @@ def test_param_list(self): assert len(param_list) == 0 def test_args2kwargs(self): - kwargs = [{ - "mean": 0.1, - "sigma": 0.1, - "xi": 0.2, - }] + kwargs = [ + { + "mean": 0.1, + "sigma": 0.1, + "xi": 0.2, + } + ] - kwargs_gauss = [{ - "mean": 0.1, - "sigma": 0.1, - }] + kwargs_gauss = [ + { + "mean": 0.1, + "sigma": 0.1, + } + ] args = self._param.kwargs2args(kwargs) kwargs_new, i = self._param.args2kwargs(args, i=0) args_new = self._param.kwargs2args(kwargs_new) diff --git a/test/test_Sampling/test_ParamManager/test_param_manager.py b/test/test_Sampling/test_ParamManager/test_param_manager.py index 37f54568..a8967af9 100644 --- a/test/test_Sampling/test_ParamManager/test_param_manager.py +++ b/test/test_Sampling/test_ParamManager/test_param_manager.py @@ -87,7 +87,7 @@ def setup_method(self): kwargs_upper_source=kwargs_upper_source, kwargs_fixed_los=kwargs_fixed_los, kwargs_lower_los=kwargs_lower_los, - kwargs_upper_los=kwargs_upper_los + kwargs_upper_los=kwargs_upper_los, ) ) @@ -118,7 +118,7 @@ def setup_method(self): kwargs_fixed_source=None, kwargs_lower_los=kwargs_lower_los, kwargs_upper_los=kwargs_upper_los, - kwargs_fixed_los=None + kwargs_fixed_los=None, ) ) self.param_list = param_list @@ -145,7 +145,6 @@ def test_kwargs2args(self): kwargs_lens = { "lambda_mst": 1, "lambda_mst_sigma": 0, - } kwargs_los = [{"mean": 0, "sigma": 0.05}] kwargs_kin = {"a_ani": 1, "a_ani_sigma": 0.3} @@ -156,17 +155,21 @@ def test_kwargs2args(self): kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_source=kwargs_source, - kwargs_los=kwargs_los + kwargs_los=kwargs_los, ) ( kwargs_cosmo_new, kwargs_lens_new, kwargs_kin_new, kwargs_source_new, - kwargs_los_new + kwargs_los_new, ) = param.args2kwargs(args) args_new = param.kwargs2args( - kwargs_cosmo_new, kwargs_lens_new, kwargs_kin_new, kwargs_source_new, kwargs_los_new + kwargs_cosmo_new, + kwargs_lens_new, + kwargs_kin_new, + kwargs_source_new, + kwargs_los_new, ) npt.assert_almost_equal(args_new, args) From 077e42f9807789b82ad54e96a78be31173ff3e16 Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Sun, 28 Apr 2024 21:37:19 -0400 Subject: [PATCH 3/4] improved unit testing --- test/test_Likelihood/test_los_distribution.py | 12 ++++++++---- .../test_ParamManager/test_los_param.py | 12 ++++++++++++ 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/test/test_Likelihood/test_los_distribution.py b/test/test_Likelihood/test_los_distribution.py index c0df55b7..463cb55a 100644 --- a/test/test_Likelihood/test_los_distribution.py +++ b/test/test_Likelihood/test_los_distribution.py @@ -2,6 +2,7 @@ from scipy.stats import genextreme import numpy as np import numpy.testing as npt +import unittest class TestLOSDistribution(object): @@ -97,7 +98,10 @@ def test_draw_bool(self): assert bool_draw is True - - - - +class TestRaise(unittest.TestCase): + def test_raise(self): + with self.assertRaises(ValueError): + los = LOSDistribution(kappa_pdf=None, kappa_bin_edges=None, + global_los_distribution=0, + los_distributions=["BAD"]) + los.draw_los(kwargs_los=[{}]) diff --git a/test/test_Sampling/test_ParamManager/test_los_param.py b/test/test_Sampling/test_ParamManager/test_los_param.py index e86ad012..c786fef1 100644 --- a/test/test_Sampling/test_ParamManager/test_los_param.py +++ b/test/test_Sampling/test_ParamManager/test_los_param.py @@ -1,6 +1,7 @@ from hierarc.Sampling.ParamManager.los_param import LOSParam import numpy.testing as npt import pytest +import unittest class TestLOSParam(object): @@ -71,5 +72,16 @@ def test_args2kwargs(self): npt.assert_almost_equal(args_new, args) +class TestRaise(unittest.TestCase): + def test_raise(self): + with self.assertRaises(ValueError): + LOSParam( + los_sampling=True, + los_distributions=["BAD"], + kwargs_fixed=[{}], + ) + + + if __name__ == "__main__": pytest.main() From 93d44e215e5f3f727298d58e1a42a9b8513b2630 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 29 Apr 2024 01:38:15 +0000 Subject: [PATCH 4/4] [pre-commit.ci] auto fixes from pre-commit.com hooks --- test/test_Likelihood/test_los_distribution.py | 9 ++++++--- test/test_Sampling/test_ParamManager/test_los_param.py | 1 - 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_Likelihood/test_los_distribution.py b/test/test_Likelihood/test_los_distribution.py index a611aa0e..3ac5434e 100644 --- a/test/test_Likelihood/test_los_distribution.py +++ b/test/test_Likelihood/test_los_distribution.py @@ -122,7 +122,10 @@ def test_draw_bool(self): class TestRaise(unittest.TestCase): def test_raise(self): with self.assertRaises(ValueError): - los = LOSDistribution(kappa_pdf=None, kappa_bin_edges=None, - global_los_distribution=0, - los_distributions=["BAD"]) + los = LOSDistribution( + kappa_pdf=None, + kappa_bin_edges=None, + global_los_distribution=0, + los_distributions=["BAD"], + ) los.draw_los(kwargs_los=[{}]) diff --git a/test/test_Sampling/test_ParamManager/test_los_param.py b/test/test_Sampling/test_ParamManager/test_los_param.py index dfb280c6..3715c4dc 100644 --- a/test/test_Sampling/test_ParamManager/test_los_param.py +++ b/test/test_Sampling/test_ParamManager/test_los_param.py @@ -88,6 +88,5 @@ def test_raise(self): ) - if __name__ == "__main__": pytest.main()