From a38dd3b97cd0aa469eb014201e602a51b83ea0e0 Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Wed, 26 Jun 2024 21:55:37 -0400 Subject: [PATCH 01/14] major restructuring of parameter distribution and design of the interpolation grid for the kinematic scaling --- hierarc/LensPosterior/base_config.py | 12 +- hierarc/LensPosterior/ddt_kin_constraints.py | 5 +- .../ddt_kin_gauss_constraints.py | 5 +- hierarc/LensPosterior/kin_constraints.py | 20 +- .../kin_constraints_composite.py | 50 ++- ...otropy_config.py => kin_scaling_config.py} | 38 +- hierarc/Likelihood/cosmo_likelihood.py | 86 +--- hierarc/Likelihood/hierarchy_likelihood.py | 334 +++++--------- hierarc/Likelihood/kin_scaling.py | 138 ++++++ hierarc/Likelihood/lens_sample_likelihood.py | 50 ++- hierarc/Likelihood/parameter_scaling.py | 227 ---------- hierarc/Sampling/Distributions/__init__.py | 0 .../Distributions/anisotropy_distributions.py | 78 ++++ .../Distributions/lens_distribution.py | 189 ++++++++ .../Distributions}/los_distributions.py | 0 .../Sampling/ParamManager/param_manager.py | 2 +- .../test_anisotropy_config.py | 24 +- .../test_kin_constraints_composite.py | 52 +-- test/test_Likelihood/test_cosmo_likelihood.py | 27 +- .../test_hierarchy_likelihood.py | 110 +++-- test/test_Likelihood/test_kin_scaling.py | 269 +++++++++++ .../test_lens_sample_likelihood.py | 7 +- test/test_Likelihood/test_los_distribution.py | 2 +- .../test_Likelihood/test_parameter_scaling.py | 417 ------------------ .../test_Distributions/__init__.py | 0 .../test_anisotropy_distribution.py | 36 ++ .../test_lens_distribution.py | 51 +++ test/test_Sampling/test_mcmc_sampling.py | 13 +- 28 files changed, 1129 insertions(+), 1113 deletions(-) rename hierarc/LensPosterior/{anisotropy_config.py => kin_scaling_config.py} (65%) create mode 100644 hierarc/Likelihood/kin_scaling.py delete mode 100644 hierarc/Likelihood/parameter_scaling.py create mode 100644 hierarc/Sampling/Distributions/__init__.py create mode 100644 hierarc/Sampling/Distributions/anisotropy_distributions.py create mode 100644 hierarc/Sampling/Distributions/lens_distribution.py rename hierarc/{Likelihood => Sampling/Distributions}/los_distributions.py (100%) create mode 100644 test/test_Likelihood/test_kin_scaling.py delete mode 100644 test/test_Likelihood/test_parameter_scaling.py create mode 100644 test/test_Sampling/test_Distributions/__init__.py create mode 100644 test/test_Sampling/test_Distributions/test_anisotropy_distribution.py create mode 100644 test/test_Sampling/test_Distributions/test_lens_distribution.py diff --git a/hierarc/LensPosterior/base_config.py b/hierarc/LensPosterior/base_config.py index 9fef09ec..258c12f7 100644 --- a/hierarc/LensPosterior/base_config.py +++ b/hierarc/LensPosterior/base_config.py @@ -1,9 +1,9 @@ from lenstronomy.Analysis.td_cosmography import TDCosmography from hierarc.LensPosterior.imaging_constraints import ImageModelPosterior -from hierarc.LensPosterior.anisotropy_config import AnisotropyConfig +from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig -class BaseLensConfig(TDCosmography, ImageModelPosterior, AnisotropyConfig): +class BaseLensConfig(TDCosmography, ImageModelPosterior, KinScalingConfig): """This class contains and manages the base configurations of the lens posteriors and makes sure that they are universally applied consistently through the different likelihood definitions.""" @@ -33,6 +33,8 @@ def __init__( num_kin_sampling=1000, multi_observations=False, cosmo_fiducial=None, + gamma_in_scaling=None, + log_m2l_scaling=None, ): """ @@ -62,6 +64,8 @@ def __init__( light profile :param cosmo_fiducial: astropy.cosmology instance, if None, uses astropy's default cosmology + :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None) + :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None) """ self._z_lens, self._z_source = z_lens, z_source @@ -105,4 +109,6 @@ def __init__( ImageModelPosterior.__init__( self, theta_E, theta_E_error, gamma, gamma_error, r_eff, r_eff_error ) - AnisotropyConfig.__init__(self, anisotropy_model, r_eff) + KinScalingConfig.__init__(self, anisotropy_model, r_eff, + gamma_in_scaling=gamma_in_scaling, + log_m2l_scaling=log_m2l_scaling) diff --git a/hierarc/LensPosterior/ddt_kin_constraints.py b/hierarc/LensPosterior/ddt_kin_constraints.py index 68b23726..3c825e2a 100644 --- a/hierarc/LensPosterior/ddt_kin_constraints.py +++ b/hierarc/LensPosterior/ddt_kin_constraints.py @@ -124,7 +124,8 @@ def hierarchy_configuration(self, num_sample_model=20): "j_model": j_model_list, "error_cov_measurement": error_cov_measurement, "error_cov_j_sqrt": error_cov_j_sqrt, - "ani_param_array": self.ani_param_array, - "ani_scaling_array_list": ani_scaling_array_list, + "kin_scaling_param_list": self.param_name_list, + "j_kin_scaling_param_axes": self.kin_scaling_param_array, + "j_kin_scaling_grid_list": ani_scaling_array_list, } return kwargs_likelihood diff --git a/hierarc/LensPosterior/ddt_kin_gauss_constraints.py b/hierarc/LensPosterior/ddt_kin_gauss_constraints.py index 4112c74a..6507d485 100644 --- a/hierarc/LensPosterior/ddt_kin_gauss_constraints.py +++ b/hierarc/LensPosterior/ddt_kin_gauss_constraints.py @@ -123,7 +123,8 @@ def hierarchy_configuration(self, num_sample_model=20): "j_model": j_model_list, "error_cov_measurement": error_cov_measurement, "error_cov_j_sqrt": error_cov_j_sqrt, - "ani_param_array": self.ani_param_array, - "ani_scaling_array_list": ani_scaling_array_list, + "kin_scaling_param_list": self.param_name_list, + "j_kin_scaling_param_axes": self.kin_scaling_param_array, + "j_kin_scaling_grid_list": ani_scaling_array_list, } return kwargs_likelihood diff --git a/hierarc/LensPosterior/kin_constraints.py b/hierarc/LensPosterior/kin_constraints.py index 74147ff1..7bae1ea4 100644 --- a/hierarc/LensPosterior/kin_constraints.py +++ b/hierarc/LensPosterior/kin_constraints.py @@ -36,6 +36,8 @@ def __init__( num_kin_sampling=1000, multi_observations=False, cosmo_fiducial=None, + gamma_in_scaling=None, + log_m2l_scaling=None, ): """ @@ -77,6 +79,8 @@ def __init__( kwargs_seeing as lists of multiple observations :param cosmo_fiducial: astropy.cosmology instance, if None, uses astropy's default + :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None) + :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None) """ self._sigma_v_measured = np.array(sigma_v_measured) self._sigma_v_error_independent = np.array(sigma_v_error_independent) @@ -111,6 +115,8 @@ def __init__( num_kin_sampling=num_kin_sampling, multi_observations=multi_observations, cosmo_fiducial=cosmo_fiducial, + gamma_in_scaling=gamma_in_scaling, + log_m2l_scaling=log_m2l_scaling, ) def j_kin_draw(self, kwargs_anisotropy, no_error=False): @@ -171,8 +177,10 @@ def hierarchy_configuration(self, num_sample_model=20): "j_model": j_model_list, "error_cov_measurement": error_cov_measurement, "error_cov_j_sqrt": error_cov_j_sqrt, - "ani_param_array": self.ani_param_array, - "ani_scaling_array_list": ani_scaling_array_list, + "kin_scaling_param_list": self.param_name_list, + "j_kin_scaling_param_axes": self.kin_scaling_param_array, + "j_kin_scaling_grid_list": ani_scaling_array_list, + } return kwargs_likelihood @@ -245,11 +253,11 @@ def _anisotropy_scaling_relative(self, j_ani_0): if self._anisotropy_model == "GOM": ani_scaling_array_list = [ - np.zeros((len(self.ani_param_array[0]), len(self.ani_param_array[1]))) + np.zeros((len(self.kin_scaling_param_array[0]), len(self.kin_scaling_param_array[1]))) for _ in range(num_data) ] - for i, a_ani in enumerate(self.ani_param_array[0]): - for j, beta_inf in enumerate(self.ani_param_array[1]): + for i, a_ani in enumerate(self.kin_scaling_param_array[0]): + for j, beta_inf in enumerate(self.kin_scaling_param_array[1]): kwargs_anisotropy = self.anisotropy_kwargs( a_ani=a_ani, beta_inf=beta_inf ) @@ -260,7 +268,7 @@ def _anisotropy_scaling_relative(self, j_ani_0): ) # perhaps change the order elif self._anisotropy_model in ["OM", "const"]: ani_scaling_array_list = [[] for _ in range(num_data)] - for a_ani in self.ani_param_array: + for a_ani in self.kin_scaling_param_array[0]: kwargs_anisotropy = self.anisotropy_kwargs(a_ani) j_kin_ani = self.j_kin_draw(kwargs_anisotropy, no_error=True) for i, j_kin in enumerate(j_kin_ani): diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py index fdc2cf0e..99e87118 100644 --- a/hierarc/LensPosterior/kin_constraints_composite.py +++ b/hierarc/LensPosterior/kin_constraints_composite.py @@ -115,6 +115,12 @@ def __init__( lens_model_list = ["GNFW", "MULTI_GAUSSIAN_KAPPA"] + # log_m2l is interpolated when sampled on the population level, otherwise marginalized + if is_m2l_population_level: + log_m2l_scaling = log_m2l_array + else: + log_m2l_scaling = None + super(KinConstraintsComposite, self).__init__( z_lens, z_source, @@ -142,6 +148,8 @@ def __init__( num_psf_sampling=num_psf_sampling, num_kin_sampling=num_kin_sampling, multi_observations=multi_observations, + gamma_in_scaling=gamma_in_array, + log_m2l_scaling=log_m2l_scaling, ) if self._check_arrays(kappa_s_array, r_s_angle_array): @@ -392,16 +400,22 @@ def hierarchy_configuration(self, num_sample_model=20): "j_model": j_model_list, "error_cov_measurement": error_cov_measurement, "error_cov_j_sqrt": error_cov_j_sqrt, - "ani_param_array": self.ani_param_array, - "gamma_in_array": self.gamma_in_array, - "log_m2l_array": self.log_m2l_array, - "param_scaling_grid_list": ani_scaling_grid_list, + + #"ani_param_array": self.kin_scaling_param_array, + #"gamma_in_array": self.gamma_in_array, + #"log_m2l_array": self.log_m2l_array, + #"param_scaling_grid_list": ani_scaling_grid_list, + "gamma_in_prior_mean": self._gamma_in_prior_mean, "gamma_in_prior_std": self._gamma_in_prior_std, + + "kin_scaling_param_list": self.param_name_list, + "j_kin_scaling_param_axes": self.kin_scaling_param_array, + "j_kin_scaling_grid_list": ani_scaling_grid_list, } - if not self._is_m2l_population_level: - kwargs_likelihood["log_m2l_array"] = None + #if not self._is_m2l_population_level: + # kwargs_likelihood["log_m2l_array"] = None return kwargs_likelihood def anisotropy_scaling(self): @@ -438,16 +452,16 @@ def _anisotropy_scaling_relative(self, j_ani_0): ani_scaling_grid_list = [ np.zeros( ( - len(self.ani_param_array[0]), - len(self.ani_param_array[1]), + len(self.kin_scaling_param_array[0]), + len(self.kin_scaling_param_array[1]), len(self.gamma_in_array), len(self.log_m2l_array), ) ) for _ in range(num_data) ] - for i, a_ani in enumerate(self.ani_param_array[0]): - for j, beta_inf in enumerate(self.ani_param_array[1]): + for i, a_ani in enumerate(self.kin_scaling_param_array[0]): + for j, beta_inf in enumerate(self.kin_scaling_param_array[1]): for k, g_in in enumerate(self.gamma_in_array): for l, log_m2l in enumerate(self.log_m2l_array): kwargs_anisotropy = self.anisotropy_kwargs( @@ -466,14 +480,14 @@ def _anisotropy_scaling_relative(self, j_ani_0): ani_scaling_grid_list = [ np.zeros( ( - len(self.ani_param_array), + len(self.kin_scaling_param_array[0]), len(self.gamma_in_array), len(self.log_m2l_array), ) ) for _ in range(num_data) ] - for i, a_ani in enumerate(self.ani_param_array): + for i, a_ani in enumerate(self.kin_scaling_param_array[0]): for k, g_in in enumerate(self.gamma_in_array): for l, log_m2l in enumerate(self.log_m2l_array): kwargs_anisotropy = self.anisotropy_kwargs(a_ani) @@ -500,15 +514,15 @@ def _anisotropy_scaling_relative_m2l(self, j_ani_0): ani_scaling_grid_list = [ np.zeros( ( - len(self.ani_param_array[0]), - len(self.ani_param_array[1]), + len(self.kin_scaling_param_array[0]), + len(self.kin_scaling_param_array[1]), len(self.gamma_in_array), ) ) for _ in range(num_data) ] - for i, a_ani in enumerate(self.ani_param_array[0]): - for j, beta_inf in enumerate(self.ani_param_array[1]): + for i, a_ani in enumerate(self.kin_scaling_param_array[0]): + for j, beta_inf in enumerate(self.kin_scaling_param_array[1]): for k, g_in in enumerate(self.gamma_in_array): kwargs_anisotropy = self.anisotropy_kwargs( a_ani=a_ani, beta_inf=beta_inf @@ -524,13 +538,13 @@ def _anisotropy_scaling_relative_m2l(self, j_ani_0): ani_scaling_grid_list = [ np.zeros( ( - len(self.ani_param_array), + len(self.kin_scaling_param_array[0]), len(self.gamma_in_array), ) ) for _ in range(num_data) ] - for i, a_ani in enumerate(self.ani_param_array): + for i, a_ani in enumerate(self.kin_scaling_param_array[0]): for k, g_in in enumerate(self.gamma_in_array): kwargs_anisotropy = self.anisotropy_kwargs(a_ani) j_kin_ani = self.j_kin_draw_composite_m2l( diff --git a/hierarc/LensPosterior/anisotropy_config.py b/hierarc/LensPosterior/kin_scaling_config.py similarity index 65% rename from hierarc/LensPosterior/anisotropy_config.py rename to hierarc/LensPosterior/kin_scaling_config.py index d827d408..8a51b2ad 100644 --- a/hierarc/LensPosterior/anisotropy_config.py +++ b/hierarc/LensPosterior/kin_scaling_config.py @@ -1,36 +1,47 @@ import numpy as np -class AnisotropyConfig(object): +class KinScalingConfig(object): """Class to manage the anisotropy model and parameters for the Posterior processing.""" - def __init__(self, anisotropy_model, r_eff): + def __init__(self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scaling=None): """ :param anisotropy_model: type of stellar anisotropy model. Supported are 'OM' and 'GOM' or 'const', see details in lenstronomy.Galkin module :param r_eff: half-light radius of the deflector galaxy + :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None) + :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None) """ self._r_eff = r_eff self._anisotropy_model = anisotropy_model + self._param_name_list = [] if self._anisotropy_model == "OM": - self._ani_param_array = np.array( - [0.1, 0.2, 0.5, 1, 2, 5] - ) # used for r_ani OsipkovMerritt anisotropy description + self._ani_param_array = [np.array([0.1, 0.2, 0.5, 1, 2, 5])] + # used for r_ani OsipkovMerritt anisotropy description + self._param_name_list = ["a_ani"] elif self._anisotropy_model == "GOM": self._ani_param_array = [ np.array([0.1, 0.2, 0.5, 1, 2, 5]), np.array([0, 0.5, 0.8, 1]), ] + self._param_name_list = ["a_ani", "beta_inf"] elif self._anisotropy_model == "const": - self._ani_param_array = np.linspace( - -0.49, 1, 7 - ) # used for constant anisotropy description + self._ani_param_array = [np.linspace(-0.49, 1, 7)] # used for constant anisotropy description + self._param_name_list = ["a_ani"] + elif self._anisotropy_model == "NONE": + self._param_name_list =[] else: raise ValueError( "anisotropy model %s not supported." % self._anisotropy_model ) + if gamma_in_scaling is not None: + self._param_name_list.append("gamma_in") + self._ani_param_array.append(np.array(gamma_in_scaling)) + if log_m2l_scaling is not None: + self._param_name_list.append("log_m2l") + self._ani_param_array.append(np.array(log_m2l_scaling)) @property def kwargs_anisotropy_base(self): @@ -57,13 +68,22 @@ def kwargs_anisotropy_base(self): return kwargs_anisotropy_0 @property - def ani_param_array(self): + def kin_scaling_param_array(self): """ :return: numpy array of anisotropy parameter values to be explored """ return self._ani_param_array + @property + def param_name_list(self): + """ + list of parameters in same order as interpolated + + :return: + """ + return self._param_name_list + def anisotropy_kwargs(self, a_ani, beta_inf=None): """ diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index d5d22a38..f900cbfa 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -14,35 +14,13 @@ def __init__( self, kwargs_likelihood_list, cosmology, + kwargs_model, kwargs_bounds, sne_likelihood=None, kwargs_sne_likelihood=None, KDE_likelihood_chain=None, kwargs_kde_likelihood=None, - ppn_sampling=False, - lambda_mst_sampling=False, - lambda_mst_distribution="delta", - anisotropy_sampling=False, - gamma_in_sampling=False, - gamma_in_distribution="NONE", - log_m2l_sampling=False, - log_m2l_distribution="NONE", - los_sampling=False, - los_distributions=None, - alpha_lambda_sampling=False, - beta_lambda_sampling=False, - alpha_gamma_in_sampling=False, - alpha_log_m2l_sampling=False, - lambda_ifu_sampling=False, - lambda_ifu_distribution="NONE", - sigma_v_systematics=False, - sne_apparent_m_sampling=False, - sne_distribution="GAUSSIAN", - z_apparent_m_anchor=0.1, - log_scatter=False, normalized=False, - anisotropy_model="OM", - anisotropy_distribution="NONE", custom_prior=None, interpolate_cosmo=True, num_redshift_interp=100, @@ -52,6 +30,8 @@ def __init__( :param kwargs_likelihood_list: keyword argument list specifying the arguments of the LensLikelihood class :param cosmology: string describing cosmological model + :param kwargs_model: model settings for ParamManager() class + :type kwargs_model: dict :param kwargs_bounds: keyword arguments of the lower and upper bounds and parameters that are held fixed. Includes: 'kwargs_lower_lens', 'kwargs_upper_lens', 'kwargs_fixed_lens', @@ -63,38 +43,6 @@ def __init__( :param sne_likelihood: (string), optional. Sampling supernovae relative expansion history likelihood, see SneLikelihood module for options :param kwargs_sne_likelihood: keyword argument for the SNe likelihood, see SneLikelihood module for options - :param ppn_sampling:post-newtonian parameter sampling - :param lambda_mst_sampling: bool, if True adds a global mass-sheet transform parameter in the sampling - :param lambda_mst_distribution: string, defines the distribution function of lambda_mst - :param lambda_ifu_sampling: bool, if True samples a separate lambda_mst for a second (e.g. IFU) data set - independently - :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter - :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly - 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 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 - :param anisotropy_distribution: string, distribution of the anisotropy parameters - :param gamma_in_sampling: bool, if True samples gNFW inner slope parameter - :param gamma_in_distribution: string, distribution function of the gamma_in parameter - :param log_m2l_sampling: bool, if True samples a global mass-to-light ratio parameter in logarithmic scale - :param log_m2l_distribution: string, distribution function of the log_m2l parameter - :param sigma_v_systematics: bool, if True samples paramaters relative to systematics in the velocity dispersion - measurement - :param sne_apparent_m_sampling: boolean, if True, samples/queries SNe unlensed magnitude distribution - (not intrinsic magnitudes but apparent!) - :param sne_distribution: string, apparent non-lensed brightness distribution (in linear space). - Currently supports: - 'GAUSSIAN': Gaussian distribution - :param z_apparent_m_anchor: redshift of pivot/anchor at which the apparent SNe brightness is defined relative to - :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space - (and thus flat prior in log) :param custom_prior: None or a definition that takes the keywords from the CosmoParam conventions and returns a log likelihood value (e.g. prior) :param interpolate_cosmo: bool, if True, uses interpolated comoving distance in the calculation for speed-up @@ -105,38 +53,16 @@ def __init__( """ self._cosmology = cosmology self._kwargs_lens_list = kwargs_likelihood_list - if sigma_v_systematics is True: + if kwargs_model.get("sigma_v_systematics", False) is True: normalized = True self._likelihoodLensSample = LensSampleLikelihood( kwargs_likelihood_list, normalized=normalized, - los_distributions=los_distributions, + kwargs_global_model=kwargs_model, ) self.param = ParamManager( cosmology, - ppn_sampling=ppn_sampling, - lambda_mst_sampling=lambda_mst_sampling, - lambda_mst_distribution=lambda_mst_distribution, - lambda_ifu_sampling=lambda_ifu_sampling, - lambda_ifu_distribution=lambda_ifu_distribution, - alpha_lambda_sampling=alpha_lambda_sampling, - beta_lambda_sampling=beta_lambda_sampling, - gamma_in_sampling=gamma_in_sampling, - gamma_in_distribution=gamma_in_distribution, - log_m2l_sampling=log_m2l_sampling, - log_m2l_distribution=log_m2l_distribution, - alpha_gamma_in_sampling=alpha_gamma_in_sampling, - alpha_log_m2l_sampling=alpha_log_m2l_sampling, - sne_apparent_m_sampling=sne_apparent_m_sampling, - sne_distribution=sne_distribution, - z_apparent_m_anchor=z_apparent_m_anchor, - sigma_v_systematics=sigma_v_systematics, - los_sampling=los_sampling, - los_distributions=los_distributions, - anisotropy_sampling=anisotropy_sampling, - anisotropy_model=anisotropy_model, - anisotropy_distribution=anisotropy_distribution, - log_scatter=log_scatter, + **kwargs_model, **kwargs_bounds ) self._lower_limit, self._upper_limit = self.param.param_bounds diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py index 5ff14904..2b5aaf7b 100644 --- a/hierarc/Likelihood/hierarchy_likelihood.py +++ b/hierarc/Likelihood/hierarchy_likelihood.py @@ -1,40 +1,59 @@ 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.Likelihood.los_distributions import LOSDistribution +from hierarc.Likelihood.kin_scaling import KinScaling +from hierarc.Sampling.Distributions.los_distributions import LOSDistribution +from hierarc.Sampling.Distributions.anisotropy_distributions import AnisotropyDistribution +from hierarc.Sampling.Distributions.lens_distribution import LensDistribution import numpy as np import copy -class LensLikelihood(TransformedCosmography, LensLikelihoodBase, ParameterScalingIFU): +class LensLikelihood(TransformedCosmography, LensLikelihoodBase, KinScaling): """Master class containing the likelihood definitions of different analysis for a single lens.""" def __init__( self, + # properties of the lens z_lens, z_source, name="name", likelihood_type="TDKin", + lambda_scaling_property=0, + lambda_scaling_property_beta=0, + kwargs_lens_properties=None, + # specific distribution settings for individual lenses + global_los_distribution=False, + mst_ifu=False, + # global distributions anisotropy_model="NONE", - ani_param_array=None, - ani_scaling_array=None, - ani_scaling_array_list=None, - gamma_in_array=None, - log_m2l_array=None, - param_scaling_grid_list=None, + anisotropy_sampling=False, + anisotroy_distribution_function="NONE", # make sure input is provided + los_distributions=None, + lambda_mst_distribution="NONE", + gamma_in_sampling=False, + gamma_in_distribution="NONE", + log_m2l_sampling=False, + log_m2l_distribution="NONE", + alpha_lambda_sampling=False, + beta_lambda_sampling=False, + alpha_gamma_in_sampling=False, + alpha_log_m2l_sampling=False, + log_scatter=False, + # kinematic model quantities + kin_scaling_param_list=None, + j_kin_scaling_param_axes=None, + j_kin_scaling_grid_list=None, + # likelihood evaluation quantities num_distribution_draws=50, - global_los_distribution=False, + normalized=True, + # kappa quantities kappa_pdf=None, kappa_bin_edges=None, - mst_ifu=False, - lambda_scaling_property=0, - lambda_scaling_property_beta=0, - normalized=True, - kwargs_lens_properties=None, - gamma_in_prior_mean=None, + # priors + gamma_in_prior_mean=None, # TODO: make a separate prior class with inputs gamma_in_prior_std=None, - los_distributions=None, + # specifics for each lens **kwargs_likelihood ): """ @@ -43,14 +62,10 @@ def __init__( :param z_source: source redshift :param name: string (optional) to name the specific lens :param likelihood_type: string to specify the likelihood type - :param ani_param_array: array of anisotropy parameter values for which the kinematics are predicted - :param ani_scaling_array: velocity dispersion sigma**2 scaling (also J scaling) of anisotropy parameter relative - to default prediction. The scaling corresponds to the ani_param_array parameter spacing - (to generate an interpolation function). A value =1 in ani_scaling_array results in the value stored in the - provided J() predictions. - :param param_scaling_grid_list: list of N-dimensional arrays with the - scalings of J() for each IFU. Needed when simultaneously scaling - anisotropy, gamma_in, and log_m2l. In that case, gamma_in_array and log_m2l_array need to be provided. + :param j_kin_scaling_param_axes: array of parameter values for each axes of j_kin_scaling_grid + :param j_kin_scaling_grid_list: list of array with the scalings of J() for each IFU + :param j_kin_scaling_param_name_list: list of strings for the parameters as they are interpolated in the same + order as j_kin_scaling_grid :param num_distribution_draws: int, number of distribution draws from the likelihood that are being averaged over :param global_los_distribution: if integer, will draw from the global kappa distribution specified in that @@ -74,46 +89,16 @@ def __init__( 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 + :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens + kinematic prediction """ 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: - ani_scaling_array_list = [ani_scaling_array] - - # AnisotropyScalingIFU.__init__( - # self, - # anisotropy_model=anisotropy_model, - # ani_param_array=ani_param_array, - # ani_scaling_array_list=ani_scaling_array_list, - # ) - if gamma_in_array is not None and log_m2l_array is not None: - if isinstance(ani_param_array, list): - param_arrays = ani_param_array + [gamma_in_array, log_m2l_array] - else: - param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] - ParameterScalingIFU.__init__( - self, - anisotropy_model=anisotropy_model, - param_arrays=param_arrays, - scaling_grid_list=param_scaling_grid_list, - ) - elif gamma_in_array is not None and log_m2l_array is None: - if isinstance(ani_param_array, list): - param_arrays = ani_param_array + [gamma_in_array] - else: - param_arrays = [ani_param_array, gamma_in_array] - ParameterScalingIFU.__init__( - self, - anisotropy_model=anisotropy_model, - param_arrays=param_arrays, - scaling_grid_list=param_scaling_grid_list, - ) - else: - ParameterScalingIFU.__init__( - self, - anisotropy_model=anisotropy_model, - param_arrays=ani_param_array, - scaling_grid_list=ani_scaling_array_list, - ) + + KinScaling.__init__(self, + j_kin_scaling_param_axes=j_kin_scaling_param_axes, + j_kin_scaling_grid_list=j_kin_scaling_grid_list, + j_kin_scaling_param_name_list=kin_scaling_param_list + ) LensLikelihoodBase.__init__( self, @@ -133,11 +118,29 @@ def __init__( 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 - self._gamma_in_array = gamma_in_array - self._log_m2l_array = log_m2l_array + kwargs_min, kwargs_max = self.param_bounds_interpol() + self._lens_distribution = LensDistribution(lambda_mst_sampling=False, + lambda_mst_distribution=lambda_mst_distribution, + gamma_in_sampling=gamma_in_sampling, + gamma_in_distribution=gamma_in_distribution, + log_m2l_sampling=log_m2l_sampling, + log_m2l_distribution=log_m2l_distribution, + alpha_lambda_sampling=alpha_lambda_sampling, + beta_lambda_sampling=beta_lambda_sampling, + alpha_gamma_in_sampling=alpha_gamma_in_sampling, + alpha_log_m2l_sampling=alpha_log_m2l_sampling, + log_scatter=log_scatter, + mst_ifu=mst_ifu, + lambda_scaling_property=lambda_scaling_property, + lambda_scaling_property_beta=lambda_scaling_property_beta, + kwargs_min=kwargs_min, + kwargs_max=kwargs_max,) + + self._aniso_distribution = AnisotropyDistribution(anisotropy_model=anisotropy_model, + anisotropy_sampling=anisotropy_sampling, + distribution_function=anisotroy_distribution_function, + kwargs_anisotropy_min=kwargs_min, + kwargs_anisotropy_max=kwargs_max) self._gamma_in_prior_mean = gamma_in_prior_mean self._gamma_in_prior_std = gamma_in_prior_std @@ -265,17 +268,19 @@ def log_likelihood_single( :param ddt: time-delay distance :param dd: angular diameter distance to the deflector :param delta_lum_dist: relative luminosity distance to pivot redshift - :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_lens: keywords of the hyperparameters of the lens model + :param kwargs_kin: keyword arguments of the kinematic model hyperparameters :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 + realization of the hyperparameter distribution """ - lambda_mst, gamma_ppn = self.draw_lens(**kwargs_lens) + kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] 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( @@ -286,13 +291,9 @@ def log_likelihood_single( kappa_ext=kappa_ext, mag_source=mag_source, ) - try: - scaling_param_array = self.draw_scaling_params( - kwargs_lens=kwargs_lens, **kwargs_kin - ) - except ValueError: - return np.nan_to_num(-np.inf) - kin_scaling = self.param_scaling(scaling_param_array) + kwargs_kin_draw = self._aniso_distribution.draw_anisotropy(**kwargs_kin) + kwargs_param = {**kwargs_lens_draw, **kwargs_kin_draw} + kin_scaling = self.kin_scaling(kwargs_param) lnlikelihood = self.log_likelihood( ddt_, @@ -305,97 +306,15 @@ def log_likelihood_single( if ( self._gamma_in_prior_mean is not None and self._gamma_in_prior_std is not None + and "gamma_in" in kwargs_lens_draw ): - if self._gamma_in_array is not None and self._log_m2l_array is not None: - lnlikelihood -= ( - self._gamma_in_prior_mean - scaling_param_array[-2] - ) ** 2 / (2 * self._gamma_in_prior_std**2) - elif self._gamma_in_array is not None and self._log_m2l_array is None: - lnlikelihood -= ( - self._gamma_in_prior_mean - scaling_param_array[-1] - ) ** 2 / (2 * self._gamma_in_prior_std**2) + gamma_in = kwargs_lens_draw["gamma_in"] + lnlikelihood -= ( + self._gamma_in_prior_mean - gamma_in + ) ** 2 / (2 * self._gamma_in_prior_std**2) return np.nan_to_num(lnlikelihood) - def draw_scaling_params(self, kwargs_lens=None, **kwargs_kin): - """Draws a realization of the anisotropy parameter scaling from the distribution - function. - - :return: array of anisotropy parameter scaling - """ - ani_param = self.draw_anisotropy(**kwargs_kin) - if self._gamma_in_array is not None and self._log_m2l_array is not None: - gamma_in, log_m2l = self.draw_lens_scaling_params(**kwargs_lens) - return np.concatenate([ani_param, [gamma_in, log_m2l]]) - elif self._gamma_in_array is not None and self._log_m2l_array is None: - gamma_in = self.draw_lens_scaling_params(**kwargs_lens) - return np.concatenate([ani_param, [gamma_in]]) - else: - return ani_param - - def draw_lens_scaling_params( - 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, - alpha_lambda=0, - beta_lambda=0, - gamma_in=1, - gamma_in_sigma=0, - alpha_gamma_in=0, - log_m2l=1, - log_m2l_sigma=0, - alpha_log_m2l=0, - ): - """Draws a realization of the anisotropy parameter scaling from the - distribution. - - :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 - :param lambda_ifu_sigma: secondary lambda_mst_sigma parameter for subset of - lenses specified for - :param alpha_lambda: float, linear slope of the lambda_int scaling relation with - lens quantity self._lambda_scaling_property - :param beta_lambda: float, a second linear slope of the lambda_int scaling - relation with lens quantity self._lambda_scaling_property_beta - :param gamma_in: inner slope of the NFW profile - :param gamma_in_sigma: spread in the distribution - :param alpha_gamma_in: float, linear slope of the gamma_in scaling relation with - lens quantity self._lambda_scaling_property - :param log_m2l: log(mass-to-light ratio) - :param log_m2l_sigma: spread in the distribution - :param alpha_log_m2l: float, linear slope of the log(m2l) scaling relation with - lens quantity self._lambda_scaling_property - :return: draw from the distributions - """ - if self._gamma_in_array is not None and self._log_m2l_array is not None: - gamma_in_draw, log_m2l_draw = self.draw_lens_parameters( - gamma_in + alpha_gamma_in * self._lambda_scaling_property, - gamma_in_sigma, - log_m2l + alpha_log_m2l * self._lambda_scaling_property, - log_m2l_sigma, - ) - return gamma_in_draw, log_m2l_draw - - elif self._gamma_in_array is not None and self._log_m2l_array is None: - gamma_in_draw = self.draw_lens_parameters( - gamma_in + alpha_gamma_in * self._lambda_scaling_property, - gamma_in_sigma, - ) - return gamma_in_draw - - else: - return None - def angular_diameter_distances(self, cosmo): """Time-delay distance Ddt, angular diameter distance to the lens (dd) @@ -441,12 +360,12 @@ def luminosity_distance_modulus(self, cosmo, z_apparent_m_anchor): 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. + hyperparameters 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 + :param kwargs_lens: lens model hyperparameter keywords + :param kwargs_kin: kinematic model hyperparameter keywords + :param kwargs_source: source brightness hyperparameter keywords + :param kwargs_los: list of dictionaries for line of sight hyperparameters :return: bool, True if delta function, else False """ lambda_mst_sigma = kwargs_lens.get("lambda_mst_sigma", 0) # scatter in MST @@ -464,62 +383,6 @@ def check_dist(self, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los): return True return False - def draw_lens( - self, - lambda_mst=1, - lambda_mst_sigma=0, - gamma_ppn=1, - lambda_ifu=1, - lambda_ifu_sigma=0, - alpha_lambda=0, - beta_lambda=0, - gamma_in=1, - gamma_in_sigma=0, - alpha_gamma_in=0, - log_m2l=1, - log_m2l_sigma=0, - alpha_log_m2l=0, - ): - """Draws a realization of a specific model from the hyper-parameter - distribution. - - :param lambda_mst: MST transform - :param lambda_mst_sigma: spread in the distribution - :param gamma_ppn: Post-Newtonian parameter - :param lambda_ifu: secondary lambda_mst parameter for subset of lenses specified - for - :param lambda_ifu_sigma: secondary lambda_mst_sigma parameter for subset of - lenses specified for - :param alpha_lambda: float, linear slope of the lambda_int scaling relation with - lens quantity self._lambda_scaling_property - :param beta_lambda: float, a second linear slope of the lambda_int scaling - relation with lens quantity self._lambda_scaling_property_beta - :param gamma_in: inner slope of the NFW profile - :param gamma_in_sigma: spread in the distribution - :param alpha_gamma_in: float, linear slope of the gamma_in scaling relation with - lens quantity self._lambda_scaling_property - :param log_m2l: log(mass-to-light ratio) - :param log_m2l_sigma: spread in the distribution - :param alpha_log_m2l: float, linear slope of the log(m2l) scaling relation with - lens quantity self._lambda_scaling_property - :return: draw from the distributions - """ - if self._mst_ifu is True: - lambda_lens = ( - lambda_ifu - + alpha_lambda * self._lambda_scaling_property - + beta_lambda * self._lambda_scaling_property_beta - ) - lambda_mst_draw = np.random.normal(lambda_lens, lambda_ifu_sigma) - else: - lambda_lens = ( - lambda_mst - + alpha_lambda * self._lambda_scaling_property - + beta_lambda * self._lambda_scaling_property_beta - ) - lambda_mst_draw = np.random.normal(lambda_lens, lambda_mst_sigma) - return lambda_mst_draw, gamma_ppn - @staticmethod def draw_source(mu_sne=1, sigma_sne=0, lum_dist=0, **kwargs): """Draws a source magnitude from a distribution specified by population @@ -546,8 +409,8 @@ def sigma_v_measured_vs_predict( 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_lens: keywords of the hyperparameters of the lens model + :param kwargs_kin: keyword arguments of the kinematic model hyperparameters :param kwargs_los: line of sight parapers :return: sigma_v_measurement, cov_error_measurement, sigma_v_predict_mean, cov_error_predict @@ -569,15 +432,15 @@ def sigma_v_measured_vs_predict( 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, gamma_ppn = self.draw_lens(**kwargs_lens) + kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] 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 ) - scaling_param_array = self.draw_scaling_params( - kwargs_lens=kwargs_lens, **kwargs_kin_copy - ) - kin_scaling = self.param_scaling(scaling_param_array) + kwargs_kin_draw = self._aniso_distribution.draw_anisotropy(**kwargs_kin) + kwargs_param = {**kwargs_lens_draw, **kwargs_kin_draw} + kin_scaling = self.kin_scaling(kwargs_param) sigma_v_predict_i, cov_error_predict_i = self.sigma_v_prediction( ddt_, dd_, kin_scaling=kin_scaling ) @@ -614,7 +477,8 @@ def ddt_dd_model_prediction(self, cosmo, kwargs_lens=None, kwargs_los=None): ddt_draws = [] dd_draws = [] for i in range(self._num_distribution_draws): - lambda_mst, gamma_ppn = self.draw_lens(**kwargs_lens) + kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) + lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] 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/kin_scaling.py b/hierarc/Likelihood/kin_scaling.py new file mode 100644 index 00000000..f016adf8 --- /dev/null +++ b/hierarc/Likelihood/kin_scaling.py @@ -0,0 +1,138 @@ +__author__ = "sibirrer", "ajshajib" + +from scipy.interpolate import interp1d +from scipy.interpolate import interp2d +from scipy.interpolate import RegularGridInterpolator +import numpy as np + + +class ParameterScalingSingleMeasurement(object): + """Class to manage anisotropy scaling for single slit observation.""" + + def __init__(self, param_grid_axes, j_kin_scaling_grid): + """ + + :param param_grid_axes: list of arrays of interpolated parameter values + :param j_kin_scaling_grid: array with the scaling of J() for single measurement bin in same dimensions as the + param_arrays + """ + self._evalute_scaling = False + # check if param arrays is 1d list or 2d list + if param_grid_axes is not None and j_kin_scaling_grid is not None: + if isinstance(param_grid_axes, list): + self._dim_scaling = len(param_grid_axes) + else: + self._dim_scaling = 1 + param_grid_axes = [param_grid_axes] + + if self._dim_scaling == 1: + self._f_ani = interp1d(param_grid_axes[0], j_kin_scaling_grid, kind="linear") + elif self._dim_scaling == 2: + self._f_ani = interp2d(param_grid_axes[0], param_grid_axes[1], j_kin_scaling_grid.T) + else: + self._f_ani = RegularGridInterpolator( + tuple(param_grid_axes), + j_kin_scaling_grid, + ) + self._evalute_scaling = True + + def j_scaling(self, param_array): + """ + + :param param_array: sorted list of parameters for the interpolation function + :return: scaling J(a_ani) for single slit + """ + + if self._evalute_scaling is not True or len(param_array) == 0: + return 1 + if self._dim_scaling == 1: + return self._f_ani(param_array[0]) + elif self._dim_scaling == 2: + return self._f_ani(param_array[0], param_array[1])[0] + else: + return self._f_ani(param_array)[0] + + +class KinScaling(object): + """Class to manage model parameter and anisotropy scalings for IFU data.""" + + def __init__( + self, j_kin_scaling_param_axes=None, j_kin_scaling_grid_list=None, j_kin_scaling_param_name_list=None + ): + """ + + :param j_kin_scaling_param_axes: array of parameter values for each axes of j_kin_scaling_grid + :param j_kin_scaling_grid_list: list of array with the scalings of J() for each IFU + :param j_kin_scaling_param_name_list: list of strings for the parameters as they are interpolated in the same + order as j_kin_scaling_grid + """ + if j_kin_scaling_param_name_list is None: + self._param_list = [] + else: + self._param_list = j_kin_scaling_param_name_list + self._param_arrays = j_kin_scaling_param_axes + if not isinstance(j_kin_scaling_param_axes, list) and j_kin_scaling_param_name_list is not None: + self._param_arrays = [j_kin_scaling_param_axes] + self._evaluate_scaling = False + self._is_log_m2l_population_level = False + if ( + j_kin_scaling_param_axes is not None + and j_kin_scaling_grid_list is not None + and j_kin_scaling_param_name_list is not None + ): + self._evaluate_scaling = True + self._j_scaling_ifu = [] + self._f_ani_list = [] + for scaling_grid in j_kin_scaling_grid_list: + self._j_scaling_ifu.append( + ParameterScalingSingleMeasurement(j_kin_scaling_param_axes, scaling_grid) + ) + + if isinstance(j_kin_scaling_param_axes, list): + self._dim_scaling = len(j_kin_scaling_param_axes) + else: + self._dim_scaling = 1 + + def _kwargs2param_array(self, kwargs): + """ + converts dictionary to sorted array in same order as interpolation grid + + :param kwargs: dictionary of all model components, must include the one that are interpolated + :return: sorted list of parameters to interpolate + """ + param_array = [] + for param in self._param_list: + if param not in kwargs: + raise ValueError("key %s not in parameters and hence kinematic scaling not possible" % param) + param_array.append(kwargs.get(param)) + return param_array + + def param_bounds_interpol(self): + """ + minimum and maximum bounds of parameters that are being used to call interpolation function + + :return: dictionaries of minimum and maximum bounds + """ + kwargs_min, kwargs_max = {}, {} + if self._evaluate_scaling is True: + for i, key in enumerate(self._param_list): + kwargs_min[key] = min(self._param_arrays[i]) + kwargs_max[key] = max(self._param_arrays[i]) + return kwargs_min, kwargs_max + + def kin_scaling(self, kwargs_param): + """ + + :param kwargs_param: dictionary of parameters for scaling the kinematics + :return: scaling J(a_ani) for the IFU's + """ + if kwargs_param is None: + return np.ones(self._dim_scaling) + param_array = self._kwargs2param_array(kwargs_param) + if self._evaluate_scaling is not True or len(param_array) == 0: + return np.ones(self._dim_scaling) + scaling_list = [] + for scaling_class in self._j_scaling_ifu: + scaling = scaling_class.j_scaling(param_array) + scaling_list.append(scaling) + return np.array(scaling_list) diff --git a/hierarc/Likelihood/lens_sample_likelihood.py b/hierarc/Likelihood/lens_sample_likelihood.py index dedf24e9..d168b097 100644 --- a/hierarc/Likelihood/lens_sample_likelihood.py +++ b/hierarc/Likelihood/lens_sample_likelihood.py @@ -9,15 +9,17 @@ 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, los_distributions=None): + def __init__(self, kwargs_lens_list, normalized=False, kwargs_global_model=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 + :param kwargs_global_model: arguments of global distribution parameters for initialization of + ParamManager() class """ + if kwargs_global_model is None: + kwargs_global_model = {} self._lens_list = [] for kwargs_lens in kwargs_lens_list: if kwargs_lens["likelihood_type"] == "DSPL": @@ -27,12 +29,10 @@ def __init__(self, kwargs_lens_list, normalized=False, los_distributions=None): DSPLikelihood(normalized=normalized, **_kwargs_lens) ) else: + kwargs_lens_ = self._merge_global2local_settings(kwargs_global_model=kwargs_global_model, + kwargs_lens=kwargs_lens) self._lens_list.append( - LensLikelihood( - normalized=normalized, - los_distributions=los_distributions, - **kwargs_lens - ) + LensLikelihood(**kwargs_lens_) ) def log_likelihood( @@ -72,3 +72,37 @@ def num_data(self): for lens in self._lens_list: num += lens.num_data() return num + + @staticmethod + def _merge_global2local_settings(kwargs_global_model, kwargs_lens): + """ + + :param kwargs_global_model: dictionary of global model settings and distribution functions + :param kwargs_lens: specific settings of an individual lens + :return: joint dictionary that overwrites global with local parameters (if needed) and only keeps the relevant + arguments that an individual lens likelihood needs + :rtype: dict + """ + kwargs_global_model_ = {} + for key in _input_param_list: + if key in kwargs_global_model: + kwargs_global_model_[key] = kwargs_global_model[key] + kwargs_global_model_subset = copy.deepcopy(kwargs_global_model_) + return {**kwargs_global_model_subset, **kwargs_lens} + + +_input_param_list = ["anisotropy_model", + "anisotropy_sampling", + "anisotroy_distribution_function", + "los_distributions", + "lambda_mst_distribution", + "gamma_in_sampling", + "gamma_in_distribution", + "log_m2l_sampling", + "log_m2l_distribution", + "alpha_lambda_sampling", + "beta_lambda_sampling", + "alpha_gamma_in_sampling", + "alpha_log_m2l_sampling", + "log_scatter" + ] diff --git a/hierarc/Likelihood/parameter_scaling.py b/hierarc/Likelihood/parameter_scaling.py deleted file mode 100644 index 4739ce57..00000000 --- a/hierarc/Likelihood/parameter_scaling.py +++ /dev/null @@ -1,227 +0,0 @@ -__author__ = "sibirrer", "ajshajib" - -from scipy.interpolate import interp1d -from scipy.interpolate import interp2d -from scipy.interpolate import RegularGridInterpolator -import numpy as np - - -class ParameterScalingSingleAperture(object): - """Class to manage anisotropy scaling for single slit observation.""" - - def __init__(self, param_arrays, scaling_grid): - """ - - :param param_arrays: list of arrays of interpolated parameter values - :param scaling_grid: array with the scaling of J() for single slit - """ - self._evalute_scaling = False - # check if param arrays is 1d list or 2d list - if param_arrays is not None and scaling_grid is not None: - if isinstance(param_arrays, list): - self._dim_scaling = len(param_arrays) - else: - self._dim_scaling = 1 - - if self._dim_scaling == 1: - self._f_ani = interp1d(param_arrays, scaling_grid, kind="linear") - elif self._dim_scaling == 2: - self._f_ani = interp2d(param_arrays[0], param_arrays[1], scaling_grid.T) - else: - self._f_ani = RegularGridInterpolator( - tuple(param_arrays), - scaling_grid, - ) - self._evalute_scaling = True - - def param_scaling(self, param_array): - """ - - :param param_array: anisotropy parameter array - :return: scaling J(a_ani) for single slit - """ - if self._evalute_scaling is not True or param_array is None: - return 1 - if self._dim_scaling == 1: - return self._f_ani(param_array[0]) - elif self._dim_scaling == 2: - return self._f_ani(param_array[0], param_array[1])[0] - else: - return self._f_ani(param_array)[0] - - -class ParameterScalingIFU(object): - """Class to manage model parameter and anisotropy scalings for IFU data.""" - - def __init__( - self, anisotropy_model="NONE", param_arrays=None, scaling_grid_list=None - ): - """ - - :param anisotropy_model: string, either 'NONE', 'OM' or 'GOM' - :param param_arrays: array of parameter values - :param scaling_grid_list: list of array with the scalings of J() for each IFU - """ - self._anisotropy_model = anisotropy_model - self._evalute_ani = False - self._is_log_m2l_population_level = False - if ( - param_arrays is not None - and scaling_grid_list is not None - and self._anisotropy_model != "NONE" - ): - self._evalute_ani = True - self._anisotropy_scaling_list = [] - self._f_ani_list = [] - for scaling_grid in scaling_grid_list: - self._anisotropy_scaling_list.append( - ParameterScalingSingleAperture(param_arrays, scaling_grid) - ) - - if isinstance(param_arrays, list): - self._dim_scaling = len(param_arrays) - else: - self._dim_scaling = 1 - - if anisotropy_model in ["OM", "const"]: - if self._dim_scaling == 1: - self._ani_param_min = np.min(param_arrays) - self._ani_param_max = np.max(param_arrays) - else: - self._ani_param_min = np.min(param_arrays[0]) - self._ani_param_max = np.max(param_arrays[0]) - - if self._dim_scaling > 1: - self._gamma_in_min = np.min(param_arrays[1]) - self._gamma_in_max = np.max(param_arrays[1]) - if self._dim_scaling > 2: - self._log_m2l_min = np.min(param_arrays[2]) - self._log_m2l_max = np.max(param_arrays[2]) - self._is_log_m2l_population_level = True - - elif anisotropy_model == "GOM": - self._ani_param_min = [min(param_arrays[0]), min(param_arrays[1])] - self._ani_param_max = [max(param_arrays[0]), max(param_arrays[1])] - - if self._dim_scaling > 2: - self._gamma_in_min = np.min(param_arrays[2]) - self._gamma_in_max = np.max(param_arrays[2]) - if self._dim_scaling > 3: - self._log_m2l_min = np.min(param_arrays[3]) - self._log_m2l_max = np.max(param_arrays[3]) - self._is_log_m2l_population_level = True - else: - raise ValueError( - f"Anisotropy model {anisotropy_model} is not recognized!" - ) - - def param_scaling(self, param_array): - """ - - :param param_array: parameter array for scaling - :return: scaling J(a_ani) for the IFU's - """ - if self._evalute_ani is not True or param_array is None: - return [1] - scaling_list = [] - for scaling_class in self._anisotropy_scaling_list: - scaling = scaling_class.param_scaling(param_array) - scaling_list.append(scaling) - return np.array(scaling_list) - - def draw_anisotropy( - self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0 - ): - """Draw Gaussian distribution and re-sample if outside bounds. - - :param a_ani: mean of the distribution - :param a_ani_sigma: std of the distribution - :param beta_inf: anisotropy at infinity (relevant for GOM model) - :param beta_inf_sigma: std of beta_inf distribution - :return: random draw from the distribution - """ - if self._anisotropy_model in ["OM", "const"]: - if a_ani < self._ani_param_min or a_ani > self._ani_param_max: - raise ValueError( - "anisotropy parameter is out of bounds of the interpolated range!" - ) - # we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM - if self._anisotropy_model == "OM": - a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani) - else: - a_ani_draw = np.random.normal(a_ani, a_ani_sigma) - if a_ani_draw < self._ani_param_min or a_ani_draw > self._ani_param_max: - return self.draw_anisotropy(a_ani, a_ani_sigma) - return np.array([a_ani_draw]) - elif self._anisotropy_model in ["GOM"]: - if ( - a_ani < self._ani_param_min[0] - or a_ani > self._ani_param_max[0] - or beta_inf < self._ani_param_min[1] - or beta_inf > self._ani_param_max[1] - ): - raise ValueError( - "anisotropy parameter is out of bounds of the interpolated range!" - ) - a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani) - beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma) - if ( - a_ani_draw < self._ani_param_min[0] - or a_ani_draw > self._ani_param_max[0] - or beta_inf_draw < self._ani_param_min[1] - or beta_inf_draw > self._ani_param_max[1] - ): - return self.draw_anisotropy( - a_ani, a_ani_sigma, beta_inf, beta_inf_sigma - ) - return np.array([a_ani_draw, beta_inf_draw]) - return None - - def draw_lens_parameters( - self, gamma_in=None, gamma_in_sigma=0, log_m2l=None, log_m2l_sigma=0 - ): - """Draw Gaussian distribution and re-sample if outside bounds. - - :param gamma_in: mean of the distribution - :param gamma_in_sigma: std of the distribution - :param log_m2l: mean of the distribution - :param log_m2l_sigma: std of the distribution - :return: random draw from the distribution - """ - if self._is_log_m2l_population_level: - if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max: - raise ValueError( - "gamma_in parameter is out of bounds of the interpolated range!" - ) - if log_m2l < self._log_m2l_min or log_m2l > self._log_m2l_max: - raise ValueError( - "m2l parameter is out of bounds of the interpolated range!" - ) - - gamma_in_draw = np.random.normal(gamma_in, gamma_in_sigma) - log_m2l_draw = np.random.normal(log_m2l, log_m2l_sigma) - - if ( - gamma_in_draw < self._gamma_in_min - or gamma_in_draw > self._gamma_in_max - or log_m2l_draw < self._log_m2l_min - or log_m2l_draw > self._log_m2l_max - ): - return self.draw_lens_parameters( - gamma_in, gamma_in_sigma, log_m2l, log_m2l_sigma - ) - - return gamma_in_draw, log_m2l_draw - - else: - if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max: - raise ValueError( - "gamma_in parameter is out of bounds of the interpolated range!" - ) - - gamma_in_draw = np.random.normal(gamma_in, gamma_in_sigma) - - if gamma_in_draw < self._gamma_in_min or gamma_in_draw > self._gamma_in_max: - return self.draw_lens_parameters(gamma_in, gamma_in_sigma) - - return gamma_in_draw diff --git a/hierarc/Sampling/Distributions/__init__.py b/hierarc/Sampling/Distributions/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py new file mode 100644 index 00000000..c5f92a0a --- /dev/null +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -0,0 +1,78 @@ +import numpy as np + + +class AnisotropyDistribution(object): + """ + class to draw anisotropy parameters from hyperparameter distributions + """ + def __init__(self, anisotropy_model, anisotropy_sampling, + distribution_function, kwargs_anisotropy_min, kwargs_anisotropy_max): + """ + + :param anisotropy_model: string, name of anisotropy model to consider + :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens + kinematic prediction + :param distribution_function: string, 'NONE', 'GAUSSIAN', description of the distribution function of the + anisotropy model parameters + """ + + self._anisotropy_model = anisotropy_model + self._anisotropy_sampling = anisotropy_sampling + self._distribution_function = distribution_function + if kwargs_anisotropy_min is None: + kwargs_anisotropy_min = {} + if kwargs_anisotropy_max is None: + kwargs_anisotropy_max = {} + self._kwargs_min = kwargs_anisotropy_min + self._kwargs_max = kwargs_anisotropy_max + self._a_ani_min, self._a_ani_max = self._kwargs_min.get("a_ani", -np.inf), self._kwargs_max.get("a_ani", np.inf) + self._beta_inf_min = self._kwargs_min.get("beta_inf", -np.inf) + self._beta_inf_max = self._kwargs_max.get("beta_inf", np.inf) + + def draw_anisotropy(self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0): + """Draw Gaussian distribution and re-sample if outside bounds. + + :param a_ani: mean of the distribution + :param a_ani_sigma: std of the distribution + :param beta_inf: anisotropy at infinity (relevant for GOM model) + :param beta_inf_sigma: std of beta_inf distribution + :return: random draw from the distribution + """ + kwargs_return = {} + if not self._anisotropy_sampling: + if a_ani is not None: + kwargs_return["a_ani"] = a_ani + if beta_inf is not None: + kwargs_return["beta_inf"] = beta_inf + return kwargs_return + if self._anisotropy_model in ["OM", "const", "GOM"]: + if a_ani < self._a_ani_min or a_ani > self._a_ani_max: + raise ValueError( + "anisotropy parameter is out of bounds of the interpolated range!" + ) + # we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM + if self._distribution_function in ["GAUSSIAN"]: + if self._anisotropy_model == "OM": + a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani) + else: + a_ani_draw = np.random.normal(a_ani, a_ani_sigma) + if a_ani_draw < self._a_ani_min or a_ani_draw > self._a_ani_max: + return self.draw_anisotropy(a_ani, a_ani_sigma) + kwargs_return["a_ani"] = a_ani_draw + else: + kwargs_return["a_ani"] = a_ani + if self._anisotropy_model in ["GOM"]: + if beta_inf < self._beta_inf_min or beta_inf > self._beta_inf_max: + raise ValueError( + "anisotropy parameter is out of bounds of the interpolated range!" + ) + if self._distribution_function in ["GAUSSIAN"]: + beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma) + else: + beta_inf_draw = beta_inf + if beta_inf_draw < self._beta_inf_min or beta_inf_draw > self._beta_inf_max: + return self.draw_anisotropy( + a_ani, a_ani_sigma, beta_inf, beta_inf_sigma + ) + kwargs_return["beta_inf"] = beta_inf_draw + return kwargs_return diff --git a/hierarc/Sampling/Distributions/lens_distribution.py b/hierarc/Sampling/Distributions/lens_distribution.py new file mode 100644 index 00000000..dd1c2738 --- /dev/null +++ b/hierarc/Sampling/Distributions/lens_distribution.py @@ -0,0 +1,189 @@ +import numpy as np + + +class LensDistribution(object): + """ + class to draw lens parameters of individual lens from distributions + """ + + def __init__( + self, + lambda_mst_sampling=False, + lambda_mst_distribution="NONE", + gamma_in_sampling=False, + gamma_in_distribution="NONE", + log_m2l_sampling=False, + log_m2l_distribution="NONE", + alpha_lambda_sampling=False, + beta_lambda_sampling=False, + alpha_gamma_in_sampling=False, + alpha_log_m2l_sampling=False, + log_scatter=False, + mst_ifu=False, + lambda_scaling_property=0, + lambda_scaling_property_beta=0, + kwargs_min=None, + kwargs_max=None, + ): + """ + + :param lambda_mst_sampling: bool, if True adds a global mass-sheet transform parameter in the sampling + :param lambda_mst_distribution: string, distribution function of the MST transform + :param gamma_in_sampling: bool, if True samples the inner slope of the GNFW profile + :param gamma_in_distribution: string, distribution function of the inner + slope of the GNFW profile + :param log_m2l_sampling: bool, if True samples the mass to light ratio of + the stars in logarithmic scale + :param log_m2l_distribution: string, distribution function of the logarithm of mass to + light ratio of the lens + :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly + according to a predefined quantity of the lens + :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly + according to a predefined quantity of the lens + :param alpha_gamma_in_sampling: bool, if True samples a parameter alpha_gamma_in, which scales gamma_in linearly + :param alpha_log_m2l_sampling: bool, if True samples a parameter alpha_log_m2l, which scales log_m2l linearly + :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space + (and thus flat prior in log) + :param mst_ifu: bool, if True replaces the lambda_mst parameter by the lambda_ifu parameter (and distribution) + in sampling this lens. + :param lambda_scaling_property: float (optional), scaling of + lambda_mst = lambda_mst_global + alpha * lambda_scaling_property + :param lambda_scaling_property_beta: float (optional), scaling of + lambda_mst = lambda_mst_global + beta * lambda_scaling_property_beta + :param kwargs_min: minimum arguments of parameters supported by each lens + :type kwargs_min: dict or None + :param kwargs_max: maximum arguments of parameters supported by each lens + :type kwargs_max: dict or None + """ + self._lambda_mst_sampling = lambda_mst_sampling + self._lambda_mst_distribution = lambda_mst_distribution + self._gamma_in_sampling = gamma_in_sampling + self._gamma_in_distribution = gamma_in_distribution + self._log_m2l_sampling = log_m2l_sampling + self._log_m2l_distribution = log_m2l_distribution + self._alpha_lambda_sampling = alpha_lambda_sampling + self._beta_lambda_sampling = beta_lambda_sampling + self._alpha_gamma_in_sampling = alpha_gamma_in_sampling + self._alpha_log_m2l_sampling = alpha_log_m2l_sampling + self._mst_ifu = mst_ifu + self._lambda_scaling_property = lambda_scaling_property + self._lambda_scaling_property_beta = lambda_scaling_property_beta + + self._log_scatter = log_scatter + if kwargs_max is None: + kwargs_max = {} + if kwargs_min is None: + kwargs_min = {} + self._gamma_in_min, self._gamma_in_max = kwargs_min.get("gamma_in", -np.inf), kwargs_max.get("gamma_in", np.inf) + self._log_m2l_min, self._log_m2l_max = kwargs_min.get("log_m2l", -np.inf), kwargs_max.get("log_m2l", np.inf) + + def draw_lens( + self, + lambda_mst=1, + lambda_mst_sigma=0, + gamma_ppn=1, + lambda_ifu=1, + lambda_ifu_sigma=0, + alpha_lambda=0, + beta_lambda=0, + gamma_in=1, + gamma_in_sigma=0, + alpha_gamma_in=0, + log_m2l=1, + log_m2l_sigma=0, + alpha_log_m2l=0, + ): + """Draws a realization of a specific model from the hyperparameter + distribution. + + :param lambda_mst: MST transform + :param lambda_mst_sigma: spread in the distribution + :param gamma_ppn: Post-Newtonian parameter + :param lambda_ifu: secondary lambda_mst parameter for subset of lenses specified + for + :param lambda_ifu_sigma: secondary lambda_mst_sigma parameter for subset of + lenses specified for + :param alpha_lambda: float, linear slope of the lambda_int scaling relation with + lens quantity self._lambda_scaling_property + :param beta_lambda: float, a second linear slope of the lambda_int scaling + relation with lens quantity self._lambda_scaling_property_beta + :param gamma_in: inner slope of the NFW profile + :param gamma_in_sigma: spread in the distribution + :param alpha_gamma_in: float, linear slope of the gamma_in scaling relation with + lens quantity self._lambda_scaling_property + :param log_m2l: log(mass-to-light ratio) + :param log_m2l_sigma: spread in the distribution + :param alpha_log_m2l: float, linear slope of the log(m2l) scaling relation with + lens quantity self._lambda_scaling_property + :return: draw from the distributions + """ + kwargs_return = {} + + if self._mst_ifu is True: + lambda_mst_mean_lens = lambda_ifu + else: + lambda_mst_mean_lens = lambda_mst + + lambda_lens = (lambda_mst_mean_lens + + alpha_lambda * self._lambda_scaling_property + + beta_lambda * self._lambda_scaling_property_beta + ) + lambda_mst_draw = lambda_lens + if self._lambda_mst_sampling: + if self._lambda_mst_distribution in ["GAUSSIAN"]: + lambda_mst_draw = np.random.normal(lambda_lens, lambda_ifu_sigma) + + kwargs_return["lambda_mst"] = lambda_mst_draw + kwargs_return["gamma_ppn"] = gamma_ppn + + if self._gamma_in_sampling: + if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max: + raise ValueError( + "gamma_in parameter is out of bounds of the interpolated range!" + ) + if self._gamma_in_distribution in ["GAUSSIAN"]: + gamma_in_lens = gamma_in + alpha_gamma_in * self._lambda_scaling_property + else: + gamma_in_lens = gamma_in + gamma_in_draw = np.random.normal(gamma_in_lens, gamma_in_sigma) + if gamma_in_draw < self._gamma_in_min or gamma_in_draw > self._gamma_in_max: + return self.draw_lens(lambda_mst=lambda_mst, + lambda_mst_sigma=lambda_mst_sigma, + gamma_ppn=gamma_ppn, + lambda_ifu=lambda_ifu, + lambda_ifu_sigma=lambda_ifu_sigma, + alpha_lambda=alpha_lambda, + beta_lambda=beta_lambda, + gamma_in=gamma_in, + gamma_in_sigma=gamma_in_sigma, + alpha_gamma_in=alpha_gamma_in, + log_m2l=log_m2l, + log_m2l_sigma=log_m2l_sigma, + alpha_log_m2l=alpha_log_m2l) + kwargs_return["gamma_in"] = gamma_in_draw + if self._log_m2l_sampling: + + if log_m2l < self._log_m2l_min or log_m2l > self._log_m2l_max: + raise ValueError( + "m2l parameter is out of bounds of the interpolated range!" + ) + + log_m2l_lens = log_m2l + alpha_log_m2l * self._lambda_scaling_property + log_m2l_draw = np.random.normal(log_m2l_lens, log_m2l_sigma) + + if log_m2l_draw < self._log_m2l_min or log_m2l_draw > self._log_m2l_max: + return self.draw_lens(lambda_mst=lambda_mst, + lambda_mst_sigma=lambda_mst_sigma, + gamma_ppn=gamma_ppn, + lambda_ifu=lambda_ifu, + lambda_ifu_sigma=lambda_ifu_sigma, + alpha_lambda=alpha_lambda, + beta_lambda=beta_lambda, + gamma_in=gamma_in, + gamma_in_sigma=gamma_in_sigma, + alpha_gamma_in=alpha_gamma_in, + log_m2l=log_m2l, + log_m2l_sigma=log_m2l_sigma, + alpha_log_m2l=alpha_log_m2l) + kwargs_return["log_m2l"] = log_m2l_draw + return kwargs_return diff --git a/hierarc/Likelihood/los_distributions.py b/hierarc/Sampling/Distributions/los_distributions.py similarity index 100% rename from hierarc/Likelihood/los_distributions.py rename to hierarc/Sampling/Distributions/los_distributions.py diff --git a/hierarc/Sampling/ParamManager/param_manager.py b/hierarc/Sampling/ParamManager/param_manager.py index 11413773..8b87363d 100644 --- a/hierarc/Sampling/ParamManager/param_manager.py +++ b/hierarc/Sampling/ParamManager/param_manager.py @@ -64,7 +64,7 @@ def __init__( according to a predefined quantity of the lens :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens - kinematic prediction + kinematic prediction :param anisotropy_distribution: string, indicating the distribution function of the anisotropy model :param gamma_in_sampling: bool, if True samples gNFW inner slope parameter :param gamma_in_distribution: string, distribution function of the gamma_in parameter diff --git a/test/test_LensPosterior/test_anisotropy_config.py b/test/test_LensPosterior/test_anisotropy_config.py index 040ed4c2..ed668744 100644 --- a/test/test_LensPosterior/test_anisotropy_config.py +++ b/test/test_LensPosterior/test_anisotropy_config.py @@ -1,4 +1,4 @@ -from hierarc.LensPosterior.anisotropy_config import AnisotropyConfig +from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig import unittest import pytest @@ -6,9 +6,9 @@ class TestAnisotropyConfig(object): def setup_method(self): self.r_eff = 2 - self.config_om = AnisotropyConfig(anisotropy_model="OM", r_eff=self.r_eff) - self.config_gom = AnisotropyConfig(anisotropy_model="GOM", r_eff=self.r_eff) - self.config_const = AnisotropyConfig(anisotropy_model="const", r_eff=self.r_eff) + self.config_om = KinScalingConfig(anisotropy_model="OM", r_eff=self.r_eff) + self.config_gom = KinScalingConfig(anisotropy_model="GOM", r_eff=self.r_eff) + self.config_const = KinScalingConfig(anisotropy_model="const", r_eff=self.r_eff) def test_kwargs_anisotropy_base(self): kwargs = self.config_om.kwargs_anisotropy_base @@ -22,15 +22,15 @@ def test_kwargs_anisotropy_base(self): assert kwargs["beta"] == 0.1 def test_ani_param_array(self): - ani_param_array = self.config_om.ani_param_array - assert len(ani_param_array) == 6 + ani_param_array = self.config_om.kin_scaling_param_array + assert len(ani_param_array[0]) == 6 - ani_param_array = self.config_gom.ani_param_array + ani_param_array = self.config_gom.kin_scaling_param_array assert len(ani_param_array[0]) == 6 assert len(ani_param_array[1]) == 4 - ani_param_array = self.config_const.ani_param_array - assert len(ani_param_array) == 7 + ani_param_array = self.config_const.kin_scaling_param_array + assert len(ani_param_array[0]) == 7 def test_anisotropy_kwargs(self): a_ani = 2 @@ -49,15 +49,15 @@ def test_anisotropy_kwargs(self): class TestRaise(unittest.TestCase): def test_raise(self): with self.assertRaises(ValueError): - AnisotropyConfig(anisotropy_model="BAD", r_eff=1) + KinScalingConfig(anisotropy_model="BAD", r_eff=1) with self.assertRaises(ValueError): - conf = AnisotropyConfig(anisotropy_model="OM", r_eff=1) + conf = KinScalingConfig(anisotropy_model="OM", r_eff=1) conf._anisotropy_model = "BAD" kwargs = conf.kwargs_anisotropy_base with self.assertRaises(ValueError): - conf = AnisotropyConfig(anisotropy_model="OM", r_eff=1) + conf = KinScalingConfig(anisotropy_model="OM", r_eff=1) conf._anisotropy_model = "BAD" kwargs = conf.anisotropy_kwargs(a_ani=1, beta_inf=1) diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py index 0ff86961..f99c6168 100644 --- a/test/test_LensPosterior/test_kin_constraints_composite.py +++ b/test/test_LensPosterior/test_kin_constraints_composite.py @@ -104,9 +104,9 @@ def test_likelihoodconfiguration_om(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(**kwargs_likelihood) + ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood) kwargs_kin = {"a_ani": 1} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin) kwargs_lens_light_test = [{"amp": [1, 1], "sigma": [1, 2]}] lens_light_model_list_test = ["MULTI_GAUSSIAN"] @@ -135,8 +135,8 @@ def test_likelihoodconfiguration_om(self): **kwargs_kin_api_settings ) - kappa_s_array = 10 ** np.random.normal(8, 0, 100) / 1e6 - r_s_angle_array = np.random.normal(0.1, 0, 100) + kappa_s_array = 10 ** np.random.normal(8, 0, 10) / 1e6 + r_s_angle_array = np.random.normal(0.1, 0, 10) kin_constraints_kappa = KinConstraintsComposite( z_lens=z_lens, @@ -224,8 +224,8 @@ def test_likelihoodconfiguration_gom(self): gamma_in_array = np.linspace(0.1, 2.9, 5) log_m2l_array = np.linspace(0.1, 1, 5) - rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0, 100) + rho0_array = 10 ** np.random.normal(8, 0, 10) / 1e6 + r_s_array = np.random.normal(0.1, 0, 10) # compute likelihood kin_constraints = KinConstraintsComposite( @@ -256,9 +256,9 @@ def test_likelihoodconfiguration_gom(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(**kwargs_likelihood) + ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood) kwargs_kin = {"a_ani": 1, "beta_inf": 0.5} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin) class TestKinConstraintsCompositeM2l(object): @@ -324,9 +324,9 @@ def test_likelihoodconfiguration_om(self): gamma_in_array = np.linspace(0.1, 2.9, 5) - log_m2l_array = np.random.uniform(0.1, 1, 100) - rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0, 100) + log_m2l_array = np.random.uniform(0.1, 1, 5) + rho0_array = 10 ** np.random.normal(8, 0, 5) / 1e6 + r_s_array = np.random.normal(0.1, 0, 5) # compute likelihood kin_constraints = KinConstraintsComposite( @@ -360,9 +360,9 @@ def test_likelihoodconfiguration_om(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(**kwargs_likelihood) + ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=False, **kwargs_likelihood) kwargs_kin = {"a_ani": 1} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 0}, kwargs_kin=kwargs_kin) def test_likelihoodconfiguration_gom(self): anisotropy_model = "GOM" @@ -423,9 +423,9 @@ def test_likelihoodconfiguration_gom(self): gamma_in_array = np.linspace(0.1, 2.9, 5) - log_m2l_array = np.random.uniform(0.1, 1, 100) - rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0, 100) + log_m2l_array = np.random.uniform(0.1, 1, 10) + rho0_array = 10 ** np.random.normal(8, 0, 10) / 1e6 + r_s_array = np.random.normal(0.1, 0, 10) # compute likelihood kin_constraints = KinConstraintsComposite( @@ -455,9 +455,9 @@ def test_likelihoodconfiguration_gom(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(**kwargs_likelihood) + ln_class = LensLikelihood(gamma_in_sampling=True, **kwargs_likelihood) kwargs_kin = {"a_ani": 1, "beta_inf": 0.5} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2}, kwargs_kin=kwargs_kin) class TestRaise(unittest.TestCase): @@ -510,8 +510,8 @@ def test_raise(self): gamma_in_array = np.linspace(0.1, 2.9, 5) log_m2l_array = np.linspace(0.1, 1, 5) - rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0.01, 100) + rho0_array = 10 ** np.random.normal(8, 0.2, 10) / 1e6 + r_s_array = np.random.normal(0.1, 0.01, 10) kin_constraints = KinConstraintsComposite( z_lens=z_lens, @@ -737,9 +737,9 @@ def test_raise_m2l(self): kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}] gamma_in_array = np.linspace(0.1, 2.9, 5) - log_m2l_array = np.random.uniform(0.1, 1, 100) - rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0.01, 100) + log_m2l_array = np.random.uniform(0.1, 1, 10) + rho0_array = 10 ** np.random.normal(8, 0.2, 10) / 1e6 + r_s_array = np.random.normal(0.1, 0.01, 10) kin_constraints = KinConstraintsComposite( z_lens=z_lens, @@ -890,9 +890,9 @@ def test_raise_m2l(self): kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}] gamma_in_array = np.linspace(0.1, 2.9, 5) - log_m2l_array = np.linspace(0.1, 1, 100) - rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6 - r_s_array = np.random.normal(0.1, 0.01, 100) + log_m2l_array = np.linspace(0.1, 1, 10) + rho0_array = 10 ** np.random.normal(8, 0.2, 10) / 1e6 + r_s_array = np.random.normal(0.1, 0.01, 10) kin_constraints = KinConstraintsComposite( z_lens=z_lens, diff --git a/test/test_Likelihood/test_cosmo_likelihood.py b/test/test_Likelihood/test_cosmo_likelihood.py index 6035cbae..aa863637 100644 --- a/test/test_Likelihood/test_cosmo_likelihood.py +++ b/test/test_Likelihood/test_cosmo_likelihood.py @@ -57,6 +57,11 @@ def setup_method(self): "kwargs_lower_cosmo": kwargs_lower_cosmo, "kwargs_upper_cosmo": kwargs_upper_cosmo, } + self.kwargs_model = {"ppn_sampling": False, + "lambda_mst_sampling": False, + "lambda_mst_distribution": "delta", + "anisotropy_sampling": False, + "anisotropy_model": "OM",} # self.kwargs_likelihood_list = [{'z_lens': self.z_L, 'z_source': self.z_S, 'likelihood_type': 'TDKinKDE', # 'dd_sample': self.dd_samples, 'ddt_sample': self.ddt_samples, @@ -66,12 +71,8 @@ def test_log_likelihood(self): cosmoL = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, - ppn_sampling=False, - lambda_mst_sampling=False, - lambda_mst_distribution="delta", - anisotropy_sampling=False, - anisotropy_model="OM", custom_prior=None, interpolate_cosmo=True, num_redshift_interp=100, @@ -84,12 +85,8 @@ def custom_prior(kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_lo cosmoL_prior = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, - ppn_sampling=False, - lambda_mst_sampling=False, - lambda_mst_distribution="delta", - anisotropy_sampling=False, - anisotropy_model="OM", custom_prior=custom_prior, interpolate_cosmo=True, num_redshift_interp=100, @@ -128,6 +125,7 @@ def test_cosmo_instance(self): cosmoL = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, interpolate_cosmo=False, cosmo_fixed=None, @@ -137,6 +135,7 @@ def test_cosmo_instance(self): cosmoL = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, interpolate_cosmo=True, num_redshift_interp=100, @@ -147,6 +146,7 @@ def test_cosmo_instance(self): cosmoL = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, interpolate_cosmo=True, num_redshift_interp=100, @@ -158,6 +158,7 @@ def test_cosmo_instance(self): cosmoL = CosmoLikelihood( self.kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, interpolate_cosmo=False, num_redshift_interp=100, @@ -209,6 +210,7 @@ def test_oLCDM_init(self): cosmoL = CosmoLikelihood( kwargs_likelihood_list, self.cosmology, + self.kwargs_model, self.kwargs_bounds, interpolate_cosmo=False, cosmo_fixed=None, @@ -218,6 +220,7 @@ def test_sne_likelihood_integration(self): cosmoL = CosmoLikelihood( [], self.cosmology, + self.kwargs_model, self.kwargs_bounds, sne_likelihood="Pantheon_binned", interpolate_cosmo=True, @@ -238,8 +241,8 @@ def test_kde_likelihood_integration(self): self.cosmology, rescale=True, ) - cosmoL = CosmoLikelihood( - [], "FLCDM", self.kwargs_bounds, KDE_likelihood_chain=chain + cosmoL = CosmoLikelihood([], "FLCDM", self.kwargs_model, + self.kwargs_bounds, KDE_likelihood_chain=chain ) kwargs_cosmo = {"h0": self.H0_true, "om": self.omega_m_true, "ok": 0} args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs_cosmo) diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py index a150d70c..f17c828a 100644 --- a/test/test_Likelihood/test_hierarchy_likelihood.py +++ b/test/test_Likelihood/test_hierarchy_likelihood.py @@ -24,56 +24,67 @@ def setup_method(self): "dd_mean": dd, "dd_sigma": dd / 10, } + kwargs_model = {"anisotropy_model": "OM", + "anisotropy_sampling": True, + "anisotroy_distribution_function": "GAUSSIAN", + "lambda_mst_distribution": "GAUSSIAN", + } + # "gamma_in_sampling" = False, + gamma_in_distribution = "NONE", + log_m2l_sampling = False, + log_m2l_distribution = "NONE", + + self.likelihood = LensLikelihood( z_lens, z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - ani_scaling_array=ani_scaling_array, + kin_scaling_param_list=["a_ani"], + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[ani_scaling_array], num_distribution_draws=200, los_distributions=["GAUSSIAN"], global_los_distribution=0, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=True, - **kwargs_likelihood + **kwargs_likelihood, + **kwargs_model ) self.likelihood_single = LensLikelihood( z_lens, z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - ani_scaling_array=ani_scaling_array, + kin_scaling_param_list=["a_ani"], + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[ani_scaling_array], num_distribution_draws=200, los_distributions=["GAUSSIAN"], global_los_distribution=0, # testing previously set to =False kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, - **kwargs_likelihood + **kwargs_likelihood, + **kwargs_model ) self.likelihood_zero_dist = LensLikelihood( z_lens, z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - ani_scaling_array=ani_scaling_array, + kin_scaling_param_list=["a_ani"], + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[ani_scaling_array], num_distribution_draws=0, los_distributions=["GAUSSIAN"], global_los_distribution=0, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=True, - **kwargs_likelihood + **kwargs_likelihood, + **kwargs_model ) kappa_posterior = np.random.normal(loc=0, scale=0.03, size=100000) @@ -83,17 +94,17 @@ def setup_method(self): z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - ani_scaling_array=ani_scaling_array, + kin_scaling_param_list=["a_ani"], + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[ani_scaling_array], num_distribution_draws=200, # los_distributions=["GAUSSIAN"], global_los_distribution=False, kappa_pdf=kappa_pdf, kappa_bin_edges=kappa_bin_edges, mst_ifu=False, - **kwargs_likelihood + **kwargs_likelihood, + **kwargs_model ) gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) @@ -102,40 +113,49 @@ def setup_method(self): np.ones_like(ani_param_array), np.outer(np.ones_like(gamma_in_array), np.ones_like(log_m2l_array)), ) + j_kin_scaling_param_axes = [ani_param_array, gamma_in_array, log_m2l_array] self.likelihood_gamma_in_log_m2l_list_ani = LensLikelihood( z_lens, z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=[ani_param_array], - ani_scaling_array_list=None, - param_scaling_grid_list=[param_scaling_array], - gamma_in_array=gamma_in_array, - log_m2l_array=log_m2l_array, + kin_scaling_param_list=["a_ani", "gamma_in", "log_m2l"], + j_kin_scaling_param_axes=j_kin_scaling_param_axes, + j_kin_scaling_grid_list=[param_scaling_array], num_distribution_draws=200, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, - **kwargs_likelihood + gamma_in_sampling=False, + gamma_in_distribution="GAUSSIAN", + log_m2l_sampling=True, + log_m2l_distribution="GAUSSIAN", + **kwargs_likelihood, + **kwargs_model ) + param_scaling_array = np.outer(np.ones_like(gamma_in_array), np.ones_like(log_m2l_array)) + + j_kin_scaling_param_axes = [gamma_in_array, log_m2l_array] + self.likelihood_gamma_in_log_m2l = LensLikelihood( z_lens, z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - param_scaling_grid_list=[param_scaling_array], - gamma_in_array=gamma_in_array, - log_m2l_array=log_m2l_array, + kin_scaling_param_list=["gamma_in", "log_m2l"], + j_kin_scaling_param_axes=j_kin_scaling_param_axes, + j_kin_scaling_grid_list=[param_scaling_array], num_distribution_draws=200, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, - **kwargs_likelihood + gamma_in_sampling=True, + gamma_in_distribution="GAUSSIAN", + log_m2l_sampling=True, + log_m2l_distribution="GAUSSIAN", + **kwargs_likelihood, + **kwargs_model # TODO: remove anisotropy sampling in that scenario? ) self.likelihood_gamma_in_fail_case = LensLikelihood( @@ -143,18 +163,16 @@ def setup_method(self): z_source, name="name", likelihood_type="DdtDdGaussian", - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=None, - param_scaling_grid_list=[param_scaling_array], - gamma_in_array=gamma_in_array, - log_m2l_array=log_m2l_array, + kin_scaling_param_list=["a_ani"], + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[ani_scaling_array], num_distribution_draws=200, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False, lambda_scaling_property=100, - **kwargs_likelihood + **kwargs_likelihood, + **kwargs_model ) def test_lens_log_likelihood(self): @@ -221,8 +239,6 @@ def test_lens_log_likelihood(self): kwargs_test = self.likelihood._kwargs_init(kwargs=None) assert type(kwargs_test) is dict - gamma_in_draw = self.likelihood.draw_lens_scaling_params() - assert gamma_in_draw is None kwargs_lens = { "gamma_in": 1, @@ -260,11 +276,11 @@ def test_lens_log_likelihood(self): dds = self.cosmo.angular_diameter_distance_z1z2(z1=z_lens, z2=z_source).value 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, kwargs_los - ) + #ln_likelihood = self.likelihood_gamma_in_fail_case.log_likelihood_single( + # ddt, dd, delta_lum_dist, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los + #) - assert ln_likelihood < -10000000 + #assert ln_likelihood < -10000000 if __name__ == "__main__": diff --git a/test/test_Likelihood/test_kin_scaling.py b/test/test_Likelihood/test_kin_scaling.py new file mode 100644 index 00000000..f66aa32f --- /dev/null +++ b/test/test_Likelihood/test_kin_scaling.py @@ -0,0 +1,269 @@ +import numpy as np +import numpy.testing as npt +import pytest + +from hierarc.Likelihood.kin_scaling import KinScaling, ParameterScalingSingleMeasurement + + +class TestKinScaling(object): + + def test_single_param(self): + param_arrays = np.linspace(0, 1, 11) + scaling_grid_list = [param_arrays**2] + param_list = ["a"] + kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list) + kwargs_param = {"a": 0.5} + j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) + npt.assert_almost_equal(j_scaling, 0.5**2, decimal=2) + kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + assert kwargs_min["a"] == 0 + + param_arrays = np.linspace(0, 1, 11) + scaling_grid_list = [param_arrays ** 2] + param_list = ["a"] + kin_scaling = KinScaling(j_kin_scaling_param_axes=[param_arrays], + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list) + kwargs_param = {"a": 0.5} + j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) + npt.assert_almost_equal(j_scaling, 0.5 ** 2, decimal=2) + kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + assert kwargs_min["a"] == 0 + + def test_two_parameters(self): + param_arrays = [np.linspace(0, 1, 11), np.linspace(0, 2, 21)] + xy, uv = np.meshgrid(param_arrays[0], param_arrays[1]) + scaling_grid_list = [xy.T * uv.T, xy.T, uv.T] + shape_scaling = np.shape(scaling_grid_list[0]) + assert shape_scaling[0] == 11 + assert shape_scaling[1] == 21 + # assert 1 == 0 + param_list = ["a", "b"] + kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list) + kwargs_param = {"a": 0.5, "b": 0.3} + j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) + print(j_scaling) + npt.assert_almost_equal(j_scaling[0], 0.5 * 0.3, decimal=2) + npt.assert_almost_equal(j_scaling[1], 0.5, decimal=2) + npt.assert_almost_equal(j_scaling[2], 0.3, decimal=2) + kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + assert kwargs_min["a"] == 0 + + def test__kwargs2param_array(self): + param_arrays = [np.linspace(0, 1, 11), np.linspace(0, 2, 11)] + xy, uv = np.meshgrid(param_arrays[0], param_arrays[1]) + scaling_grid_list = [xy * uv, xy, uv] + + param_list = ["a", "b"] + kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list) + kwargs_param = {"a": 0.5, "b": 0.3} + param_array = kin_scaling._kwargs2param_array(kwargs_param) + assert param_array[0] == kwargs_param["a"] + assert param_array[1] == kwargs_param["b"] + kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + assert kwargs_min["a"] == 0 + assert kwargs_min["b"] == 0 + assert kwargs_max["b"] == 2 + assert kwargs_max["a"] == 1 + + def test_empty(self): + kin_scaling = KinScaling(j_kin_scaling_param_axes=None, j_kin_scaling_grid_list=None, j_kin_scaling_param_name_list=None) + output = kin_scaling.kin_scaling(kwargs_param=None) + assert output == 1 + kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + +class TestParameterScalingSingleAperture(object): + def setup_method(self): + ani_param_array = np.linspace(start=0, stop=1, num=10) + param_scaling_array = ani_param_array * 2 + self.scaling = ParameterScalingSingleMeasurement( + ani_param_array, param_scaling_array + ) + + ani_param_array = [ + np.linspace(start=0, stop=1, num=10), + np.linspace(start=1, stop=2, num=5), + ] + param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) + print(np.shape(param_scaling_array), 'test shape') + self.scaling_2d = ParameterScalingSingleMeasurement( + ani_param_array, param_scaling_array + ) + + ani_param_array = np.linspace(start=0, stop=1, num=10) + gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) + log_m2l_array = np.linspace(start=0.1, stop=1, num=10) + + param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.outer(gamma_in_array, log_m2l_array), + ) + self.scaling_nfw = ParameterScalingSingleMeasurement( + param_arrays, param_scaling_array + ) + + gom_param_array = [ + np.linspace(start=0, stop=1, num=10), + np.linspace(start=1, stop=2, num=5), + ] + param_arrays = [ + gom_param_array[0], + gom_param_array[1], + gamma_in_array, + log_m2l_array, + ] + param_scaling_array = np.multiply.outer( + gom_param_array[0], + np.multiply.outer( + gom_param_array[1], + np.outer(gamma_in_array, log_m2l_array), + ), + ) + + self.scaling_nfw_2d = ParameterScalingSingleMeasurement( + param_arrays, param_scaling_array + ) + + def test_param_scaling(self): + scaling = self.scaling.j_scaling(param_array=[1]) + assert scaling == np.array([2]) + + scaling = self.scaling.j_scaling(param_array=[]) + assert scaling == 1 + + scaling = self.scaling_2d.j_scaling(param_array=[1, 2]) + assert scaling == 2 + + scaling = self.scaling_nfw.j_scaling(param_array=[1, 2.9, 0.5]) + assert scaling == 1 * 2.9 * 0.5 + + scaling = self.scaling_nfw_2d.j_scaling(param_array=[1, 2, 2.9, 0.5]) + assert scaling == 1 * 2 * 2.9 * 0.5 + + +class TestParameterScalingIFU(object): + def setup_method(self): + ani_param_array = np.linspace(start=0, stop=1, num=10) + param_scaling_array = ani_param_array * 2 + self.scaling = KinScaling( + j_kin_scaling_param_axes=ani_param_array, j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a"] + ) + + ani_param_array = [ + np.linspace(start=0, stop=1, num=10), + np.linspace(start=1, stop=2, num=5), + ] + param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) + self.scaling_2d = KinScaling(j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b"] + ) + + ani_param_array = np.linspace(start=0, stop=1, num=10) + gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) + log_m2l_array = np.linspace(start=0.1, stop=1, num=10) + + param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.outer(gamma_in_array, log_m2l_array), + ) + self.scaling_nfw = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c"] + ) + + + gom_param_array = [ + np.linspace(start=0, stop=1, num=10), + np.linspace(start=1, stop=2, num=5), + ] + param_arrays = [ + gom_param_array[0], + gom_param_array[1], + gamma_in_array, + log_m2l_array, + ] + param_scaling_array = np.multiply.outer( + gom_param_array[0], + np.multiply.outer( + gom_param_array[1], + np.outer(gamma_in_array, log_m2l_array), + ), + ) + self.scaling_nfw_2d = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c", "d"] + ) + + + param_arrays = [ani_param_array, gamma_in_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + gamma_in_array, + ) + self.scaling_nfw_no_m2l = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b"] + ) + + gom_param_array = [ + np.linspace(start=0, stop=1, num=10), + np.linspace(start=1, stop=2, num=5), + ] + param_arrays = [ + gom_param_array[0], + gom_param_array[1], + gamma_in_array, + ] + param_scaling_array = np.multiply.outer( + gom_param_array[0], + np.multiply.outer( + gom_param_array[1], + gamma_in_array, + ), + ) + self.scaling_nfw_2d_no_m2l = KinScaling(j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c"] + ) + + def test_kin_scaling(self): + + scaling = self.scaling.kin_scaling(kwargs_param=None) + assert scaling[0] == 1 + + scaling = self.scaling.kin_scaling(kwargs_param={"a":1}) + assert scaling[0] == 2 + + kwargs_param = {"a": 1, "b": 2} + scaling = self.scaling_2d.kin_scaling(kwargs_param=kwargs_param) + assert scaling[0] == 2 + + kwargs_param = {"a": 1, "b": 2.9, "c": 0.5} + scaling = self.scaling_nfw.kin_scaling(kwargs_param=kwargs_param) + assert scaling[0] == 1 * 2.9 * 0.5 + + kwargs_param = {"a": 1, "b": 2., "c": 2.9, "d": 0.5} + scaling = self.scaling_nfw_2d.kin_scaling(kwargs_param=kwargs_param) + assert scaling[0] == 1 * 2 * 2.9 * 0.5 + + kwargs_param = {"a": 1, "b": 2.9} + scaling = self.scaling_nfw_no_m2l.kin_scaling(kwargs_param=kwargs_param) + assert scaling[0] == 1 * 2.9 + + kwargs_param = {"a": 1, "b": 2, "c": 2.9} + scaling = self.scaling_nfw_2d_no_m2l.kin_scaling(kwargs_param=kwargs_param) + assert scaling[0] == 1 * 2 * 2.9 + + +if __name__ == "__main__": + pytest.main() diff --git a/test/test_Likelihood/test_lens_sample_likelihood.py b/test/test_Likelihood/test_lens_sample_likelihood.py index 8bd3acc0..f0d3b487 100644 --- a/test/test_Likelihood/test_lens_sample_likelihood.py +++ b/test/test_Likelihood/test_lens_sample_likelihood.py @@ -47,10 +47,13 @@ def setup_method(self): "likelihood_type": "DsDdsGaussian", "ds_dds_mean": lensCosmo.ds / lensCosmo.dds, "ds_dds_sigma": 1, - "ani_param_array": ani_param_array, - "ani_scaling_array": ani_scaling_array, + "kin_scaling_param_list": ["a_ani"], + "j_kin_scaling_param_axes": ani_param_array, + "j_kin_scaling_grid_list": [ani_scaling_array], }, ] + kwargs_global_model = {"anisotropy_sampling": True} + self.likelihood = LensSampleLikelihood(kwargs_lens_list=self.kwargs_lens_list) def test_log_likelihood(self): diff --git a/test/test_Likelihood/test_los_distribution.py b/test/test_Likelihood/test_los_distribution.py index 3ac5434e..97b25d08 100644 --- a/test/test_Likelihood/test_los_distribution.py +++ b/test/test_Likelihood/test_los_distribution.py @@ -1,4 +1,4 @@ -from hierarc.Likelihood.los_distributions import LOSDistribution +from hierarc.Sampling.Distributions.los_distributions import LOSDistribution from scipy.stats import genextreme import numpy as np import numpy.testing as npt diff --git a/test/test_Likelihood/test_parameter_scaling.py b/test/test_Likelihood/test_parameter_scaling.py deleted file mode 100644 index d96695fd..00000000 --- a/test/test_Likelihood/test_parameter_scaling.py +++ /dev/null @@ -1,417 +0,0 @@ -import pytest -import numpy as np -import unittest -from hierarc.Likelihood.parameter_scaling import ( - ParameterScalingSingleAperture, - ParameterScalingIFU, -) - - -class TestParameterScalingSingleAperture(object): - def setup_method(self): - ani_param_array = np.linspace(start=0, stop=1, num=10) - param_scaling_array = ani_param_array * 2 - self.scaling = ParameterScalingSingleAperture( - ani_param_array, param_scaling_array - ) - - ani_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - self.scaling_2d = ParameterScalingSingleAperture( - ani_param_array, param_scaling_array - ) - - ani_param_array = np.linspace(start=0, stop=1, num=10) - gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) - log_m2l_array = np.linspace(start=0.1, stop=1, num=10) - - param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] - param_scaling_array = np.multiply.outer( - ani_param_array, - np.outer(gamma_in_array, log_m2l_array), - ) - self.scaling_nfw = ParameterScalingSingleAperture( - param_arrays, param_scaling_array - ) - - gom_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - param_arrays = [ - gom_param_array[0], - gom_param_array[1], - gamma_in_array, - log_m2l_array, - ] - param_scaling_array = np.multiply.outer( - gom_param_array[0], - np.multiply.outer( - gom_param_array[1], - np.outer(gamma_in_array, log_m2l_array), - ), - ) - - self.scaling_nfw_2d = ParameterScalingSingleAperture( - param_arrays, param_scaling_array - ) - - def test_param_scaling(self): - scaling = self.scaling.param_scaling(param_array=[1]) - assert scaling == np.array([2]) - - scaling = self.scaling.param_scaling(param_array=None) - assert scaling == 1 - - scaling = self.scaling_2d.param_scaling(param_array=[1, 2]) - assert scaling == 2 - - scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 0.5]) - assert scaling == 1 * 2.9 * 0.5 - - scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 0.5]) - assert scaling == 1 * 2 * 2.9 * 0.5 - - -class TestParameterScalingIFU(object): - def setup_method(self): - ani_param_array = np.linspace(start=0, stop=1, num=10) - param_scaling_array = ani_param_array * 2 - self.scaling = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=ani_param_array, - scaling_grid_list=[param_scaling_array], - ) - - ani_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - self.scaling_2d = ParameterScalingIFU( - anisotropy_model="GOM", - param_arrays=ani_param_array, - scaling_grid_list=[param_scaling_array], - ) - - ani_param_array = np.linspace(start=0, stop=1, num=10) - gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) - log_m2l_array = np.linspace(start=0.1, stop=1, num=10) - - param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] - param_scaling_array = np.multiply.outer( - ani_param_array, - np.outer(gamma_in_array, log_m2l_array), - ) - self.scaling_nfw = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - - gom_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - param_arrays = [ - gom_param_array[0], - gom_param_array[1], - gamma_in_array, - log_m2l_array, - ] - param_scaling_array = np.multiply.outer( - gom_param_array[0], - np.multiply.outer( - gom_param_array[1], - np.outer(gamma_in_array, log_m2l_array), - ), - ) - self.scaling_nfw_2d = ParameterScalingIFU( - anisotropy_model="GOM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - - param_arrays = [ani_param_array, gamma_in_array] - param_scaling_array = np.multiply.outer( - ani_param_array, - gamma_in_array, - ) - self.scaling_nfw_no_m2l = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - - gom_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - param_arrays = [ - gom_param_array[0], - gom_param_array[1], - gamma_in_array, - ] - param_scaling_array = np.multiply.outer( - gom_param_array[0], - np.multiply.outer( - gom_param_array[1], - gamma_in_array, - ), - ) - self.scaling_nfw_2d_no_m2l = ParameterScalingIFU( - anisotropy_model="GOM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - - def test_param_scaling(self): - scaling = self.scaling.param_scaling(param_array=None) - assert scaling[0] == 1 - - scaling = self.scaling.param_scaling(param_array=[1]) - assert scaling[0] == 2 - - scaling = self.scaling_2d.param_scaling(param_array=[1, 2]) - assert scaling[0] == 2 - - scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 0.5]) - assert scaling[0] == 1 * 2.9 * 0.5 - - scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 0.5]) - assert scaling[0] == 1 * 2 * 2.9 * 0.5 - - scaling = self.scaling_nfw_no_m2l.param_scaling(param_array=[1, 2.9]) - assert scaling[0] == 1 * 2.9 - - scaling = self.scaling_nfw_2d_no_m2l.param_scaling(param_array=[1, 2, 2.9]) - assert scaling[0] == 1 * 2 * 2.9 - - def test_draw_anisotropy(self): - a_ani = 1 - beta_inf = 1.5 - param_draw = self.scaling.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - for i in range(100): - param_draw = self.scaling.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - self.scaling._anisotropy_model = "const" - param_draw = self.scaling.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == 1 - - param_draw = self.scaling_2d.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - assert param_draw[1] == beta_inf - for i in range(100): - param_draw = self.scaling_2d.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - - param_draw = self.scaling_nfw.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - for i in range(100): - param_draw = self.scaling_nfw.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - - param_draw = self.scaling_nfw_2d.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - assert param_draw[1] == beta_inf - for i in range(100): - param_draw = self.scaling_nfw_2d.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - - param_draw = self.scaling_nfw_no_m2l.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - for i in range(100): - param_draw = self.scaling_nfw_no_m2l.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - - param_draw = self.scaling_nfw_2d_no_m2l.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw[0] == a_ani - assert param_draw[1] == beta_inf - for i in range(100): - param_draw = self.scaling_nfw_2d_no_m2l.draw_anisotropy( - a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1 - ) - - scaling = ParameterScalingIFU(anisotropy_model="NONE") - param_draw = scaling.draw_anisotropy( - a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0 - ) - assert param_draw is None - - def test_draw_lens_parameters(self): - param_draw = self.scaling_nfw.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0 - ) - assert param_draw[0] == 1 - assert param_draw[1] == 0.5 - for i in range(100): - param_draw = self.scaling_nfw.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3 - ) - - param_draw = self.scaling_nfw_2d.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0 - ) - assert param_draw[0] == 1 - assert param_draw[1] == 0.5 - - for i in range(100): - param_draw = self.scaling_nfw_2d.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3 - ) - - param_draw = self.scaling_nfw_no_m2l.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0 - ) - assert param_draw == 1 - for i in range(100): - param_draw = self.scaling_nfw_no_m2l.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3 - ) - - param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0 - ) - assert param_draw == 1 - for i in range(100): - param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3 - ) - - -class TestRaise(unittest.TestCase): - def test_raise(self): - with self.assertRaises(ValueError): - ParameterScalingIFU( - anisotropy_model="blabla", - param_arrays=np.array([0, 1]), - scaling_grid_list=[np.array([0, 1])], - ) - - ani_param_array = np.linspace(start=0, stop=1, num=10) - ani_scaling_array = ani_param_array * 2 - scaling = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=ani_param_array, - scaling_grid_list=[ani_scaling_array], - ) - with self.assertRaises(ValueError): - scaling.draw_anisotropy( - a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0 - ) - - ani_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - scaling = ParameterScalingIFU( - anisotropy_model="GOM", - param_arrays=ani_param_array, - scaling_grid_list=[ani_scaling_array], - ) - with self.assertRaises(ValueError): - scaling.draw_anisotropy( - a_ani=0.5, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0 - ) - - ani_param_array = np.linspace(start=0, stop=1, num=10) - gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) - log_m2l_array = np.linspace(start=0.1, stop=1, num=10) - - param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] - param_scaling_array = np.multiply.outer( - ani_param_array, - np.multiply.outer(gamma_in_array, log_m2l_array), - ) - print(param_scaling_array.shape) - scaling = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - with self.assertRaises(ValueError): - scaling.draw_lens_parameters( - gamma_in=-1, gamma_in_sigma=0.1, log_m2l=0.5, log_m2l_sigma=0.1 - ) - with self.assertRaises(ValueError): - scaling.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0.1, log_m2l=-0.5, log_m2l_sigma=0.1 - ) - - ani_param_array = np.linspace(start=0, stop=1, num=10) - gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) - log_m2l_array = np.linspace(start=0.1, stop=1, num=10) - - param_arrays = [ani_param_array, gamma_in_array, log_m2l_array] - param_scaling_array = np.multiply.outer( - ani_param_array, - np.multiply.outer(gamma_in_array, log_m2l_array), - ) - - scaling = ParameterScalingIFU( - anisotropy_model="OM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - with self.assertRaises(ValueError): - scaling.draw_lens_parameters( - gamma_in=1, gamma_in_sigma=0, log_m2l=0, log_m2l_sigma=0 - ) - - gom_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - ] - - gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) - - param_arrays = [ - gom_param_array[0], - gom_param_array[1], - gamma_in_array, - ] - param_scaling_array = np.multiply.outer( - gom_param_array[0], - np.multiply.outer( - gom_param_array[1], - gamma_in_array, - ), - ) - self.scaling_nfw_2d_no_m2l = ParameterScalingIFU( - anisotropy_model="GOM", - param_arrays=param_arrays, - scaling_grid_list=[param_scaling_array], - ) - - with self.assertRaises(ValueError): - param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters( - gamma_in=-1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3 - ) - - -if __name__ == "__main__": - pytest.main() diff --git a/test/test_Sampling/test_Distributions/__init__.py b/test/test_Sampling/test_Distributions/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py new file mode 100644 index 00000000..6dd65ce8 --- /dev/null +++ b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py @@ -0,0 +1,36 @@ +from hierarc.Sampling.Distributions.anisotropy_distributions import AnisotropyDistribution + + +class TestAnisotropyDistribution(object): + + def setup_method(self): + anisotropy_model = "GOM" + distribution_function = "GAUSSIAN" + kwargs_anisotropy_min = {"a_ani": 0, "beta_inf": 0.1} + kwargs_anisotropy_max = {"a_ani": 5, "beta_inf": 1} + + self._ani_dist = AnisotropyDistribution(anisotropy_model=anisotropy_model, + anisotropy_sampling=True, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max) + + ani_dist = AnisotropyDistribution(anisotropy_model=anisotropy_model, + anisotropy_sampling=False, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max) + + def test_draw_anisotropy(self): + kwargs_anisotropy = {"a_ani": 1, "beta_inf": 0.8, "a_ani_sigma": 0.1, "beta_inf_sigma": 0.2} + kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) + assert "a_ani" in kwargs_drawn + assert "beta_inf" in kwargs_drawn + + ani_dist = AnisotropyDistribution(anisotropy_model="NONE", + anisotropy_sampling=False, + distribution_function="NONE", + kwargs_anisotropy_min=None, + kwargs_anisotropy_max=None) + kwargs_drawn = ani_dist.draw_anisotropy() + assert "a_ani" not in kwargs_drawn diff --git a/test/test_Sampling/test_Distributions/test_lens_distribution.py b/test/test_Sampling/test_Distributions/test_lens_distribution.py new file mode 100644 index 00000000..6f71712f --- /dev/null +++ b/test/test_Sampling/test_Distributions/test_lens_distribution.py @@ -0,0 +1,51 @@ +import copy + +from hierarc.Sampling.Distributions.lens_distribution import LensDistribution + + +class TestLensDistribution(object): + + def setup_method(self): + self.kwargs_sampling = {"lambda_mst_distribution": "GAUSSIAN", + "gamma_in_sampling": True, + "gamma_in_distribution": "GAUSSIAN", + "log_m2l_sampling": True, + "log_m2l_distribution" : "GAUSSIAN", + "alpha_lambda_sampling": True, + "beta_lambda_sampling": True, + "alpha_gamma_in_sampling": True, + "alpha_log_m2l_sampling": True, + "log_scatter": False, # change for different tests + "mst_ifu": False, # change for different tests + "lambda_scaling_property": 0.1, + "lambda_scaling_property_beta": 0.2, + "kwargs_min": {"gamma_in": 0, "log_m2l": -3}, + "kwargs_max": {"gamma_in": 3, "log_m2l": 3}} + + self.kwargs_lens = {"lambda_mst": 1.1, + "lambda_mst_sigma": 0.1, + "gamma_ppn": 0.9, + "lambda_ifu": 0.5, + "lambda_ifu_sigma": 0.2, + "alpha_lambda": -0.2, + "beta_lambda": 0.3, + "gamma_in": 1.5, + "gamma_in_sigma": 0.2, + "alpha_gamma_in": 0.2, + "log_m2l": 0.6, + "log_m2l_sigma": 0.2, + "alpha_log_m2l": -0.1} + + def test_draw_lens(self): + lens_dist = LensDistribution(**self.kwargs_sampling) + kwargs_return = lens_dist.draw_lens(**self.kwargs_lens) + + assert "lambda_mst" in kwargs_return + + kwargs_sampling = copy.deepcopy(self.kwargs_sampling) + kwargs_sampling["log_scatter"] = True + kwargs_sampling["lambda_ifu"] = True + lens_dist = LensDistribution(kwargs_sampling) + kwargs_return = lens_dist.draw_lens(**self.kwargs_lens) + + assert "lambda_mst" in kwargs_return diff --git a/test/test_Sampling/test_mcmc_sampling.py b/test/test_Sampling/test_mcmc_sampling.py index 8926e01d..e27550ec 100644 --- a/test/test_Sampling/test_mcmc_sampling.py +++ b/test/test_Sampling/test_mcmc_sampling.py @@ -67,16 +67,19 @@ def test_mcmc_emcee(self): backend = emcee.backends.HDFBackend(backup_filename) kwargs_emcee = {"backend": backend} + kwargs_model = {"ppn_sampling": False, + "lambda_mst_sampling": False, + "lambda_mst_distribution": "delta", + "anisotropy_sampling": False, + "anisotropy_model": "OM", + + } mcmc_sampler = MCMCSampler( kwargs_likelihood_list=kwargs_likelihood_list, cosmology=cosmology, kwargs_bounds=kwargs_bounds, - ppn_sampling=False, - lambda_mst_sampling=False, - lambda_mst_distribution="delta", - anisotropy_sampling=False, - anisotropy_model="OM", + kwargs_model=kwargs_model, custom_prior=None, interpolate_cosmo=True, num_redshift_interp=100, From 4bb09ba9f96c261a633a4ac2b4f7529dbcf549b9 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jun 2024 01:57:35 +0000 Subject: [PATCH 02/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- hierarc/LensPosterior/base_config.py | 10 +- hierarc/LensPosterior/kin_constraints.py | 8 +- .../kin_constraints_composite.py | 13 +-- hierarc/LensPosterior/kin_scaling_config.py | 13 ++- hierarc/Likelihood/cosmo_likelihood.py | 6 +- hierarc/Likelihood/hierarchy_likelihood.py | 84 +++++++++------ hierarc/Likelihood/kin_scaling.py | 37 +++++-- hierarc/Likelihood/lens_sample_likelihood.py | 38 +++---- .../Distributions/anisotropy_distributions.py | 23 ++-- .../Distributions/lens_distribution.py | 84 ++++++++------- .../test_kin_constraints_composite.py | 28 +++-- test/test_Likelihood/test_cosmo_likelihood.py | 24 +++-- .../test_hierarchy_likelihood.py | 29 ++--- test/test_Likelihood/test_kin_scaling.py | 101 ++++++++++-------- .../test_anisotropy_distribution.py | 47 +++++--- .../test_lens_distribution.py | 62 ++++++----- test/test_Sampling/test_mcmc_sampling.py | 14 +-- 17 files changed, 365 insertions(+), 256 deletions(-) diff --git a/hierarc/LensPosterior/base_config.py b/hierarc/LensPosterior/base_config.py index 258c12f7..ee39e438 100644 --- a/hierarc/LensPosterior/base_config.py +++ b/hierarc/LensPosterior/base_config.py @@ -109,6 +109,10 @@ def __init__( ImageModelPosterior.__init__( self, theta_E, theta_E_error, gamma, gamma_error, r_eff, r_eff_error ) - KinScalingConfig.__init__(self, anisotropy_model, r_eff, - gamma_in_scaling=gamma_in_scaling, - log_m2l_scaling=log_m2l_scaling) + KinScalingConfig.__init__( + self, + anisotropy_model, + r_eff, + gamma_in_scaling=gamma_in_scaling, + log_m2l_scaling=log_m2l_scaling, + ) diff --git a/hierarc/LensPosterior/kin_constraints.py b/hierarc/LensPosterior/kin_constraints.py index 7bae1ea4..15cf2a9a 100644 --- a/hierarc/LensPosterior/kin_constraints.py +++ b/hierarc/LensPosterior/kin_constraints.py @@ -180,7 +180,6 @@ def hierarchy_configuration(self, num_sample_model=20): "kin_scaling_param_list": self.param_name_list, "j_kin_scaling_param_axes": self.kin_scaling_param_array, "j_kin_scaling_grid_list": ani_scaling_array_list, - } return kwargs_likelihood @@ -253,7 +252,12 @@ def _anisotropy_scaling_relative(self, j_ani_0): if self._anisotropy_model == "GOM": ani_scaling_array_list = [ - np.zeros((len(self.kin_scaling_param_array[0]), len(self.kin_scaling_param_array[1]))) + np.zeros( + ( + len(self.kin_scaling_param_array[0]), + len(self.kin_scaling_param_array[1]), + ) + ) for _ in range(num_data) ] for i, a_ani in enumerate(self.kin_scaling_param_array[0]): diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py index 99e87118..60bcc145 100644 --- a/hierarc/LensPosterior/kin_constraints_composite.py +++ b/hierarc/LensPosterior/kin_constraints_composite.py @@ -400,21 +400,18 @@ def hierarchy_configuration(self, num_sample_model=20): "j_model": j_model_list, "error_cov_measurement": error_cov_measurement, "error_cov_j_sqrt": error_cov_j_sqrt, - - #"ani_param_array": self.kin_scaling_param_array, - #"gamma_in_array": self.gamma_in_array, - #"log_m2l_array": self.log_m2l_array, - #"param_scaling_grid_list": ani_scaling_grid_list, - + # "ani_param_array": self.kin_scaling_param_array, + # "gamma_in_array": self.gamma_in_array, + # "log_m2l_array": self.log_m2l_array, + # "param_scaling_grid_list": ani_scaling_grid_list, "gamma_in_prior_mean": self._gamma_in_prior_mean, "gamma_in_prior_std": self._gamma_in_prior_std, - "kin_scaling_param_list": self.param_name_list, "j_kin_scaling_param_axes": self.kin_scaling_param_array, "j_kin_scaling_grid_list": ani_scaling_grid_list, } - #if not self._is_m2l_population_level: + # if not self._is_m2l_population_level: # kwargs_likelihood["log_m2l_array"] = None return kwargs_likelihood diff --git a/hierarc/LensPosterior/kin_scaling_config.py b/hierarc/LensPosterior/kin_scaling_config.py index 8a51b2ad..d6979f3d 100644 --- a/hierarc/LensPosterior/kin_scaling_config.py +++ b/hierarc/LensPosterior/kin_scaling_config.py @@ -5,7 +5,9 @@ class KinScalingConfig(object): """Class to manage the anisotropy model and parameters for the Posterior processing.""" - def __init__(self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scaling=None): + def __init__( + self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scaling=None + ): """ :param anisotropy_model: type of stellar anisotropy model. Supported are 'OM' and 'GOM' or 'const', see details in lenstronomy.Galkin module @@ -28,10 +30,12 @@ def __init__(self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scali ] self._param_name_list = ["a_ani", "beta_inf"] elif self._anisotropy_model == "const": - self._ani_param_array = [np.linspace(-0.49, 1, 7)] # used for constant anisotropy description + self._ani_param_array = [ + np.linspace(-0.49, 1, 7) + ] # used for constant anisotropy description self._param_name_list = ["a_ani"] elif self._anisotropy_model == "NONE": - self._param_name_list =[] + self._param_name_list = [] else: raise ValueError( "anisotropy model %s not supported." % self._anisotropy_model @@ -77,8 +81,7 @@ def kin_scaling_param_array(self): @property def param_name_list(self): - """ - list of parameters in same order as interpolated + """List of parameters in same order as interpolated. :return: """ diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index f900cbfa..8ad2a06f 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -60,11 +60,7 @@ def __init__( normalized=normalized, kwargs_global_model=kwargs_model, ) - self.param = ParamManager( - cosmology, - **kwargs_model, - **kwargs_bounds - ) + self.param = ParamManager(cosmology, **kwargs_model, **kwargs_bounds) self._lower_limit, self._upper_limit = self.param.param_bounds self._prior_add = False if custom_prior is not None: diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py index 2b5aaf7b..16ba54fb 100644 --- a/hierarc/Likelihood/hierarchy_likelihood.py +++ b/hierarc/Likelihood/hierarchy_likelihood.py @@ -2,7 +2,9 @@ from hierarc.Likelihood.LensLikelihood.base_lens_likelihood import LensLikelihoodBase from hierarc.Likelihood.kin_scaling import KinScaling from hierarc.Sampling.Distributions.los_distributions import LOSDistribution -from hierarc.Sampling.Distributions.anisotropy_distributions import AnisotropyDistribution +from hierarc.Sampling.Distributions.anisotropy_distributions import ( + AnisotropyDistribution, +) from hierarc.Sampling.Distributions.lens_distribution import LensDistribution import numpy as np import copy @@ -94,11 +96,12 @@ def __init__( """ TransformedCosmography.__init__(self, z_lens=z_lens, z_source=z_source) - KinScaling.__init__(self, - j_kin_scaling_param_axes=j_kin_scaling_param_axes, - j_kin_scaling_grid_list=j_kin_scaling_grid_list, - j_kin_scaling_param_name_list=kin_scaling_param_list - ) + KinScaling.__init__( + self, + j_kin_scaling_param_axes=j_kin_scaling_param_axes, + j_kin_scaling_grid_list=j_kin_scaling_grid_list, + j_kin_scaling_param_name_list=kin_scaling_param_list, + ) LensLikelihoodBase.__init__( self, @@ -119,28 +122,32 @@ def __init__( los_distributions=los_distributions, ) kwargs_min, kwargs_max = self.param_bounds_interpol() - self._lens_distribution = LensDistribution(lambda_mst_sampling=False, - lambda_mst_distribution=lambda_mst_distribution, - gamma_in_sampling=gamma_in_sampling, - gamma_in_distribution=gamma_in_distribution, - log_m2l_sampling=log_m2l_sampling, - log_m2l_distribution=log_m2l_distribution, - alpha_lambda_sampling=alpha_lambda_sampling, - beta_lambda_sampling=beta_lambda_sampling, - alpha_gamma_in_sampling=alpha_gamma_in_sampling, - alpha_log_m2l_sampling=alpha_log_m2l_sampling, - log_scatter=log_scatter, - mst_ifu=mst_ifu, - lambda_scaling_property=lambda_scaling_property, - lambda_scaling_property_beta=lambda_scaling_property_beta, - kwargs_min=kwargs_min, - kwargs_max=kwargs_max,) - - self._aniso_distribution = AnisotropyDistribution(anisotropy_model=anisotropy_model, - anisotropy_sampling=anisotropy_sampling, - distribution_function=anisotroy_distribution_function, - kwargs_anisotropy_min=kwargs_min, - kwargs_anisotropy_max=kwargs_max) + self._lens_distribution = LensDistribution( + lambda_mst_sampling=False, + lambda_mst_distribution=lambda_mst_distribution, + gamma_in_sampling=gamma_in_sampling, + gamma_in_distribution=gamma_in_distribution, + log_m2l_sampling=log_m2l_sampling, + log_m2l_distribution=log_m2l_distribution, + alpha_lambda_sampling=alpha_lambda_sampling, + beta_lambda_sampling=beta_lambda_sampling, + alpha_gamma_in_sampling=alpha_gamma_in_sampling, + alpha_log_m2l_sampling=alpha_log_m2l_sampling, + log_scatter=log_scatter, + mst_ifu=mst_ifu, + lambda_scaling_property=lambda_scaling_property, + lambda_scaling_property_beta=lambda_scaling_property_beta, + kwargs_min=kwargs_min, + kwargs_max=kwargs_max, + ) + + self._aniso_distribution = AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=anisotropy_sampling, + distribution_function=anisotroy_distribution_function, + kwargs_anisotropy_min=kwargs_min, + kwargs_anisotropy_max=kwargs_max, + ) self._gamma_in_prior_mean = gamma_in_prior_mean self._gamma_in_prior_std = gamma_in_prior_std @@ -278,7 +285,10 @@ def log_likelihood_single( realization of the hyperparameter distribution """ kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) - lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] + lambda_mst, gamma_ppn = ( + kwargs_lens_draw["lambda_mst"], + kwargs_lens_draw["gamma_ppn"], + ) kappa_ext = self._los.draw_los(kwargs_los) # draw intrinsic source magnitude @@ -309,9 +319,9 @@ def log_likelihood_single( and "gamma_in" in kwargs_lens_draw ): gamma_in = kwargs_lens_draw["gamma_in"] - lnlikelihood -= ( - self._gamma_in_prior_mean - gamma_in - ) ** 2 / (2 * self._gamma_in_prior_std**2) + lnlikelihood -= (self._gamma_in_prior_mean - gamma_in) ** 2 / ( + 2 * self._gamma_in_prior_std**2 + ) return np.nan_to_num(lnlikelihood) @@ -433,7 +443,10 @@ def sigma_v_measured_vs_predict( cov_error_predict = np.zeros_like(cov_error_measurement) for i in range(self._num_distribution_draws): kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) - lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] + lambda_mst, gamma_ppn = ( + kwargs_lens_draw["lambda_mst"], + kwargs_lens_draw["gamma_ppn"], + ) 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 @@ -478,7 +491,10 @@ def ddt_dd_model_prediction(self, cosmo, kwargs_lens=None, kwargs_los=None): dd_draws = [] for i in range(self._num_distribution_draws): kwargs_lens_draw = self._lens_distribution.draw_lens(**kwargs_lens) - lambda_mst, gamma_ppn = kwargs_lens_draw["lambda_mst"], kwargs_lens_draw["gamma_ppn"] + lambda_mst, gamma_ppn = ( + kwargs_lens_draw["lambda_mst"], + kwargs_lens_draw["gamma_ppn"], + ) 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/kin_scaling.py b/hierarc/Likelihood/kin_scaling.py index f016adf8..40e75f58 100644 --- a/hierarc/Likelihood/kin_scaling.py +++ b/hierarc/Likelihood/kin_scaling.py @@ -26,9 +26,13 @@ def __init__(self, param_grid_axes, j_kin_scaling_grid): param_grid_axes = [param_grid_axes] if self._dim_scaling == 1: - self._f_ani = interp1d(param_grid_axes[0], j_kin_scaling_grid, kind="linear") + self._f_ani = interp1d( + param_grid_axes[0], j_kin_scaling_grid, kind="linear" + ) elif self._dim_scaling == 2: - self._f_ani = interp2d(param_grid_axes[0], param_grid_axes[1], j_kin_scaling_grid.T) + self._f_ani = interp2d( + param_grid_axes[0], param_grid_axes[1], j_kin_scaling_grid.T + ) else: self._f_ani = RegularGridInterpolator( tuple(param_grid_axes), @@ -57,7 +61,10 @@ class KinScaling(object): """Class to manage model parameter and anisotropy scalings for IFU data.""" def __init__( - self, j_kin_scaling_param_axes=None, j_kin_scaling_grid_list=None, j_kin_scaling_param_name_list=None + self, + j_kin_scaling_param_axes=None, + j_kin_scaling_grid_list=None, + j_kin_scaling_param_name_list=None, ): """ @@ -71,7 +78,10 @@ def __init__( else: self._param_list = j_kin_scaling_param_name_list self._param_arrays = j_kin_scaling_param_axes - if not isinstance(j_kin_scaling_param_axes, list) and j_kin_scaling_param_name_list is not None: + if ( + not isinstance(j_kin_scaling_param_axes, list) + and j_kin_scaling_param_name_list is not None + ): self._param_arrays = [j_kin_scaling_param_axes] self._evaluate_scaling = False self._is_log_m2l_population_level = False @@ -85,7 +95,9 @@ def __init__( self._f_ani_list = [] for scaling_grid in j_kin_scaling_grid_list: self._j_scaling_ifu.append( - ParameterScalingSingleMeasurement(j_kin_scaling_param_axes, scaling_grid) + ParameterScalingSingleMeasurement( + j_kin_scaling_param_axes, scaling_grid + ) ) if isinstance(j_kin_scaling_param_axes, list): @@ -94,22 +106,25 @@ def __init__( self._dim_scaling = 1 def _kwargs2param_array(self, kwargs): - """ - converts dictionary to sorted array in same order as interpolation grid + """Converts dictionary to sorted array in same order as interpolation grid. - :param kwargs: dictionary of all model components, must include the one that are interpolated + :param kwargs: dictionary of all model components, must include the one that are + interpolated :return: sorted list of parameters to interpolate """ param_array = [] for param in self._param_list: if param not in kwargs: - raise ValueError("key %s not in parameters and hence kinematic scaling not possible" % param) + raise ValueError( + "key %s not in parameters and hence kinematic scaling not possible" + % param + ) param_array.append(kwargs.get(param)) return param_array def param_bounds_interpol(self): - """ - minimum and maximum bounds of parameters that are being used to call interpolation function + """Minimum and maximum bounds of parameters that are being used to call + interpolation function. :return: dictionaries of minimum and maximum bounds """ diff --git a/hierarc/Likelihood/lens_sample_likelihood.py b/hierarc/Likelihood/lens_sample_likelihood.py index d168b097..a8e1431d 100644 --- a/hierarc/Likelihood/lens_sample_likelihood.py +++ b/hierarc/Likelihood/lens_sample_likelihood.py @@ -29,11 +29,10 @@ def __init__(self, kwargs_lens_list, normalized=False, kwargs_global_model=None) DSPLikelihood(normalized=normalized, **_kwargs_lens) ) else: - kwargs_lens_ = self._merge_global2local_settings(kwargs_global_model=kwargs_global_model, - kwargs_lens=kwargs_lens) - self._lens_list.append( - LensLikelihood(**kwargs_lens_) + kwargs_lens_ = self._merge_global2local_settings( + kwargs_global_model=kwargs_global_model, kwargs_lens=kwargs_lens ) + self._lens_list.append(LensLikelihood(**kwargs_lens_)) def log_likelihood( self, @@ -91,18 +90,19 @@ def _merge_global2local_settings(kwargs_global_model, kwargs_lens): return {**kwargs_global_model_subset, **kwargs_lens} -_input_param_list = ["anisotropy_model", - "anisotropy_sampling", - "anisotroy_distribution_function", - "los_distributions", - "lambda_mst_distribution", - "gamma_in_sampling", - "gamma_in_distribution", - "log_m2l_sampling", - "log_m2l_distribution", - "alpha_lambda_sampling", - "beta_lambda_sampling", - "alpha_gamma_in_sampling", - "alpha_log_m2l_sampling", - "log_scatter" - ] +_input_param_list = [ + "anisotropy_model", + "anisotropy_sampling", + "anisotroy_distribution_function", + "los_distributions", + "lambda_mst_distribution", + "gamma_in_sampling", + "gamma_in_distribution", + "log_m2l_sampling", + "log_m2l_distribution", + "alpha_lambda_sampling", + "beta_lambda_sampling", + "alpha_gamma_in_sampling", + "alpha_log_m2l_sampling", + "log_scatter", +] diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py index c5f92a0a..1050a1fb 100644 --- a/hierarc/Sampling/Distributions/anisotropy_distributions.py +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -2,11 +2,16 @@ class AnisotropyDistribution(object): - """ - class to draw anisotropy parameters from hyperparameter distributions - """ - def __init__(self, anisotropy_model, anisotropy_sampling, - distribution_function, kwargs_anisotropy_min, kwargs_anisotropy_max): + """Class to draw anisotropy parameters from hyperparameter distributions.""" + + def __init__( + self, + anisotropy_model, + anisotropy_sampling, + distribution_function, + kwargs_anisotropy_min, + kwargs_anisotropy_max, + ): """ :param anisotropy_model: string, name of anisotropy model to consider @@ -25,11 +30,15 @@ def __init__(self, anisotropy_model, anisotropy_sampling, kwargs_anisotropy_max = {} self._kwargs_min = kwargs_anisotropy_min self._kwargs_max = kwargs_anisotropy_max - self._a_ani_min, self._a_ani_max = self._kwargs_min.get("a_ani", -np.inf), self._kwargs_max.get("a_ani", np.inf) + self._a_ani_min, self._a_ani_max = self._kwargs_min.get( + "a_ani", -np.inf + ), self._kwargs_max.get("a_ani", np.inf) self._beta_inf_min = self._kwargs_min.get("beta_inf", -np.inf) self._beta_inf_max = self._kwargs_max.get("beta_inf", np.inf) - def draw_anisotropy(self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0): + def draw_anisotropy( + self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0 + ): """Draw Gaussian distribution and re-sample if outside bounds. :param a_ani: mean of the distribution diff --git a/hierarc/Sampling/Distributions/lens_distribution.py b/hierarc/Sampling/Distributions/lens_distribution.py index dd1c2738..0c98b168 100644 --- a/hierarc/Sampling/Distributions/lens_distribution.py +++ b/hierarc/Sampling/Distributions/lens_distribution.py @@ -2,9 +2,7 @@ class LensDistribution(object): - """ - class to draw lens parameters of individual lens from distributions - """ + """Class to draw lens parameters of individual lens from distributions.""" def __init__( self, @@ -74,8 +72,12 @@ def __init__( kwargs_max = {} if kwargs_min is None: kwargs_min = {} - self._gamma_in_min, self._gamma_in_max = kwargs_min.get("gamma_in", -np.inf), kwargs_max.get("gamma_in", np.inf) - self._log_m2l_min, self._log_m2l_max = kwargs_min.get("log_m2l", -np.inf), kwargs_max.get("log_m2l", np.inf) + self._gamma_in_min, self._gamma_in_max = kwargs_min.get( + "gamma_in", -np.inf + ), kwargs_max.get("gamma_in", np.inf) + self._log_m2l_min, self._log_m2l_max = kwargs_min.get( + "log_m2l", -np.inf + ), kwargs_max.get("log_m2l", np.inf) def draw_lens( self, @@ -93,8 +95,7 @@ def draw_lens( log_m2l_sigma=0, alpha_log_m2l=0, ): - """Draws a realization of a specific model from the hyperparameter - distribution. + """Draws a realization of a specific model from the hyperparameter distribution. :param lambda_mst: MST transform :param lambda_mst_sigma: spread in the distribution @@ -124,10 +125,11 @@ def draw_lens( else: lambda_mst_mean_lens = lambda_mst - lambda_lens = (lambda_mst_mean_lens - + alpha_lambda * self._lambda_scaling_property - + beta_lambda * self._lambda_scaling_property_beta - ) + lambda_lens = ( + lambda_mst_mean_lens + + alpha_lambda * self._lambda_scaling_property + + beta_lambda * self._lambda_scaling_property_beta + ) lambda_mst_draw = lambda_lens if self._lambda_mst_sampling: if self._lambda_mst_distribution in ["GAUSSIAN"]: @@ -142,24 +144,28 @@ def draw_lens( "gamma_in parameter is out of bounds of the interpolated range!" ) if self._gamma_in_distribution in ["GAUSSIAN"]: - gamma_in_lens = gamma_in + alpha_gamma_in * self._lambda_scaling_property + gamma_in_lens = ( + gamma_in + alpha_gamma_in * self._lambda_scaling_property + ) else: gamma_in_lens = gamma_in gamma_in_draw = np.random.normal(gamma_in_lens, gamma_in_sigma) if gamma_in_draw < self._gamma_in_min or gamma_in_draw > self._gamma_in_max: - return self.draw_lens(lambda_mst=lambda_mst, - lambda_mst_sigma=lambda_mst_sigma, - gamma_ppn=gamma_ppn, - lambda_ifu=lambda_ifu, - lambda_ifu_sigma=lambda_ifu_sigma, - alpha_lambda=alpha_lambda, - beta_lambda=beta_lambda, - gamma_in=gamma_in, - gamma_in_sigma=gamma_in_sigma, - alpha_gamma_in=alpha_gamma_in, - log_m2l=log_m2l, - log_m2l_sigma=log_m2l_sigma, - alpha_log_m2l=alpha_log_m2l) + return self.draw_lens( + lambda_mst=lambda_mst, + lambda_mst_sigma=lambda_mst_sigma, + gamma_ppn=gamma_ppn, + lambda_ifu=lambda_ifu, + lambda_ifu_sigma=lambda_ifu_sigma, + alpha_lambda=alpha_lambda, + beta_lambda=beta_lambda, + gamma_in=gamma_in, + gamma_in_sigma=gamma_in_sigma, + alpha_gamma_in=alpha_gamma_in, + log_m2l=log_m2l, + log_m2l_sigma=log_m2l_sigma, + alpha_log_m2l=alpha_log_m2l, + ) kwargs_return["gamma_in"] = gamma_in_draw if self._log_m2l_sampling: @@ -172,18 +178,20 @@ def draw_lens( log_m2l_draw = np.random.normal(log_m2l_lens, log_m2l_sigma) if log_m2l_draw < self._log_m2l_min or log_m2l_draw > self._log_m2l_max: - return self.draw_lens(lambda_mst=lambda_mst, - lambda_mst_sigma=lambda_mst_sigma, - gamma_ppn=gamma_ppn, - lambda_ifu=lambda_ifu, - lambda_ifu_sigma=lambda_ifu_sigma, - alpha_lambda=alpha_lambda, - beta_lambda=beta_lambda, - gamma_in=gamma_in, - gamma_in_sigma=gamma_in_sigma, - alpha_gamma_in=alpha_gamma_in, - log_m2l=log_m2l, - log_m2l_sigma=log_m2l_sigma, - alpha_log_m2l=alpha_log_m2l) + return self.draw_lens( + lambda_mst=lambda_mst, + lambda_mst_sigma=lambda_mst_sigma, + gamma_ppn=gamma_ppn, + lambda_ifu=lambda_ifu, + lambda_ifu_sigma=lambda_ifu_sigma, + alpha_lambda=alpha_lambda, + beta_lambda=beta_lambda, + gamma_in=gamma_in, + gamma_in_sigma=gamma_in_sigma, + alpha_gamma_in=alpha_gamma_in, + log_m2l=log_m2l, + log_m2l_sigma=log_m2l_sigma, + alpha_log_m2l=alpha_log_m2l, + ) kwargs_return["log_m2l"] = log_m2l_draw return kwargs_return diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py index f99c6168..5b7b7791 100644 --- a/test/test_LensPosterior/test_kin_constraints_composite.py +++ b/test/test_LensPosterior/test_kin_constraints_composite.py @@ -104,9 +104,13 @@ def test_likelihoodconfiguration_om(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood) + ln_class = LensLikelihood( + gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood + ) kwargs_kin = {"a_ani": 1} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood( + cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin + ) kwargs_lens_light_test = [{"amp": [1, 1], "sigma": [1, 2]}] lens_light_model_list_test = ["MULTI_GAUSSIAN"] @@ -256,9 +260,13 @@ def test_likelihoodconfiguration_gom(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood) + ln_class = LensLikelihood( + gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood + ) kwargs_kin = {"a_ani": 1, "beta_inf": 0.5} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood( + cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin + ) class TestKinConstraintsCompositeM2l(object): @@ -360,9 +368,13 @@ def test_likelihoodconfiguration_om(self): kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) kwargs_likelihood["normalized"] = False - ln_class = LensLikelihood(gamma_in_sampling=True, log_m2l_sampling=False, **kwargs_likelihood) + ln_class = LensLikelihood( + gamma_in_sampling=True, log_m2l_sampling=False, **kwargs_likelihood + ) kwargs_kin = {"a_ani": 1} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 0}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood( + cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 0}, kwargs_kin=kwargs_kin + ) def test_likelihoodconfiguration_gom(self): anisotropy_model = "GOM" @@ -457,7 +469,9 @@ def test_likelihoodconfiguration_gom(self): kwargs_likelihood["normalized"] = False ln_class = LensLikelihood(gamma_in_sampling=True, **kwargs_likelihood) kwargs_kin = {"a_ani": 1, "beta_inf": 0.5} - ln_class.lens_log_likelihood(cosmo, kwargs_lens={"gamma_in": 2}, kwargs_kin=kwargs_kin) + ln_class.lens_log_likelihood( + cosmo, kwargs_lens={"gamma_in": 2}, kwargs_kin=kwargs_kin + ) class TestRaise(unittest.TestCase): diff --git a/test/test_Likelihood/test_cosmo_likelihood.py b/test/test_Likelihood/test_cosmo_likelihood.py index aa863637..fcbd2fc9 100644 --- a/test/test_Likelihood/test_cosmo_likelihood.py +++ b/test/test_Likelihood/test_cosmo_likelihood.py @@ -57,11 +57,13 @@ def setup_method(self): "kwargs_lower_cosmo": kwargs_lower_cosmo, "kwargs_upper_cosmo": kwargs_upper_cosmo, } - self.kwargs_model = {"ppn_sampling": False, - "lambda_mst_sampling": False, - "lambda_mst_distribution": "delta", - "anisotropy_sampling": False, - "anisotropy_model": "OM",} + self.kwargs_model = { + "ppn_sampling": False, + "lambda_mst_sampling": False, + "lambda_mst_distribution": "delta", + "anisotropy_sampling": False, + "anisotropy_model": "OM", + } # self.kwargs_likelihood_list = [{'z_lens': self.z_L, 'z_source': self.z_S, 'likelihood_type': 'TDKinKDE', # 'dd_sample': self.dd_samples, 'ddt_sample': self.ddt_samples, @@ -79,7 +81,9 @@ def test_log_likelihood(self): cosmo_fixed=None, ) - def custom_prior(kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los): + def custom_prior( + kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los + ): return -1 cosmoL_prior = CosmoLikelihood( @@ -241,8 +245,12 @@ def test_kde_likelihood_integration(self): self.cosmology, rescale=True, ) - cosmoL = CosmoLikelihood([], "FLCDM", self.kwargs_model, - self.kwargs_bounds, KDE_likelihood_chain=chain + cosmoL = CosmoLikelihood( + [], + "FLCDM", + self.kwargs_model, + self.kwargs_bounds, + KDE_likelihood_chain=chain, ) kwargs_cosmo = {"h0": self.H0_true, "om": self.omega_m_true, "ok": 0} args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs_cosmo) diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py index f17c828a..f89e8a2d 100644 --- a/test/test_Likelihood/test_hierarchy_likelihood.py +++ b/test/test_Likelihood/test_hierarchy_likelihood.py @@ -24,16 +24,16 @@ def setup_method(self): "dd_mean": dd, "dd_sigma": dd / 10, } - kwargs_model = {"anisotropy_model": "OM", - "anisotropy_sampling": True, - "anisotroy_distribution_function": "GAUSSIAN", - "lambda_mst_distribution": "GAUSSIAN", - } + kwargs_model = { + "anisotropy_model": "OM", + "anisotropy_sampling": True, + "anisotroy_distribution_function": "GAUSSIAN", + "lambda_mst_distribution": "GAUSSIAN", + } # "gamma_in_sampling" = False, - gamma_in_distribution = "NONE", - log_m2l_sampling = False, - log_m2l_distribution = "NONE", - + gamma_in_distribution = ("NONE",) + log_m2l_sampling = (False,) + log_m2l_distribution = ("NONE",) self.likelihood = LensLikelihood( z_lens, @@ -134,7 +134,9 @@ def setup_method(self): **kwargs_model ) - param_scaling_array = np.outer(np.ones_like(gamma_in_array), np.ones_like(log_m2l_array)) + param_scaling_array = np.outer( + np.ones_like(gamma_in_array), np.ones_like(log_m2l_array) + ) j_kin_scaling_param_axes = [gamma_in_array, log_m2l_array] @@ -239,7 +241,6 @@ def test_lens_log_likelihood(self): kwargs_test = self.likelihood._kwargs_init(kwargs=None) assert type(kwargs_test) is dict - kwargs_lens = { "gamma_in": 1, "gamma_in_sigma": 0, @@ -276,11 +277,11 @@ def test_lens_log_likelihood(self): dds = self.cosmo.angular_diameter_distance_z1z2(z1=z_lens, z2=z_source).value ddt = (1.0 + z_lens) * dd * ds / dds - #ln_likelihood = self.likelihood_gamma_in_fail_case.log_likelihood_single( + # ln_likelihood = self.likelihood_gamma_in_fail_case.log_likelihood_single( # ddt, dd, delta_lum_dist, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los - #) + # ) - #assert ln_likelihood < -10000000 + # assert ln_likelihood < -10000000 if __name__ == "__main__": diff --git a/test/test_Likelihood/test_kin_scaling.py b/test/test_Likelihood/test_kin_scaling.py index f66aa32f..f1201660 100644 --- a/test/test_Likelihood/test_kin_scaling.py +++ b/test/test_Likelihood/test_kin_scaling.py @@ -11,9 +11,11 @@ def test_single_param(self): param_arrays = np.linspace(0, 1, 11) scaling_grid_list = [param_arrays**2] param_list = ["a"] - kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=scaling_grid_list, - j_kin_scaling_param_name_list=param_list) + kin_scaling = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list, + ) kwargs_param = {"a": 0.5} j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) npt.assert_almost_equal(j_scaling, 0.5**2, decimal=2) @@ -21,14 +23,16 @@ def test_single_param(self): assert kwargs_min["a"] == 0 param_arrays = np.linspace(0, 1, 11) - scaling_grid_list = [param_arrays ** 2] + scaling_grid_list = [param_arrays**2] param_list = ["a"] - kin_scaling = KinScaling(j_kin_scaling_param_axes=[param_arrays], - j_kin_scaling_grid_list=scaling_grid_list, - j_kin_scaling_param_name_list=param_list) + kin_scaling = KinScaling( + j_kin_scaling_param_axes=[param_arrays], + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list, + ) kwargs_param = {"a": 0.5} j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) - npt.assert_almost_equal(j_scaling, 0.5 ** 2, decimal=2) + npt.assert_almost_equal(j_scaling, 0.5**2, decimal=2) kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() assert kwargs_min["a"] == 0 @@ -41,9 +45,11 @@ def test_two_parameters(self): assert shape_scaling[1] == 21 # assert 1 == 0 param_list = ["a", "b"] - kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=scaling_grid_list, - j_kin_scaling_param_name_list=param_list) + kin_scaling = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list, + ) kwargs_param = {"a": 0.5, "b": 0.3} j_scaling = kin_scaling.kin_scaling(kwargs_param=kwargs_param) print(j_scaling) @@ -59,9 +65,11 @@ def test__kwargs2param_array(self): scaling_grid_list = [xy * uv, xy, uv] param_list = ["a", "b"] - kin_scaling = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=scaling_grid_list, - j_kin_scaling_param_name_list=param_list) + kin_scaling = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=scaling_grid_list, + j_kin_scaling_param_name_list=param_list, + ) kwargs_param = {"a": 0.5, "b": 0.3} param_array = kin_scaling._kwargs2param_array(kwargs_param) assert param_array[0] == kwargs_param["a"] @@ -73,11 +81,16 @@ def test__kwargs2param_array(self): assert kwargs_max["a"] == 1 def test_empty(self): - kin_scaling = KinScaling(j_kin_scaling_param_axes=None, j_kin_scaling_grid_list=None, j_kin_scaling_param_name_list=None) + kin_scaling = KinScaling( + j_kin_scaling_param_axes=None, + j_kin_scaling_grid_list=None, + j_kin_scaling_param_name_list=None, + ) output = kin_scaling.kin_scaling(kwargs_param=None) assert output == 1 kwargs_min, kwargs_max = kin_scaling.param_bounds_interpol() + class TestParameterScalingSingleAperture(object): def setup_method(self): ani_param_array = np.linspace(start=0, stop=1, num=10) @@ -91,7 +104,7 @@ def setup_method(self): np.linspace(start=1, stop=2, num=5), ] param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - print(np.shape(param_scaling_array), 'test shape') + print(np.shape(param_scaling_array), "test shape") self.scaling_2d = ParameterScalingSingleMeasurement( ani_param_array, param_scaling_array ) @@ -153,8 +166,9 @@ def setup_method(self): ani_param_array = np.linspace(start=0, stop=1, num=10) param_scaling_array = ani_param_array * 2 self.scaling = KinScaling( - j_kin_scaling_param_axes=ani_param_array, j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a"] + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a"], ) ani_param_array = [ @@ -162,10 +176,11 @@ def setup_method(self): np.linspace(start=1, stop=2, num=5), ] param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - self.scaling_2d = KinScaling(j_kin_scaling_param_axes=ani_param_array, - j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a", "b"] - ) + self.scaling_2d = KinScaling( + j_kin_scaling_param_axes=ani_param_array, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b"], + ) ani_param_array = np.linspace(start=0, stop=1, num=10) gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) @@ -176,11 +191,11 @@ def setup_method(self): ani_param_array, np.outer(gamma_in_array, log_m2l_array), ) - self.scaling_nfw = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a", "b", "c"] - ) - + self.scaling_nfw = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c"], + ) gom_param_array = [ np.linspace(start=0, stop=1, num=10), @@ -199,21 +214,22 @@ def setup_method(self): np.outer(gamma_in_array, log_m2l_array), ), ) - self.scaling_nfw_2d = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a", "b", "c", "d"] - ) - + self.scaling_nfw_2d = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c", "d"], + ) param_arrays = [ani_param_array, gamma_in_array] param_scaling_array = np.multiply.outer( ani_param_array, gamma_in_array, ) - self.scaling_nfw_no_m2l = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a", "b"] - ) + self.scaling_nfw_no_m2l = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b"], + ) gom_param_array = [ np.linspace(start=0, stop=1, num=10), @@ -231,17 +247,18 @@ def setup_method(self): gamma_in_array, ), ) - self.scaling_nfw_2d_no_m2l = KinScaling(j_kin_scaling_param_axes=param_arrays, - j_kin_scaling_grid_list=[param_scaling_array], - j_kin_scaling_param_name_list=["a", "b", "c"] - ) + self.scaling_nfw_2d_no_m2l = KinScaling( + j_kin_scaling_param_axes=param_arrays, + j_kin_scaling_grid_list=[param_scaling_array], + j_kin_scaling_param_name_list=["a", "b", "c"], + ) def test_kin_scaling(self): scaling = self.scaling.kin_scaling(kwargs_param=None) assert scaling[0] == 1 - scaling = self.scaling.kin_scaling(kwargs_param={"a":1}) + scaling = self.scaling.kin_scaling(kwargs_param={"a": 1}) assert scaling[0] == 2 kwargs_param = {"a": 1, "b": 2} @@ -252,7 +269,7 @@ def test_kin_scaling(self): scaling = self.scaling_nfw.kin_scaling(kwargs_param=kwargs_param) assert scaling[0] == 1 * 2.9 * 0.5 - kwargs_param = {"a": 1, "b": 2., "c": 2.9, "d": 0.5} + kwargs_param = {"a": 1, "b": 2.0, "c": 2.9, "d": 0.5} scaling = self.scaling_nfw_2d.kin_scaling(kwargs_param=kwargs_param) assert scaling[0] == 1 * 2 * 2.9 * 0.5 diff --git a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py index 6dd65ce8..55b954f5 100644 --- a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py +++ b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py @@ -1,4 +1,6 @@ -from hierarc.Sampling.Distributions.anisotropy_distributions import AnisotropyDistribution +from hierarc.Sampling.Distributions.anisotropy_distributions import ( + AnisotropyDistribution, +) class TestAnisotropyDistribution(object): @@ -9,28 +11,39 @@ def setup_method(self): kwargs_anisotropy_min = {"a_ani": 0, "beta_inf": 0.1} kwargs_anisotropy_max = {"a_ani": 5, "beta_inf": 1} - self._ani_dist = AnisotropyDistribution(anisotropy_model=anisotropy_model, - anisotropy_sampling=True, - distribution_function=distribution_function, - kwargs_anisotropy_min=kwargs_anisotropy_min, - kwargs_anisotropy_max=kwargs_anisotropy_max) + self._ani_dist = AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=True, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max, + ) - ani_dist = AnisotropyDistribution(anisotropy_model=anisotropy_model, - anisotropy_sampling=False, - distribution_function=distribution_function, - kwargs_anisotropy_min=kwargs_anisotropy_min, - kwargs_anisotropy_max=kwargs_anisotropy_max) + ani_dist = AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=False, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max, + ) def test_draw_anisotropy(self): - kwargs_anisotropy = {"a_ani": 1, "beta_inf": 0.8, "a_ani_sigma": 0.1, "beta_inf_sigma": 0.2} + kwargs_anisotropy = { + "a_ani": 1, + "beta_inf": 0.8, + "a_ani_sigma": 0.1, + "beta_inf_sigma": 0.2, + } kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) assert "a_ani" in kwargs_drawn assert "beta_inf" in kwargs_drawn - ani_dist = AnisotropyDistribution(anisotropy_model="NONE", - anisotropy_sampling=False, - distribution_function="NONE", - kwargs_anisotropy_min=None, - kwargs_anisotropy_max=None) + ani_dist = AnisotropyDistribution( + anisotropy_model="NONE", + anisotropy_sampling=False, + distribution_function="NONE", + kwargs_anisotropy_min=None, + kwargs_anisotropy_max=None, + ) kwargs_drawn = ani_dist.draw_anisotropy() assert "a_ani" not in kwargs_drawn diff --git a/test/test_Sampling/test_Distributions/test_lens_distribution.py b/test/test_Sampling/test_Distributions/test_lens_distribution.py index 6f71712f..f5a30a27 100644 --- a/test/test_Sampling/test_Distributions/test_lens_distribution.py +++ b/test/test_Sampling/test_Distributions/test_lens_distribution.py @@ -6,35 +6,39 @@ class TestLensDistribution(object): def setup_method(self): - self.kwargs_sampling = {"lambda_mst_distribution": "GAUSSIAN", - "gamma_in_sampling": True, - "gamma_in_distribution": "GAUSSIAN", - "log_m2l_sampling": True, - "log_m2l_distribution" : "GAUSSIAN", - "alpha_lambda_sampling": True, - "beta_lambda_sampling": True, - "alpha_gamma_in_sampling": True, - "alpha_log_m2l_sampling": True, - "log_scatter": False, # change for different tests - "mst_ifu": False, # change for different tests - "lambda_scaling_property": 0.1, - "lambda_scaling_property_beta": 0.2, - "kwargs_min": {"gamma_in": 0, "log_m2l": -3}, - "kwargs_max": {"gamma_in": 3, "log_m2l": 3}} - - self.kwargs_lens = {"lambda_mst": 1.1, - "lambda_mst_sigma": 0.1, - "gamma_ppn": 0.9, - "lambda_ifu": 0.5, - "lambda_ifu_sigma": 0.2, - "alpha_lambda": -0.2, - "beta_lambda": 0.3, - "gamma_in": 1.5, - "gamma_in_sigma": 0.2, - "alpha_gamma_in": 0.2, - "log_m2l": 0.6, - "log_m2l_sigma": 0.2, - "alpha_log_m2l": -0.1} + self.kwargs_sampling = { + "lambda_mst_distribution": "GAUSSIAN", + "gamma_in_sampling": True, + "gamma_in_distribution": "GAUSSIAN", + "log_m2l_sampling": True, + "log_m2l_distribution": "GAUSSIAN", + "alpha_lambda_sampling": True, + "beta_lambda_sampling": True, + "alpha_gamma_in_sampling": True, + "alpha_log_m2l_sampling": True, + "log_scatter": False, # change for different tests + "mst_ifu": False, # change for different tests + "lambda_scaling_property": 0.1, + "lambda_scaling_property_beta": 0.2, + "kwargs_min": {"gamma_in": 0, "log_m2l": -3}, + "kwargs_max": {"gamma_in": 3, "log_m2l": 3}, + } + + self.kwargs_lens = { + "lambda_mst": 1.1, + "lambda_mst_sigma": 0.1, + "gamma_ppn": 0.9, + "lambda_ifu": 0.5, + "lambda_ifu_sigma": 0.2, + "alpha_lambda": -0.2, + "beta_lambda": 0.3, + "gamma_in": 1.5, + "gamma_in_sigma": 0.2, + "alpha_gamma_in": 0.2, + "log_m2l": 0.6, + "log_m2l_sigma": 0.2, + "alpha_log_m2l": -0.1, + } def test_draw_lens(self): lens_dist = LensDistribution(**self.kwargs_sampling) diff --git a/test/test_Sampling/test_mcmc_sampling.py b/test/test_Sampling/test_mcmc_sampling.py index e27550ec..22287009 100644 --- a/test/test_Sampling/test_mcmc_sampling.py +++ b/test/test_Sampling/test_mcmc_sampling.py @@ -67,13 +67,13 @@ def test_mcmc_emcee(self): backend = emcee.backends.HDFBackend(backup_filename) kwargs_emcee = {"backend": backend} - kwargs_model = {"ppn_sampling": False, - "lambda_mst_sampling": False, - "lambda_mst_distribution": "delta", - "anisotropy_sampling": False, - "anisotropy_model": "OM", - - } + kwargs_model = { + "ppn_sampling": False, + "lambda_mst_sampling": False, + "lambda_mst_distribution": "delta", + "anisotropy_sampling": False, + "anisotropy_model": "OM", + } mcmc_sampler = MCMCSampler( kwargs_likelihood_list=kwargs_likelihood_list, From 58d9bd4a1ded80cd6455a79ab5581d9ac00bcb14 Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Thu, 27 Jun 2024 14:29:48 -0400 Subject: [PATCH 03/14] minor change in plotting if there is only a single lens --- hierarc/Diagnostics/goodness_of_fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hierarc/Diagnostics/goodness_of_fit.py b/hierarc/Diagnostics/goodness_of_fit.py index c65158d0..79c0b01a 100644 --- a/hierarc/Diagnostics/goodness_of_fit.py +++ b/hierarc/Diagnostics/goodness_of_fit.py @@ -222,7 +222,7 @@ def plot_kin_fit( sigma_v_model_error_list, ) = 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)) + f, ax = plt.subplots(1, 1, figsize=(max(int(len(sigma_v_name_list) / 2), 1), 4)) ax.errorbar( np.arange(len(sigma_v_name_list)), sigma_v_measurement_list, From c0ed4a6578cc3b10c29b688e08c3e2bf24475c5c Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Thu, 27 Jun 2024 14:44:49 -0400 Subject: [PATCH 04/14] improved tests in lens distribution draws --- .../test_lens_distribution.py | 34 ++++++++++++++----- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/test/test_Sampling/test_Distributions/test_lens_distribution.py b/test/test_Sampling/test_Distributions/test_lens_distribution.py index 6f71712f..f1442ca6 100644 --- a/test/test_Sampling/test_Distributions/test_lens_distribution.py +++ b/test/test_Sampling/test_Distributions/test_lens_distribution.py @@ -1,4 +1,5 @@ import copy +import numpy.testing as npt from hierarc.Sampling.Distributions.lens_distribution import LensDistribution @@ -6,11 +7,13 @@ class TestLensDistribution(object): def setup_method(self): - self.kwargs_sampling = {"lambda_mst_distribution": "GAUSSIAN", + self.kwargs_sampling = { + "lambda_mst_sampling": True, + "lambda_mst_distribution": "GAUSSIAN", "gamma_in_sampling": True, "gamma_in_distribution": "GAUSSIAN", "log_m2l_sampling": True, - "log_m2l_distribution" : "GAUSSIAN", + "log_m2l_distribution": "GAUSSIAN", "alpha_lambda_sampling": True, "beta_lambda_sampling": True, "alpha_gamma_in_sampling": True, @@ -19,8 +22,8 @@ def setup_method(self): "mst_ifu": False, # change for different tests "lambda_scaling_property": 0.1, "lambda_scaling_property_beta": 0.2, - "kwargs_min": {"gamma_in": 0, "log_m2l": -3}, - "kwargs_max": {"gamma_in": 3, "log_m2l": 3}} + "kwargs_min": {"gamma_in": 1, "log_m2l": 0}, + "kwargs_max": {"gamma_in": 2, "log_m2l": 1}} self.kwargs_lens = {"lambda_mst": 1.1, "lambda_mst_sigma": 0.1, @@ -30,10 +33,10 @@ def setup_method(self): "alpha_lambda": -0.2, "beta_lambda": 0.3, "gamma_in": 1.5, - "gamma_in_sigma": 0.2, + "gamma_in_sigma": 1, "alpha_gamma_in": 0.2, "log_m2l": 0.6, - "log_m2l_sigma": 0.2, + "log_m2l_sigma": 1, "alpha_log_m2l": -0.1} def test_draw_lens(self): @@ -46,6 +49,21 @@ def test_draw_lens(self): kwargs_sampling["log_scatter"] = True kwargs_sampling["lambda_ifu"] = True lens_dist = LensDistribution(kwargs_sampling) - kwargs_return = lens_dist.draw_lens(**self.kwargs_lens) + for i in range(100): + kwargs_return = lens_dist.draw_lens(**self.kwargs_lens) - assert "lambda_mst" in kwargs_return + assert "lambda_mst" in kwargs_return + + def test_raises(self): + + with npt.assert_raises(ValueError): + lens_dist = LensDistribution(**self.kwargs_sampling) + kwargs_lens = copy.deepcopy(self.kwargs_lens) + kwargs_lens["gamma_in"] = -10 + kwargs_return = lens_dist.draw_lens(**kwargs_lens) + + with npt.assert_raises(ValueError): + lens_dist = LensDistribution(**self.kwargs_sampling) + kwargs_lens = copy.deepcopy(self.kwargs_lens) + kwargs_lens["log_m2l"] = -100 + kwargs_return = lens_dist.draw_lens(**kwargs_lens) From 567b0d458364777eb9d983611c559179fe870447 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jun 2024 18:45:59 +0000 Subject: [PATCH 05/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- .../test_lens_distribution.py | 61 ++++++++++--------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/test/test_Sampling/test_Distributions/test_lens_distribution.py b/test/test_Sampling/test_Distributions/test_lens_distribution.py index f1442ca6..bb9513aa 100644 --- a/test/test_Sampling/test_Distributions/test_lens_distribution.py +++ b/test/test_Sampling/test_Distributions/test_lens_distribution.py @@ -8,36 +8,39 @@ class TestLensDistribution(object): def setup_method(self): self.kwargs_sampling = { - "lambda_mst_sampling": True, - "lambda_mst_distribution": "GAUSSIAN", - "gamma_in_sampling": True, - "gamma_in_distribution": "GAUSSIAN", - "log_m2l_sampling": True, - "log_m2l_distribution": "GAUSSIAN", - "alpha_lambda_sampling": True, - "beta_lambda_sampling": True, - "alpha_gamma_in_sampling": True, - "alpha_log_m2l_sampling": True, - "log_scatter": False, # change for different tests - "mst_ifu": False, # change for different tests - "lambda_scaling_property": 0.1, - "lambda_scaling_property_beta": 0.2, - "kwargs_min": {"gamma_in": 1, "log_m2l": 0}, - "kwargs_max": {"gamma_in": 2, "log_m2l": 1}} + "lambda_mst_sampling": True, + "lambda_mst_distribution": "GAUSSIAN", + "gamma_in_sampling": True, + "gamma_in_distribution": "GAUSSIAN", + "log_m2l_sampling": True, + "log_m2l_distribution": "GAUSSIAN", + "alpha_lambda_sampling": True, + "beta_lambda_sampling": True, + "alpha_gamma_in_sampling": True, + "alpha_log_m2l_sampling": True, + "log_scatter": False, # change for different tests + "mst_ifu": False, # change for different tests + "lambda_scaling_property": 0.1, + "lambda_scaling_property_beta": 0.2, + "kwargs_min": {"gamma_in": 1, "log_m2l": 0}, + "kwargs_max": {"gamma_in": 2, "log_m2l": 1}, + } - self.kwargs_lens = {"lambda_mst": 1.1, - "lambda_mst_sigma": 0.1, - "gamma_ppn": 0.9, - "lambda_ifu": 0.5, - "lambda_ifu_sigma": 0.2, - "alpha_lambda": -0.2, - "beta_lambda": 0.3, - "gamma_in": 1.5, - "gamma_in_sigma": 1, - "alpha_gamma_in": 0.2, - "log_m2l": 0.6, - "log_m2l_sigma": 1, - "alpha_log_m2l": -0.1} + self.kwargs_lens = { + "lambda_mst": 1.1, + "lambda_mst_sigma": 0.1, + "gamma_ppn": 0.9, + "lambda_ifu": 0.5, + "lambda_ifu_sigma": 0.2, + "alpha_lambda": -0.2, + "beta_lambda": 0.3, + "gamma_in": 1.5, + "gamma_in_sigma": 1, + "alpha_gamma_in": 0.2, + "log_m2l": 0.6, + "log_m2l_sigma": 1, + "alpha_log_m2l": -0.1, + } def test_draw_lens(self): lens_dist = LensDistribution(**self.kwargs_sampling) From caa13ded9678e8b290e67a028f2cbcdbfde4b24b Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Thu, 27 Jun 2024 21:50:29 -0400 Subject: [PATCH 06/14] improve testing and stability --- .../Distributions/anisotropy_distributions.py | 2 +- .../test_anisotropy_distribution.py | 42 +++++++++++++++++++ 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py index 1050a1fb..0b66627b 100644 --- a/hierarc/Sampling/Distributions/anisotropy_distributions.py +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -66,7 +66,7 @@ def draw_anisotropy( else: a_ani_draw = np.random.normal(a_ani, a_ani_sigma) if a_ani_draw < self._a_ani_min or a_ani_draw > self._a_ani_max: - return self.draw_anisotropy(a_ani, a_ani_sigma) + return self.draw_anisotropy(a_ani, a_ani_sigma, beta_inf, beta_inf_sigma) kwargs_return["a_ani"] = a_ani_draw else: kwargs_return["a_ani"] = a_ani diff --git a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py index 55b954f5..70e9324e 100644 --- a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py +++ b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py @@ -1,3 +1,4 @@ +import numpy.testing as npt from hierarc.Sampling.Distributions.anisotropy_distributions import ( AnisotropyDistribution, ) @@ -47,3 +48,44 @@ def test_draw_anisotropy(self): ) kwargs_drawn = ani_dist.draw_anisotropy() assert "a_ani" not in kwargs_drawn + + ani_dist = AnisotropyDistribution( + anisotropy_model="GOM", + anisotropy_sampling=True, + distribution_function="NONE", + kwargs_anisotropy_min=None, + kwargs_anisotropy_max=None, + ) + kwargs_drawn = ani_dist.draw_anisotropy(a_ani=1, beta_inf=0.9) + assert kwargs_drawn["a_ani"] == 1 + assert kwargs_drawn["beta_inf"] == 0.9 + + kwargs_anisotropy = { + "a_ani": 1, + "beta_inf": 0.8, + "a_ani_sigma": 2, + "beta_inf_sigma": 2, + } + + for i in range(100): + kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) + + def test_raises(self): + + with npt.assert_raises(ValueError): + kwargs_anisotropy = { + "a_ani": -1, + "beta_inf": 0.8, + "a_ani_sigma": 0.1, + "beta_inf_sigma": 0.2, + } + kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) + + with npt.assert_raises(ValueError): + kwargs_anisotropy = { + "a_ani": 1, + "beta_inf": -1, + "a_ani_sigma": 0.1, + "beta_inf_sigma": 0.2, + } + kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) From f5c5b13bd7abacf2917b98fc80d1d8897d20b848 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 01:50:42 +0000 Subject: [PATCH 07/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- hierarc/Sampling/Distributions/anisotropy_distributions.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py index 0b66627b..231c9f76 100644 --- a/hierarc/Sampling/Distributions/anisotropy_distributions.py +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -66,7 +66,9 @@ def draw_anisotropy( else: a_ani_draw = np.random.normal(a_ani, a_ani_sigma) if a_ani_draw < self._a_ani_min or a_ani_draw > self._a_ani_max: - return self.draw_anisotropy(a_ani, a_ani_sigma, beta_inf, beta_inf_sigma) + return self.draw_anisotropy( + a_ani, a_ani_sigma, beta_inf, beta_inf_sigma + ) kwargs_return["a_ani"] = a_ani_draw else: kwargs_return["a_ani"] = a_ani From 4758eab4d6e54caa63a4d5f0e8de2b7f6adf2660 Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Fri, 28 Jun 2024 18:15:21 -0400 Subject: [PATCH 08/14] re-definition of distribution functions for OM and GOM model for a_ani with scaled function separately --- .../Distributions/anisotropy_distributions.py | 29 ++++++++--- hierarc/Sampling/ParamManager/kin_param.py | 4 +- .../test_anisotropy_distribution.py | 50 ++++++++++++++++++- 3 files changed, 71 insertions(+), 12 deletions(-) diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py index 0b66627b..0aec213d 100644 --- a/hierarc/Sampling/Distributions/anisotropy_distributions.py +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -1,5 +1,8 @@ import numpy as np +_SUPPORTED_DISTRIBUTIONS = ["GAUSSIAN", "GAUSSIAN_SCALED", "NONE"] +_SUPPORTED_MODELS = ["OM", "GOM", "const", "NONE"] + class AnisotropyDistribution(object): """Class to draw anisotropy parameters from hyperparameter distributions.""" @@ -17,12 +20,20 @@ def __init__( :param anisotropy_model: string, name of anisotropy model to consider :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens kinematic prediction - :param distribution_function: string, 'NONE', 'GAUSSIAN', description of the distribution function of the - anisotropy model parameters + :param distribution_function: string, 'NONE', 'GAUSSIAN', 'GAUSSIAN_SCALED', + description of the distribution function of the anisotropy model parameters """ - + if anisotropy_model not in _SUPPORTED_MODELS: + raise ValueError("Anisotropy model %s not supported. Chose among %s." + % (anisotropy_model, _SUPPORTED_MODELS)) self._anisotropy_model = anisotropy_model self._anisotropy_sampling = anisotropy_sampling + if distribution_function not in _SUPPORTED_DISTRIBUTIONS: + raise ValueError("Anisotropy distribution function %s not supported. Chose among %s." + % (distribution_function, _SUPPORTED_DISTRIBUTIONS)) + if anisotropy_model not in ["OM", "GOM"] and distribution_function == "GAUSSIAN_SCALED": + raise ValueError("GAUSSIAN_SCALED distribution only supported for 'OM' and 'GOM' models, not for %s." + % anisotropy_model) self._distribution_function = distribution_function if kwargs_anisotropy_min is None: kwargs_anisotropy_min = {} @@ -60,22 +71,24 @@ def draw_anisotropy( "anisotropy parameter is out of bounds of the interpolated range!" ) # we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM - if self._distribution_function in ["GAUSSIAN"]: - if self._anisotropy_model == "OM": - a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani) - else: + if self._distribution_function in ["GAUSSIAN", "GAUSSIAN_SCALED"]: + if self._distribution_function in ["GAUSSIAN"]: a_ani_draw = np.random.normal(a_ani, a_ani_sigma) + else: + a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani) + if a_ani_draw < self._a_ani_min or a_ani_draw > self._a_ani_max: return self.draw_anisotropy(a_ani, a_ani_sigma, beta_inf, beta_inf_sigma) kwargs_return["a_ani"] = a_ani_draw else: kwargs_return["a_ani"] = a_ani + if self._anisotropy_model in ["GOM"]: if beta_inf < self._beta_inf_min or beta_inf > self._beta_inf_max: raise ValueError( "anisotropy parameter is out of bounds of the interpolated range!" ) - if self._distribution_function in ["GAUSSIAN"]: + if self._distribution_function in ["GAUSSIAN", "GAUSSIAN_SCALED"]: beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma) else: beta_inf_draw = beta_inf diff --git a/hierarc/Sampling/ParamManager/kin_param.py b/hierarc/Sampling/ParamManager/kin_param.py index d8023c76..155e2a6b 100644 --- a/hierarc/Sampling/ParamManager/kin_param.py +++ b/hierarc/Sampling/ParamManager/kin_param.py @@ -49,7 +49,7 @@ def param_list(self, latex_style=False): list.append(r"$\langle a_{\rm ani}\rangle$") else: list.append("a_ani") - if self._distribution_function in ["GAUSSIAN"]: + if self._distribution_function in ["GAUSSIAN", "GAUSSIAN_SCALED"]: if "a_ani_sigma" not in self._kwargs_fixed: if latex_style is True: if self._log_scatter is True: @@ -64,7 +64,7 @@ def param_list(self, latex_style=False): list.append(r"$\beta_{\infty}$") else: list.append("beta_inf") - if self._distribution_function in ["GAUSSIAN"]: + if self._distribution_function in ["GAUSSIAN", "GAUSSIAN_SCALED"]: if "beta_inf_sigma" not in self._kwargs_fixed: if latex_style is True: if self._log_scatter is True: diff --git a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py index 70e9324e..218bb52d 100644 --- a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py +++ b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py @@ -20,10 +20,10 @@ def setup_method(self): kwargs_anisotropy_max=kwargs_anisotropy_max, ) - ani_dist = AnisotropyDistribution( + self._ani_dist_scaled = AnisotropyDistribution( anisotropy_model=anisotropy_model, anisotropy_sampling=False, - distribution_function=distribution_function, + distribution_function="GAUSSIAN_SCALED", kwargs_anisotropy_min=kwargs_anisotropy_min, kwargs_anisotropy_max=kwargs_anisotropy_max, ) @@ -39,6 +39,10 @@ def test_draw_anisotropy(self): assert "a_ani" in kwargs_drawn assert "beta_inf" in kwargs_drawn + kwargs_drawn = self._ani_dist_scaled.draw_anisotropy(**kwargs_anisotropy) + assert "a_ani" in kwargs_drawn + assert "beta_inf" in kwargs_drawn + ani_dist = AnisotropyDistribution( anisotropy_model="NONE", anisotropy_sampling=False, @@ -89,3 +93,45 @@ def test_raises(self): "beta_inf_sigma": 0.2, } kwargs_drawn = self._ani_dist.draw_anisotropy(**kwargs_anisotropy) + + with npt.assert_raises(ValueError): + anisotropy_model = "const" + distribution_function = "GAUSSIAN_SCALED" + kwargs_anisotropy_min = {"a_ani": 0, "beta_inf": 0.1} + kwargs_anisotropy_max = {"a_ani": 5, "beta_inf": 1} + + AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=True, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max, + ) + + with npt.assert_raises(ValueError): + anisotropy_model = "const" + distribution_function = "INVALID" + kwargs_anisotropy_min = {"a_ani": 0, "beta_inf": 0.1} + kwargs_anisotropy_max = {"a_ani": 5, "beta_inf": 1} + + AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=True, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max, + ) + + with npt.assert_raises(ValueError): + anisotropy_model = "INVALID" + distribution_function = "GAUSSIAN" + kwargs_anisotropy_min = {"a_ani": 0, "beta_inf": 0.1} + kwargs_anisotropy_max = {"a_ani": 5, "beta_inf": 1} + + AnisotropyDistribution( + anisotropy_model=anisotropy_model, + anisotropy_sampling=True, + distribution_function=distribution_function, + kwargs_anisotropy_min=kwargs_anisotropy_min, + kwargs_anisotropy_max=kwargs_anisotropy_max, + ) From 73a3464a666edfeb12b2d10f4a093cf77c6f2dd1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 22:15:59 +0000 Subject: [PATCH 09/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- .../Distributions/anisotropy_distributions.py | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/hierarc/Sampling/Distributions/anisotropy_distributions.py b/hierarc/Sampling/Distributions/anisotropy_distributions.py index cf83256b..abc97388 100644 --- a/hierarc/Sampling/Distributions/anisotropy_distributions.py +++ b/hierarc/Sampling/Distributions/anisotropy_distributions.py @@ -24,16 +24,25 @@ def __init__( description of the distribution function of the anisotropy model parameters """ if anisotropy_model not in _SUPPORTED_MODELS: - raise ValueError("Anisotropy model %s not supported. Chose among %s." - % (anisotropy_model, _SUPPORTED_MODELS)) + raise ValueError( + "Anisotropy model %s not supported. Chose among %s." + % (anisotropy_model, _SUPPORTED_MODELS) + ) self._anisotropy_model = anisotropy_model self._anisotropy_sampling = anisotropy_sampling if distribution_function not in _SUPPORTED_DISTRIBUTIONS: - raise ValueError("Anisotropy distribution function %s not supported. Chose among %s." - % (distribution_function, _SUPPORTED_DISTRIBUTIONS)) - if anisotropy_model not in ["OM", "GOM"] and distribution_function == "GAUSSIAN_SCALED": - raise ValueError("GAUSSIAN_SCALED distribution only supported for 'OM' and 'GOM' models, not for %s." - % anisotropy_model) + raise ValueError( + "Anisotropy distribution function %s not supported. Chose among %s." + % (distribution_function, _SUPPORTED_DISTRIBUTIONS) + ) + if ( + anisotropy_model not in ["OM", "GOM"] + and distribution_function == "GAUSSIAN_SCALED" + ): + raise ValueError( + "GAUSSIAN_SCALED distribution only supported for 'OM' and 'GOM' models, not for %s." + % anisotropy_model + ) self._distribution_function = distribution_function if kwargs_anisotropy_min is None: kwargs_anisotropy_min = {} From 506df698aa21c517b9e6feee8e04f4733739cb0b Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Fri, 28 Jun 2024 18:30:11 -0400 Subject: [PATCH 10/14] minor testing changes --- test/test_LensPosterior/test_kin_scaling_config.py | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 test/test_LensPosterior/test_kin_scaling_config.py diff --git a/test/test_LensPosterior/test_kin_scaling_config.py b/test/test_LensPosterior/test_kin_scaling_config.py new file mode 100644 index 00000000..aed0bcf7 --- /dev/null +++ b/test/test_LensPosterior/test_kin_scaling_config.py @@ -0,0 +1,9 @@ +from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig + +class TestKinScalingConfig(object): + + def setup_method(self): + pass + + def test_init(self): + KinScalingConfig(anisotropy_model="NONE", r_eff=None, gamma_in_scaling=None, log_m2l_scaling=None) From 7ce952a31080248cfd925366efeebe4250c911a7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 22:30:34 +0000 Subject: [PATCH 11/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- test/test_LensPosterior/test_kin_scaling_config.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/test_LensPosterior/test_kin_scaling_config.py b/test/test_LensPosterior/test_kin_scaling_config.py index aed0bcf7..48a2df92 100644 --- a/test/test_LensPosterior/test_kin_scaling_config.py +++ b/test/test_LensPosterior/test_kin_scaling_config.py @@ -1,9 +1,15 @@ from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig + class TestKinScalingConfig(object): def setup_method(self): pass def test_init(self): - KinScalingConfig(anisotropy_model="NONE", r_eff=None, gamma_in_scaling=None, log_m2l_scaling=None) + KinScalingConfig( + anisotropy_model="NONE", + r_eff=None, + gamma_in_scaling=None, + log_m2l_scaling=None, + ) From 36c98e42cd7598c7aa290c445775009158712a1b Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Fri, 28 Jun 2024 18:33:05 -0400 Subject: [PATCH 12/14] minor testing changes --- hierarc/LensPosterior/kin_scaling_config.py | 2 +- test/test_LensPosterior/test_kin_scaling_config.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/hierarc/LensPosterior/kin_scaling_config.py b/hierarc/LensPosterior/kin_scaling_config.py index d6979f3d..efb927f7 100644 --- a/hierarc/LensPosterior/kin_scaling_config.py +++ b/hierarc/LensPosterior/kin_scaling_config.py @@ -87,7 +87,7 @@ def param_name_list(self): """ return self._param_name_list - def anisotropy_kwargs(self, a_ani, beta_inf=None): + def anisotropy_kwargs(self, a_ani=None, beta_inf=None): """ :param a_ani: anisotropy parameter diff --git a/test/test_LensPosterior/test_kin_scaling_config.py b/test/test_LensPosterior/test_kin_scaling_config.py index aed0bcf7..01a3f773 100644 --- a/test/test_LensPosterior/test_kin_scaling_config.py +++ b/test/test_LensPosterior/test_kin_scaling_config.py @@ -1,4 +1,5 @@ from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig +import numpy.testing as npt class TestKinScalingConfig(object): @@ -6,4 +7,8 @@ def setup_method(self): pass def test_init(self): - KinScalingConfig(anisotropy_model="NONE", r_eff=None, gamma_in_scaling=None, log_m2l_scaling=None) + kin_scaling = KinScalingConfig(anisotropy_model="NONE", r_eff=None, gamma_in_scaling=None, log_m2l_scaling=None) + kin_scaling._anisotropy_model = "BAD" + + with npt.assert_raises(ValueError): + kin_scaling.anisotropy_kwargs() From 3e6a1e1fb4f096a1056cb8c75004846f1dbef9a5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 22:33:45 +0000 Subject: [PATCH 13/14] [pre-commit.ci] auto fixes from pre-commit.com hooks --- test/test_LensPosterior/test_kin_scaling_config.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/test/test_LensPosterior/test_kin_scaling_config.py b/test/test_LensPosterior/test_kin_scaling_config.py index e8e39853..20677b8c 100644 --- a/test/test_LensPosterior/test_kin_scaling_config.py +++ b/test/test_LensPosterior/test_kin_scaling_config.py @@ -8,7 +8,12 @@ def setup_method(self): pass def test_init(self): - kin_scaling = KinScalingConfig(anisotropy_model="NONE", r_eff=None, gamma_in_scaling=None, log_m2l_scaling=None) + kin_scaling = KinScalingConfig( + anisotropy_model="NONE", + r_eff=None, + gamma_in_scaling=None, + log_m2l_scaling=None, + ) kin_scaling._anisotropy_model = "BAD" with npt.assert_raises(ValueError): From a68ea3e54f7a2ed1a267546415b90e845de2a54b Mon Sep 17 00:00:00 2001 From: Simon Birrer Date: Fri, 28 Jun 2024 21:10:07 -0400 Subject: [PATCH 14/14] minor testing changes --- test/test_Likelihood/test_kin_scaling.py | 4 ++++ .../test_Distributions/test_anisotropy_distribution.py | 2 +- .../test_Distributions/test_lens_distribution.py | 1 + 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/test/test_Likelihood/test_kin_scaling.py b/test/test_Likelihood/test_kin_scaling.py index f1201660..9b8d1c63 100644 --- a/test/test_Likelihood/test_kin_scaling.py +++ b/test/test_Likelihood/test_kin_scaling.py @@ -80,6 +80,10 @@ def test__kwargs2param_array(self): assert kwargs_max["b"] == 2 assert kwargs_max["a"] == 1 + with npt.assert_raises(ValueError): + kwargs_param = {"a": 0.5} # remove parameter "b" and expect a raise + param_array = kin_scaling._kwargs2param_array(kwargs_param) + def test_empty(self): kin_scaling = KinScaling( j_kin_scaling_param_axes=None, diff --git a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py index 218bb52d..51683fca 100644 --- a/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py +++ b/test/test_Sampling/test_Distributions/test_anisotropy_distribution.py @@ -22,7 +22,7 @@ def setup_method(self): self._ani_dist_scaled = AnisotropyDistribution( anisotropy_model=anisotropy_model, - anisotropy_sampling=False, + anisotropy_sampling=True, distribution_function="GAUSSIAN_SCALED", kwargs_anisotropy_min=kwargs_anisotropy_min, kwargs_anisotropy_max=kwargs_anisotropy_max, diff --git a/test/test_Sampling/test_Distributions/test_lens_distribution.py b/test/test_Sampling/test_Distributions/test_lens_distribution.py index bb9513aa..db250ea6 100644 --- a/test/test_Sampling/test_Distributions/test_lens_distribution.py +++ b/test/test_Sampling/test_Distributions/test_lens_distribution.py @@ -51,6 +51,7 @@ def test_draw_lens(self): kwargs_sampling = copy.deepcopy(self.kwargs_sampling) kwargs_sampling["log_scatter"] = True kwargs_sampling["lambda_ifu"] = True + kwargs_sampling["gamma_in_sampling"] = False lens_dist = LensDistribution(kwargs_sampling) for i in range(100): kwargs_return = lens_dist.draw_lens(**self.kwargs_lens)