Skip to content

Commit

Permalink
changed galaxycluser.py for quadrupole
Browse files Browse the repository at this point in the history
  • Loading branch information
@tae-h-shin
tae-h-shin committed Jul 17, 2024
1 parent 8525ba7 commit 40358a7
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 41 deletions.
4 changes: 3 additions & 1 deletion clmm/dataops/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,7 @@ def compute_tangential_and_cross_components(
_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
Expand All @@ -221,12 +222,13 @@ def compute_tangential_and_cross_components(
# If the is_deltasigma flag is True, multiply the results by Sigma_crit.
if sigma_c is not None:
four_theta_comp *= _sigma_c_arr
const_comp *= _sigma_c_arr
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
):
Expand Down
109 changes: 69 additions & 40 deletions clmm/galaxycluster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,25 +36,25 @@ class GalaxyCluster:
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 `is_quadrupole==True`.
Only needed when `include_quadrupole==True`.
ra_mem : array_like
RAs of member galaxies (in degrees).
Only needed when `is_quadrupole==True`.
Only needed when `include_quadrupole==True`.
dec_mem : array_like
DECs of member galaxies (in degrees).
Only needed when `is_quadrupole==True`.
Only needed when `include_quadrupole==True`.
weight_mem : array_like
weights of member galaxies (to calculate major axis).
Only needed when `is_quadrupole==True`.
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
is_quadrupole: bool
include_quadrupole: bool
If quadrupole WL be calculated.
"""

def __init__(self, *args, validate_input=True, is_quadrupole=False, **kwargs):
def __init__(self, *args, validate_input=True, include_quadrupole=False, **kwargs):
self.unique_id = None
self.ra = None
self.dec = None
Expand All @@ -65,7 +65,7 @@ def __init__(self, *args, validate_input=True, is_quadrupole=False, **kwargs):
self.dec_mem = None
self.weight_mem = None
self.validate_input = validate_input
self.is_quadrupole = is_quadrupole
self.include_quadrupole = include_quadrupole
if len(args) > 0 or len(kwargs) > 0:
self._add_values(*args, **kwargs)
self._check_types()
Expand All @@ -79,7 +79,7 @@ def _add_values(self, unique_id: str, ra: float, dec: float, z: float, galcat: G
self.dec = dec
self.z = z
self.galcat = galcat
if self.is_quadrupole:
if self.include_quadrupole:
self.phi_major = phi_major
self.ra_mem = ra_mem
self.dec_mem = dec_mem
Expand All @@ -96,7 +96,7 @@ def _check_types(self):
self.ra = float(self.ra)
self.dec = float(self.dec)
self.z = float(self.z)
if self.is_quadrupole:
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)
Expand Down Expand Up @@ -264,8 +264,8 @@ def compute_tangential_and_cross_components(
shear2: `galcat` shape_component2
geometry: `input` geometry
is_deltasigma: `input` is_deltasigma
is_quadrupole: `input` is_quadrupole
[if is_quadrupole:]
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
Expand Down Expand Up @@ -311,12 +311,11 @@ def compute_tangential_and_cross_components(
-------
angsep: array_like
Angular separation between lens and each source galaxy in radians
[if not self.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 self.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 self.include_quadrupole:]
quad_4theta_component: array_like
4theta quadrupole shear (or assimilated quantity) for each source galaxy
quad_const_component: array_like
Expand All @@ -335,14 +334,14 @@ def compute_tangential_and_cross_components(
cols = self._get_input_galdata(col_dict)

# compute shears
if self.is_quadrupole:
angsep, four_theta_comp, const_comp = compute_tangential_and_cross_components(
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,
is_quadupole=self.is_quarupole,
include_quadrupole=self.include_quadrupole,
phi_major=self.phi_major,
ra_mem=self.ra_mem,
dec_mem=self.dec_mem,
Expand All @@ -351,13 +350,17 @@ def compute_tangential_and_cross_components(
)
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, four_theta_comp, const_comp
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,
Expand Down Expand Up @@ -579,7 +582,7 @@ 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 is_quadrupole
(quad_4theta_component_in, quad_const_component_in, z) IF include_quadrupole
angsep: `galcat` theta
angsep_units: 'radians'
bin_units: `input` bin_units
Expand Down Expand Up @@ -671,8 +674,15 @@ def make_radial_profile(
"""
# Too many local variables (19/15)
# pylint: disable=R0914

if self.is_quadrupole:
if not all(
t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
):
raise TypeError(
"Shear or ellipticity information is missing. Galaxy catalog must have tangential"
"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")
):
Expand All @@ -681,21 +691,16 @@ def make_radial_profile(
"and constant quadrupole shears (g4theta, gconst) or ellipticities (e4theta, econst). "
"Run compute_tangential_and_cross_components first."
)
else:
if not all(
t_ in self.galcat.columns for t_ in (tan_component_in, cross_component_in, "theta")
):
raise TypeError(
"Shear or ellipticity information is missing. Galaxy catalog must have tangential"
"and cross shears (gt, gx) or ellipticities (et, ex). "
"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
if self.is_quadrupole:
if self.include_quadrupole:
profile_table, binnumber = make_radial_profile(
[self.galcat[n].data for n in (quad_4theta_component_in, quad_const_component_in, "z")],
[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,
Expand All @@ -708,12 +713,16 @@ def make_radial_profile(
validate_input=self.validate_input,
components_error=[
None if n is None else self.galcat[n].data
for n in (quad_4theta_component_in_err, quad_const_component_in_err, None)
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([quad_4theta_component_out, quad_const_component_out, "z"]):
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:
Expand Down Expand Up @@ -823,14 +832,33 @@ 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)
if self.is_quadrupole:
for col in (quad_4theta_component, quad_const_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 (quad_4theta_component_error, quad_const_component_error):
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],
Expand All @@ -848,6 +876,7 @@ def plot_profiles(
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:
Expand Down

0 comments on commit 40358a7

Please sign in to comment.