diff --git a/clmm/__init__.py b/clmm/__init__.py index ee4a5e9c3..7d98303e4 100644 --- a/clmm/__init__.py +++ b/clmm/__init__.py @@ -2,7 +2,7 @@ from .gcdata import GCData from .galaxycluster import GalaxyCluster from .clusterensemble import ClusterEnsemble -from .dataops import compute_tangential_and_cross_components, make_radial_profile +from .dataops import compute_tangential_and_cross_components, make_radial_profile, make_binned_estimators_triaxiality from .utils import compute_radial_averages, make_bins, convert_units from .theory import ( compute_reduced_shear_from_convergence, diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py index a93cba312..cd04cd60b 100644 --- a/clmm/dataops/__init__.py +++ b/clmm/dataops/__init__.py @@ -2,6 +2,7 @@ import warnings import numpy as np import scipy +from scipy.stats import binned_statistic from astropy.coordinates import SkyCoord from astropy import units as u from ..gcdata import GCData @@ -32,7 +33,7 @@ def compute_tangential_and_cross_components( geometry="curve", is_deltasigma=False, sigma_c=None, - is_quadrupole=False, + include_quadrupole=False, phi_major=None, ra_mem=None, dec_mem=None, @@ -117,21 +118,21 @@ def compute_tangential_and_cross_components( sigma_c : None, array_like Critical (effective) surface density in units of :math:`M_\odot\ Mpc^{-2}`. Used only when is_deltasigma=True. - is_quadrupole: bool + include_quadrupole: bool If `True`, the quadrupole shear components (g_4theta, g_const; Shin+2018) are calculated instead of g_t and g_x phi_major : float, optional the direction of the major axis of the input cluster in the unit of radian. - only needed when `is_quadupole` is `True`. + only needed when `include_quadrupole` is `True`. Users could choose to provide ra_mem, dec_mem and weight_mem instead of this quantity. ra_mem : array, optional right ascentions of the member galaxies of the input cluster, to calculate the direction of the major axis of the cluster. - Only needed when `is_quadrupole` is `True` and `phi_major` is not provided. + Only needed when `include_quadrupole` is `True` and `phi_major` is not provided. dec_mem : array, optional declinations of the member galaxies of the input cluster, to calculate the direction of the major axis of the cluster. - Only needed when `is_quadrupole` is `True` and `phi_major` is not provided. + Only needed when `include_quadrupole` is `True` and `phi_major` is not provided. weight_mem : array, optional weights for the member galaxies to be applied when calcualting the major axis of the input cluster out of ra_mem and dec_mem. @@ -143,16 +144,15 @@ def compute_tangential_and_cross_components( ------- angsep: array_like Angular separation between lens and each source galaxy in radians - IF is_quadrupole: + tangential_component: array_like + Tangential shear (or assimilated quantity) for each source galaxy + cross_component: array_like + Cross shear (or assimilated quantity) for each source galaxy + IF include_quadrupole: 4theta_component: array_like 4theta shear component (or assimilated quantity) for each source galaxy const_component: array_like constant shear component (or assimilated quantity) for each source galaxy - ELSE: - tangential_component: array_like - Tangential shear (or assimilated quantity) for each source galaxy - cross_component: array_like - Cross shear (or assimilated quantity) for each source galaxy """ # pylint: disable-msg=too-many-locals # Note: we make these quantities to be np.array so that a name is not passed from astropy @@ -166,7 +166,7 @@ def compute_tangential_and_cross_components( validate_argument(locals(), "shear2", "float_array") validate_argument(locals(), "geometry", str) validate_argument(locals(), "is_deltasigma", bool) - validate_argument(locals(), "is_quadrupole", bool) + validate_argument(locals(), "include_quadrupole", bool) validate_argument(locals(), "sigma_c", "float_array", none_ok=True) ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency( [ra_source, dec_source, shear1, shear2], @@ -174,6 +174,12 @@ def compute_tangential_and_cross_components( prefix="Tangential- and Cross- shape components sources", ) _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c) + if phi_major is not None: + validate_argument(locals(), "phi_major", float) + else: + validate_argument(locals(), "ra_mem", "float_array") + validate_argument(locals(), "dec_mem", "float_array") + validate_argument(locals(), "weight_mem", "float_array") elif np.iterable(ra_source): ra_source_, dec_source_, shear1_, shear2_ = ( np.array(col) for col in [ra_source, dec_source, shear1, shear2] @@ -192,7 +198,16 @@ def compute_tangential_and_cross_components( angsep, phi = _compute_lensing_angles_astropy(ra_lens, dec_lens, ra_source_, dec_source_) else: raise NotImplementedError(f"Sky geometry {geometry} is not currently supported") - if is_quadupole: + # Compute the tangential and cross shears + tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi) + cross_comp = _compute_cross_shear(shear1_, shear2_, phi) + # If the is_deltasigma flag is True, multiply the results by Sigma_crit. + if sigma_c is not None: + _sigma_c_arr = np.array(sigma_c) + tangential_comp *= _sigma_c_arr + cross_comp *= _sigma_c_arr + + if include_quadrupole: if phi_major is not None: phi_major_ = phi_major else: @@ -206,23 +221,14 @@ def compute_tangential_and_cross_components( const_comp = rotated_shear1 # If the is_deltasigma flag is True, multiply the results by Sigma_crit. if sigma_c is not None: - _sigma_c_arr = np.array(sigma_c) four_theta_comp *= _sigma_c_arr - const_comp *= _sigma_c_arr - else: - # Compute the tangential and cross shears - tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi) - cross_comp = _compute_cross_shear(shear1_, shear2_, phi) - # If the is_deltasigma flag is True, multiply the results by Sigma_crit. - if sigma_c is not None: - _sigma_c_arr = np.array(sigma_c) - tangential_comp *= _sigma_c_arr - cross_comp *= _sigma_c_arr - if is_quadrupole: - return angsep, four_theta_comp, const_comp + const_comp *= _sigma_c_arr + if include_quadrupole: + return angsep, tangential_comp, cross_comp, four_theta_comp, const_comp else: return angsep, tangential_comp, cross_comp + def compute_background_probability( z_lens, z_src=None, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True ): @@ -689,7 +695,7 @@ def make_stacked_radial_profile(angsep, weights, components): Parameters ---------- angsep: 2d array - Transvesal distances corresponding to each object with shape `n_obj, n_rad_bins`. + Transversal distances corresponding to each object with shape `n_obj, n_rad_bins`. weights: 2d array Weights corresponding to each objects with shape `n_obj, n_rad_bins`. components: list of 2d arrays @@ -708,3 +714,189 @@ def make_stacked_radial_profile(angsep, weights, components): np.average(component, axis=0, weights=weights) for component in components ] return staked_angsep, stacked_components + +def measure_Delta_Sigma_const_triaxiality(w, gamma1, Sigma_crit) : + """Measure Delta sigma const quadrupole value for individual galaxies. + + Parameters + ---------- + w: float or array + Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366) + gamma1: float or array + Shear component gamma 1 + Sigma_crit: float + Critical surface mass density in M_{sun}/Mpc^{2} + + Returns + ------- + dsconst_data: The measured "4theta" quadrupole component from the input shear measurements for each galaxy + """ + dsconst_data = Sigma_crit * gamma1 + return dsconst_data + +def measure_Delta_Sigma_4theta_triaxiality(w1, w2, gamma1, gamma2, theta, Sigma_crit) : + """Measure Delta sigma 4theta quadrupole value for individual galaxies. + + Parameters + ---------- + w1: float or array + Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366) + w2: float or array + Weight computed for each galaxy by the descrition as given in Shin et al. (https://doi.org/10.1093/mnras/stx3366) + gamma1: float or array + Shear component gamma 1 + gamma2: float or array + Shear component gamma 2 + theta: float or array + Position Angle of galaxy measured counter-clockwise from the gamma1 direction + Sigma_crit: float + Critical surface mass density in M_{sun}/Mpc^{2} + + Returns + ------- + ds4theta_data: The measured "const" quadrupole component from the input shear measurements for each galaxy + """ + ds4theta_data = Sigma_crit * (w1*gamma1/np.cos(4*theta) + w2*gamma2/np.sin(4*theta)) / (w1 + w2) + return ds4theta_data + +def measure_weights_triaxiality(Sigma_crit, theta, Sigma_shape=0.0001, Sigma_meas=0) : + """Measure weight values for quadrupole measurement of individual galaxies. + Using Equations 35, 32, 33 from Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366) + + Parameters + ---------- + Sigma_crit: float + Critical surface mass density in M_{sun}/Mpc^{2} + theta: float or array + Position Angle of galaxy measured counter-clockwise from the gamma1 direction + Sigma_shape: Float, optional + Scatter in the shape measurement from shape noise + Sigma_meas: Float, optional + Measurement error for shapes + + Returns + ------- + w, w1, w2: float or array + Weights for each galaxy + """ + + w = 1 / (Sigma_crit**2 * (Sigma_shape**2 + Sigma_meas**2)) + w1 = np.cos(4*theta)**2 * w + w2 = np.sin(4*theta)**2 * w + return w, w1, w2 + +def make_binned_estimators_triaxiality(gamma1, gamma2, x_arcsec, y_arcsec, z_cluster, z_source, cosmo, monopole_bins=15, + quadrupole_bins=15, r_min=0.3, r_max=2.5): + """ + Makes radially binned profiles for Monopole, Quadrupole (4theta and const) as described in Shin et al. 2018 + Parameters + ---------- + gamma1: float or array + Shear component gamma 1 + gamma2: float or array + Shear component gamma 2 + x_arcsec: float or array + X position in arcsec when cluster center is at 0.0 and x axis is assumed to be along the major axis of the cluster + y_arcsec: float or array + Y position in arcsec when cluster center is at 0.0 + z_cluster: float + Redshift of lens cluster + cosmo: clmm.cosmology.Cosmology object + CLMM Cosmology object + monopole_bins: float, optional + Number of bins for monopole lensing profile + quadrupole_bins: float, optional + Number of bins for quadrupole lensing profiles (4theta, const) + r_min: Float, optional + Minimum radial distance to measure monopole and quadrupole shear profiles in Mpc (Default is 0.3 Mpc) + r_max: Float, optional + Maximum radial distance to measure monopole and quadrupole shear profiles in Mpc (Default is 2.5 Mpc) + + Returns + ------- + ds_mono: float or array + Binned monopole lensing profile + ds_mono_err: float or array + Uncertainty for the binned monopole lensing shear profile + r_mono: float or array + Bin centers for the monopole lensing shear profile + ds_4theta: float or array + Binned binned quadrupole 4-theta lensing shear profile + ds_4theta_err: float or array + Uncertainty for the binned quadrupole 4-theta lensing shear profile + ds_const: float or array + Binned binned quadrupole const lensing shear profile + ds_const_err: float or array + Uncertainty for the binned quadrupole const lensing shear profile + r_quad: float or array + Bin centers for the quadrupole lensing shear profiles (4theta and const) + """ + + sigma_crit = cosmo.eval_sigma_crit(z_cluster,z_source) + r = np.sqrt((x_arcsec**2 + y_arcsec**2)) #In arcsecs + theta = np.arctan2(y_arcsec, x_arcsec) + r_mpc = r*cosmo.eval_da(z_cluster) * np.pi/180.0 * 1/3600 #In Mpc + + w, w1, w2 = measure_weights_triaxiality(sigma_crit, theta) + DS4theta = measure_Delta_Sigma_4theta_triaxiality(w1, w2, gamma1, gamma2, theta, sigma_crit) + DSconst = measure_Delta_Sigma_const_triaxiality(w, gamma1, sigma_crit) + + #Binning for Quadrupole measurements DS4theta and DSconst + bins=quadrupole_bins + r_min = r_min + r_max = r_max + + bin_edges = np.logspace(np.log10(r_min), np.log10(r_max), bins) + N_i = [] + for i in np.arange(bins-1): + N_i.append(len(r_mpc[(r_mpc > bin_edges[i]) & (r_mpc < bin_edges[i+1])])) + N_i=np.array(N_i) + + + result = binned_statistic(r_mpc, gamma1, statistic='mean', bins=bin_edges) + gamma1_i = result.statistic + res = binned_statistic(r_mpc, gamma2, statistic='mean', bins=bin_edges) + gamma2_i = res.statistic + res = binned_statistic(r_mpc, DS4theta, statistic='mean', bins=bin_edges) + DS4theta_i_err = binned_statistic(r_mpc, DS4theta, statistic='std', bins=bin_edges).statistic/np.sqrt(N_i) + DS4theta_i = res.statistic + res = binned_statistic(r_mpc, DSconst, statistic='mean', bins=bin_edges) + DSconst_i = res.statistic + DSconst_i_err = binned_statistic(r_mpc, DSconst, statistic='std', bins=bin_edges).statistic/np.sqrt(N_i) + r_i = bin_edges + + #Binning for Monopole Measurements: + bins_mono=monopole_bins + bin_edges = np.logspace(np.log10(r_min), np.log10(r_max), bins_mono) + N_i = [] + for i in np.arange(bins_mono-1): + N_i.append(len(r_mpc[(r_mpc > bin_edges[i]) & (r_mpc < bin_edges[i+1])])) + N_i=np.array(N_i) + + r_mono = bin_edges + res = binned_statistic(r_mpc, -gamma1*np.cos(2*theta)-gamma2*np.sin(2*theta), statistic='mean', bins=bin_edges) + gammat_mono = res.statistic + ds_mono_err = binned_statistic(r_mpc, -gamma1*np.cos(2*theta)-gamma2*np.sin(2*theta), statistic='std', + bins=bin_edges).statistic/np.sqrt(N_i)*sigma_crit + ds_mono = gammat_mono*sigma_crit + + # SAFEGUARD AGAINST BINS WITH NANs and 0.0s + + ind = np.invert(np.isnan(ds_mono) | np.isnan(ds_mono_err)) + ds_mono = ds_mono[ind] + ds_mono_err = ds_mono_err[ind] + r_mono = np.sqrt(r_mono[:-1]*r_mono[1:])[ind] + ind = (ds_mono!= 0.0) & (ds_mono_err!= 0.0) + ds_mono = ds_mono[ind] + ds_mono_err = np.abs(ds_mono_err[ind]) + r_mono = r_mono[ind] + + ind = np.invert(np.isnan(DS4theta_i) | np.isnan(DS4theta_i_err) | np.isnan(DSconst_i) | np.isnan(DSconst_i_err)) + ds_4theta = DS4theta_i[ind] + ds_4theta_err = np.abs(DS4theta_i_err[ind]) + ds_const = DSconst_i[ind] + ds_const_err = np.abs(DSconst_i_err[ind]) + r_quad = np.sqrt(r_i[:-1]*r_i[1:])[ind] + + return ds_mono,ds_mono_err,r_mono,ds_4theta,ds_4theta_err,ds_const,ds_const_err,r_quad + diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py index b4db5446c..35a109c65 100644 --- a/clmm/galaxycluster.py +++ b/clmm/galaxycluster.py @@ -33,31 +33,57 @@ class GalaxyCluster: Declination of galaxy cluster center (in degrees) z : float Redshift of galaxy cluster center + phi_major : float + Direction of the cluster major axis (in radian). + Users may provide `ra_mem`, `dec_mem` and `weight_mem` instead of this parameter. + Only needed when `include_quadrupole==True`. + ra_mem : array_like + RAs of member galaxies (in degrees). + Only needed when `include_quadrupole==True`. + dec_mem : array_like + DECs of member galaxies (in degrees). + Only needed when `include_quadrupole==True`. + weight_mem : array_like + weights of member galaxies (to calculate major axis). + Only needed when `include_quadrupole==True`. galcat : GCData Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z validate_input: bool Validade each input argument + include_quadrupole: bool + If quadrupole WL be calculated. """ - def __init__(self, *args, validate_input=True, **kwargs): + def __init__(self, *args, validate_input=True, include_quadrupole=False, **kwargs): self.unique_id = None self.ra = None self.dec = None self.z = None self.galcat = None + self.phi_major = None + self.ra_mem = None + self.dec_mem = None + self.weight_mem = None self.validate_input = validate_input + self.include_quadrupole = include_quadrupole if len(args) > 0 or len(kwargs) > 0: self._add_values(*args, **kwargs) self._check_types() self.set_ra_lower(ra_low=0) - def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: GCData): + def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: GCData, + phi_major=None, ra_mem=None, dec_mem=None, weight_mem=None): """Add values for all attributes""" self.unique_id = unique_id self.ra = ra self.dec = dec self.z = z self.galcat = galcat + if self.include_quadrupole: + self.phi_major = phi_major + self.ra_mem = ra_mem + self.dec_mem = dec_mem + self.weight_mem = weight_mem def _check_types(self): """Check types of all attributes""" @@ -70,6 +96,17 @@ def _check_types(self): self.ra = float(self.ra) self.dec = float(self.dec) self.z = float(self.z) + if self.include_quadrupole: + if self.phi_major is not None: + validate_argument(vars(self), "phi_major", float) + self.phi_major = float(self.phi_major) + else: + validate_argument(vars(self), "ra_mem", "float_array") + validate_argument(vars(self), "dec_mem", "float_array") + validate_argument(vars(self), "weight_mem", "float_array") + self.ra_mem = np.array(self.ra_mem, dtype=float) + self.dec_mem = np.array(self.dec_mem, dtype=float) + self.weight_mem = np.array(self.weight_mem, dtype=float) def save(self, filename, **kwargs): """Saves GalaxyCluster object to filename using Pickle""" @@ -208,6 +245,8 @@ def compute_tangential_and_cross_components( shape_component2="e2", tan_component="et", cross_component="ex", + quad_4theta_component="e4theta", + quad_const_component="econst", geometry="curve", is_deltasigma=False, use_pdz=False, @@ -225,6 +264,13 @@ def compute_tangential_and_cross_components( shear2: `galcat` shape_component2 geometry: `input` geometry is_deltasigma: `input` is_deltasigma + include_quadrupole: `input` include_quadrupole + [if include_quadrupole:] + phi_major: `cluster` major axis direction (in radian with respect to +x) + (OR) + ra_mem: `cluster` RAs of member galaxies for calculating major axis of a given cluster + dec_mem: `cluster` DECs of member galaxies for calculating major axis of a given cluster + weight_mem: `cluster` weights to be applied to member galaxies for calculating major axis Parameters ---------- @@ -242,6 +288,14 @@ def compute_tangential_and_cross_components( Name of the column to be added to the `galcat` astropy table that will contain the cross component computed from columns `shape_component1` and `shape_component2`. Default: `ex` + quad_4theta_component: string, optional + Name of the column to be added to the `galcat` astropy table that will contain the + 4theta quarupole component computed from columns `shape_component1` and `shape_component2`. + Default: `e4theta` + quad_const_component: string, optional + Name of the column to be added to the `galcat` astropy table that will contain the + constant quarupole component computed from columns `shape_component1` and `shape_component2`. + Default: `econst` geometry: str, optional Sky geometry to compute angular separation. Options are curve (uses astropy) or flat. @@ -261,6 +315,11 @@ def compute_tangential_and_cross_components( Tangential shear (or assimilated quantity) for each source galaxy cross_component: array_like Cross shear (or assimilated quantity) for each source galaxy + [if self.include_quadrupole:] + quad_4theta_component: array_like + 4theta quadrupole shear (or assimilated quantity) for each source galaxy + quad_const_component: array_like + constatnt quadrupole shear (or assimilated quantity) for each source galaxy """ # Check is all the required data is available col_dict = { @@ -275,23 +334,51 @@ def compute_tangential_and_cross_components( cols = self._get_input_galdata(col_dict) # compute shears - angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components( - is_deltasigma=is_deltasigma, - ra_lens=self.ra, - dec_lens=self.dec, - geometry=geometry, - validate_input=self.validate_input, - **cols, - ) - if add: - self.galcat["theta"] = angsep - self.galcat[tan_component] = tangential_comp - self.galcat[cross_component] = cross_comp - if is_deltasigma: - sigmac_type = "effective" if use_pdz else "standard" - self.galcat.meta[f"{tan_component}_sigmac_type"] = sigmac_type - self.galcat.meta[f"{cross_component}_sigmac_type"] = sigmac_type - return angsep, tangential_comp, cross_comp + if self.include_quadrupole: + angsep, tangential_comp, cross_comp, four_theta_comp, const_comp = compute_tangential_and_cross_components( + is_deltasigma=is_deltasigma, + ra_lens=self.ra, + dec_lens=self.dec, + geometry=geometry, + validate_input=self.validate_input, + include_quadrupole=self.include_quadrupole, + phi_major=self.phi_major, + ra_mem=self.ra_mem, + dec_mem=self.dec_mem, + weight_mem=self.weight_mem, + **cols, + ) + if add: + self.galcat["theta"] = angsep + self.galcat[tan_component] = tangential_comp + self.galcat[cross_component] = cross_comp + self.galcat[quad_4theta_component] = four_theta_comp + self.galcat[quad_const_component] = const_comp + if is_deltasigma: + sigmac_type = "effective" if use_pdz else "standard" + self.galcat.meta[f"{tan_component}_sigmac_type"] = sigmac_type + self.galcat.meta[f"{cross_component}_sigmac_type"] = sigmac_type + self.galcat.meta[f"{quad_4theta_component}_sigmac_type"] = sigmac_type + self.galcat.meta[f"{quad_const_component}_sigmac_type"] = sigmac_type + return angsep, tangential_comp, cross_comp, four_theta_comp, const_comp + else: + angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components( + is_deltasigma=is_deltasigma, + ra_lens=self.ra, + dec_lens=self.dec, + geometry=geometry, + validate_input=self.validate_input, + **cols, + ) + if add: + self.galcat["theta"] = angsep + self.galcat[tan_component] = tangential_comp + self.galcat[cross_component] = cross_comp + if is_deltasigma: + sigmac_type = "effective" if use_pdz else "standard" + self.galcat.meta[f"{tan_component}_sigmac_type"] = sigmac_type + self.galcat.meta[f"{cross_component}_sigmac_type"] = sigmac_type + return angsep, tangential_comp, cross_comp def compute_background_probability( self, use_pdz=False, add=True, p_background_name="p_background" @@ -468,10 +555,16 @@ def make_radial_profile( cosmo=None, tan_component_in="et", cross_component_in="ex", + quad_4theta_component_in="e4theta", + quad_const_component_in="econst", tan_component_out="gt", cross_component_out="gx", + quad_4theta_component_out="g4theta", + quad_const_component_out="gconst", tan_component_in_err=None, cross_component_in_err=None, + quad_4theta_component_in_err=None, + quad_const_component_in_err=None, include_empty_bins=False, gal_ids_in_bins=False, add=True, @@ -488,6 +581,8 @@ def make_radial_profile( Calls `clmm.dataops.make_radial_profile` with the following arguments: components: `galcat` components (tan_component_in, cross_component_in, z) + OR + (quad_4theta_component_in, quad_const_component_in, z) IF include_quadrupole angsep: `galcat` theta angsep_units: 'radians' bin_units: `input` bin_units @@ -521,18 +616,36 @@ def make_radial_profile( cross_component_in: string, optional Name of the cross component column in `galcat` to be binned. Default: 'ex' + quad_4theta_component_in: string, optional + Name of the 4theta quadrupole component column in `galcat` to be binned. + Default: 'e4theta' + quad_const_component_in: string, optional + Name of the constant quadrupole component column in `galcat` to be binned. + Default: 'econst' tan_component_out: string, optional Name of the tangetial component binned column to be added in profile table. Default: 'gt' cross_component_out: string, optional Name of the cross component binned profile column to be added in profile table. Default: 'gx' + quad_4theta_component_out: string, optional + Name of the 4theta quadrupole component binned column to be added in profile table. + Default: 'g4theta' + quad_const_component_out: string, optional + Name of the constant quadrupole component binned profile column to be added in profile table. + Default: 'gconst' tan_component_in_err: string, None, optional Name of the tangential component error column in `galcat` to be binned. Default: None cross_component_in_err: string, None, optional Name of the cross component error column in `galcat` to be binned. Default: None + quad_4theta_component_in_err: string, None, optional + Name of the 4theta quadrupole component error column in `galcat` to be binned. + Default: None + quad_const_component_in_err: string, None, optional + Name of the constant quadrupole component error column in `galcat` to be binned. + Default: None include_empty_bins: bool, optional Also include empty bins in the returned table gal_ids_in_bins: bool, optional @@ -561,7 +674,6 @@ def make_radial_profile( """ # Too many local variables (19/15) # pylint: disable=R0914 - if not all( t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta") ): @@ -570,31 +682,72 @@ def make_radial_profile( "and cross shears (gt, gx) or ellipticities (et, ex). " "Run compute_tangential_and_cross_components first." ) + if self.include_quadrupole: + if not all( + t_ in self.galcat.columns for t_ in (quad_4theta_component_in, quad_const_component_in, "theta") + ): + raise TypeError( + "Shear or ellipticity information is missing. Galaxy catalog must have 4theta" + "and constant quadrupole shears (g4theta, gconst) or ellipticities (e4theta, econst). " + "Run compute_tangential_and_cross_components first." + ) if "z" not in self.galcat.columns: raise TypeError("Missing galaxy redshifts!") # Compute the binned averages and associated errors - profile_table, binnumber = make_radial_profile( - [self.galcat[n].data for n in (tan_component_in, cross_component_in, "z")], - angsep=self.galcat["theta"], - angsep_units="radians", - bin_units=bin_units, - bins=bins, - error_model=error_model, - include_empty_bins=include_empty_bins, - return_binnumber=True, - cosmo=cosmo, - z_lens=self.z, - validate_input=self.validate_input, - components_error=[ - None if n is None else self.galcat[n].data - for n in (tan_component_in_err, cross_component_in_err, None) - ], - weights=self.galcat[weights_in].data if use_weights else None, - ) - # Reaname table columns - for i, name in enumerate([tan_component_out, cross_component_out, "z"]): - profile_table.rename_column(f"p_{i}", name) - profile_table.rename_column(f"p_{i}_err", f"{name}_err") + if self.include_quadrupole: + profile_table, binnumber = make_radial_profile( + [self.galcat[n].data for n in + (tan_component_in, cross_component_in, + quad_4theta_component_in, quad_const_component_in, "z" + ) + ], + angsep=self.galcat["theta"], + angsep_units="radians", + bin_units=bin_units, + bins=bins, + error_model=error_model, + include_empty_bins=include_empty_bins, + return_binnumber=True, + cosmo=cosmo, + z_lens=self.z, + validate_input=self.validate_input, + components_error=[ + None if n is None else self.galcat[n].data + for n in (tan_component_in_err, cross_component_in_err, + quad_4theta_component_in_err, quad_const_component_in_err, None + ) + ], + weights=self.galcat[weights_in].data if use_weights else None, + ) + # Reaname table columns + for i, name in enumerate([tan_component_out, cross_component_out, + quad_4theta_component_out, quad_const_component_out, "z"] + ): + profile_table.rename_column(f"p_{i}", name) + profile_table.rename_column(f"p_{i}_err", f"{name}_err") + else: + profile_table, binnumber = make_radial_profile( + [self.galcat[n].data for n in (tan_component_in, cross_component_in, "z")], + angsep=self.galcat["theta"], + angsep_units="radians", + bin_units=bin_units, + bins=bins, + error_model=error_model, + include_empty_bins=include_empty_bins, + return_binnumber=True, + cosmo=cosmo, + z_lens=self.z, + validate_input=self.validate_input, + components_error=[ + None if n is None else self.galcat[n].data + for n in (tan_component_in_err, cross_component_in_err, None) + ], + weights=self.galcat[weights_in].data if use_weights else None, + ) + # Reaname table columns + for i, name in enumerate([tan_component_out, cross_component_out, "z"]): + profile_table.rename_column(f"p_{i}", name) + profile_table.rename_column(f"p_{i}_err", f"{name}_err") # Reaname weights columns profile_table.rename_column("weights_sum", weights_out) # add galaxy IDs @@ -626,6 +779,10 @@ def plot_profiles( tangential_component_error="gt_err", cross_component="gx", cross_component_error="gx_err", + quad_4theta_component="g4theta", + quad_4theta_component_error="g4theta_err", + quad_const_component="gconst", + quad_const_component_error="gconst_err", table_name="profile", xscale="linear", yscale="linear", @@ -646,6 +803,18 @@ def plot_profiles( cross_component_error: str, optional Name of the column in the galcat Table corresponding to the uncertainty in the cross component of the shear or reduced shear. Default: 'gx_err' + quad_4theta_component: str, optional + Name of the column in the galcat Table corresponding to the 4theta quadrupole component of + the shear or reduced shear (Delta Sigma not yet implemented). Default: 'g4theta' + quad_4theta_component_error: str, optional + Name of the column in the galcat Table corresponding to the uncertainty in 4theta quadrupole + component of the shear or reduced shear. Default: 'g4theta_err' + quad_const_component: str, optional + Name of the column in the galcat Table corresponding to the constant quadrupole component of the + shear or reduced shear. Default: 'gconst' + quad_const_component_error: str, optional + Name of the column in the galcat Table corresponding to the uncertainty in the constant quadrupole + component of the shear or reduced shear. Default: 'gconst_er' table_name: str, optional Name of the GalaxyCluster() `.profile` attribute. Default: 'profile' xscale: @@ -663,30 +832,76 @@ def plot_profiles( if not hasattr(self, table_name): raise ValueError(f"GalaxyClusters does not have a '{table_name}' table.") profile = getattr(self, table_name) - for col in (tangential_component, cross_component): - if col not in profile.columns: - raise ValueError(f"Column for plotting '{col}' does not exist.") - for col in (tangential_component_error, cross_component_error): - if col not in profile.columns: - warnings.warn(f"Column for plotting '{col}' does not exist.") - return plot_profiles( - rbins=profile["radius"], - r_units=profile.meta["bin_units"], - tangential_component=profile[tangential_component], - tangential_component_error=( - profile[tangential_component_error] - if tangential_component_error in profile.columns - else None - ), - cross_component=profile[cross_component], - cross_component_error=( - profile[cross_component_error] if cross_component_error in profile.columns else None - ), - xscale=xscale, - yscale=yscale, - tangential_component_label=tangential_component, - cross_component_label=cross_component, - ) + if self.include_quadrupole: + for col in (tangential_component, cross_component, quad_4theta_component, quad_const_component): + if col not in profile.columns: + raise ValueError(f"Column for plotting '{col}' does not exist.") + for col in (tangential_component_error, cross_component_error, + quad_4theta_component_error, quad_const_component_error + ): + if col not in profile.columns: + warnings.warn(f"Column for plotting '{col}' does not exist.") + return plot_profiles( + rbins=profile["radius"], + r_units=profile.meta["bin_units"], + tangential_component=profile[tangential_component], + tangential_component_error=( + profile[tangential_component_error] + if tangential_component_error in profile.columns + else None + ), + cross_component=profile[cross_component], + cross_component_error=( + profile[cross_component_error] if cross_component_error in profile.columns else None + ), + xscale=xscale, + yscale=yscale, + tangential_component_label=tangential_component, + cross_component_label=cross_component, + ), plot_profiles( + rbins=profile["radius"], + r_units=profile.meta["bin_units"], + tangential_component=profile[quad_4theta_component], + tangential_component_error=( + profile[quad_4theta_component_error] + if quad_4theta_component_error in profile.columns + else None + ), + cross_component=profile[quad_const_component], + cross_component_error=( + profile[quad_const_component_error] if quad_const_component_error in profile.columns else None + ), + xscale=xscale, + yscale=yscale, + tangential_component_label=quad_4theta_component, + cross_component_label=quad_const_component, + ) + + else: + for col in (tangential_component, cross_component): + if col not in profile.columns: + raise ValueError(f"Column for plotting '{col}' does not exist.") + for col in (tangential_component_error, cross_component_error): + if col not in profile.columns: + warnings.warn(f"Column for plotting '{col}' does not exist.") + return plot_profiles( + rbins=profile["radius"], + r_units=profile.meta["bin_units"], + tangential_component=profile[tangential_component], + tangential_component_error=( + profile[tangential_component_error] + if tangential_component_error in profile.columns + else None + ), + cross_component=profile[cross_component], + cross_component_error=( + profile[cross_component_error] if cross_component_error in profile.columns else None + ), + xscale=xscale, + yscale=yscale, + tangential_component_label=tangential_component, + cross_component_label=cross_component, + ) def set_ra_lower(self, ra_low=0): """ diff --git a/clmm/theory/func_layer.py b/clmm/theory/func_layer.py index d69e629a0..699d462e0 100644 --- a/clmm/theory/func_layer.py +++ b/clmm/theory/func_layer.py @@ -1308,7 +1308,7 @@ def compute_delta_sigma_4theta_triaxiality( verbose=False, validate_input=True, ): - r"""Compute the "4-theta"-quadrupole value + r"""Compute the "4-theta" component of quadrupole shear as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366) Parameters ---------- @@ -1351,7 +1351,7 @@ def compute_delta_sigma_4theta_triaxiality( Returns ------- ds4theta: array - Delta sigma 4-theta value at a position for the elliptical halo specified + Delta sigma 4-theta component value at a position for the elliptical halo specified """ ### DEFINING INTEGRALS: r_arr = np.linspace(0.01, 3*np.max(r_source), sample_N) @@ -1390,7 +1390,7 @@ def compute_delta_sigma_const_triaxiality( verbose=False, validate_input=True, ): - r"""Compute the excess surface density lensing profile for "const"-quadrupole value + r"""Compute the "const" component of quadrupole shear as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366) Parameters ---------- @@ -1474,7 +1474,8 @@ def compute_delta_sigma_excess_triaxiality( verbose=False, validate_input=True, ): - r"""Compute the excess surface density lensing profile for the monopole component along with second order expansion term (e**2) + r"""Compute the excess surface density lensing profile for the monopole component along + with second order expansion term (e**2) as given in Shin et al. 2018 (https://doi.org/10.1093/mnras/stx3366) Parameters ---------- diff --git a/examples/demo_mock_ensemble.ipynb b/examples/demo_mock_ensemble.ipynb index 84354b111..a14eb088f 100644 --- a/examples/demo_mock_ensemble.ipynb +++ b/examples/demo_mock_ensemble.ipynb @@ -618,9 +618,9 @@ ], "metadata": { "kernelspec": { - "display_name": "wrk", + "display_name": "Python [conda env:tf-gpu]", "language": "python", - "name": "wrk" + "name": "conda-env-tf-gpu-py" }, "language_info": { "codemirror_mode": { @@ -632,7 +632,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.9.18" } }, "nbformat": 4,