Coverage for python/lsst/analysis/tools/actions/vector/calcRhoStatistics.py: 95%
90 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-02 02:22 -0700
1# This file is part of analysis_tools.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = (
25 "CalcRhoStatistics",
26 "TreecorrConfig",
27)
29import logging
30from collections.abc import Mapping
31from typing import TYPE_CHECKING, Any, cast
33import numpy as np
34import treecorr # type: ignore[import]
35from deprecated.sphinx import deprecated
37from lsst.meas.algorithms.treecorrUtils import TreecorrConfig as TreecorrConfigNew
38from lsst.pex.config import ChoiceField, ConfigField, Field
40from ...interfaces import KeyedData, KeyedDataAction, Vector
41from .calcMomentSize import CalcMomentSize
42from .ellipticity import CalcE, CalcEDiff
43from .mathActions import FractionalDifference
45if TYPE_CHECKING:
46 from treecorr import GGCorrelation, KKCorrelation
48 from ...interfaces import KeyedDataSchema
50_LOG = logging.getLogger(__name__)
52# We need to explicitly turn off multiprocessing in treecorr.
53treecorr.set_max_omp_threads(1)
56# TO DO: Remove TreecorrConfig in here for next major release (DM-47072)
57@deprecated(
58 reason=(
59 "TreecorrConfig is no longer a part of analysis_tools (DM-45899). "
60 "Please use lsst.meas.algorithms.treecorrUtils.TreecorrConfig instead."
61 ),
62 version="v28.0",
63 category=FutureWarning,
64)
65class TreecorrConfig(TreecorrConfigNew):
66 pass
69class CalcRhoStatistics(KeyedDataAction):
70 r"""Calculate rho statistics.
72 Rho statistics refer to a collection of correlation functions involving
73 PSF ellipticity and size residuals. They quantify the contribution from PSF
74 leakage due to errors in PSF modeling to the weak lensing shear correlation
75 functions.
77 .. _rho_definitions:
79 The exact definitions of rho statistics as defined in [1]_ are given below.
81 .. math::
83 \rho_1(\theta) &= \left\langle
84 \delta e^*_{PSF}(x)
85 \delta e_{PSF}(x+\theta)
86 \right\rangle
88 \rho_2(\theta) &= \left\langle
89 e^*_{PSF}(x)
90 \delta e_{PSF}(x+\theta
91 \right\rangle
93 \rho_3(\theta) &= \left\langle
94 (e^*_{PSF}\frac{\delta T_{PSF}}{T_{PSF}}(x))
95 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)
96 \right\rangle
98 \rho_4(\theta) &= \left\langle
99 \delta e^*_{PSF}(x)
100 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)
101 \right\rangle
103 \rho_5(\theta) &= \left\langle
104 e^*_{PSF}(x)
105 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)
106 \right\rangle
109 In addition to these five, we also compute the auto-correlation function of
110 the fractional size residuals and call it as the :math:`\rho'_3( \theta )`,
111 as referred to in Melchior et al. (2015) [2]_.
113 .. math::
115 \rho'_3(\theta) = \left\langle\frac{\delta T_{PSF}}{T_{PSF}}(x)
116 \frac{\delta T_{PSF}}{T_{PSF}}(x+\theta)
117 \right\rangle
120 The definition of ellipticity used in [1]_ correspond to shear-type,
121 which is typically smaller by a factor of 4 than using distortion-type.
123 References
124 ----------
126 .. [1] Jarvis, M., Sheldon, E., Zuntz, J., Kacprzak, T., Bridle, S. L.,
127 et. al (2016).
128 The DES Science Verification weak lensing shear catalogues
129 MNRAS, 460, 2245–2281.
130 https://doi.org/10.1093/mnras/stw990;
131 https://arxiv.org/abs/1507.05603
132 .. [2] Melchior, P., et. al (2015)
133 Mass and galaxy distributions of four massive galaxy clusters from
134 Dark Energy Survey Science Verification data
135 MNRAS, 449, no. 3, pp. 2219–2238.
136 https://doi:10.1093/mnras/stv398
137 https://arxiv.org/abs/1405.4285
138 """
140 colRa = Field[str](doc="RA column", default="coord_ra")
142 colDec = Field[str](doc="Dec column", default="coord_dec")
144 colXx = Field[str](doc="The column name to get the xx shape component from.", default="{band}_ixx")
146 colYy = Field[str](doc="The column name to get the yy shape component from.", default="{band}_iyy")
148 colXy = Field[str](doc="The column name to get the xy shape component from.", default="{band}_ixy")
150 colPsfXx = Field[str](
151 doc="The column name to get the PSF xx shape component from.", default="{band}_ixxPSF"
152 )
154 colPsfYy = Field[str](
155 doc="The column name to get the PSF yy shape component from.", default="{band}_iyyPSF"
156 )
158 colPsfXy = Field[str](
159 doc="The column name to get the PSF xy shape component from.", default="{band}_ixyPSF"
160 )
162 ellipticityType = ChoiceField[str](
163 doc="The type of ellipticity to calculate",
164 optional=False,
165 allowed={
166 "distortion": r"Distortion, measured as :math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy})`",
167 "shear": (
168 r"Shear, measured as :math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^2})`"
169 ),
170 },
171 default="distortion",
172 )
174 sizeType = ChoiceField[str](
175 doc="The type of size to calculate",
176 default="trace",
177 allowed={
178 "trace": "trace radius",
179 "determinant": "determinant radius",
180 },
181 )
183 treecorr = ConfigField[TreecorrConfigNew](
184 doc="TreeCorr configuration",
185 )
187 def setDefaults(self):
188 super().setDefaults()
189 self.treecorr = TreecorrConfigNew()
190 self.treecorr.sep_units = "arcmin"
191 self.treecorr.max_sep = 100.0
193 def getInputSchema(self) -> KeyedDataSchema:
194 return (
195 (self.colRa, Vector),
196 (self.colDec, Vector),
197 (self.colXx, Vector),
198 (self.colYy, Vector),
199 (self.colXy, Vector),
200 (self.colPsfXx, Vector),
201 (self.colPsfYy, Vector),
202 (self.colPsfXy, Vector),
203 )
205 def __call__(self, data: KeyedData, **kwargs) -> KeyedData:
206 calcEMeas = CalcE(
207 colXx=self.colXx,
208 colYy=self.colYy,
209 colXy=self.colXy,
210 ellipticityType=self.ellipticityType,
211 )
212 calcEpsf = CalcE(
213 colXx=self.colPsfXx,
214 colYy=self.colPsfYy,
215 colXy=self.colPsfXy,
216 ellipticityType=self.ellipticityType,
217 )
219 calcEDiff = CalcEDiff(colA=calcEMeas, colB=calcEpsf)
221 calcSizeResidual = FractionalDifference(
222 actionA=CalcMomentSize(
223 colXx=self.colPsfXx,
224 colYy=self.colPsfYy,
225 colXy=self.colPsfXy,
226 sizeType=self.sizeType,
227 ),
228 actionB=CalcMomentSize(
229 colXx=self.colXx,
230 colYy=self.colYy,
231 colXy=self.colXy,
232 sizeType=self.sizeType,
233 ),
234 )
236 # distortion-type ellipticity has a shear response of 2, so we need to
237 # divide by 2 so that the rho-stats do not depend on the
238 # ellipticity-type.
239 # Note: For distortion, the responsitivity is 2(1 - e^2_{rms}),
240 # where e_rms is the root mean square ellipticity per component.
241 # This is expected to be small and we ignore it here.
242 # This definition of responsitivity is consistent with the definions
243 # used in the rho-statistics calculations for the HSC shear catalog
244 # papers (Mandelbaum et al. 2018, Li et al., 2022).
245 responsitivity = 2.0 if self.ellipticityType == "distortion" else 1.0
247 # Call the actions on the data.
248 eMEAS = calcEMeas(data, **kwargs)
249 if self.ellipticityType == "distortion": 249 ↛ 251line 249 didn't jump to line 251 because the condition on line 249 was always true
250 _LOG.debug("Correction value of responsitivity would be %f", 2 - np.mean(np.abs(eMEAS) ** 2))
251 eMEAS /= responsitivity # type: ignore
252 e1, e2 = np.real(eMEAS), np.imag(eMEAS)
253 eRes = calcEDiff(data, **kwargs)
254 eRes /= responsitivity # type: ignore
255 e1Res, e2Res = np.real(eRes), np.imag(eRes)
256 sizeRes = -1 * calcSizeResidual(data, **kwargs) # sign flip residual to T_psf - T_model
258 # Scale the sizeRes by ellipticities
259 e1SizeRes = e1 * sizeRes
260 e2SizeRes = e2 * sizeRes
262 # Package the arguments to capture auto-/cross-correlations for the
263 # Rho statistics.
264 args = {
265 0: (sizeRes, None),
266 1: (e1Res, e2Res, None, None),
267 2: (e1, e2, e1Res, e2Res),
268 3: (e1SizeRes, e2SizeRes, None, None),
269 4: (e1Res, e2Res, e1SizeRes, e2SizeRes),
270 5: (e1, e2, e1SizeRes, e2SizeRes),
271 }
273 ra: Vector = data[self.colRa] # type: ignore
274 dec: Vector = data[self.colDec] # type: ignore
276 treecorr_config_dict = self.treecorr.toDict()
278 # Swap rng_seed with an rng instance in treecorr config.
279 rng = np.random.RandomState(treecorr_config_dict.pop("rng_seed"))
280 treecorr_config_dict["rng"] = rng
282 # Pass the appropriate arguments to the correlator and build a dict
283 rhoStats: Mapping[str, treecorr.BinnedCorr2] = {}
284 for rhoIndex in range(1, 6):
285 _LOG.info("Calculating rho-%d", rhoIndex)
286 rhoStats[f"rho{rhoIndex}"] = self._corrSpin2( # type: ignore[index]
287 ra,
288 dec,
289 *(args[rhoIndex]),
290 treecorr_config_dict=treecorr_config_dict,
291 )
293 _LOG.info("Calculating rho3alt")
294 rhoStats["rho3alt"] = self._corrSpin0( # type: ignore[index]
295 ra,
296 dec,
297 *(args[0]),
298 treecorr_config_dict=treecorr_config_dict,
299 )
300 return cast(KeyedData, rhoStats)
302 @classmethod
303 def _corrSpin0(
304 cls,
305 ra: Vector,
306 dec: Vector,
307 k1: Vector,
308 k2: Vector | None = None,
309 raUnits: str = "degrees",
310 decUnits: str = "degrees",
311 treecorr_config_dict: Mapping[str, Any] | None = None,
312 ) -> KKCorrelation:
313 """Function to compute correlations between at most two scalar fields.
315 This is used to compute rho3alt statistics, given the appropriate
316 spin-0 (scalar) fields, usually fractional size residuals.
318 Parameters
319 ----------
320 ra : `numpy.array`
321 The right ascension values of entries in the catalog.
322 dec : `numpy.array`
323 The declination values of entries in the catalog.
324 k1 : `numpy.array`
325 The primary scalar field.
326 k2 : `numpy.array`, optional
327 The secondary scalar field.
328 Autocorrelation of the primary field is computed if `None`.
329 raUnits : `str`, optional
330 Unit of the right ascension values. Valid options are
331 "degrees", "arcmin", "arcsec", "hours" or "radians".
332 decUnits : `str`, optional
333 Unit of the declination values. Valid options are
334 "degrees", "arcmin", "arcsec", "hours" or "radians".
335 treecorr_config_dict: `dict`, optional
336 Config dictionary to be passed to `treecorr`
337 (`treecorr.KKCorrelation` or `treecorr.Catalog`).
339 Returns
340 -------
341 xy : `treecorr.KKCorrelation`
342 A `treecorr.KKCorrelation` object containing the correlation
343 function.
344 """
345 _LOG.debug(
346 "No. of entries: %d. The number of pairs in the resulting KKCorrelation cannot exceed %d",
347 len(ra),
348 len(ra) * (len(ra) - 1) / 2,
349 )
350 xy = treecorr.KKCorrelation(config=treecorr_config_dict)
351 catA = treecorr.Catalog(
352 config=treecorr_config_dict,
353 ra=ra,
354 dec=dec,
355 k=k1,
356 ra_units=raUnits,
357 dec_units=decUnits,
358 logger=_LOG,
359 )
360 if k2 is None: 360 ↛ 364line 360 didn't jump to line 364 because the condition on line 360 was always true
361 # Calculate the auto-correlation
362 xy.process(catA)
363 else:
364 catB = treecorr.Catalog(
365 config=treecorr_config_dict,
366 ra=ra,
367 dec=dec,
368 k=k2,
369 ra_units=raUnits,
370 dec_units=decUnits,
371 logger=_LOG,
372 patch_centers=catA.patch_centers,
373 )
374 # Calculate the cross-correlation
375 xy.process(catA, catB)
377 _LOG.debug("Correlated %d pairs based on the config set.", sum(xy.npairs))
378 return xy
380 @classmethod
381 def _corrSpin2(
382 cls,
383 ra: Vector,
384 dec: Vector,
385 g1a: Vector,
386 g2a: Vector,
387 g1b: Vector | None = None,
388 g2b: Vector | None = None,
389 raUnits: str = "degrees",
390 decUnits: str = "degrees",
391 treecorr_config_dict: Mapping[str, Any] | None = None,
392 ) -> GGCorrelation:
393 """Function to compute correlations between shear-like fields.
395 This is used to compute Rho statistics, given the appropriate spin-2
396 (shear-like) fields.
398 Parameters
399 ----------
400 ra : `numpy.array`
401 The right ascension values of entries in the catalog.
402 dec : `numpy.array`
403 The declination values of entries in the catalog.
404 g1a : `numpy.array`
405 The first component of the primary shear-like field.
406 g2a : `numpy.array`
407 The second component of the primary shear-like field.
408 g1b : `numpy.array`, optional
409 The first component of the secondary shear-like field.
410 Autocorrelation of the primary field is computed if `None`.
411 g2b : `numpy.array`, optional
412 The second component of the secondary shear-like field.
413 Autocorrelation of the primary field is computed if `None`.
414 raUnits : `str`, optional
415 Unit of the right ascension values. Valid options are
416 "degrees", "arcmin", "arcsec", "hours" or "radians".
417 decUnits : `str`, optional
418 Unit of the declination values. Valid options are
419 "degrees", "arcmin", "arcsec", "hours" or "radians".
420 treecorr_config_dict : `dict`, optional
421 Config dictionary to be passed to `treecorr`
422 (`treecorr.GGCorrelation` or `treecorr.Catalog`).
424 Returns
425 -------
426 xy : `treecorr.GGCorrelation`
427 A `treecorr.GGCorrelation` object containing the correlation
428 function.
429 """
430 _LOG.debug(
431 "No. of entries: %d. The number of pairs in the resulting GGCorrelation cannot exceed %d",
432 len(ra),
433 len(ra) * (len(ra) - 1) / 2,
434 )
435 xy = treecorr.GGCorrelation(config=treecorr_config_dict)
436 catA = treecorr.Catalog(
437 config=treecorr_config_dict,
438 ra=ra,
439 dec=dec,
440 g1=g1a,
441 g2=g2a,
442 ra_units=raUnits,
443 dec_units=decUnits,
444 logger=_LOG,
445 )
446 if g1b is None or g2b is None:
447 # Calculate the auto-correlation
448 xy.process(catA)
449 else:
450 catB = treecorr.Catalog(
451 config=treecorr_config_dict,
452 ra=ra,
453 dec=dec,
454 g1=g1b,
455 g2=g2b,
456 ra_units=raUnits,
457 dec_units=decUnits,
458 logger=_LOG,
459 patch_centers=catA.patch_centers,
460 )
461 # Calculate the cross-correlation
462 xy.process(catA, catB)
464 _LOG.debug("Correlated %d pairs based on the config set.", sum(xy.npairs))
465 return xy