# This file is part of pipe_tasks.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""Cluster lensing shear profile measurement and multiplicative bias estimation.
Provides :class:`HaloMcBiasMultibandPipe`, a pipeline task that matches
detected sources to truth catalogs around NFW halos, measures tangential
shear profiles, and estimates multiplicative bias from the ratio of
measured-to-true reduced shear.
"""
__all__ = [
"HaloMcBiasMultibandPipeConfig",
"HaloMcBiasMultibandPipe",
"HaloMcBiasMultibandPipeConnections",
]
from typing import Any
import astropy.units as u
import lsst.pipe.base.connectionTypes as cT
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import SkyCoord
from lsst.pex.config import Field
from lsst.pipe.base import (
PipelineTask,
PipelineTaskConfig,
PipelineTaskConnections,
Struct,
)
from lsst.skymap import BaseSkyMap
from lsst.utils.logging import LsstLogAdapter
from scipy.interpolate import interp1d
from scipy.spatial import cKDTree
[docs]
class HaloMcBiasMultibandPipeConnections(
PipelineTaskConnections,
dimensions=("skymap", "band"),
defaultTemplates={
"coaddName": "deep",
"dataType": "",
"version": "",
},
):
"""Butler connections for :class:`HaloMcBiasMultibandPipe`."""
[docs]
skymap = cT.Input(
doc="SkyMap to use in processing",
name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
storageClass="SkyMap",
dimensions=("skymap",),
)
[docs]
src00List = cT.Input(
doc="Source catalog with all the measurement generated in this task",
name="{coaddName}_rot0_coadd_anacal_{dataType}",
dimensions=(
"skymap",
"tract",
"patch",
),
storageClass="ArrowAstropy",
multiple=True,
deferLoad=True,
)
[docs]
src01List = cT.Input(
doc="Source catalog with all the measurement generated in this task",
name="{coaddName}_rot1_coadd_anacal_{dataType}",
dimensions=(
"skymap",
"tract",
"patch",
),
storageClass="ArrowAstropy",
multiple=True,
deferLoad=True,
)
[docs]
truth00List = cT.Input(
doc="input truth catalog",
name="{coaddName}_rot0_coadd_truthCatalog",
storageClass="ArrowAstropy",
dimensions=(
"skymap",
"tract",
),
multiple=True,
deferLoad=True,
)
[docs]
truth01List = cT.Input(
doc="input truth catalog",
name="{coaddName}_rot1_coadd_truthCatalog",
storageClass="ArrowAstropy",
dimensions=(
"skymap",
"tract",
),
multiple=True,
deferLoad=True,
)
[docs]
outputSummary = cT.Output(
doc="Summary statistics",
name="{coaddName}_halo_mc_{dataType}_summary_stats",
storageClass="ArrowAstropy",
dimensions=("skymap",),
)
[docs]
summaryPlot = cT.Output(
doc="simple plot of summary stats",
storageClass="Plot",
name="{coaddName}_halo_mc_summary_{dataType}_plot",
dimensions=("skymap",),
)
def __init__(self, *, config=None):
super().__init__(config=config)
[docs]
class HaloMcBiasMultibandPipeConfig(
PipelineTaskConfig,
pipelineConnections=HaloMcBiasMultibandPipeConnections,
):
"""Configuration for :class:`HaloMcBiasMultibandPipe`."""
[docs]
ename = Field[str](
doc="ellipticity column name",
default="e",
)
[docs]
xname = Field[str](
doc="detection coordinate row name",
default="x1_det",
)
[docs]
yname = Field[str](
doc="detection coordinate column name",
default="x2_det",
)
[docs]
wname = Field[str](
doc="weight column name",
default="wsel",
)
[docs]
mass = Field[float](
doc="halo mass",
default=5e-14,
)
[docs]
conc = Field[float](
doc="halo concertration",
default=1.0,
)
[docs]
z_lens = Field[float](
doc="halo redshift",
default=1.0,
)
[docs]
z_source = Field[float](
doc="source redshift",
default=None,
)
[docs]
def validate(self):
super().validate()
if len(self.connections.dataType) == 0:
raise ValueError("connections.dataTape missing")
[docs]
class HaloMcBiasMultibandPipe(PipelineTask):
"""Measure tangential shear profiles around NFW halos and estimate
multiplicative shear bias.
Matches detected source catalogs to truth catalogs, computes
tangential and cross shear in radial bins, and compares measured
shear to the true input shear to derive the multiplicative bias.
"""
[docs]
_DefaultName = "HaloMcTask"
[docs]
ConfigClass = HaloMcBiasMultibandPipeConfig
def __init__(
self,
*,
config: HaloMcBiasMultibandPipeConfig | None = None,
log: LsstLogAdapter | None = None,
initInputs: dict[str, Any] | None = None,
**kwargs: Any,
):
super().__init__(
config=config,
log=log,
initInputs=initInputs,
**kwargs,
)
assert isinstance(
self.config,
HaloMcBiasMultibandPipeConfig,
)
[docs]
self.ename = self.config.ename
[docs]
self.egname = lambda x, y: (
"fpfs_d" + self.ename + str(x) + "_" + "dg" + str(y)
)
[docs]
self.xname = self.config.xname
[docs]
self.yname = self.config.yname
[docs]
self.wname = self.config.wname
[docs]
self.wgname = lambda x: "d" + self.wname + "_dg" + str(x)
return
[docs]
def runQuantum(self, butlerQC, inputRefs, outputRefs):
assert isinstance(
self.config,
HaloMcBiasMultibandPipeConfig,
)
# Retrieve the filename of the input exposure
inputs = butlerQC.get(inputRefs)
tractList = []
for ref in inputRefs.src00List:
tract_id = ref.dataId["tract"]
tractList.append(tract_id)
outputs = self.run(
**inputs, tractList=tractList,
)
butlerQC.put(outputs, outputRefs)
return
@staticmethod
[docs]
def _rotate_spin_2_vec(e1, e2, angle):
"""
Rotate a spin-2 field by an array of angles (one per e1, e2 pair)
"""
# Ensure e1, e2, and angle are numpy arrays of the same length
e1 = np.asarray(e1)
e2 = np.asarray(e2)
angle = np.asarray(angle)
assert (
e1.shape == e2.shape == angle.shape
), "e1, e2, and angle must have the same shape"
# Create an empty output array for the rotated values
output = np.zeros((2, len(e1)))
# Compute cos(2*angle) and sin(2*angle) for each angle
cos_2angle = np.cos(2 * angle)
sin_2angle = np.sin(2 * angle)
# Apply the rotation for each e1, e2 pair
# invert the sign of output so the tangential shear is positive
output[0] = (cos_2angle * e1 - sin_2angle * e2) * (-1) # Rotated e1
output[1] = (sin_2angle * e1 + cos_2angle * e2) * (-1) # Rotated e2
return output
@staticmethod
[docs]
def _rotate_spin_2_matrix(R11, R22, angle):
# off diagonal terms are assumed to be zero
# R is active rotation matrix
output = np.zeros((2, len(R11)))
output[0] = np.cos(2 * angle) ** 2 * R11 + np.sin(2 * angle) ** 2 * R22
output[1] = np.sin(2 * angle) ** 2 * R11 + np.cos(2 * angle) ** 2 * R22
return output
@staticmethod
[docs]
def _get_response_from_w_and_der(e1, e2, w, e1_g1, e2_g2, w_g1, w_g2):
r11 = e1_g1 * w + e1 * w_g1
r22 = e2_g2 * w + e2 * w_g2
return r11, r22
@staticmethod
[docs]
def _get_eT_eX_rT_rX_sum(
eT,
eX,
w,
rT,
rX,
gT_true_matched,
gX_true_matched,
gT_true,
gX_true,
kappa_true,
dist,
lensed_shift,
radial_lensed_shift,
radial_bin_edges,
match_dist,
m00,
m20,
all_true_gT,
all_true_gX,
all_true_dist,
):
"""calculate the sum of eT, eX, and rT in each radial bin for a single
halo
Args:
eT (array): tangential shape
eX (array): cross shape
w (weight): ancal weight
rT (array): tangential response,
rX (array): cross response,
gT_true_matched (array): true tangential shear of matched objects
gX_true_matched (array): true cross shear of matched objects
gT_true (array): true tangential shear calculated with distance
gX_true (array): true cross shear calculated with distance
kappa_true (array): true convergence
dist (array): angular distance from the halo center
lensed_shift (array):
distance between the lensed and prelensed position
radial_lensed_shift (array):
distance between the lensed and prelensed position in radial
direction
radial_bin_edges (array): radial bin edges in pixel
match_dist (array):
the distance between detection and matched input
m00 (array): fpfs shapelet mode m00
m20 (array): fpfs shapelet mode m20
all_true_gT (array): true tangential shear of all input objects
all_true_gX (array): true cross shear of matched objects
all_true_dist (array):
angular distance from the halo center for all input objects
Returns:
eT(array): sum of eT in each radial bin
eX(array): sum of eX in each radial bin
rT(array): sum of tangential resposne in each radial bin
rX(array): sum of radial response in each radial bin
gT_true_matched(array):
sum of true tangential shear of matches in each radial bin
gX_true_matched(array):
sum of true cross shear of matches in each radial bin
gT_true (array):
sum of true tan shear in each radial bin
gX_true (array):
sum of true cross shear in each radial bin
all_true_gT(array):
mean of true tan shear of all input objects in each radial bin
all_true_gX(array):
mean of true cross shear of all input objects in each radial bin
kappa_true(array): sum of true convergence in each radial bin
lensed_shift(array): mean of lensed shift in each radial bin
radial_lensed_shift(array):
mean of lensed shift in each radial bin, projected on the
radial direction
r_weighted_gT_matched(array):
sum of tangential shear weighted by rT in each radial bin
r_weighted_gX_matched(array):
sum of cross shear weighted by rX in each radial bin
r_weighted_gT(array):
sum of tangential shear weighted by rT in each radial bin
r_weighted_gX(array):
sum of cross shear weighted by rX in each radial bin
ngal_in_bin(array): number of galaxies in each radial bin
eT_std_list(array):
per galaxy standard deviation of eT in each radial bin
eX_std_list(array):
per galaxy standard deviation of eX in each radial bin
mean_dist(array): mean distance of galaxies of the halo center
median_match_dist_list(array):
median of the match distance in each radial bin, expected to be
around 0.5
match_failure_rate_list(array):
fraction of match distance larger than 2 in each radial bin
m00_list(array): mean of shapelet mode m00 in each radial bin
m20_list(array): mean of shapelet mode m20 in each radial bin
"""
n_bins = len(radial_bin_edges) - 1
# this list stores results in each radial bin
eT_list = []
eX_list = []
rT_list = []
rX_list = []
gT_true_matched_list = []
gX_true_matched_list = []
gT_true_list = []
gX_true_list = []
kappa_true_list = []
r_weighted_gT_matched_list = []
r_weighted_gX_matched_list = []
r_weighted_gT_list = []
r_weighted_gX_list = []
ngal_in_bin = []
eT_std_list = []
eX_std_list = []
lensed_shift_list = []
radial_lensed_shift_list = []
mean_dist = []
median_matched_dist_list = []
match_failure_rate_list = []
m00_list = []
m20_list = []
all_true_gT_list = []
all_true_gX_list = []
for i_bin in range(n_bins):
mask = (dist >= radial_bin_edges[i_bin]) & (
dist < radial_bin_edges[i_bin + 1]
)
all_mask = (all_true_dist >= radial_bin_edges[i_bin]) & (
all_true_dist < radial_bin_edges[i_bin + 1]
)
eT_sum = np.sum(eT[mask] * w[mask])
eX_sum = np.sum(eX[mask] * w[mask])
rT_sum = np.sum(rT[mask])
rX_sum = np.sum(rX[mask])
gT_true_matched_list.append(np.sum(gT_true_matched[mask]))
gX_true_matched_list.append(np.sum(gX_true_matched[mask]))
gT_true_list.append(np.sum(gT_true[mask]))
gX_true_list.append(np.sum(gX_true[mask]))
all_true_gT_list.append(np.mean(all_true_gT[all_mask]))
all_true_gX_list.append(np.mean(all_true_gX[all_mask]))
kappa_true_list.append(np.sum(kappa_true[mask]))
r_weighted_gT_matched_list.append(
np.sum(gT_true_matched[mask] * rT[mask])
)
r_weighted_gX_matched_list.append(
np.sum(gX_true_matched[mask] * rX[mask])
)
r_weighted_gT_list.append(np.sum(gT_true[mask] * rT[mask]))
r_weighted_gX_list.append(np.sum(gX_true[mask] * rX[mask]))
ngal_in_bin.append(np.sum(mask))
eT_list.append(eT_sum)
eX_list.append(eX_sum)
rT_list.append(rT_sum)
rX_list.append(rX_sum)
eT_std_list.append(np.std(eT[mask]))
eX_std_list.append(np.std(eX[mask]))
lensed_shift_list.append(np.mean(lensed_shift[mask]))
radial_lensed_shift_list.append(np.mean(radial_lensed_shift[mask]))
mean_dist.append(np.mean(dist[mask]))
median_matched_dist_list.append(np.median(match_dist[mask]))
match_failure_rate_list.append(
np.sum(match_dist[mask] > 2) / np.sum(mask)
)
m00_list.append(np.mean(m00[mask]))
m20_list.append(np.mean(m20[mask]))
return (
np.array(eT_list),
np.array(eX_list),
np.array(rT_list),
np.array(rX_list),
np.array(gT_true_matched_list),
np.array(gX_true_matched_list),
np.array(gT_true_list),
np.array(gX_true_list),
np.array(all_true_gT_list),
np.array(all_true_gX_list),
np.array(kappa_true_list),
np.array(lensed_shift_list),
np.array(radial_lensed_shift_list),
np.array(r_weighted_gT_matched_list),
np.array(r_weighted_gX_matched_list),
np.array(r_weighted_gT_list),
np.array(r_weighted_gX_list),
np.array(ngal_in_bin),
np.array(eT_std_list),
np.array(eX_std_list),
np.array(mean_dist),
np.array(median_matched_dist_list),
np.array(match_failure_rate_list),
np.array(m00_list),
np.array(m20_list),
)
@staticmethod
@staticmethod
[docs]
def get_summary_struct(n_halos, n_bins):
dt = [
(
"angular_bin_left",
f"({n_bins},)f8",
),
(
"angular_bin_right",
f"({n_bins},)f8",
),
("ngal_in_bin", f"({n_bins},)i4"),
("eT", f"({n_bins},)f8"),
("eT_std", f"({n_bins},)f8"),
("eX", f"({n_bins},)f8"),
("eX_std", f"({n_bins},)f8"),
("rT", f"({n_bins},)f8"),
("rT_std", f"({n_bins},)f8"),
("rX", f"({n_bins},)f8"),
("rX_std", f"({n_bins},)f8"),
("gT_true_matched", f"({n_bins},)f8"),
("gX_true_matched", f"({n_bins},)f8"),
("gT_true", f"({n_bins},)f8"),
("gX_true", f"({n_bins},)f8"),
("all_true_gT", f"({n_bins},)f8"),
("all_true_gX", f"({n_bins},)f8"),
("kappa_true", f"({n_bins},)f8"),
("lensed_shift", f"({n_bins},)f8"),
(
"radial_lensed_shift",
f"({n_bins},)f8",
),
("r_weighted_gT_matched", f"({n_bins},)f8"),
("r_weighted_gX_matched", f"({n_bins},)f8"),
("r_weighted_gT", f"({n_bins},)f8"),
("r_weighted_gX", f"({n_bins},)f8"),
(
"median_match_dist",
f"({n_bins},)f8",
),
("mean_dist", f"({n_bins},)f8"),
(
"match_failure_rate",
f"({n_bins},)f8",
),
("mean_m00", f"({n_bins},)f8"),
("mean_m20", f"({n_bins},)f8"),
]
return np.zeros(n_halos, dtype=dt)
@staticmethod
[docs]
def get_per_galaxy_struct(n_gal):
dt = [
("rT", "f8"),
("gT", "f8"),
("x", "f8"),
("y", "f8"),
("lensed_x", "f8"),
("lensed_y", "f8"),
("dist", "f8"),
]
return np.zeros(n_gal, dtype=dt)
@staticmethod
[docs]
def angsep(ra, dec, ra2=200.0, dec2=0.0):
c1 = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs")
c2 = SkyCoord(ra=ra2 * u.deg, dec=dec2 * u.deg, frame="icrs")
return c1.separation(c2).arcsec
@staticmethod
[docs]
def position_angle_ccw_from_east(ra1, dec1, ra2, dec2):
"""
Compute the PA of (ra2,dec2) relative to (ra1,dec1),
measured from +RA (east=0°), increasing CCW towards +Dec (north).
Returns an Angle in [0,360)°.
"""
c1 = SkyCoord(ra1 * u.deg, dec1 * u.deg, frame="icrs")
c2 = SkyCoord(ra2 * u.deg, dec2 * u.deg, frame="icrs")
# pa_north = angle east of north (0° at north, + toward east)
pa_north = c1.position_angle(c2)
# shift zero to east and make it CCW (north positive)
pa_ccw = (90 * u.deg - pa_north).wrap_at(360 * u.deg)
return pa_ccw
@staticmethod
[docs]
def generate_summary_plot(summary_table):
area = np.mean(
np.pi
* (
(summary_table["angular_bin_right"] / 3600 * np.pi / 180.0) ** 2
- (summary_table["angular_bin_left"] / 3600 * np.pi / 180.0)
** 2
)
* (60 * 180 / np.pi) ** 2,
axis=0,
)
fig, axes = plt.subplots(6, 2, figsize=(12, 13))
angular_bin_mid = (
(
np.mean(
summary_table["angular_bin_left"],
axis=0,
)
+ np.mean(
summary_table["angular_bin_right"],
axis=0,
)
)
/ 2
/ 60
)
start = 0
mean_eX = np.mean(summary_table["eX"], axis=0)
mean_rX = np.mean(summary_table["rX"], axis=0)
sigma_eX = np.std(summary_table["eX"], axis=0) / np.sqrt(
len(summary_table)
)
sigma_rX = np.std(summary_table["rX"], axis=0) / np.sqrt(
len(summary_table)
)
sigma_gX = np.sqrt(
(sigma_eX / mean_rX) ** 2 + (mean_eX * sigma_rX / mean_rX**2) ** 2
)
true_gX_matched = np.mean(
summary_table["r_weighted_gX_matched"], axis=0
) / np.mean(summary_table["rX"], axis=0)
true_gX = np.mean(summary_table["r_weighted_gX"], axis=0) / np.mean(
summary_table["rX"], axis=0
)
axes[0, 0].errorbar(
angular_bin_mid[start:],
(mean_eX / mean_rX)[start:],
label="measured gX",
yerr=sigma_gX[start:],
)
axes[0, 0].plot(
angular_bin_mid[start:],
true_gX_matched[start:],
label="r weighted matched true gX",
)
axes[0, 0].plot(
angular_bin_mid[start:],
true_gX[start:],
label="r weighted true gX",
)
# axes[0, 0].plot(
# angular_bin_mid[start:],
# np.mean(summary_table["all_true_gX"], axis=0)[start:],
# label="input gX",
# )
axes[0, 0].set_ylim(-0.01, 0.01)
axes[0, 0].set_ylabel(r"$g^X$")
axes[0, 0].set_xlabel("arcmin")
axes[0, 0].legend()
axes[1, 0].errorbar(
angular_bin_mid[start:],
(mean_eX / mean_rX - true_gX)[start:],
yerr=sigma_gX[start:],
)
axes[1, 0].set_ylabel(r"$g^X_m - g^X_t$")
mean_eT = np.mean(summary_table["eT"], axis=0)
mean_rT = np.mean(summary_table["rT"], axis=0)
sigma_eT = np.std(summary_table["eT"], axis=0) / np.sqrt(
len(summary_table)
)
sigma_rT = np.std(summary_table["rT"], axis=0) / np.sqrt(
len(summary_table)
)
sigma_gT = np.sqrt(
(sigma_eT / mean_rT) ** 2 + (mean_eT * sigma_rT / mean_rT**2) ** 2
)
true_gT_matched = np.mean(
summary_table["r_weighted_gT_matched"], axis=0
) / np.mean(summary_table["rT"], axis=0)
true_gT = np.mean(summary_table["r_weighted_gT"], axis=0) / np.mean(
summary_table["rT"], axis=0
)
axes[0, 1].errorbar(
angular_bin_mid[start:],
(mean_eT / mean_rT)[start:],
label="measured gT",
yerr=sigma_gT[start:],
)
axes[0, 1].plot(
angular_bin_mid[start:],
true_gT_matched[start:],
label="r weighted matched gT",
)
axes[0, 1].plot(
angular_bin_mid[start:],
true_gT[start:],
label="r weighted true gT",
)
axes[0, 1].legend()
axes[0, 1].set_xlabel("arcmin")
axes[0, 1].set_ylabel(r"$g^T$")
axes[1, 1].errorbar(
angular_bin_mid[start:],
(mean_eT / mean_rT / true_gT - 1)[start:],
label="true gT",
yerr=(sigma_gT / np.abs(true_gT))[start:],
)
axes[1, 1].set_ylabel(r"$g^T_{m}/g^T_{t} - 1$")
axes[3, 0].errorbar(
angular_bin_mid,
np.mean(
summary_table["kappa_true"],
axis=0,
)
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
yerr=np.std(
summary_table["kappa_true"],
axis=0,
)
/ np.sqrt(len(summary_table))
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
)
axes[3, 0].set_ylabel(r"$\langle \kappa \rangle$")
axes[3, 0].set_xlabel("arcmin")
axes[2, 0].errorbar(
angular_bin_mid,
np.mean(summary_table["rX"], axis=0)
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
yerr=np.std(summary_table["rX"], axis=0)
/ np.sqrt(len(summary_table))
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
)
axes[2, 0].set_xlabel("arcmin")
axes[2, 0].set_ylabel(r"$\langle r^X \rangle$")
# axes[1, 0].set_ylim(0.16, 0.185)
axes[2, 1].errorbar(
angular_bin_mid,
np.mean(summary_table["rT"], axis=0)
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
yerr=np.std(summary_table["rT"], axis=0)
/ np.sqrt(len(summary_table))
/ np.mean(
summary_table["ngal_in_bin"],
axis=0,
),
)
axes[2, 1].set_xlabel("arcmin")
axes[2, 1].set_ylabel(r"$\langle r^T \rangle$")
# axes[1, 1].set_ylim(0.16, 0.185)
axes[3, 1].errorbar(
angular_bin_mid,
np.mean(
summary_table["ngal_in_bin"],
axis=0,
)
/ area,
yerr=np.std(
summary_table["ngal_in_bin"],
axis=0,
)
/ np.sqrt(len(summary_table))
/ area,
label=f"n_halo={len(summary_table)}",
)
axes[3, 1].set_xlabel("arcmin")
axes[3, 1].set_ylabel(r"Detection Number density [arcmin$^{-2}$]")
axes[3, 1].legend()
axes[4, 0].errorbar(
angular_bin_mid,
np.mean(
summary_table["lensed_shift"],
axis=0,
),
yerr=np.std(
summary_table["lensed_shift"],
axis=0,
)
/ np.sqrt(len(summary_table)),
label="shift",
)
axes[4, 0].errorbar(
angular_bin_mid,
np.mean(
summary_table["radial_lensed_shift"],
axis=0,
),
yerr=np.std(
summary_table["radial_lensed_shift"],
axis=0,
)
/ np.sqrt(len(summary_table)),
label="radial_shift",
)
axes[4, 0].set_ylabel("Lens shift [arcsec]")
axes[4, 0].legend()
axes[4, 1].plot(
angular_bin_mid,
np.mean(
summary_table["median_match_dist"],
axis=0,
),
label="match failure rate",
)
axes[4, 1].set_ylabel("Median match distance pixel")
axes[4, 1].set_ylim(0, 0.5)
axes[5, 0].plot(
angular_bin_mid,
np.mean(
summary_table["match_failure_rate"],
axis=0,
),
label="match failure rate",
)
axes[5, 0].set_ylabel("Match failure rate (> 2 pixel)")
axes[5, 0].set_ylim(0, 0.6)
axes[5, 1].errorbar(
angular_bin_mid,
np.mean(
summary_table["mean_m00"] + summary_table["mean_m20"],
axis=0,
),
yerr=(
np.std(
summary_table["mean_m00"],
axis=0,
)
+ np.std(
summary_table["mean_m20"],
axis=0,
)
)
/ np.sqrt(len(summary_table)),
)
axes[5, 1].set_ylabel(r"$M_{00} + M_{20}$")
axes[5, 1].legend()
for ax in axes.flatten():
ax.set_xlim(0, 10)
ax.set_xlabel("Angular separation [arcmin]")
plt.subplots_adjust(hspace=0.2, wspace=0.2)
return fig
[docs]
def run(
self,
*,
skymap,
src00List,
src01List,
truth00List,
truth01List,
tractList,
**kwargs,
):
self.log.info("load truth list")
self.log.info(f"len truth00List: {len(truth00List)}")
self.log.info(f"len truth01List: {len(truth01List)}")
pixel_scale = skymap.config.pixelScale # arcsec per pixel
image_dim = skymap.config.patchInnerDimensions[0] # in pixels
max_pixel = (image_dim - 64) / 2
self.log.info("image dim", image_dim)
self.log.info("pixel scale", pixel_scale)
self.log.info("max pixel", max_pixel)
self.log.info(
"max pixel in arcsec",
max_pixel * pixel_scale,
)
n_bins = 15
# starting from 15 arcsec
pixel_bin_edges = np.linspace(15.0 / 0.2, max_pixel, n_bins + 1)
angular_bin_edges = pixel_bin_edges * pixel_scale
en = self.ename
e1n = "fpfs_" + en + "1"
e2n = "fpfs_" + en + "2"
e1g1n = self.egname(1, 1)
e2g2n = self.egname(2, 2)
xn = self.xname
yn = self.yname
wn = self.wname
wg1n = self.wgname(1)
wg2n = self.wgname(2)
self.log.info(
"The length of source list is",
len(src00List),
len(src01List),
)
n_realization = len(src00List)
# binned arrays
rT_ensemble = np.empty((len(src00List), n_bins))
rX_ensemble = np.empty((len(src00List), n_bins))
eT_ensemble = np.empty((len(src00List), n_bins))
eX_ensemble = np.empty((len(src00List), n_bins))
gT_true_matched_ensemble = np.empty((len(src00List), n_bins))
gX_true_matched_ensemble = np.empty((len(src00List), n_bins))
gT_true_ensemble = np.empty((len(src00List), n_bins))
gX_true_ensemble = np.empty((len(src00List), n_bins))
all_true_gT_ensemble = np.empty((len(src00List), n_bins))
all_true_gX_ensemble = np.empty((len(src00List), n_bins))
r_weighted_gT_matched_ensemble = np.empty((len(src00List), n_bins))
r_weighted_gX_matched_ensemble = np.empty((len(src00List), n_bins))
r_weighted_gT_ensemble = np.empty((len(src00List), n_bins))
r_weighted_gX_ensemble = np.empty((len(src00List), n_bins))
kappa_true_ensemble = np.empty((len(src00List), n_bins))
ngal_in_bin_ensemble = np.empty((len(src00List), n_bins))
eT_std_ensemble = np.empty((len(src00List), n_bins))
eX_std_ensemble = np.empty((len(src00List), n_bins))
lensed_shift_ensemble = np.empty((len(src00List), n_bins))
radial_lensed_shift_ensemble = np.empty((len(src00List), n_bins))
mean_dist_ensemble = np.empty((len(src00List), n_bins))
median_match_dist_ensemble = np.empty((len(src00List), n_bins))
match_failure_rate_ensemble = np.empty((len(src00List), n_bins))
m00_ensemble = np.empty((len(src00List), n_bins))
m20_ensemble = np.empty((len(src00List), n_bins))
# detection length arrays
ind_rT_list = []
ind_gT_true_list = []
ind_x_list = []
ind_y_list = []
ind_dist_list = []
ind_lensed_x_list = []
ind_lensed_y_list = []
for i, cats in enumerate(
zip(
src00List,
src01List,
truth00List,
truth01List,
tractList,
)
):
src00, src01, truth00, truth01, tract_id = (
cats[0],
cats[1],
cats[2],
cats[3],
cats[4],
)
sr_00_res = src00.get()
sr_01_res = src01.get()
truth_00_res = truth00.get()
truth_01_res = truth01.get()
wcs = skymap[tract_id].getWcs()
true_x_00, true_y_00 = wcs.skyToPixelArray(
truth_00_res["ra"], truth_00_res["dec"], degrees=True,
)
true_x_01, true_y_01 = wcs.skyToPixelArray(
truth_01_res["ra"], truth_01_res["dec"], degrees=True,
)
idx_00, match_dist_00 = self._match_input_to_det(
true_x_00,
true_y_00,
sr_00_res[xn] / pixel_scale,
sr_00_res[yn] / pixel_scale,
)
idx_01, match_dist_01 = self._match_input_to_det(
true_x_01,
true_y_01,
sr_01_res[xn] / pixel_scale,
sr_01_res[yn] / pixel_scale,
)
match_dist = np.concatenate([match_dist_00, match_dist_01])
gamma1_true = np.concatenate(
[
truth_00_res["gamma1"][idx_00],
truth_01_res["gamma1"][idx_01],
]
)
gamma2_true = np.concatenate(
[
truth_00_res["gamma2"][idx_00],
truth_01_res["gamma2"][idx_01],
]
)
kappa_true = np.concatenate(
[
truth_00_res["kappa"][idx_00],
truth_01_res["kappa"][idx_01],
]
)
g1_true = gamma1_true / (1 - kappa_true)
g2_true = gamma2_true / (1 - kappa_true)
e1 = np.concatenate([sr_00_res[e1n], sr_01_res[e1n]])
e2 = np.concatenate([sr_00_res[e2n], sr_01_res[e2n]])
print(f"i: {i}, e1: {e1.shape}, e2: {e2.shape}")
det_x = np.concatenate(
[
sr_00_res[xn] / pixel_scale,
sr_01_res[xn] / pixel_scale,
]
)
det_y = np.concatenate(
[
sr_00_res[yn] / pixel_scale,
sr_01_res[yn] / pixel_scale,
]
)
det_ra, det_dec = wcs.pixelToSkyArray(det_x, det_y, degrees=True)
e1_g1 = np.concatenate([sr_00_res[e1g1n], sr_01_res[e1g1n]])
e2_g2 = np.concatenate([sr_00_res[e2g2n], sr_01_res[e2g2n]])
w = np.concatenate([sr_00_res[wn], sr_01_res[wn]])
w_g1 = np.concatenate([sr_00_res[wg1n], sr_01_res[wg1n]])
w_g2 = np.concatenate([sr_00_res[wg2n], sr_01_res[wg2n]])
m00 = np.concatenate([sr_00_res["fpfs_m0"], sr_01_res["fpfs_m0"]])
m20 = np.concatenate([sr_00_res["fpfs_m2"], sr_01_res["fpfs_m2"]])
self.log.info(f"i: {i}, e1: {e1.shape}, e2: {e2.shape}")
e1_g1 = np.concatenate(
[
sr_00_res[e1g1n],
sr_01_res[e1g1n],
]
)
e2_g2 = np.concatenate(
[
sr_00_res[e2g2n],
sr_01_res[e2g2n],
]
)
w = np.concatenate([sr_00_res[wn], sr_01_res[wn]])
w_g1 = np.concatenate(
[
sr_00_res[wg1n],
sr_01_res[wg1n],
]
)
w_g2 = np.concatenate(
[
sr_00_res[wg2n],
sr_01_res[wg2n],
]
)
m00 = np.concatenate(
[
sr_00_res["fpfs_m0"],
sr_01_res["fpfs_m0"],
]
)
m20 = np.concatenate(
[
sr_00_res["fpfs_m2"],
sr_01_res["fpfs_m2"],
]
)
# Convert prelensed ra/dec to pixel positions via WCS
pre_ra = np.concatenate(
[
truth_00_res["prelensed_ra"][idx_00],
truth_01_res["prelensed_ra"][idx_01],
]
)
pre_dec = np.concatenate(
[
truth_00_res["prelensed_dec"][idx_00],
truth_01_res["prelensed_dec"][idx_01],
]
)
x, y = wcs.skyToPixelArray(pre_ra, pre_dec, degrees=True)
# Convert lensed ra/dec to pixel positions via WCS
lensed_ra = np.concatenate(
[
truth_00_res["ra"][idx_00],
truth_01_res["ra"][idx_01],
]
)
lensed_dec = np.concatenate(
[
truth_00_res["dec"][idx_00],
truth_01_res["dec"][idx_01],
]
)
lensed_x, lensed_y = wcs.skyToPixelArray(
lensed_ra, lensed_dec, degrees=True,
)
all_true_ra = np.concatenate(
[truth_00_res["ra"], truth_01_res["ra"]]
)
all_true_dec = np.concatenate(
[truth_00_res["dec"], truth_01_res["dec"]]
)
all_true_gamma1 = np.concatenate(
[truth_00_res["gamma1"], truth_01_res["gamma1"]]
)
all_true_gamma2 = np.concatenate(
[truth_00_res["gamma2"], truth_01_res["gamma2"]]
)
all_true_kappa = np.concatenate(
[truth_00_res["kappa"], truth_01_res["kappa"]]
)
all_true_g1 = all_true_gamma1 / (1 - all_true_kappa)
all_true_g2 = all_true_gamma2 / (1 - all_true_kappa)
ind_x_list.append(x)
ind_y_list.append(y)
ind_lensed_x_list.append(lensed_x)
ind_lensed_y_list.append(lensed_y)
pix_origin = wcs.getPixelOrigin()
lensed_shift = (
np.sqrt((lensed_x - x) ** 2 + (lensed_y - y) ** 2) * pixel_scale
)
radial_dist_lensed = np.sqrt(
(lensed_x - pix_origin.x) ** 2 + (lensed_y - pix_origin.y) ** 2
)
radial_dist = np.sqrt(
(x - pix_origin.x) ** 2 + (y - pix_origin.y) ** 2
)
radial_lensed_shift = (
radial_dist_lensed - radial_dist
) * pixel_scale
sky_origin = wcs.getSkyOrigin()
sky_ra = sky_origin.getRa().asDegrees()
sky_dec = sky_origin.getDec().asDegrees()
angle = self.position_angle_ccw_from_east(
sky_ra, sky_dec, det_ra, det_dec
).rad
all_true_angle = self.position_angle_ccw_from_east(
sky_ra,
sky_dec,
all_true_ra,
all_true_dec,
).rad
# negative since we are rotating axes
eT, eX = self._rotate_spin_2_vec(e1, e2, angle)
gT_true_matched, gX_true_matched = self._rotate_spin_2_vec(
g1_true, g2_true, angle
)
all_true_gT, all_true_gX = self._rotate_spin_2_vec(
all_true_g1, all_true_g2, all_true_angle
)
ind_gT_true_list.append(gT_true_matched)
dist = self.angsep(det_ra, det_dec)
all_true_dist = self.angsep(all_true_ra, all_true_dec)
ind_dist_list.append(dist)
dist_sort_idx = np.argsort(all_true_dist)
dist_sorted = all_true_dist[dist_sort_idx]
gT_sorted = all_true_gT[dist_sort_idx]
gX_sorted = all_true_gX[dist_sort_idx]
gT_interp = interp1d(
dist_sorted, gT_sorted, bounds_error=False, fill_value=np.nan
)
gX_interp = interp1d(
dist_sorted, gX_sorted, bounds_error=False, fill_value=np.nan
)
gT_true = gT_interp(dist)
gX_true = gX_interp(dist)
r11, r22 = self._get_response_from_w_and_der(
e1,
e2,
w,
e1_g1,
e2_g2,
w_g1,
w_g2,
)
rT, rX = self._rotate_spin_2_matrix(r11, r22, angle)
ind_rT_list.append(rT)
(
eT_list,
eX_list,
rT_list,
rX_list,
gT_true_matched_list,
gX_true_matched_list,
gT_true_list,
gX_true_list,
all_true_gT_list,
all_true_gX_list,
kappa_true_list,
lensed_shift_list,
radial_lensed_shift_list,
r_weighted_gT_matched_list,
r_weighted_gX_matched_list,
r_weighted_gT_list,
r_weighted_gX_list,
ngal_in_bin,
eT_std_list,
eX_std_list,
mean_dist_list,
median_match_dist,
match_failure_rate,
m00_list,
m20_list,
) = self._get_eT_eX_rT_rX_sum(
eT,
eX,
w,
rT,
rX,
gT_true_matched,
gX_true_matched,
gT_true,
gX_true,
kappa_true,
dist,
lensed_shift,
radial_lensed_shift,
angular_bin_edges,
match_dist,
m00,
m20,
all_true_gT,
all_true_gX,
all_true_dist,
)
rT_ensemble[i, :] = rT_list
rX_ensemble[i, :] = rX_list
eT_ensemble[i, :] = eT_list
eX_ensemble[i, :] = eX_list
gT_true_matched_ensemble[i, :] = gT_true_list
gX_true_matched_ensemble[i, :] = gX_true_list
gT_true_ensemble[i, :] = gT_true_list
gX_true_ensemble[i, :] = gX_true_list
all_true_gT_ensemble[i, :] = all_true_gT_list
all_true_gX_ensemble[i, :] = all_true_gX_list
r_weighted_gT_matched_ensemble[i, :] = r_weighted_gT_matched_list
r_weighted_gX_matched_ensemble[i, :] = r_weighted_gX_matched_list
r_weighted_gT_ensemble[i, :] = r_weighted_gT_list
r_weighted_gX_ensemble[i, :] = r_weighted_gX_list
kappa_true_ensemble[i, :] = kappa_true_list
lensed_shift_ensemble[i, :] = lensed_shift_list
radial_lensed_shift_ensemble[i, :] = radial_lensed_shift_list
ngal_in_bin_ensemble[i, :] = ngal_in_bin
eT_std_ensemble[i, :] = eT_std_list / np.sqrt(ngal_in_bin)
eX_std_ensemble[i, :] = eX_std_list / np.sqrt(ngal_in_bin)
mean_dist_ensemble[i, :] = mean_dist_list
median_match_dist_ensemble[i, :] = median_match_dist
match_failure_rate_ensemble[i, :] = match_failure_rate
m00_ensemble[i, :] = m00_list
m20_ensemble[i, :] = m20_list
rT_ind = np.concatenate(ind_rT_list)
ind_gT_true = np.concatenate(ind_gT_true_list)
x_ind = np.concatenate(ind_x_list)
y_ind = np.concatenate(ind_y_list)
lensed_x_ind = np.concatenate(ind_lensed_x_list)
lensed_y_ind = np.concatenate(ind_lensed_y_list)
dist_ind = np.concatenate(ind_dist_list)
per_galaxy_struct = self.get_per_galaxy_struct(len(rT_ind))
per_galaxy_struct["rT"] = rT_ind
per_galaxy_struct["gT"] = ind_gT_true
per_galaxy_struct["x"] = x_ind
per_galaxy_struct["y"] = y_ind
per_galaxy_struct["lensed_x"] = lensed_x_ind
per_galaxy_struct["lensed_y"] = lensed_y_ind
per_galaxy_struct["dist"] = dist_ind
summary_stats = self.get_summary_struct(
n_realization,
len(angular_bin_edges) - 1,
)
# Populate the structured array directly with the ensemble variables
summary_stats["angular_bin_left"] = np.tile(
angular_bin_edges[:-1],
(n_realization, 1),
)
summary_stats["angular_bin_right"] = np.tile(
angular_bin_edges[1:],
(n_realization, 1),
)
summary_stats["ngal_in_bin"] = (
ngal_in_bin_ensemble # Shape (n_halos, n_bins)
)
summary_stats["eT"] = eT_ensemble # Shape (n_halos, n_bins)
summary_stats["eT_std"] = eT_std_ensemble # Shape (n_halos, n_bins)
summary_stats["eX"] = eX_ensemble # Shape (n_halos, n_bins)
summary_stats["eX_std"] = eX_std_ensemble # Shape (n_halos, n_bins)
summary_stats["rT"] = rT_ensemble # Shape (n_halos, n_bins)
summary_stats["rX"] = rX_ensemble # Shape (n_halos, n_bins)
summary_stats["gT_true_matched"] = (
gT_true_matched_ensemble # Shape (n_halos, n_bins)
)
summary_stats["gX_true_matched"] = (
gX_true_matched_ensemble # Shape (n_halos, n_bins)
)
summary_stats["gT_true"] = gT_true_ensemble # Shape (n_halos, n_bins)
summary_stats["gX_true"] = gX_true_ensemble # Shape (n_halos, n_bins)
summary_stats["all_true_gT"] = (
all_true_gT_ensemble # Shape (n_halos, n_bins)
)
summary_stats["all_true_gX"] = all_true_gX_ensemble # Shape
summary_stats["kappa_true"] = (
kappa_true_ensemble # Shape (n_halos, n_bins)
)
summary_stats["lensed_shift"] = (
lensed_shift_ensemble # Shape (n_halos, n_bins)
)
summary_stats["radial_lensed_shift"] = (
radial_lensed_shift_ensemble # Shape (n_halos, n_bins)
)
summary_stats["r_weighted_gT_matched"] = (
r_weighted_gT_matched_ensemble # Shape (n_halos, n_bins)
)
summary_stats["r_weighted_gX_matched"] = (
r_weighted_gX_matched_ensemble # Shape (n_halos, n_bins)
)
summary_stats["r_weighted_gT"] = (
r_weighted_gT_ensemble # Shape (n_halos, n_bins)
)
summary_stats["r_weighted_gX"] = (
r_weighted_gX_ensemble # Shape (n_halos, n_bins)
)
summary_stats["radial_lensed_shift"] = (
radial_lensed_shift_ensemble # Shape (n_halos, n_bins)
)
summary_stats["mean_dist"] = mean_dist_ensemble
summary_stats["median_match_dist"] = (
median_match_dist_ensemble # Shape (n_halos, n_bins)
)
summary_stats["match_failure_rate"] = (
match_failure_rate_ensemble # Shape (n_halos, n_bins)
)
summary_stats["mean_m00"] = m00_ensemble
summary_stats["mean_m20"] = m20_ensemble
summary_plot = self.generate_summary_plot(summary_stats)
return Struct(
outputSummary=summary_stats,
summaryPlot=summary_plot,
perGalaxy=per_galaxy_struct,
)