Merge branch 'issue/591/triaxiality' of https://github.com/LSSTDESC/CLMM

 into issue/591/triaxiality

merge in dataops updates

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,

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,14 +166,20 @@ 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],
             names=("Ra", "Dec", "Shear1", "Shear2"),
             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
+