diff --git a/hierarc/Diagnostics/goodness_of_fit.py b/hierarc/Diagnostics/goodness_of_fit.py index 8f14cb1..c65158d 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 b742183..c0b5fa1 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 46c9c2a..c627dbe 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 3294023..5ff1490 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 3b0bc53..dedf24e 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 9a7bd8a..2e02b3e 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 3ed09df..e229221 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 34c00b5..1141377 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 82ec674..0c9be7d 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 9dc3e67..a150d70 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 463cb55..a611aa0 100644 --- a/test/test_Likelihood/test_los_distribution.py +++ b/test/test_Likelihood/test_los_distribution.py @@ -28,37 +28,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 @@ -76,24 +86,35 @@ 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 c786fef..dfb280c 100644 --- a/test/test_Sampling/test_ParamManager/test_los_param.py +++ b/test/test_Sampling/test_ParamManager/test_los_param.py @@ -18,11 +18,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"], @@ -46,16 +48,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_gauss = [{ - "mean": 0.1, - "sigma": 0.1, - }] + 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) diff --git a/test/test_Sampling/test_ParamManager/test_param_manager.py b/test/test_Sampling/test_ParamManager/test_param_manager.py index 37f5456..a8967af 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)