Coverage for python/lsst/analysis/tools/actions/vector/calcRhoStatistics.py: 38%

1# This file is part of analysis_tools.

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.

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__)

53# TO DO: Remove TreecorrConfig in here for next major release (DM-47072)

54@deprecated(

55 reason=(

56 "TreecorrConfig is no longer a part of analysis_tools (DM-45899). "

57 "Please use lsst.meas.algorithms.treecorrUtils.TreecorrConfig instead."

58 ),

59 version="v28.0",

60 category=FutureWarning,

61)

62class TreecorrConfig(TreecorrConfigNew):

63 pass

66class CalcRhoStatistics(KeyedDataAction):

67 r"""Calculate rho statistics.

69 Rho statistics refer to a collection of correlation functions involving

70 PSF ellipticity and size residuals. They quantify the contribution from PSF

71 leakage due to errors in PSF modeling to the weak lensing shear correlation

72 functions.

74 .. _rho_definitions:

76 The exact definitions of rho statistics as defined in [1]_ are given below.

78 .. math::

80 \rho_1(\theta) &= \left\langle

81 \delta e^*_{PSF}(x)

82 \delta e_{PSF}(x+\theta)

83 \right\rangle

85 \rho_2(\theta) &= \left\langle

86 e^*_{PSF}(x)

87 \delta e_{PSF}(x+\theta

88 \right\rangle

90 \rho_3(\theta) &= \left\langle

91 (e^*_{PSF}\frac{\delta T_{PSF}}{T_{PSF}}(x))

92 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)

93 \right\rangle

95 \rho_4(\theta) &= \left\langle

96 \delta e^*_{PSF}(x)

97 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)

98 \right\rangle

100 \rho_5(\theta) &= \left\langle

101 e^*_{PSF}(x)

102 (e_{PSF}\frac{\delta T_{PSF}}{T_{PSF}})(x+\theta)

103 \right\rangle

104

105

106 In addition to these five, we also compute the auto-correlation function of

107 the fractional size residuals and call it as the :math:`\rho'_3( \theta )`,

108 as referred to in Melchior et al. (2015) [2]_.

109

110 .. math::

111

112 \rho'_3(\theta) = \left\langle\frac{\delta T_{PSF}}{T_{PSF}}(x)

113 \frac{\delta T_{PSF}}{T_{PSF}}(x+\theta)

114 \right\rangle

115

116

117 The definition of ellipticity used in [1]_ correspond to shear-type,

118 which is typically smaller by a factor of 4 than using distortion-type.

119

120 References

121 ----------

122

123 .. [1] Jarvis, M., Sheldon, E., Zuntz, J., Kacprzak, T., Bridle, S. L.,

124 et. al (2016).

125 The DES Science Verification weak lensing shear catalogues

126 MNRAS, 460, 2245–2281.

127 https://doi.org/10.1093/mnras/stw990;

128 https://arxiv.org/abs/1507.05603

129 .. [2] Melchior, P., et. al (2015)

130 Mass and galaxy distributions of four massive galaxy clusters from

131 Dark Energy Survey Science Verification data

132 MNRAS, 449, no. 3, pp. 2219–2238.

133 https://doi:10.1093/mnras/stv398

134 https://arxiv.org/abs/1405.4285

135 """

136

137 colRa = Field[str](doc="RA column", default="coord_ra")

138

139 colDec = Field[str](doc="Dec column", default="coord_dec")

140

141 colXx = Field[str](doc="The column name to get the xx shape component from.", default="{band}_ixx")

142

143 colYy = Field[str](doc="The column name to get the yy shape component from.", default="{band}_iyy")

144

145 colXy = Field[str](doc="The column name to get the xy shape component from.", default="{band}_ixy")

146

147 colPsfXx = Field[str](

148 doc="The column name to get the PSF xx shape component from.", default="{band}_ixxPSF"

149 )

150

151 colPsfYy = Field[str](

152 doc="The column name to get the PSF yy shape component from.", default="{band}_iyyPSF"

153 )

154

155 colPsfXy = Field[str](

156 doc="The column name to get the PSF xy shape component from.", default="{band}_ixyPSF"

157 )

158

159 ellipticityType = ChoiceField[str](

160 doc="The type of ellipticity to calculate",

161 optional=False,

162 allowed={

163 "distortion": r"Distortion, measured as :math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy})`",

164 "shear": (

165 r"Shear, measured as :math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^2})`"

166 ),

167 },

168 default="distortion",

169 )

170

171 sizeType = ChoiceField[str](

172 doc="The type of size to calculate",

173 default="trace",

174 allowed={

175 "trace": "trace radius",

176 "determinant": "determinant radius",

177 },

178 )

179

180 treecorr = ConfigField[TreecorrConfigNew](

181 doc="TreeCorr configuration",

182 )

183

184 def setDefaults(self):

185 super().setDefaults()

186 self.treecorr = TreecorrConfigNew()

187 self.treecorr.sep_units = "arcmin"

188 self.treecorr.max_sep = 100.0

189

190 def getInputSchema(self) -> KeyedDataSchema:

191 return (

192 (self.colRa, Vector),

193 (self.colDec, Vector),

194 (self.colXx, Vector),

195 (self.colYy, Vector),

196 (self.colXy, Vector),

197 (self.colPsfXx, Vector),

198 (self.colPsfYy, Vector),

199 (self.colPsfXy, Vector),

200 )

201

202 def __call__(self, data: KeyedData, **kwargs) -> KeyedData:

203 calcEMeas = CalcE(

204 colXx=self.colXx,

205 colYy=self.colYy,

206 colXy=self.colXy,

207 ellipticityType=self.ellipticityType,

208 )

209 calcEpsf = CalcE(

210 colXx=self.colPsfXx,

211 colYy=self.colPsfYy,

212 colXy=self.colPsfXy,

213 ellipticityType=self.ellipticityType,

214 )

215

216 calcEDiff = CalcEDiff(colA=calcEMeas, colB=calcEpsf)

217

218 calcSizeResidual = FractionalDifference(

219 actionA=CalcMomentSize(

220 colXx=self.colPsfXx,

221 colYy=self.colPsfYy,

222 colXy=self.colPsfXy,

223 sizeType=self.sizeType,

224 ),

225 actionB=CalcMomentSize(

226 colXx=self.colXx,

227 colYy=self.colYy,

228 colXy=self.colXy,

229 sizeType=self.sizeType,

230 ),

231 )

232

233 # distortion-type ellipticity has a shear response of 2, so we need to

234 # divide by 2 so that the rho-stats do not depend on the

235 # ellipticity-type.

236 # Note: For distortion, the responsitivity is 2(1 - e^2_{rms}),

237 # where e_rms is the root mean square ellipticity per component.

238 # This is expected to be small and we ignore it here.

239 # This definition of responsitivity is consistent with the definions

240 # used in the rho-statistics calculations for the HSC shear catalog

241 # papers (Mandelbaum et al. 2018, Li et al., 2022).

242 responsitivity = 2.0 if self.ellipticityType == "distortion" else 1.0

243

244 # Call the actions on the data.

245 eMEAS = calcEMeas(data, **kwargs)

246 if self.ellipticityType == "distortion":

247 _LOG.debug("Correction value of responsitivity would be %f", 2 - np.mean(np.abs(eMEAS) ** 2))

248 eMEAS /= responsitivity # type: ignore

249 e1, e2 = np.real(eMEAS), np.imag(eMEAS)

250 eRes = calcEDiff(data, **kwargs)

251 eRes /= responsitivity # type: ignore

252 e1Res, e2Res = np.real(eRes), np.imag(eRes)

253 sizeRes = -1 * calcSizeResidual(data, **kwargs) # sign flip residual to T_psf - T_model

254

255 # Scale the sizeRes by ellipticities

256 e1SizeRes = e1 * sizeRes

257 e2SizeRes = e2 * sizeRes

258

259 # Package the arguments to capture auto-/cross-correlations for the

260 # Rho statistics.

261 args = {

262 0: (sizeRes, None),

263 1: (e1Res, e2Res, None, None),

264 2: (e1, e2, e1Res, e2Res),

265 3: (e1SizeRes, e2SizeRes, None, None),

266 4: (e1Res, e2Res, e1SizeRes, e2SizeRes),

267 5: (e1, e2, e1SizeRes, e2SizeRes),

268 }

269

270 ra: Vector = data[self.colRa] # type: ignore

271 dec: Vector = data[self.colDec] # type: ignore

272

273 treecorr_config_dict = self.treecorr.toDict()

274

275 # Swap rng_seed with an rng instance in treecorr config.

276 rng = np.random.RandomState(treecorr_config_dict.pop("rng_seed"))

277 treecorr_config_dict["rng"] = rng

278

279 # Pass the appropriate arguments to the correlator and build a dict

280 rhoStats: Mapping[str, treecorr.BinnedCorr2] = {}

281 for rhoIndex in range(1, 6):

282 _LOG.info("Calculating rho-%d", rhoIndex)

283 rhoStats[f"rho{rhoIndex}"] = self._corrSpin2( # type: ignore[index]

284 ra,

285 dec,

286 *(args[rhoIndex]),

287 treecorr_config_dict=treecorr_config_dict,

288 )

289

290 _LOG.info("Calculating rho3alt")

291 rhoStats["rho3alt"] = self._corrSpin0( # type: ignore[index]

292 ra,

293 dec,

294 *(args[0]),

295 treecorr_config_dict=treecorr_config_dict,

296 )

297 return cast(KeyedData, rhoStats)

298

299 @classmethod

300 def _corrSpin0(

301 cls,

302 ra: Vector,

303 dec: Vector,

304 k1: Vector,

305 k2: Vector | None = None,

306 raUnits: str = "degrees",

307 decUnits: str = "degrees",

308 treecorr_config_dict: Mapping[str, Any] | None = None,

309 ) -> KKCorrelation:

310 """Function to compute correlations between at most two scalar fields.

311

312 This is used to compute rho3alt statistics, given the appropriate

313 spin-0 (scalar) fields, usually fractional size residuals.

314

315 Parameters

316 ----------

317 ra : `numpy.array`

318 The right ascension values of entries in the catalog.

319 dec : `numpy.array`

320 The declination values of entries in the catalog.

321 k1 : `numpy.array`

322 The primary scalar field.

323 k2 : `numpy.array`, optional

324 The secondary scalar field.

325 Autocorrelation of the primary field is computed if `None`.

326 raUnits : `str`, optional

327 Unit of the right ascension values. Valid options are

328 "degrees", "arcmin", "arcsec", "hours" or "radians".

329 decUnits : `str`, optional

330 Unit of the declination values. Valid options are

331 "degrees", "arcmin", "arcsec", "hours" or "radians".

332 treecorr_config_dict: `dict`, optional

333 Config dictionary to be passed to `treecorr`

334 (`treecorr.KKCorrelation` or `treecorr.Catalog`).

335

336 Returns

337 -------

338 xy : `treecorr.KKCorrelation`

339 A `treecorr.KKCorrelation` object containing the correlation

340 function.

341 """

342 _LOG.debug(

343 "No. of entries: %d. The number of pairs in the resulting KKCorrelation cannot exceed %d",

344 len(ra),

345 len(ra) * (len(ra) - 1) / 2,

346 )

347 xy = treecorr.KKCorrelation(config=treecorr_config_dict)

348 catA = treecorr.Catalog(

349 config=treecorr_config_dict,

350 ra=ra,

351 dec=dec,

352 k=k1,

353 ra_units=raUnits,

354 dec_units=decUnits,

355 logger=_LOG,

356 )

357 if k2 is None:

358 # Calculate the auto-correlation

359 xy.process(catA)

360 else:

361 catB = treecorr.Catalog(

362 config=treecorr_config_dict,

363 ra=ra,

364 dec=dec,

365 k=k2,

366 ra_units=raUnits,

367 dec_units=decUnits,

368 logger=_LOG,

369 patch_centers=catA.patch_centers,

370 )

371 # Calculate the cross-correlation

372 xy.process(catA, catB)

373

374 _LOG.debug("Correlated %d pairs based on the config set.", sum(xy.npairs))

375 return xy

376

377 @classmethod

378 def _corrSpin2(

379 cls,

380 ra: Vector,

381 dec: Vector,

382 g1a: Vector,

383 g2a: Vector,

384 g1b: Vector | None = None,

385 g2b: Vector | None = None,

386 raUnits: str = "degrees",

387 decUnits: str = "degrees",

388 treecorr_config_dict: Mapping[str, Any] | None = None,

389 ) -> GGCorrelation:

390 """Function to compute correlations between shear-like fields.

391

392 This is used to compute Rho statistics, given the appropriate spin-2

393 (shear-like) fields.

394

395 Parameters

396 ----------

397 ra : `numpy.array`

398 The right ascension values of entries in the catalog.

399 dec : `numpy.array`

400 The declination values of entries in the catalog.

401 g1a : `numpy.array`

402 The first component of the primary shear-like field.

403 g2a : `numpy.array`

404 The second component of the primary shear-like field.

405 g1b : `numpy.array`, optional

406 The first component of the secondary shear-like field.

407 Autocorrelation of the primary field is computed if `None`.

408 g2b : `numpy.array`, optional

409 The second component of the secondary shear-like field.

410 Autocorrelation of the primary field is computed if `None`.

411 raUnits : `str`, optional

412 Unit of the right ascension values. Valid options are

413 "degrees", "arcmin", "arcsec", "hours" or "radians".

414 decUnits : `str`, optional

415 Unit of the declination values. Valid options are

416 "degrees", "arcmin", "arcsec", "hours" or "radians".

417 treecorr_config_dict : `dict`, optional

418 Config dictionary to be passed to `treecorr`

419 (`treecorr.GGCorrelation` or `treecorr.Catalog`).

420

421 Returns

422 -------

423 xy : `treecorr.GGCorrelation`

424 A `treecorr.GGCorrelation` object containing the correlation

425 function.

426 """

427 _LOG.debug(

428 "No. of entries: %d. The number of pairs in the resulting GGCorrelation cannot exceed %d",

429 len(ra),

430 len(ra) * (len(ra) - 1) / 2,

431 )

432 xy = treecorr.GGCorrelation(config=treecorr_config_dict)

433 catA = treecorr.Catalog(

434 config=treecorr_config_dict,

435 ra=ra,

436 dec=dec,

437 g1=g1a,

438 g2=g2a,

439 ra_units=raUnits,

440 dec_units=decUnits,

441 logger=_LOG,

442 )

443 if g1b is None or g2b is None:

444 # Calculate the auto-correlation

445 xy.process(catA)

446 else:

447 catB = treecorr.Catalog(

448 config=treecorr_config_dict,

449 ra=ra,

450 dec=dec,

451 g1=g1b,

452 g2=g2b,

453 ra_units=raUnits,

454 dec_units=decUnits,

455 logger=_LOG,

456 patch_centers=catA.patch_centers,

457 )

458 # Calculate the cross-correlation

459 xy.process(catA, catB)

460

461 _LOG.debug("Correlated %d pairs based on the config set.", sum(xy.npairs))

462 return xy

Coverage for python / lsst / analysis / tools / actions / vector / calcRhoStatistics.py: 38%

89 statements