diff --git a/hierarc/LensPosterior/base_config.py b/hierarc/LensPosterior/base_config.py index 832cca6..9fef09e 100644 --- a/hierarc/LensPosterior/base_config.py +++ b/hierarc/LensPosterior/base_config.py @@ -22,6 +22,7 @@ def __init__( kwargs_seeing, kwargs_numerics_galkin, anisotropy_model, + lens_model_list=None, kwargs_lens_light=None, lens_light_model_list=["HERNQUIST"], MGE_light=False, @@ -31,6 +32,7 @@ def __init__( num_psf_sampling=100, num_kin_sampling=1000, multi_observations=False, + cosmo_fiducial=None, ): """ @@ -42,27 +44,43 @@ def __init__( :param gamma_error: 1-sigma uncertainty on power-law slope :param r_eff: half-light radius of the deflector (arc seconds) :param r_eff_error: uncertainty on half-light radius - :param kwargs_aperture: spectroscopic aperture keyword arguments, see lenstronomy.Galkin.aperture for options - :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf - :param kwargs_numerics_galkin: numerical settings for the integrated line-of-sight velocity dispersion - :param anisotropy_model: type of stellar anisotropy model. See details in MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy - :param multi_observations: bool, if True, interprets kwargs_aperture and kwargs_seeing as lists of multiple - observations + :param kwargs_aperture: spectroscopic aperture keyword arguments, see + lenstronomy.Galkin.aperture for options + :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds + to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf + :param kwargs_numerics_galkin: numerical settings for the integrated + line-of-sight velocity dispersion + :param anisotropy_model: type of stellar anisotropy model. See details in + MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy + :param multi_observations: bool, if True, interprets kwargs_aperture and + kwargs_seeing as lists of multiple observations + :param lens_model_list: keyword argument list of lens model (optional) :param kwargs_lens_light: keyword argument list of lens light model (optional) - :param kwargs_mge_light: keyword arguments that go into the MGE decomposition routine - :param hernquist_approx: bool, if True, uses the Hernquist approximation for the light profile + :param kwargs_mge_light: keyword arguments that go into the MGE decomposition + routine + :param hernquist_approx: bool, if True, uses the Hernquist approximation for the + light profile + :param cosmo_fiducial: astropy.cosmology instance, if None, + uses astropy's default cosmology """ self._z_lens, self._z_source = z_lens, z_source - kwargs_model = { - "lens_model_list": ["SPP"], - "lens_light_model_list": lens_light_model_list, - } + + if lens_model_list is None: + kwargs_model = { + "lens_model_list": ["SPP"], + "lens_light_model_list": lens_light_model_list, + } + else: + kwargs_model = { + "lens_model_list": lens_model_list, + "lens_light_model_list": lens_light_model_list, + } TDCosmography.__init__( self, z_lens, z_source, kwargs_model, - cosmo_fiducial=None, + cosmo_fiducial=cosmo_fiducial, lens_model_kinematics_bool=None, light_model_kinematics_bool=None, kwargs_seeing=kwargs_seeing, diff --git a/hierarc/LensPosterior/kin_constraints.py b/hierarc/LensPosterior/kin_constraints.py index fe1e28a..74147ff 100644 --- a/hierarc/LensPosterior/kin_constraints.py +++ b/hierarc/LensPosterior/kin_constraints.py @@ -27,6 +27,7 @@ def __init__( sigma_v_error_cov_matrix=None, kwargs_lens_light=None, lens_light_model_list=["HERNQUIST"], + lens_model_list=None, MGE_light=False, kwargs_mge_light=None, hernquist_approx=True, @@ -34,6 +35,7 @@ def __init__( num_psf_sampling=100, num_kin_sampling=1000, multi_observations=False, + cosmo_fiducial=None, ): """ @@ -45,23 +47,36 @@ def __init__( :param gamma_error: 1-sigma uncertainty on power-law slope :param r_eff: half-light radius of the deflector (arc seconds) :param r_eff_error: uncertainty on half-light radius - :param sigma_v_measured: numpy array of IFU velocity dispersion of the main deflector in km/s - :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity dispersion of the IFU + :param sigma_v_measured: numpy array of IFU velocity dispersion of the main + deflector in km/s + :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity + dispersion of the IFU observation independent of each other - :param sigma_v_error_covariant: covariant error in the measured kinematics shared among all IFU measurements - :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v measurements (km/s)^2 - :type sigma_v_error_cov_matrix: nxn matrix with n the length of the sigma_v_measured array - :param kwargs_aperture: spectroscopic aperture keyword arguments, see lenstronomy.Galkin.aperture for options - :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds to kwargs_psf in the GalKin + :param sigma_v_error_covariant: covariant error in the measured kinematics + shared among all IFU measurements + :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v + measurements (km/s)^2 + :type sigma_v_error_cov_matrix: nxn matrix with n the length of the + sigma_v_measured array + :param kwargs_aperture: spectroscopic aperture keyword arguments, see + lenstronomy.Galkin.aperture for options + :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds + to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf - :param kwargs_numerics_galkin: numerical settings for the integrated line-of-sight velocity dispersion - :param anisotropy_model: type of stellar anisotropy model. See details in MamonLokasAnisotropy() class of - lenstronomy.GalKin.anisotropy + :param kwargs_numerics_galkin: numerical settings for the integrated + line-of-sight velocity dispersion + :param anisotropy_model: type of stellar anisotropy model. See details in + MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy + :param lens_model_list: keyword argument list of lens model (optional) :param kwargs_lens_light: keyword argument list of lens light model (optional) - :param kwargs_mge_light: keyword arguments that go into the MGE decomposition routine - :param hernquist_approx: bool, if True, uses the Hernquist approximation for the light profile - :param multi_observations: bool, if True, interprets kwargs_aperture and kwargs_seeing as lists of multiple - observations + :param kwargs_mge_light: keyword arguments that go into the MGE decomposition + routine + :param hernquist_approx: bool, if True, uses the Hernquist approximation for the + light profile + :param multi_observations: bool, if True, interprets kwargs_aperture and + kwargs_seeing as lists of multiple observations + :param cosmo_fiducial: astropy.cosmology instance, if None, + uses astropy's default """ self._sigma_v_measured = np.array(sigma_v_measured) self._sigma_v_error_independent = np.array(sigma_v_error_independent) @@ -70,6 +85,7 @@ def __init__( self._kwargs_lens_light = kwargs_lens_light self._anisotropy_model = anisotropy_model + BaseLensConfig.__init__( self, z_lens, @@ -84,6 +100,7 @@ def __init__( kwargs_seeing, kwargs_numerics_galkin, anisotropy_model, + lens_model_list=lens_model_list, kwargs_lens_light=kwargs_lens_light, lens_light_model_list=lens_light_model_list, MGE_light=MGE_light, @@ -93,6 +110,7 @@ def __init__( num_psf_sampling=num_psf_sampling, num_kin_sampling=num_kin_sampling, multi_observations=multi_observations, + cosmo_fiducial=cosmo_fiducial, ) def j_kin_draw(self, kwargs_anisotropy, no_error=False): @@ -161,9 +179,10 @@ def hierarchy_configuration(self, num_sample_model=20): def model_marginalization(self, num_sample_model=20): """ - :param num_sample_model: number of samples drawn from the lens and light model posterior to compute the - dimensionless kinematic component J() - :return: J() as array for each measurement prediction, covariance matrix in sqrt(J) + :param num_sample_model: number of samples drawn from the lens and light model + posterior to compute the dimensionless kinematic component J() + :return: J() as array for each measurement prediction, covariance matrix in + sqrt(J) """ num_data = len(self._sigma_v_measured) j_kin_matrix = np.zeros( diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py new file mode 100644 index 0000000..f7b2e77 --- /dev/null +++ b/hierarc/LensPosterior/kin_constraints_composite.py @@ -0,0 +1,355 @@ +__author__ = "ajshjib" + +import copy + +import numpy as np +from hierarc.LensPosterior.kin_constraints import KinConstraints +from lenstronomy.Util import constants as const +from lenstronomy.Analysis.light_profile import LightProfileAnalysis +from lenstronomy.LightModel.light_model import LightModel + + +class KinConstraintsComposite(KinConstraints): + def __init__( + self, + z_lens, + z_source, + gamma_in_array, + m2l_array, + rho0_array, + r_s_array, + theta_E, + theta_E_error, + gamma, + gamma_error, + r_eff, + r_eff_error, + sigma_v_measured, + kwargs_aperture, + kwargs_seeing, + kwargs_numerics_galkin, + anisotropy_model, + sigma_v_error_independent=None, + sigma_v_error_covariant=None, + sigma_v_error_cov_matrix=None, + kwargs_lens_light=None, + lens_light_model_list=["HERNQUIST"], + kwargs_mge_light=None, + sampling_number=1000, + num_psf_sampling=100, + num_kin_sampling=1000, + multi_observations=False, + ): + """ + + :param z_lens: lens redshift + :param z_source: source redshift + :param gamma_in_array: array of power-law slopes of the mass model + :param m2l_array: array of mass-to-light ratios of the stellar component, + needs to be in the unit/scaling such that m2l * amp in the + kwargs_lens_light provides the convergence amplitude of the stars + :param rho0_array: array of halo mass normalizations in M_sun / Mpc^3 + :param r_s_array: array of halo scale radii in arc seconds + :param theta_E: Einstein radius (in arc seconds) + :param theta_E_error: 1-sigma error on Einstein radius + :param gamma: power-law slope (2 = isothermal) estimated from imaging data + :param gamma_error: 1-sigma uncertainty on power-law slope + :param r_eff: half-light radius of the deflector (arc seconds) + :param r_eff_error: uncertainty on half-light radius + :param sigma_v_measured: numpy array of IFU velocity dispersion of the main + deflector in km/s + :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity + dispersion of the IFU + observation independent of each other + :param sigma_v_error_covariant: covariant error in the measured kinematics + shared among all IFU measurements + :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v + measurements (km/s)^2 + :type sigma_v_error_cov_matrix: nxn matrix with n the length of the + sigma_v_measured array + :param kwargs_aperture: spectroscopic aperture keyword arguments, see + lenstronomy.Galkin.aperture for options + :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds + to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf + :param kwargs_numerics_galkin: numerical settings for the integrated + line-of-sight velocity dispersion + :param anisotropy_model: type of stellar anisotropy model. See details in + MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy + :param kwargs_lens_light: keyword argument list of lens light model (optional) + :param kwargs_mge_light: keyword arguments that go into the MGE decomposition + routine + :param multi_observations: bool, if True, interprets kwargs_aperture and + kwargs_seeing as lists of multiple observations + """ + self._light_profile_analysis = LightProfileAnalysis( + light_model=LightModel(light_model_list=lens_light_model_list) + ) + + ( + amps, + sigmas, + center_x, + center_y, + ) = self._light_profile_analysis.multi_gaussian_decomposition( + kwargs_lens_light, r_h=r_eff, **kwargs_mge_light + ) + + lens_light_model_list = ["MULTI_GAUSSIAN"] + kwargs_lens_light = [{"amp": amps, "sigma": sigmas}] + + lens_model_list = ["GNFW", "MULTI_GAUSSIAN_KAPPA"] + + super(KinConstraintsComposite, self).__init__( + z_lens, + z_source, + theta_E, + theta_E_error, + gamma, + gamma_error, + r_eff, + r_eff_error, + sigma_v_measured, + kwargs_aperture, + kwargs_seeing, + kwargs_numerics_galkin, + anisotropy_model, + sigma_v_error_independent=sigma_v_error_independent, + sigma_v_error_covariant=sigma_v_error_covariant, + sigma_v_error_cov_matrix=sigma_v_error_cov_matrix, + kwargs_lens_light=kwargs_lens_light, + lens_light_model_list=lens_light_model_list, + lens_model_list=lens_model_list, + MGE_light=False, # set False, as MGE is already done as default + kwargs_mge_light=None, + hernquist_approx=False, + sampling_number=sampling_number, + num_psf_sampling=num_psf_sampling, + num_kin_sampling=num_kin_sampling, + multi_observations=multi_observations, + ) + + if len(rho0_array) != len(r_s_array): + raise ValueError("rho0 and r_s must have the same length!") + + self._rho_s_array = rho0_array + self._r_scale_array = r_s_array + self._kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle( + rho0_array, r_s_array + ) + self.gamma_in_array = gamma_in_array + self.m2l_array = m2l_array + + def get_kappa_s_r_s_angle(self, rho_s, r_scale): + """Computes the surface mass density of the NFW halo at the scale radius. + + :param rho_s: halo mass normalization in M_sun / Mpc^3 + :param r_scale: halo scale radius in arc seconds + :return: surface mass density divided by the critical density + """ + r_s_angle = r_scale / self.lensCosmo.dd / const.arcsec # Rs in arcsec + kappa_s = rho_s * r_scale / self.lensCosmo.sigma_crit + + return kappa_s, r_s_angle + + def draw_lens(self, no_error=False): + """Draws a lens model from the posterior. + + :param no_error: bool, if True, does not render from the uncertainty but uses + the mean values instead + """ + if no_error is True: + return ( + np.mean(self._rho_s_array), + np.mean(self._r_scale_array), + self._r_eff, + 1, + ) + + random_index = np.random.randint(low=0, high=len(self._rho_s_array)) + kappa_s_draw = self._kappa_s_array[random_index] + r_scale_angle_draw = self._r_scale_angle_array[random_index] + + # we make sure no negative r_eff are being sampled + delta_r_eff = np.maximum( + np.random.normal(loc=1, scale=self._r_eff_error / self._r_eff), 0.001 + ) + r_eff_draw = delta_r_eff * self._r_eff + + return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff + + def model_marginalization(self, num_sample_model=20): + """ + + :param num_sample_model: number of samples drawn from the lens and light model + posterior to compute the dimensionless kinematic component J() + :return: J() as array for each measurement prediction, covariance matrix in + sqrt(J) + """ + num_data = len(self._sigma_v_measured) + j_kin_matrix = np.zeros( + (num_sample_model, num_data) + ) # matrix that contains the sampled J() distribution + for i in range(num_sample_model): + j_kin = self.j_kin_draw_composite( + self.kwargs_anisotropy_base, + np.mean(self.gamma_in_array), + np.mean(self.m2l_array), + no_error=False, + ) + j_kin_matrix[i, :] = j_kin + + error_cov_j_sqrt = np.cov(np.sqrt(j_kin_matrix.T)) + j_model_list = np.mean(j_kin_matrix, axis=0) + return j_model_list, error_cov_j_sqrt + + def j_kin_draw_composite(self, kwargs_anisotropy, gamma_in, m2l, no_error=False): + """One simple sampling realization of the dimensionless kinematics of the model. + + :param kwargs_anisotropy: keyword argument of anisotropy setting + :param gamma_in: power-law slope of the mass model + :param m2l: mass-to-light ratio of the stellar component + :param no_error: bool, if True, does not render from the uncertainty but uses + the mean values instead + :return: dimensionless kinematic component J() Birrer et al. 2016, 2019 + """ + kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff = self.draw_lens( + no_error=no_error + ) + + kwargs_lens_stars = copy.deepcopy(self._kwargs_lens_light[0]) + + kwargs_lens_stars["amp"] *= m2l + + kwargs_lens_stars["sigma"] *= delta_r_eff + + kwargs_light = copy.deepcopy(self._kwargs_lens_light) + for kwargs in kwargs_light: + kwargs["sigma"] *= delta_r_eff + + kwargs_lens = [ + { + "Rs": r_scale_angle_draw, + "gamma_in": gamma_in, + "kappa_s": kappa_s_draw, + "center_x": 0, + "center_y": 0, + }, + kwargs_lens_stars, + ] + + j_kin = self.velocity_dispersion_map_dimension_less( + kwargs_lens=kwargs_lens, + kwargs_lens_light=kwargs_light, + kwargs_anisotropy=kwargs_anisotropy, + r_eff=r_eff_draw, + theta_E=self._theta_E, # send this to avoid unnecessary recomputation + gamma=self._gamma, # send this to avoid unnecessary recomputation + ) + return j_kin + + def hierarchy_configuration(self, num_sample_model=20): + """Routine to configure the likelihood to be used in the hierarchical sampling. + In particular, a default configuration is set to compute the Gaussian + approximation of Ds/Dds by sampling the posterior and the estimate of the + variance of the sample. The anisotropy scaling is then performed. Different + anisotropy models are supported. + + :param num_sample_model: number of samples drawn from the lens and light model + posterior to compute the dimensionless kinematic component J() + :return: keyword arguments + """ + + j_model_list, error_cov_j_sqrt = self.model_marginalization(num_sample_model) + ani_scaling_grid_list = self.anisotropy_scaling() + + error_cov_measurement = self.error_cov_measurement + # configuration keyword arguments for the hierarchical sampling + kwargs_likelihood = { + "z_lens": self._z_lens, + "z_source": self._z_source, + "likelihood_type": "IFUKinCov", + "sigma_v_measurement": self._sigma_v_measured, + "anisotropy_model": self._anisotropy_model, + "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, + "m2l_array": self.m2l_array, + "param_scaling_grid_list": ani_scaling_grid_list, + } + return kwargs_likelihood + + def anisotropy_scaling(self): + """ + + :return: anisotropy scaling grid along the axes defined by ani_param_array + """ + j_ani_0 = self.j_kin_draw_composite( + self.kwargs_anisotropy_base, + np.mean(self.gamma_in_array), + np.mean(self.m2l_array), + no_error=True, + ) + return self._anisotropy_scaling_relative(j_ani_0) + + def _anisotropy_scaling_relative(self, j_ani_0): + """Anisotropy scaling relative to a default J prediction. + + :param j_ani_0: default J() prediction for default anisotropy + :return: list of arrays (for the number of measurements) according to anisotropy + scaling + """ + num_data = len(self._sigma_v_measured) + + if self._anisotropy_model == "GOM": + ani_scaling_grid_list = [ + np.zeros( + ( + len(self.ani_param_array[0]), + len(self.ani_param_array[1]), + len(self.gamma_in_array), + len(self.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 k, g_in in enumerate(self.gamma_in_array): + for l, m2l in enumerate(self.m2l_array): + kwargs_anisotropy = self.anisotropy_kwargs( + a_ani=a_ani, beta_inf=beta_inf + ) + j_kin_ani = self.j_kin_draw_composite( + kwargs_anisotropy, g_in, m2l, no_error=True + ) + + for m, j_kin in enumerate(j_kin_ani): + ani_scaling_grid_list[m][i, j, k, l] = ( + j_kin / j_ani_0[m] + ) + # perhaps change the order + elif self._anisotropy_model in ["OM", "const"]: + ani_scaling_grid_list = [ + np.zeros( + ( + len(self.ani_param_array), + len(self.gamma_in_array), + len(self.m2l_array), + ) + ) + for _ in range(num_data) + ] + for i, a_ani in enumerate(self.ani_param_array): + for k, g_in in enumerate(self.gamma_in_array): + for l, m2l in enumerate(self.m2l_array): + kwargs_anisotropy = self.anisotropy_kwargs(a_ani) + j_kin_ani = self.j_kin_draw_composite( + kwargs_anisotropy, g_in, m2l, no_error=True + ) + for m, j_kin in enumerate(j_kin_ani): + ani_scaling_grid_list[m][i, k, l] = j_kin / j_ani_0[m] + else: + raise ValueError("anisotropy model %s not valid." % self._anisotropy_model) + return ani_scaling_grid_list diff --git a/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py b/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py index 854b7f4..c618311 100644 --- a/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py @@ -154,13 +154,13 @@ def num_data(self): return self._lens_type.num_data def log_likelihood( - self, ddt, dd, aniso_scaling=None, sigma_v_sys_error=None, mu_intrinsic=None + self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None, mu_intrinsic=None ): """ :param ddt: time-delay distance [physical Mpc] :param dd: angular diameter distance to the lens [physical Mpc] - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :param sigma_v_sys_error: unaccounted uncertainty in the velocity dispersion measurement @@ -175,12 +175,12 @@ def log_likelihood( ]: return self._lens_type.log_likelihood(ddt, dd) elif self.likelihood_type in ["DdtDdKDE", "DdtDdGaussian", "DsDdsGaussian"]: - return self._lens_type.log_likelihood(ddt, dd, aniso_scaling=aniso_scaling) + return self._lens_type.log_likelihood(ddt, dd, kin_scaling=kin_scaling) elif self.likelihood_type in ["DdtHistKin", "IFUKinCov", "DdtGaussKin"]: return self._lens_type.log_likelihood( ddt, dd, - aniso_scaling=aniso_scaling, + kin_scaling=kin_scaling, sigma_v_sys_error=sigma_v_sys_error, ) elif self.likelihood_type in ["Mag"]: @@ -220,14 +220,14 @@ def sigma_v_measurement(self, sigma_v_sys_error=None): ) return None, None - def sigma_v_prediction(self, ddt, dd, aniso_scaling=None): + def sigma_v_prediction(self, ddt, dd, kin_scaling=None): """ :param ddt: ddt in physical Mpc :param dd: dd in physical Mpc - :param aniso_scaling: anisotropy scaling in J + :param kin_scaling: anisotropy scaling in J :return: model predicted velocity dispersion (vector) and model covariance matrix thereof """ if self.likelihood_type in ["DdtHistKin", "IFUKinCov", "DdtGaussKin"]: - return self._lens_type.sigma_v_prediction(ddt, dd, aniso_scaling) + return self._lens_type.sigma_v_prediction(ddt, dd, kin_scaling) return None, None diff --git a/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py index 8730d07..f825300 100644 --- a/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py @@ -25,18 +25,18 @@ def __init__(self, z_lens, z_source, ddt_mean, ddt_sigma, dd_mean, dd_sigma): ) self.num_data = 2 - def log_likelihood(self, ddt, dd, aniso_scaling=None): + def log_likelihood(self, ddt, dd, kin_scaling=None): """ :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :return: log likelihood given the single lens analysis """ - if aniso_scaling is not None: - dd_ = dd * aniso_scaling[0] + if kin_scaling is not None: + dd_ = dd * kin_scaling[0] else: dd_ = dd lnlikelihood = ( diff --git a/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py index 97ead86..9590cfb 100644 --- a/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py @@ -52,18 +52,18 @@ def __init__( self._interpol = interpol self.num_data = 2 - def log_likelihood(self, ddt, dd, aniso_scaling=None): + def log_likelihood(self, ddt, dd, kin_scaling=None): """ :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :return: log likelihood given the single lens analysis """ - if aniso_scaling is not None: - dd_ = dd * aniso_scaling[0] + if kin_scaling is not None: + dd_ = dd * kin_scaling[0] else: dd_ = dd if self._interpol is True: diff --git a/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py index f9cfbdd..d3a5822 100644 --- a/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py @@ -55,7 +55,7 @@ def log_likelihood( self, ddt, dd, - aniso_scaling=None, + kin_scaling=None, sigma_v_sys_error=None, sigma_v_sys_offset=None, ): @@ -63,7 +63,7 @@ def log_likelihood( :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature @@ -76,19 +76,19 @@ def log_likelihood( ) + self._kinlikelihood.log_likelihood( ddt, dd, - aniso_scaling, + kin_scaling, sigma_v_sys_error=sigma_v_sys_error, sigma_v_sys_offset=sigma_v_sys_offset, ) return lnlikelihood - def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): + def sigma_v_prediction(self, ddt, dd, kin_scaling=1): """Model prediction mean velocity dispersion vector and model prediction covariance matrix. :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the @@ -96,7 +96,7 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): :return: model prediction mean velocity dispersion vector and model prediction covariance matrix """ - return self._kinlikelihood.sigma_v_prediction(ddt, dd, aniso_scaling) + return self._kinlikelihood.sigma_v_prediction(ddt, dd, kin_scaling) def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None): """ diff --git a/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py index ec37141..760dadf 100644 --- a/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py @@ -61,28 +61,28 @@ def __init__( ) self.num_data = self._tdLikelihood.num_data + self._kinlikelihood.num_data - def log_likelihood(self, ddt, dd, aniso_scaling=None, sigma_v_sys_error=None): + def log_likelihood(self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None): """ :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: numpy array of anisotropy scaling on prediction of Ds/Dds + :param kin_scaling: numpy array of anisotropy scaling on prediction of Ds/Dds :return: log likelihood given the single lens analysis """ lnlikelihood = self._tdLikelihood.log_likelihood( ddt ) + self._kinlikelihood.log_likelihood( - ddt, dd, aniso_scaling, sigma_v_sys_error=sigma_v_sys_error + ddt, dd, kin_scaling, sigma_v_sys_error=sigma_v_sys_error ) return lnlikelihood - def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): + def sigma_v_prediction(self, ddt, dd, kin_scaling=1): """Model prediction mean velocity dispersion vector and model prediction covariance matrix. :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the @@ -90,7 +90,7 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): :return: model prediction mean velocity dispersion vector and model prediction covariance matrix """ - return self._kinlikelihood.sigma_v_prediction(ddt, dd, aniso_scaling) + return self._kinlikelihood.sigma_v_prediction(ddt, dd, kin_scaling) def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None): """ diff --git a/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py b/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py index 03ebbaa..870022a 100644 --- a/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py @@ -19,20 +19,20 @@ def __init__(self, z_lens, z_source, ds_dds_mean, ds_dds_sigma): self._ds_dds_sigma2 = ds_dds_sigma**2 self.num_data = 1 - def log_likelihood(self, ddt, dd, aniso_scaling=None): + def log_likelihood(self, ddt, dd, kin_scaling=None): """ Note: kinematics + imaging data can constrain Ds/Dds. The input of Ddt, Dd is transformed here to match Ds/Dds :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :return: log likelihood given the single lens analysis """ ds_dds = ddt / dd / (1 + self._z_lens) - if aniso_scaling is not None: - scaling = aniso_scaling[0] + if kin_scaling is not None: + scaling = kin_scaling[0] else: scaling = 1 ds_dds_ = ds_dds / scaling diff --git a/hierarc/Likelihood/LensLikelihood/kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/kin_likelihood.py index e2e13b3..acee56f 100644 --- a/hierarc/Likelihood/LensLikelihood/kin_likelihood.py +++ b/hierarc/Likelihood/LensLikelihood/kin_likelihood.py @@ -45,7 +45,7 @@ def log_likelihood( self, ddt, dd, - aniso_scaling=None, + kin_scaling=None, sigma_v_sys_error=None, sigma_v_sys_offset=None, ): @@ -54,7 +54,7 @@ def log_likelihood( :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature @@ -63,10 +63,10 @@ def log_likelihood( :return: log likelihood given the single lens analysis """ ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0) - if aniso_scaling is None: + if kin_scaling is None: scaling_ifu = 1 else: - scaling_ifu = aniso_scaling + scaling_ifu = kin_scaling sigma_v_predict = self.sigma_v_model(ds_dds, scaling_ifu) delta = self.sigma_v_measurement_mean(sigma_v_sys_offset) - sigma_v_predict cov_error = self.cov_error_measurement( @@ -98,26 +98,24 @@ def sigma_v_measurement_mean(self, sigma_v_sys_offset=None): else: return self._sigma_v_measured * (1 + sigma_v_sys_offset) - def sigma_v_model(self, ds_dds, aniso_scaling=1): + def sigma_v_model(self, ds_dds, kin_scaling=1): """Model predicted velocity dispersion for the IFU's. :param ds_dds: Ds/Dds - :param aniso_scaling: scaling of the anisotropy affecting sigma_v^2 + :param kin_scaling: scaling of the anisotropy affecting sigma_v^2 :return: array of predicted velocity dispersions """ - sigma_v_predict = ( - np.sqrt(self._j_model * ds_dds * aniso_scaling) * const.c / 1000 - ) + sigma_v_predict = np.sqrt(self._j_model * ds_dds * kin_scaling) * const.c / 1000 return sigma_v_predict - def cov_error_model(self, ds_dds, aniso_scaling=1): + def cov_error_model(self, ds_dds, kin_scaling=1): """ :param ds_dds: Ds/Dds - :param aniso_scaling: scaling of the anisotropy affecting sigma_v^2 + :param kin_scaling: scaling of the anisotropy affecting sigma_v^2 :return: covariance matrix of the error in the predicted model (from mass model uncertainties) """ - scaling_matix = np.outer(np.sqrt(aniso_scaling), np.sqrt(aniso_scaling)) + scaling_matix = np.outer(np.sqrt(kin_scaling), np.sqrt(kin_scaling)) return self._error_cov_j_sqrt * scaling_matix * ds_dds * (const.c / 1000) ** 2 def cov_error_measurement(self, sigma_v_sys_error=None): @@ -134,13 +132,13 @@ def cov_error_measurement(self, sigma_v_sys_error=None): else: return self._error_cov_measurement - def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): + def sigma_v_prediction(self, ddt, dd, kin_scaling=1): """Model prediction mean velocity dispersion vector and model prediction covariance matrix. :param ddt: time-delay distance :param dd: angular diameter distance to the deflector - :param aniso_scaling: array of size of the velocity dispersion measurement or + :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the @@ -149,8 +147,8 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1): covariance matrix """ ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0) - sigma_v_predict = self.sigma_v_model(ds_dds, aniso_scaling) - cov_error_predict = self.cov_error_model(ds_dds, aniso_scaling) + sigma_v_predict = self.sigma_v_model(ds_dds, kin_scaling) + cov_error_predict = self.cov_error_model(ds_dds, kin_scaling) return sigma_v_predict, cov_error_predict def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None): diff --git a/hierarc/Likelihood/anisotropy_scaling.py b/hierarc/Likelihood/anisotropy_scaling.py deleted file mode 100644 index 0afec08..0000000 --- a/hierarc/Likelihood/anisotropy_scaling.py +++ /dev/null @@ -1,156 +0,0 @@ -__author__ = "sibirrer" -from scipy.interpolate import interp1d, interp2d -import numpy as np - - -class AnisotropyScalingSingleAperture(object): - """Class to manage anisotropy scaling for single slit observation.""" - - def __init__(self, ani_param_array, ani_scaling_array): - """ - - :param ani_param_array: array of anisotropy parameter value - :param ani_scaling_array: array with the scaling of J() for single slit - """ - self._evalute_ani = False - if ani_param_array is not None and ani_scaling_array is not None: - if isinstance(ani_param_array, list): - self._dim_scaling = len(ani_param_array) - else: - self._dim_scaling = 1 - if self._dim_scaling == 1: - self._f_ani = interp1d( - ani_param_array, ani_scaling_array, kind="linear" - ) - elif self._dim_scaling == 2: - self._f_ani = interp2d( - ani_param_array[0], ani_param_array[1], ani_scaling_array.T - ) - else: - raise ValueError( - "anisotropy scaling with dimension %s not supported." - % self._dim_scaling - ) - self._evalute_ani = True - - def ani_scaling(self, aniso_param_array): - """ - - :param aniso_param_array: anisotropy parameter array - :return: scaling J(a_ani) for single slit - """ - if self._evalute_ani is not True or aniso_param_array is None: - return 1 - if self._dim_scaling == 1: - return self._f_ani(aniso_param_array[0]) - elif self._dim_scaling == 2: - return self._f_ani(aniso_param_array[0], aniso_param_array[1])[0] - - -class AnisotropyScalingIFU(object): - """Class to manage anisotropy scalings for IFU data.""" - - def __init__( - self, anisotropy_model="NONE", ani_param_array=None, ani_scaling_array_list=None - ): - """ - - :param anisotropy_model: string, either 'NONE', 'OM' or 'GOM' - :param ani_param_array: array of anisotropy parameter value (1d for 'OM' model, 2d for 'GOM' model) - :param ani_scaling_array_list: list of array with the scalings of J() for each IFU - """ - self._anisotropy_model = anisotropy_model - self._evalute_ani = False - if ( - ani_param_array is not None - and ani_scaling_array_list is not None - and self._anisotropy_model != "NONE" - ): - self._evalute_ani = True - self._anisotropy_scaling_list = [] - self._f_ani_list = [] - for ani_scaling_array in ani_scaling_array_list: - self._anisotropy_scaling_list.append( - AnisotropyScalingSingleAperture( - ani_param_array=ani_param_array, - ani_scaling_array=ani_scaling_array, - ) - ) - - if isinstance(ani_param_array, list): - self._dim_scaling = len(ani_param_array) - else: - self._dim_scaling = 1 - if self._dim_scaling == 1 and anisotropy_model in ["OM", "const"]: - self._ani_param_min = np.min(ani_param_array) - self._ani_param_max = np.max(ani_param_array) - elif self._dim_scaling == 2 and anisotropy_model == "GOM": - self._ani_param_min = [min(ani_param_array[0]), min(ani_param_array[1])] - self._ani_param_max = [max(ani_param_array[0]), max(ani_param_array[1])] - else: - raise ValueError( - "anisotropy scaling with dimension %s does not match anisotropy model %s" - % (self._dim_scaling, self._anisotropy_model) - ) - - def ani_scaling(self, aniso_param_array): - """ - - :param aniso_param_array: anisotropy parameter array - :return: scaling J(a_ani) for the IFU's - """ - if self._evalute_ani is not True or aniso_param_array is None: - return [1] - scaling_list = [] - for scaling_class in self._anisotropy_scaling_list: - scaling = scaling_class.ani_scaling(aniso_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 diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py index f1f9826..5d78870 100644 --- a/hierarc/Likelihood/cosmo_likelihood.py +++ b/hierarc/Likelihood/cosmo_likelihood.py @@ -23,10 +23,16 @@ def __init__( lambda_mst_sampling=False, lambda_mst_distribution="delta", anisotropy_sampling=False, + gamma_in_sampling=False, + gamma_in_distribution="NONE", + m2l_sampling=False, + m2l_distribution="NONE", kappa_ext_sampling=False, kappa_ext_distribution="NONE", alpha_lambda_sampling=False, beta_lambda_sampling=False, + alpha_gamma_in_sampling=False, + alpha_m2l_sampling=False, lambda_ifu_sampling=False, lambda_ifu_distribution="NONE", sigma_v_systematics=False, @@ -73,6 +79,10 @@ def __init__( 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 m2l_sampling: bool, if True samples a global mass-to-light ratio parameter + :param m2l_distribution: string, distribution function of the 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 @@ -107,6 +117,12 @@ def __init__( 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, + m2l_sampling=m2l_sampling, + m2l_distribution=m2l_distribution, + alpha_gamma_in_sampling=alpha_gamma_in_sampling, + alpha_m2l_sampling=alpha_m2l_sampling, sne_apparent_m_sampling=sne_apparent_m_sampling, sne_distribution=sne_distribution, z_apparent_m_anchor=z_apparent_m_anchor, diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py index de3a934..533f2c3 100644 --- a/hierarc/Likelihood/hierarchy_likelihood.py +++ b/hierarc/Likelihood/hierarchy_likelihood.py @@ -1,12 +1,12 @@ from hierarc.Likelihood.transformed_cosmography import TransformedCosmography from hierarc.Likelihood.LensLikelihood.base_lens_likelihood import LensLikelihoodBase -from hierarc.Likelihood.anisotropy_scaling import AnisotropyScalingIFU +from hierarc.Likelihood.parameter_scaling import ParameterScalingIFU from hierarc.Util.distribution_util import PDFSampling import numpy as np import copy -class LensLikelihood(TransformedCosmography, LensLikelihoodBase, AnisotropyScalingIFU): +class LensLikelihood(TransformedCosmography, LensLikelihoodBase, ParameterScalingIFU): """Master class containing the likelihood definitions of different analysis for s single lens.""" @@ -18,8 +18,11 @@ def __init__( likelihood_type="TDKin", anisotropy_model="NONE", ani_param_array=None, - ani_scaling_array_list=None, ani_scaling_array=None, + ani_scaling_array_list=None, + gamma_in_array=None, + m2l_array=None, + param_scaling_grid_list=None, num_distribution_draws=50, kappa_ext_bias=False, kappa_pdf=None, @@ -42,7 +45,9 @@ def __init__( 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 ani_scaling_array_list: list of array with the scalings of J() for each IFU + :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 m2l. In that case, gamma_in_array and m2l_array need to be provided. :param num_distribution_draws: int, number of distribution draws from the likelihood that are being averaged over :param kappa_ext_bias: bool, if True incorporates the global external selection function into the likelihood. @@ -65,12 +70,32 @@ def __init__( 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, - ) + + # 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 m2l_array is not None: + if isinstance(ani_param_array, list): + param_arrays = ani_param_array + [gamma_in_array, m2l_array] + else: + param_arrays = [ani_param_array, gamma_in_array, m2l_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, + ) + LensLikelihoodBase.__init__( self, z_lens=z_lens, @@ -93,6 +118,8 @@ def __init__( self._draw_kappa = False self._lambda_scaling_property = lambda_scaling_property self._lambda_scaling_property_beta = lambda_scaling_property_beta + self._gamma_in_array = gamma_in_array + self._m2l_array = m2l_array def lens_log_likelihood( self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None @@ -116,7 +143,7 @@ def lens_log_likelihood( delta_lum_dist = self.luminosity_distance_modulus(cosmo, z_apparent_m_anchor) # here we effectively change the posteriors of the lens, but rather than changing the instance of the KDE we # displace the predicted angular diameter distances in the opposite direction - return self.hyper_param_likelihood( + a = self.hyper_param_likelihood( ddt, dd, delta_lum_dist, @@ -125,6 +152,7 @@ def lens_log_likelihood( kwargs_source=kwargs_source, cosmo=cosmo, ) + return a def hyper_param_likelihood( self, @@ -221,18 +249,88 @@ def log_likelihood_single( kappa_ext=kappa_ext, mag_source=mag_source, ) - aniso_param_array = self.draw_anisotropy(**kwargs_kin) - aniso_scaling = self.ani_scaling(aniso_param_array) + scaling_param_array = self.draw_scaling_params( + kwargs_lens=kwargs_lens, **kwargs_kin + ) + kin_scaling = self.param_scaling(scaling_param_array) lnlikelihood = self.log_likelihood( ddt_, dd_, - aniso_scaling=aniso_scaling, + kin_scaling=kin_scaling, sigma_v_sys_error=sigma_v_sys_error, mu_intrinsic=mag_source_, ) 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._m2l_array is not None: + gamma_in, m2l = self.draw_lens_scaling_params(**kwargs_lens) + return np.concatenate([ani_param, [gamma_in, m2l]]) + 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, + m2l=1, + m2l_sigma=0, + alpha_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 m2l: mass-to-light ratio + :param m2l_sigma: spread in the distribution + :param alpha_m2l: float, linear slope of the 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._m2l_array is not None: + gamma_in_draw, m2l_draw = self.draw_lens_parameters( + gamma_in + alpha_gamma_in * self._lambda_scaling_property, + gamma_in_sigma, + m2l + alpha_m2l * self._lambda_scaling_property, + m2l_sigma, + ) + return gamma_in_draw, m2l_draw + else: + return None, None + def angular_diameter_distances(self, cosmo): """Time-delay distance Ddt, angular diameter distance to the lens (dd) @@ -312,6 +410,12 @@ def draw_lens( lambda_ifu_sigma=0, alpha_lambda=0, beta_lambda=0, + gamma_in=1, + gamma_in_sigma=0, + alpha_gamma_in=0, + m2l=1, + m2l_sigma=0, + alpha_m2l=0, ): """Draws a realization of a specific model from the hyper-parameter distribution. @@ -329,6 +433,14 @@ def draw_lens( 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 m2l: mass-to-light ratio + :param m2l_sigma: spread in the distribution + :param alpha_m2l: float, linear slope of the m2l scaling relation with lens + quantity self._lambda_scaling_property :return: draw from the distributions """ if self._mst_ifu is True: @@ -403,10 +515,12 @@ def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None): ddt_, dd_, _ = self.displace_prediction( ddt, dd, gamma_ppn=gamma_ppn, lambda_mst=lambda_mst, kappa_ext=kappa_ext ) - aniso_param_array = self.draw_anisotropy(**kwargs_kin_copy) - aniso_scaling = self.ani_scaling(aniso_param_array) + scaling_param_array = self.draw_scaling_params( + kwargs_lens=kwargs_lens, **kwargs_kin_copy + ) + kin_scaling = self.param_scaling(scaling_param_array) sigma_v_predict_i, cov_error_predict_i = self.sigma_v_prediction( - ddt_, dd_, aniso_scaling=aniso_scaling + ddt_, dd_, kin_scaling=kin_scaling ) sigma_v_predict_mean += sigma_v_predict_i cov_error_predict += cov_error_predict_i diff --git a/hierarc/Likelihood/parameter_scaling.py b/hierarc/Likelihood/parameter_scaling.py new file mode 100644 index 0000000..3a00f82 --- /dev/null +++ b/hierarc/Likelihood/parameter_scaling.py @@ -0,0 +1,208 @@ +__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 + 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._m2l_min = np.min(param_arrays[2]) + self._m2l_max = np.max(param_arrays[2]) + + 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._m2l_min = np.min(param_arrays[3]) + self._m2l_max = np.max(param_arrays[3]) + 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, m2l=None, 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 m2l: mean of the distribution + :param m2l_sigma: std of the distribution + :return: random draw from the distribution + """ + 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 m2l < self._m2l_min or m2l > self._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) + m2l_draw = np.random.normal(m2l, m2l_sigma) + + if ( + gamma_in_draw < self._gamma_in_min + or gamma_in_draw > self._gamma_in_max + or m2l_draw < self._m2l_min + or m2l_draw > self._m2l_max + ): + return self.draw_lens_parameters(gamma_in, gamma_in_sigma, m2l, m2l_sigma) + + return gamma_in_draw, m2l_draw diff --git a/hierarc/Sampling/ParamManager/lens_param.py b/hierarc/Sampling/ParamManager/lens_param.py index aefb3c3..60b94de 100644 --- a/hierarc/Sampling/ParamManager/lens_param.py +++ b/hierarc/Sampling/ParamManager/lens_param.py @@ -8,12 +8,18 @@ def __init__( self, lambda_mst_sampling=False, lambda_mst_distribution="NONE", - kappa_ext_sampling=False, - kappa_ext_distribution="NONE", lambda_ifu_sampling=False, lambda_ifu_distribution="NONE", + gamma_in_sampling=False, + gamma_in_distribution="NONE", + m2l_sampling=False, + m2l_distribution="NONE", + kappa_ext_sampling=False, + kappa_ext_distribution="NONE", alpha_lambda_sampling=False, beta_lambda_sampling=False, + alpha_gamma_in_sampling=False, + alpha_m2l_sampling=False, kwargs_fixed=None, log_scatter=False, ): @@ -21,15 +27,26 @@ def __init__( :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 kappa_ext_sampling: bool, if True samples a global external convergence parameter - :param kappa_ext_distribution: string, distribution function of the kappa_ext parameter - :param lambda_ifu_sampling: bool, if True samples a separate lambda_mst for a second (e.g. IFU) data set - independently + :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 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 m2l_sampling: bool, if True samples the mass to light ratio of + the stars + :param m2l_distribution: string, distribution function of the mass to + light ratio of the lens + :param kappa_ext_sampling: bool, if True samples a global external convergence + parameter + :param kappa_ext_distribution: string, distribution function of the kappa_ext + parameter :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly - according to a predefined quantity of the lens + 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 + 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_m2l_sampling: bool, if True samples a parameter alpha_m2l, which scales m2l linearly :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space (and thus flat prior in log) :param kwargs_fixed: keyword arguments that are held fixed through the sampling """ @@ -37,10 +54,17 @@ def __init__( self._lambda_mst_distribution = lambda_mst_distribution self._lambda_ifu_sampling = lambda_ifu_sampling self._lambda_ifu_distribution = lambda_ifu_distribution + self._gamma_in_sampling = gamma_in_sampling + self._gamma_in_distribution = gamma_in_distribution + self._m2l_sampling = m2l_sampling + self._m2l_distribution = m2l_distribution self._kappa_ext_sampling = kappa_ext_sampling self._kappa_ext_distribution = kappa_ext_distribution self._alpha_lambda_sampling = alpha_lambda_sampling self._beta_lambda_sampling = beta_lambda_sampling + self._alpha_gamma_in_sampling = alpha_gamma_in_sampling + self._alpha_m2l_sampling = alpha_m2l_sampling + self._log_scatter = log_scatter if kwargs_fixed is None: kwargs_fixed = {} @@ -83,6 +107,36 @@ def param_list(self, latex_style=False): list.append(r"$\sigma(\lambda_{\rm ifu})$") else: list.append("lambda_ifu_sigma") + if self._gamma_in_sampling: + if "gamma_in" not in self._kwargs_fixed: + if latex_style is True: + list.append(r"$\gamma_{\rm in}$") + else: + list.append("gamma_in") + if self._gamma_in_distribution == "GAUSSIAN": + if "gamma_in_sigma" not in self._kwargs_fixed: + if latex_style is True: + if self._log_scatter is True: + list.append(r"$\log_{10}\sigma(\gamma_{\rm in})$") + else: + list.append(r"$\sigma(\gamma_{\rm in})$") + else: + list.append("gamma_in_sigma") + if self._m2l_sampling is True: + if "m2l" not in self._kwargs_fixed: + if latex_style is True: + list.append(r"$\Upsilon_{\rm stars}$") + else: + list.append("m2l") + if self._m2l_distribution == "GAUSSIAN": + if "m2l_sigma" not in self._kwargs_fixed: + if latex_style is True: + if self._log_scatter is True: + list.append(r"$\log_{10}\sigma(\Upsilon_{\rm stars})$") + else: + list.append(r"$\sigma(\Upsilon_{\rm stars})$") + else: + list.append("m2l_sigma") if self._kappa_ext_sampling is True: if "kappa_ext" not in self._kwargs_fixed: if latex_style is True: @@ -107,6 +161,18 @@ def param_list(self, latex_style=False): list.append(r"$\beta_{\lambda}$") else: list.append("beta_lambda") + if self._alpha_gamma_in_sampling is True: + if "alpha_gamma_in" not in self._kwargs_fixed: + if latex_style is True: + list.append(r"$\alpha_{\gamma_{\rm in}}$") + else: + list.append("alpha_gamma_in") + if self._alpha_m2l_sampling is True: + if "alpha_m2l" not in self._kwargs_fixed: + if latex_style is True: + list.append(r"$\alpha_{\Upsilon_{\rm stars}}$") + else: + list.append("alpha_m2l") return list def args2kwargs(self, args, i=0): @@ -146,6 +212,36 @@ def args2kwargs(self, args, i=0): else: kwargs["lambda_ifu_sigma"] = args[i] i += 1 + if self._gamma_in_sampling is True: + if "gamma_in" in self._kwargs_fixed: + kwargs["gamma_in"] = self._kwargs_fixed["gamma_in"] + else: + kwargs["gamma_in"] = args[i] + i += 1 + if self._gamma_in_distribution == "GAUSSIAN": + if "gamma_in_sigma" in self._kwargs_fixed: + kwargs["gamma_in_sigma"] = self._kwargs_fixed["gamma_in_sigma"] + else: + if self._log_scatter is True: + kwargs["gamma_in_sigma"] = 10 ** (args[i]) + else: + kwargs["gamma_in_sigma"] = args[i] + i += 1 + if self._m2l_sampling is True: + if "m2l" in self._kwargs_fixed: + kwargs["m2l"] = self._kwargs_fixed["m2l"] + else: + kwargs["m2l"] = args[i] + i += 1 + if self._m2l_distribution == "GAUSSIAN": + if "m2l_sigma" in self._kwargs_fixed: + kwargs["m2l_sigma"] = self._kwargs_fixed["m2l_sigma"] + else: + if self._log_scatter is True: + kwargs["m2l_sigma"] = 10 ** (args[i]) + else: + kwargs["m2l_sigma"] = args[i] + i += 1 if self._kappa_ext_sampling is True: if "kappa_ext" in self._kwargs_fixed: kwargs["kappa_ext"] = self._kwargs_fixed["kappa_ext"] @@ -170,6 +266,19 @@ def args2kwargs(self, args, i=0): else: kwargs["beta_lambda"] = args[i] i += 1 + if self._alpha_gamma_in_sampling is True: + if "alpha_gamma_in" in self._kwargs_fixed: + kwargs["alpha_gamma_in"] = self._kwargs_fixed["alpha_gamma_in"] + else: + kwargs["alpha_gamma_in"] = args[i] + i += 1 + if self._alpha_m2l_sampling is True: + if "alpha_m2l" in self._kwargs_fixed: + kwargs["alpha_m2l"] = self._kwargs_fixed["alpha_m2l"] + else: + kwargs["alpha_m2l"] = args[i] + i += 1 + return kwargs, i def kwargs2args(self, kwargs): @@ -197,6 +306,24 @@ def kwargs2args(self, kwargs): args.append(np.log10(kwargs["lambda_ifu_sigma"])) else: args.append(kwargs["lambda_ifu_sigma"]) + if self._gamma_in_sampling is True: + if "gamma_in" not in self._kwargs_fixed: + args.append(kwargs["gamma_in"]) + if self._gamma_in_distribution == "GAUSSIAN": + if "gamma_in_sigma" not in self._kwargs_fixed: + if self._log_scatter is True: + args.append(np.log10(kwargs["gamma_in_sigma"])) + else: + args.append(kwargs["gamma_in_sigma"]) + if self._m2l_sampling is True: + if "m2l" not in self._kwargs_fixed: + args.append(kwargs["m2l"]) + if self._m2l_distribution == "GAUSSIAN": + if "m2l_sigma" not in self._kwargs_fixed: + if self._log_scatter is True: + args.append(np.log10(kwargs["m2l_sigma"])) + else: + args.append(kwargs["m2l_sigma"]) if self._kappa_ext_sampling is True: if "kappa_ext" not in self._kwargs_fixed: args.append(kwargs["kappa_ext"]) @@ -209,4 +336,10 @@ def kwargs2args(self, kwargs): if self._beta_lambda_sampling is True: if "beta_lambda" not in self._kwargs_fixed: args.append(kwargs["beta_lambda"]) + if self._alpha_gamma_in_sampling is True: + if "alpha_gamma_in" not in self._kwargs_fixed: + args.append(kwargs["alpha_gamma_in"]) + if self._alpha_m2l_sampling is True: + if "alpha_m2l" not in self._kwargs_fixed: + args.append(kwargs["alpha_m2l"]) return args diff --git a/hierarc/Sampling/ParamManager/param_manager.py b/hierarc/Sampling/ParamManager/param_manager.py index 093cb8a..9a4f3bd 100644 --- a/hierarc/Sampling/ParamManager/param_manager.py +++ b/hierarc/Sampling/ParamManager/param_manager.py @@ -16,12 +16,18 @@ def __init__( anisotropy_sampling=False, anisotropy_model="OM", anisotropy_distribution="NONE", + gamma_in_sampling=False, + gamma_in_distribution="NONE", + m2l_sampling=False, + m2l_distribution="NONE", kappa_ext_sampling=False, kappa_ext_distribution="NONE", lambda_ifu_sampling=False, lambda_ifu_distribution="NONE", alpha_lambda_sampling=False, beta_lambda_sampling=False, + alpha_gamma_in_sampling=False, + alpha_m2l_sampling=False, sigma_v_systematics=False, sne_apparent_m_sampling=False, sne_distribution="GAUSSIAN", @@ -56,6 +62,14 @@ def __init__( :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens 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 + :param m2l_sampling: bool, if True samples the mass-to-light ratio of the stars + :param m2l_distribution: string, distribution function of the m2l parameter + :param alpha_gamma_in_sampling: bool, if True samples a parameter alpha_gamma_in, which scales gamma_in linearly + according to a predefined quantity of the lens + :param alpha_m2l_sampling: bool, if True samples a parameter alpha_m2l, which scales m2l linearly + according to a predefined quantity of the lens :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). @@ -83,10 +97,16 @@ def __init__( lambda_mst_distribution=lambda_mst_distribution, lambda_ifu_sampling=lambda_ifu_sampling, lambda_ifu_distribution=lambda_ifu_distribution, + gamma_in_sampling=gamma_in_sampling, + gamma_in_distribution=gamma_in_distribution, + m2l_sampling=m2l_sampling, + m2l_distribution=m2l_distribution, kappa_ext_sampling=kappa_ext_sampling, kappa_ext_distribution=kappa_ext_distribution, alpha_lambda_sampling=alpha_lambda_sampling, beta_lambda_sampling=beta_lambda_sampling, + alpha_gamma_in_sampling=alpha_gamma_in_sampling, + alpha_m2l_sampling=alpha_m2l_sampling, log_scatter=log_scatter, kwargs_fixed=kwargs_fixed_lens, ) diff --git a/requirements.txt b/requirements.txt index 00e04f0..ab46f10 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,7 +7,6 @@ mpmath h5py matplotlib pandas -lenstronomy>=1.11.3 -# -e git://github.com/sibirrer/lenstronomy.git@main#egg=lenstronomy +lenstronomy @ git+https://github.com/lenstronomy/lenstronomy@main scikit-learn>=0.22 pandas diff --git a/test/test_LensPosterior/test_kin_constraints.py b/test/test_LensPosterior/test_kin_constraints.py index 0e86bc1..c01a900 100644 --- a/test/test_LensPosterior/test_kin_constraints.py +++ b/test/test_LensPosterior/test_kin_constraints.py @@ -6,7 +6,7 @@ import unittest -class TestIFUKinPosterior(object): +class TestKinConstraints(object): def setup(self): pass diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py new file mode 100644 index 0000000..f855ee2 --- /dev/null +++ b/test/test_LensPosterior/test_kin_constraints_composite.py @@ -0,0 +1,361 @@ +from hierarc.LensPosterior.kin_constraints_composite import KinConstraintsComposite +from lenstronomy.Analysis.kinematics_api import KinematicsAPI +from hierarc.Likelihood.hierarchy_likelihood import LensLikelihood +import numpy.testing as npt +import numpy as np +import pytest +import unittest +from astropy.cosmology import FlatLambdaCDM + + +class TestKinConstraintsComposite(object): + def setup(self): + pass + + def test_likelihoodconfiguration_om(self): + anisotropy_model = "OM" + kwargs_aperture = { + "aperture_type": "shell", + "r_in": 0, + "r_out": 3 / 2.0, + "center_ra": 0.0, + "center_dec": 0, + } + kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4} + + # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed) + kwargs_numerics_galkin = { + "interpol_grid_num": 50, # numerical interpolation, should converge -> infinity + "log_integration": True, + # log or linear interpolation of surface brightness and mass models + "max_integrate": 100, + "min_integrate": 0.001, + } # lower/upper bound of numerical integrals + + # redshift + z_lens = 0.5 + z_source = 1.5 + + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + theta_E = 1.0 + r_eff = 1 + gamma = 2.1 + + kwargs_mge_light = { + "grid_spacing": 0.1, + "grid_num": 10, + "n_comp": 20, + "center_x": 0, + "center_y": 0, + } + + kwargs_kin_api_settings = { + "multi_observations": False, + "kwargs_numerics_galkin": kwargs_numerics_galkin, + "kwargs_mge_light": kwargs_mge_light, + "sampling_number": 50, + "num_kin_sampling": 50, + "num_psf_sampling": 50, + } + + kwargs_lens_light = [ + { + "R_sersic": 2, + "amp": 1, + "n_sersic": 2, + "center_x": 0, + "center_y": 0, + } + ] + lens_light_model_list = ["SERSIC"] + + gamma_in_array = np.linspace(0.1, 2.9, 5) + m2l_array = np.linspace(0.1, 10, 5) + + rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6 + r_s_array = np.random.normal(0.1, 0, 100) + + # compute likelihood + kin_constraints = KinConstraintsComposite( + z_lens=z_lens, + z_source=z_source, + gamma_in_array=gamma_in_array, + m2l_array=m2l_array, + rho0_array=rho0_array, + r_s_array=r_s_array, + theta_E=theta_E, + theta_E_error=0.01, + gamma=gamma, + gamma_error=0.02, + r_eff=r_eff, + r_eff_error=0.05, + sigma_v_measured=[200], + sigma_v_error_independent=[10], + sigma_v_error_covariant=0, + kwargs_aperture=kwargs_aperture, + kwargs_seeing=kwargs_seeing, + anisotropy_model=anisotropy_model, + kwargs_lens_light=kwargs_lens_light, + lens_light_model_list=lens_light_model_list, + **kwargs_kin_api_settings + ) + + kin_constraints.draw_lens(no_error=True) + + kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) + kwargs_likelihood["normalized"] = False + ln_class = LensLikelihood(**kwargs_likelihood) + kwargs_kin = {"a_ani": 1} + ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + + def test_likelihoodconfiguration_gom(self): + anisotropy_model = "GOM" + kwargs_aperture = { + "aperture_type": "shell", + "r_in": 0, + "r_out": 3 / 2.0, + "center_ra": 0.0, + "center_dec": 0, + } + kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4} + + # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed) + kwargs_numerics_galkin = { + "interpol_grid_num": 50, # numerical interpolation, should converge -> infinity + "log_integration": True, + # log or linear interpolation of surface brightness and mass models + "max_integrate": 100, + "min_integrate": 0.001, + } # lower/upper bound of numerical integrals + + # redshift + z_lens = 0.5 + z_source = 1.5 + + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + theta_E = 1.0 + r_eff = 1 + gamma = 2.1 + + kwargs_mge_light = { + "grid_spacing": 0.1, + "grid_num": 10, + "n_comp": 20, + "center_x": 0, + "center_y": 0, + } + + kwargs_kin_api_settings = { + "multi_observations": False, + "kwargs_numerics_galkin": kwargs_numerics_galkin, + "kwargs_mge_light": kwargs_mge_light, + "sampling_number": 10, + "num_kin_sampling": 10, + "num_psf_sampling": 10, + } + + kwargs_lens_light = [ + { + "R_sersic": 2, + "amp": 1, + "n_sersic": 2, + "center_x": 0, + "center_y": 0, + } + ] + lens_light_model_list = ["SERSIC"] + + gamma_in_array = np.linspace(0.1, 2.9, 5) + m2l_array = np.linspace(0.1, 10, 5) + + rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6 + r_s_array = np.random.normal(0.1, 0, 100) + + # compute likelihood + kin_constraints = KinConstraintsComposite( + z_lens=z_lens, + z_source=z_source, + gamma_in_array=gamma_in_array, + m2l_array=m2l_array, + rho0_array=rho0_array, + r_s_array=r_s_array, + theta_E=theta_E, + theta_E_error=0.01, + gamma=gamma, + gamma_error=0.02, + r_eff=r_eff, + r_eff_error=0.05, + sigma_v_measured=[200], + sigma_v_error_independent=[10], + sigma_v_error_covariant=0, + kwargs_aperture=kwargs_aperture, + kwargs_seeing=kwargs_seeing, + anisotropy_model=anisotropy_model, + kwargs_lens_light=kwargs_lens_light, + lens_light_model_list=lens_light_model_list, + **kwargs_kin_api_settings + ) + + kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5) + kwargs_likelihood["normalized"] = False + ln_class = LensLikelihood(**kwargs_likelihood) + kwargs_kin = {"a_ani": 1, "beta_inf": 0.5} + ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin) + + +class TestRaise(unittest.TestCase): + def test_raise(self): + with self.assertRaises(ValueError): + anisotropy_model = "OM" + kwargs_aperture = { + "aperture_type": "shell", + "r_in": 0, + "r_out": 3 / 2.0, + "center_ra": 0.0, + "center_dec": 0, + } + kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4} + + # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed) + kwargs_numerics_galkin = { + "interpol_grid_num": 50, # numerical interpolation, should converge -> infinity + "log_integration": True, + # log or linear interpolation of surface brightness and mass models + "max_integrate": 100, + "min_integrate": 0.001, + } # lower/upper bound of numerical integrals + + # redshift + z_lens = 0.5 + z_source = 1.5 + theta_E = 1.0 + r_eff = 1 + gamma = 2.1 + + kwargs_mge_light = { + "grid_spacing": 0.1, + "grid_num": 10, + "n_comp": 20, + "center_x": 0, + "center_y": 0, + } + + # settings for kinematics calculation with KinematicsAPI of lenstronomy + kwargs_kin_api_settings = { + "multi_observations": False, + "kwargs_numerics_galkin": kwargs_numerics_galkin, + "kwargs_mge_light": kwargs_mge_light, + "sampling_number": 10, + "num_kin_sampling": 10, + "num_psf_sampling": 10, + } + kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}] + + gamma_in_array = np.linspace(0.1, 2.9, 5) + m2l_array = np.linspace(0.1, 10, 5) + rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6 + r_s_array = np.random.normal(0.1, 0.01, 100) + + kin_constraints = KinConstraintsComposite( + z_lens=z_lens, + z_source=z_source, + gamma_in_array=gamma_in_array, + m2l_array=m2l_array, + rho0_array=rho0_array, + r_s_array=r_s_array, + theta_E=theta_E, + theta_E_error=0.01, + gamma=gamma, + gamma_error=0.02, + r_eff=r_eff, + r_eff_error=0.05, + sigma_v_measured=[200], + sigma_v_error_independent=[10], + sigma_v_error_covariant=0, + kwargs_aperture=kwargs_aperture, + kwargs_seeing=kwargs_seeing, + anisotropy_model=anisotropy_model, + kwargs_lens_light=kwargs_lens_light, + **kwargs_kin_api_settings + ) + kin_constraints._anisotropy_model = "BAD" + kin_constraints._anisotropy_scaling_relative(j_ani_0=1) + + with self.assertRaises(ValueError): + anisotropy_model = "OM" + kwargs_aperture = { + "aperture_type": "shell", + "r_in": 0, + "r_out": 3 / 2.0, + "center_ra": 0.0, + "center_dec": 0, + } + kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4} + + # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed) + kwargs_numerics_galkin = { + "interpol_grid_num": 100, # numerical interpolation, should converge -> infinity + "log_integration": True, + # log or linear interpolation of surface brightness and mass models + "max_integrate": 100, + "min_integrate": 0.001, + } # lower/upper bound of numerical integrals + + # redshift + z_lens = 0.5 + z_source = 1.5 + theta_E = 1.0 + r_eff = 1 + gamma = 2.1 + + kwargs_mge_light = { + "grid_spacing": 0.1, + "grid_num": 10, + "n_comp": 20, + "center_x": 0, + "center_y": 0, + } + + # settings for kinematics calculation with KinematicsAPI of lenstronomy + kwargs_kin_api_settings = { + "multi_observations": False, + "kwargs_numerics_galkin": kwargs_numerics_galkin, + "kwargs_mge_light": kwargs_mge_light, + "sampling_number": 50, + "num_kin_sampling": 50, + "num_psf_sampling": 50, + } + + kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}] + gamma_in_array = np.linspace(0.1, 2.9, 5) + m2l_array = np.linspace(0.1, 10, 5) + + rho0_array = 10 ** np.random.normal(8, 0.2, 5) / 1e6 + r_s_array = np.random.normal(0.1, 0.01, 6) + + KinConstraintsComposite( + z_lens=z_lens, + z_source=z_source, + gamma_in_array=gamma_in_array, + m2l_array=m2l_array, + rho0_array=rho0_array, + r_s_array=r_s_array, + theta_E=theta_E, + theta_E_error=0.01, + gamma=gamma, + gamma_error=0.02, + r_eff=r_eff, + r_eff_error=0.05, + sigma_v_measured=[200], + sigma_v_error_independent=[10], + sigma_v_error_covariant=0, + kwargs_aperture=kwargs_aperture, + kwargs_seeing=kwargs_seeing, + anisotropy_model=anisotropy_model, + kwargs_lens_light=kwargs_lens_light, + **kwargs_kin_api_settings + ) + + +if __name__ == "__main__": + pytest.main() diff --git a/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py index 41d057a..436a685 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py @@ -101,7 +101,7 @@ def test_log_likelihood(self): ) print(likelihood_type) logl = likelihood.log_likelihood( - ddt=1, dd=1, aniso_scaling=None, sigma_v_sys_error=1, mu_intrinsic=1 + ddt=1, dd=1, kin_scaling=None, sigma_v_sys_error=1, mu_intrinsic=1 ) print(logl) assert logl > -np.inf @@ -116,7 +116,7 @@ def test_predictions_measurements(self): ) ddt_measurement = likelihood.ddt_measurement() likelihood.sigma_v_measurement(sigma_v_sys_error=0) - likelihood.sigma_v_prediction(ddt=1, dd=1, aniso_scaling=1) + likelihood.sigma_v_prediction(ddt=1, dd=1, kin_scaling=1) assert len(ddt_measurement) == 2 diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py index 5a4f634..2387d28 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py @@ -35,9 +35,9 @@ def test_log_likelihood(self): ) npt.assert_almost_equal(lnlog, -0.5, decimal=5) - aniso_scaling = [1 + self.dd_sigma / self.dd_mean] + kin_scaling = [1 + self.dd_sigma / self.dd_mean] lnlog = likelihood.log_likelihood( - ddt=self.ddt_mean, dd=self.dd_mean, aniso_scaling=aniso_scaling + ddt=self.ddt_mean, dd=self.dd_mean, kin_scaling=kin_scaling ) npt.assert_almost_equal(lnlog, -0.5, decimal=5) diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py index 19defa9..088639b 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py @@ -63,11 +63,11 @@ def setup(self): def test_log_likelihood(self): ddt, dd = 9, 0.9 lnlog_tot = self.ddt_gauss_kin_likelihood.log_likelihood( - ddt, dd, aniso_scaling=None, sigma_v_sys_error=None + ddt, dd, kin_scaling=None, sigma_v_sys_error=None ) lnlog_ddt = self.ddt_gauss_likelihood.log_likelihood(ddt, dd) lnlog_kin = self.kin_likelihood.log_likelihood( - ddt, dd, aniso_scaling=None, sigma_v_sys_error=None + ddt, dd, kin_scaling=None, sigma_v_sys_error=None ) npt.assert_almost_equal(lnlog_tot, lnlog_ddt + lnlog_kin, decimal=5) @@ -88,7 +88,7 @@ def test_sigma_v_measurement(self): sigma_v_predict, error_cov_predict, ) = self.ddt_gauss_kin_likelihood.sigma_v_prediction( - self.ddt_mean, self.dd_mean, aniso_scaling=1 + self.ddt_mean, self.dd_mean, kin_scaling=1 ) assert sigma_v_predict[0] == self.sigma_v_measurement[0] assert error_cov_predict[0, 0] == 0 diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py index 3288f9b..01f9e65 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py @@ -64,7 +64,7 @@ def setup(self): def test_log_likelihood(self): logl_max = self._ddt_kin_likelihood.log_likelihood(ddt=self._ddt, dd=self._dd) logl = self._ddt_kin_likelihood.log_likelihood( - self._ddt, self._dd / (1 + self._sigma) ** 2, aniso_scaling=None + self._ddt, self._dd / (1 + self._sigma) ** 2, kin_scaling=None ) npt.assert_almost_equal(logl - logl_max, -self._num_ifu / 2, decimal=5) @@ -85,7 +85,7 @@ def test_sigma_v_measurement(self): sigma_v_predict, error_cov_predict, ) = self._ddt_kin_likelihood.sigma_v_prediction( - self._ddt, self._dd, aniso_scaling=1 + self._ddt, self._dd, kin_scaling=1 ) assert sigma_v_predict[0] == self.sigma_v_measurement[0] assert error_cov_predict[0, 0] == 0 diff --git a/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py index 44a51dd..665c9cf 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py @@ -35,7 +35,7 @@ def test_log_likelihood(self): npt.assert_almost_equal(lnlog, -0.5, decimal=5) aniso_scaling = [1.0 / (1 + self.ds_dds_sigma / self.ds_dds_mean)] - lnlog = likelihood.log_likelihood(ddt=ddt, dd=dd, aniso_scaling=aniso_scaling) + lnlog = likelihood.log_likelihood(ddt=ddt, dd=dd, kin_scaling=aniso_scaling) npt.assert_almost_equal(lnlog, -0.5, decimal=5) diff --git a/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py index a52bffe..8fc5108 100644 --- a/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py +++ b/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py @@ -33,12 +33,12 @@ def test_log_likelihood(self): error_cov_j_sqrt, normalized=False, ) - logl = kin_likelihood.log_likelihood(ddt, dd, aniso_scaling=None) + logl = kin_likelihood.log_likelihood(ddt, dd, kin_scaling=None) npt.assert_almost_equal(logl, 0, decimal=5) - logl = kin_likelihood.log_likelihood(ddt * (0.9**2), dd, aniso_scaling=None) + logl = kin_likelihood.log_likelihood(ddt * (0.9**2), dd, kin_scaling=None) npt.assert_almost_equal(logl, -num_ifu / 2, decimal=5) logl = kin_likelihood.log_likelihood( - ddt * (1 - np.sqrt(0.1**2 + 0.1**2)) ** 2, dd, aniso_scaling=None + ddt * (1 - np.sqrt(0.1**2 + 0.1**2)) ** 2, dd, kin_scaling=None ) npt.assert_almost_equal(logl, -num_ifu, decimal=5) @@ -50,7 +50,7 @@ def test_log_likelihood(self): assert error_cov_measurement_[0, 0] == error_cov_measurement[0, 0] sigma_v_predict, error_cov_predict = kin_likelihood.sigma_v_prediction( - ddt, dd, aniso_scaling=1 + ddt, dd, kin_scaling=1 ) assert sigma_v_predict[0] == sigma_v_measurement[0] assert error_cov_predict[0, 0] == 0 @@ -82,7 +82,7 @@ def test_error_systematic(self): logl = ifu_likelihood.log_likelihood( ddt * (1 - np.sqrt(0.1**2)) ** 2, dd, - aniso_scaling=1, + kin_scaling=1, sigma_v_sys_error=0.1, ) npt.assert_almost_equal(logl, -1 / 2.0, decimal=5) @@ -112,10 +112,10 @@ def test_log_likelihood_marg(self): normalized=True, sigma_sys_error_include=True, ) - logl = ifu_likelihood.log_likelihood(ddt, dd, aniso_scaling=1) + logl = ifu_likelihood.log_likelihood(ddt, dd, kin_scaling=1) aniso_scaling = 0.8 logl_ani = ifu_likelihood.log_likelihood( - ddt, dd * aniso_scaling, aniso_scaling=aniso_scaling + ddt, dd * aniso_scaling, kin_scaling=aniso_scaling ) npt.assert_almost_equal(logl_ani - logl, 0, decimal=1) @@ -123,14 +123,14 @@ def test_log_likelihood_marg(self): # analytical Gaussian calculation in the covariance matrix normalization sigma_v_sys_error = 0.1 logl = ifu_likelihood.log_likelihood( - ddt, dd, aniso_scaling=1, sigma_v_sys_error=sigma_v_sys_error + ddt, dd, kin_scaling=1, sigma_v_sys_error=sigma_v_sys_error ) l_sum = 0 num_sample = 10000 for i in range(num_sample): sigma_v_pert = np.random.normal(loc=0, scale=sigma_v_sys_error) logl_i = ifu_likelihood.log_likelihood( - ddt, dd, aniso_scaling=1, sigma_v_sys_offset=sigma_v_pert + ddt, dd, kin_scaling=1, sigma_v_sys_offset=sigma_v_pert ) l_sum += np.exp(logl_i) logl_average = np.log(l_sum / num_sample) @@ -149,7 +149,7 @@ def test_log_likelihood_invert_error(self): normalized=True, sigma_sys_error_include=True, ) - logl = ifu_likelihood.log_likelihood(ddt=1, dd=1, aniso_scaling=1) + logl = ifu_likelihood.log_likelihood(ddt=1, dd=1, kin_scaling=1) assert logl == -np.inf diff --git a/test/test_Likelihood/test_anisotropy_scaling.py b/test/test_Likelihood/test_anisotropy_scaling.py deleted file mode 100644 index df1e7bd..0000000 --- a/test/test_Likelihood/test_anisotropy_scaling.py +++ /dev/null @@ -1,147 +0,0 @@ -import pytest -import numpy as np -import unittest -from hierarc.Likelihood.anisotropy_scaling import ( - AnisotropyScalingSingleAperture, - AnisotropyScalingIFU, -) - - -class TestAnisotropyScalingSingleAperture(object): - def setup(self): - ani_param_array = np.linspace(start=0, stop=1, num=10) - ani_scaling_array = ani_param_array * 2 - self.scaling = AnisotropyScalingSingleAperture( - ani_param_array, ani_scaling_array - ) - - 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]) - self.scaling_2d = AnisotropyScalingSingleAperture( - ani_param_array, ani_scaling_array - ) - - def test_ani_scaling(self): - scaling = self.scaling.ani_scaling(aniso_param_array=[1]) - assert scaling == 2 - - scaling = self.scaling.ani_scaling(aniso_param_array=None) - assert scaling == 1 - - scaling = self.scaling_2d.ani_scaling(aniso_param_array=[1, 2]) - assert scaling == 2 - - -class TestAnisotropyScalingIFU(object): - def setup(self): - ani_param_array = np.linspace(start=0, stop=1, num=10) - ani_scaling_array = ani_param_array * 2 - self.scaling = AnisotropyScalingIFU( - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=[ani_scaling_array], - ) - - 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]) - self.scaling_2d = AnisotropyScalingIFU( - anisotropy_model="GOM", - ani_param_array=ani_param_array, - ani_scaling_array_list=[ani_scaling_array], - ) - - def test_ani_scaling(self): - scaling = self.scaling.ani_scaling(aniso_param_array=[1]) - assert scaling[0] == 2 - - scaling = self.scaling.ani_scaling(aniso_param_array=None) - assert scaling[0] == 1 - - scaling = self.scaling_2d.ani_scaling(aniso_param_array=[1, 2]) - assert scaling[0] == 2 - - 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 - ) - - 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 - ) - - scaling = AnisotropyScalingIFU(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 - - -class TestRaise(unittest.TestCase): - def test_raise(self): - with self.assertRaises(ValueError): - ani_param_array = [ - np.linspace(start=0, stop=1, num=10), - np.linspace(start=1, stop=2, num=5), - 1, - ] - ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1]) - self.scaling_2d = AnisotropyScalingSingleAperture( - ani_param_array, ani_scaling_array - ) - - with self.assertRaises(ValueError): - AnisotropyScalingIFU( - anisotropy_model="blabla", - ani_param_array=np.array([0, 1]), - ani_scaling_array_list=[np.array([0, 1])], - ) - - with self.assertRaises(ValueError): - ani_param_array = np.linspace(start=0, stop=1, num=10) - ani_scaling_array = ani_param_array * 2 - scaling = AnisotropyScalingIFU( - anisotropy_model="OM", - ani_param_array=ani_param_array, - ani_scaling_array_list=[ani_scaling_array], - ) - scaling.draw_anisotropy( - a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0 - ) - - with self.assertRaises(ValueError): - 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 = AnisotropyScalingIFU( - anisotropy_model="GOM", - ani_param_array=ani_param_array, - ani_scaling_array_list=[ani_scaling_array], - ) - scaling.draw_anisotropy( - a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0 - ) - - -if __name__ == "__main__": - pytest.main() diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py index 319c3a5..a0ce60e 100644 --- a/test/test_Likelihood/test_hierarchy_likelihood.py +++ b/test/test_Likelihood/test_hierarchy_likelihood.py @@ -92,6 +92,50 @@ def setup(self): **kwargs_likelihood ) + gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5) + m2l_array = np.linspace(start=1, stop=10, num=10) + param_scaling_array = np.multiply.outer( + np.ones_like(ani_param_array), + np.outer(np.ones_like(gamma_in_array), np.ones_like(m2l_array)), + ) + self.likelihood_gamma_in_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, + m2l_array=m2l_array, + num_distribution_draws=200, + kappa_ext_bias=False, + kappa_pdf=None, + kappa_bin_edges=None, + mst_ifu=False, + **kwargs_likelihood + ) + + self.likelihood_gamma_in_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, + m2l_array=m2l_array, + num_distribution_draws=200, + kappa_ext_bias=False, + kappa_pdf=None, + kappa_bin_edges=None, + mst_ifu=False, + **kwargs_likelihood + ) + def test_lens_log_likelihood(self): np.random.seed(42) kwargs_lens = { @@ -132,9 +176,36 @@ def test_lens_log_likelihood(self): ) assert ln_inf < -10000000 + ln_inf = self.likelihood_single.lens_log_likelihood( + self.cosmo, kwargs_lens=None, kwargs_kin=kwargs_kin + ) + npt.assert_almost_equal(ln_inf, 0.0, decimal=1) + + ln_inf = self.likelihood_single.sigma_v_measured_vs_predict( + self.cosmo, kwargs_lens=None, kwargs_kin=None + ) + assert np.all([x is None for x in ln_inf]) + kwargs_test = self.likelihood._kwargs_init(kwargs=None) assert type(kwargs_test) is dict + gamma_in_draw, m2l_draw = self.likelihood.draw_lens_scaling_params() + assert gamma_in_draw is None + assert m2l_draw is None + + kwargs_lens = { + "gamma_in": 1, + "gamma_in_sigma": 0, + "alpha_gamma_in": 0, + "m2l": 1, + "m2l_sigma": 0, + "alpha_m2l": 0, + } + ln_likelihood = self.likelihood_gamma_in_m2l.lens_log_likelihood( + self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin + ) + npt.assert_almost_equal(ln_likelihood, -0.0, decimal=1) + if __name__ == "__main__": pytest.main() diff --git a/test/test_Likelihood/test_parameter_scaling.py b/test/test_Likelihood/test_parameter_scaling.py new file mode 100644 index 0000000..dad4f0c --- /dev/null +++ b/test/test_Likelihood/test_parameter_scaling.py @@ -0,0 +1,310 @@ +import pytest +import numpy as np +import unittest +from hierarc.Likelihood.parameter_scaling import ( + ParameterScalingSingleAperture, + ParameterScalingIFU, +) + + +class TestParameterScalingSingleAperture(object): + def setup(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) + m2l_array = np.linspace(start=1, stop=10, num=10) + + param_arrays = [ani_param_array, gamma_in_array, m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.outer(gamma_in_array, 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, + m2l_array, + ] + param_scaling_array = np.multiply.outer( + gom_param_array[0], + np.multiply.outer( + gom_param_array[1], + np.outer(gamma_in_array, 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, 10]) + assert scaling == 1 * 2.9 * 10 + + scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 10]) + assert scaling == 1 * 2 * 2.9 * 10 + + +class TestParameterScalingIFU(object): + def setup(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) + m2l_array = np.linspace(start=1, stop=10, num=10) + + param_arrays = [ani_param_array, gamma_in_array, m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.outer(gamma_in_array, 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, + m2l_array, + ] + param_scaling_array = np.multiply.outer( + gom_param_array[0], + np.multiply.outer( + gom_param_array[1], + np.outer(gamma_in_array, m2l_array), + ), + ) + self.scaling_nfw_2d = 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, 10]) + assert scaling[0] == 1 * 2.9 * 10 + + scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 10]) + assert scaling[0] == 1 * 2 * 2.9 * 10 + + 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 + ) + + 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, m2l=5, m2l_sigma=0 + ) + assert param_draw[0] == 1 + assert param_draw[1] == 5 + for i in range(100): + param_draw = self.scaling_nfw.draw_lens_parameters( + gamma_in=1, gamma_in_sigma=1, m2l=5, m2l_sigma=3 + ) + + param_draw = self.scaling_nfw_2d.draw_lens_parameters( + gamma_in=1, gamma_in_sigma=0, m2l=5, m2l_sigma=0 + ) + assert param_draw[0] == 1 + assert param_draw[1] == 5 + for i in range(100): + param_draw = self.scaling_nfw_2d.draw_lens_parameters( + gamma_in=1, gamma_in_sigma=1, m2l=5, 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) + m2l_array = np.linspace(start=1, stop=10, num=10) + + param_arrays = [ani_param_array, gamma_in_array, m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.multiply.outer(gamma_in_array, 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, m2l=5, m2l_sigma=0.1 + ) + with self.assertRaises(ValueError): + scaling.draw_lens_parameters( + gamma_in=1, gamma_in_sigma=0.1, m2l=-5, 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) + m2l_array = np.linspace(start=1, stop=10, num=10) + + param_arrays = [ani_param_array, gamma_in_array, m2l_array] + param_scaling_array = np.multiply.outer( + ani_param_array, + np.multiply.outer(gamma_in_array, 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, m2l=0, m2l_sigma=0 + ) + + +if __name__ == "__main__": + pytest.main() diff --git a/test/test_Sampling/test_ParamManager/test_lens_param.py b/test/test_Sampling/test_ParamManager/test_lens_param.py index 577438e..5b09406 100644 --- a/test/test_Sampling/test_ParamManager/test_lens_param.py +++ b/test/test_Sampling/test_ParamManager/test_lens_param.py @@ -91,3 +91,86 @@ def test_args2kwargs(self): kwargs_new, i = self._param_fixed.args2kwargs(args, i=0) args_new = self._param_fixed.kwargs2args(kwargs_new) npt.assert_almost_equal(args_new, args) + + +class TestLensParamGammaInnerM2l(object): + def setup(self): + self._param = LensParam( + gamma_in_sampling=True, + gamma_in_distribution="GAUSSIAN", + m2l_sampling=True, + m2l_distribution="GAUSSIAN", + alpha_gamma_in_sampling=True, + alpha_m2l_sampling=True, + kwargs_fixed={}, + log_scatter=False, + ) + + kwargs_fixed = { + "gamma_in": 1, + "gamma_in_sigma": 0.1, + "m2l": 5, + "m2l_sigma": 1, + "alpha_gamma_in": 0.1, + "alpha_m2l": -0.5, + } + self._param_fixed = LensParam( + gamma_in_sampling=True, + gamma_in_distribution="GAUSSIAN", + m2l_sampling=True, + m2l_distribution="GAUSSIAN", + alpha_gamma_in_sampling=True, + alpha_m2l_sampling=True, + kwargs_fixed=kwargs_fixed, + log_scatter=False, + ) + self._param_log_scatter = LensParam( + gamma_in_sampling=True, + gamma_in_distribution="GAUSSIAN", + m2l_sampling=True, + m2l_distribution="GAUSSIAN", + alpha_gamma_in_sampling=True, + alpha_m2l_sampling=True, + kwargs_fixed={}, + log_scatter=True, + ) + + def test_param_list(self): + param_list = self._param.param_list(latex_style=False) + assert len(param_list) == 6 + param_list = self._param.param_list(latex_style=True) + assert len(param_list) == 6 + + param_list = self._param_log_scatter.param_list(latex_style=False) + assert len(param_list) == 6 + param_list = self._param_log_scatter.param_list(latex_style=True) + assert len(param_list) == 6 + + param_list = self._param_fixed.param_list(latex_style=False) + assert len(param_list) == 0 + param_list = self._param_fixed.param_list(latex_style=True) + assert len(param_list) == 0 + + def test_args2kwargs(self): + kwargs = { + "gamma_in": 1, + "gamma_in_sigma": 0.1, + "m2l": 5, + "m2l_sigma": 1, + "alpha_gamma_in": 0.1, + "alpha_m2l": -0.5, + } + args = self._param.kwargs2args(kwargs) + kwargs_new, i = self._param.args2kwargs(args, i=0) + args_new = self._param.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args) + + args = self._param_log_scatter.kwargs2args(kwargs) + kwargs_new, i = self._param_log_scatter.args2kwargs(args, i=0) + args_new = self._param_log_scatter.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args) + + args = self._param_fixed.kwargs2args(kwargs) + kwargs_new, i = self._param_fixed.args2kwargs(args, i=0) + args_new = self._param_fixed.kwargs2args(kwargs_new) + npt.assert_almost_equal(args_new, args)