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 typing import TYPE_CHECKING, Any, Mapping, cast

32import numpy as np

33import treecorr # type: ignore[import]

34from deprecated.sphinx import deprecated

35from lsst.meas.algorithms.treecorrUtils import TreecorrConfig as TreecorrConfigNew

36from lsst.pex.config import ChoiceField, ConfigField, Field

38from ...interfaces import KeyedData, KeyedDataAction, Vector

39from .calcMomentSize import CalcMomentSize

40from .ellipticity import CalcE, CalcEDiff

41from .mathActions import FractionalDifference

43if TYPE_CHECKING:

44 from treecorr import GGCorrelation, KKCorrelation

46 from ...interfaces import KeyedDataSchema

48_LOG = logging.getLogger(__name__)

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

52@deprecated(

53 reason=(

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

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

56 ),

57 version="v28.0",

58 category=FutureWarning,

59)

60class TreecorrConfig(TreecorrConfigNew):

61 pass

64class CalcRhoStatistics(KeyedDataAction):

65 r"""Calculate rho statistics.

67 Rho statistics refer to a collection of correlation functions involving

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

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

70 functions.

72 .. _rho_definitions:

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

76 .. math::

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

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

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

81 \right\rangle

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

84 e^*_{PSF}(x)

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

86 \right\rangle

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

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

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

91 \right\rangle

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

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

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

96 \right\rangle

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

99 e^*_{PSF}(x)

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

101 \right\rangle

102

103

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

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

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

107

108 .. math::

109

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

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

112 \right\rangle

113

114

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

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

117

118 References

119 ----------

120

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

122 et. al (2016).

123 The DES Science Verification weak lensing shear catalogues

124 MNRAS, 460, 2245–2281.

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

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

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

128 Mass and galaxy distributions of four massive galaxy clusters from

129 Dark Energy Survey Science Verification data

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

131 https://doi:10.1093/mnras/stv398

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

133 """

134

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

136

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

138

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

140

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

142

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

144

145 colPsfXx = Field[str](

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

147 )

148

149 colPsfYy = Field[str](

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

151 )

152

153 colPsfXy = Field[str](

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

155 )

156

157 ellipticityType = ChoiceField[str](

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

159 optional=False,

160 allowed={

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

162 "shear": (

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

164 ),

165 },

166 default="distortion",

167 )

168

169 sizeType = ChoiceField[str](

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

171 default="trace",

172 allowed={

173 "trace": "trace radius",

174 "determinant": "determinant radius",

175 },

176 )

177

178 treecorr = ConfigField[TreecorrConfigNew](

179 doc="TreeCorr configuration",

180 )

181

182 def setDefaults(self):

183 super().setDefaults()

184 self.treecorr = TreecorrConfigNew()

185 self.treecorr.sep_units = "arcmin"

186 self.treecorr.max_sep = 100.0

187

188 def getInputSchema(self) -> KeyedDataSchema:

189 return (

190 (self.colRa, Vector),

191 (self.colDec, Vector),

192 (self.colXx, Vector),

193 (self.colYy, Vector),

194 (self.colXy, Vector),

195 (self.colPsfXx, Vector),

196 (self.colPsfYy, Vector),

197 (self.colPsfXy, Vector),

198 )

199

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

201 calcEMeas = CalcE(

202 colXx=self.colXx,

203 colYy=self.colYy,

204 colXy=self.colXy,

205 ellipticityType=self.ellipticityType,

206 )

207 calcEpsf = CalcE(

208 colXx=self.colPsfXx,

209 colYy=self.colPsfYy,

210 colXy=self.colPsfXy,

211 ellipticityType=self.ellipticityType,

212 )

213

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

215

216 calcSizeResidual = FractionalDifference(

217 actionA=CalcMomentSize(

218 colXx=self.colPsfXx,

219 colYy=self.colPsfYy,

220 colXy=self.colPsfXy,

221 sizeType=self.sizeType,

222 ),

223 actionB=CalcMomentSize(

224 colXx=self.colXx,

225 colYy=self.colYy,

226 colXy=self.colXy,

227 sizeType=self.sizeType,

228 ),

229 )

230

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

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

233 # ellipticity-type.

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

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

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

237 # This definition of responsitivity is consistent with the definions

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

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

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

241

242 # Call the actions on the data.

243 eMEAS = calcEMeas(data, **kwargs)

244 if self.ellipticityType == "distortion":

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

246 eMEAS /= responsitivity # type: ignore

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

248 eRes = calcEDiff(data, **kwargs)

249 eRes /= responsitivity # type: ignore

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

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

252

253 # Scale the sizeRes by ellipticities

254 e1SizeRes = e1 * sizeRes

255 e2SizeRes = e2 * sizeRes

256

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

258 # Rho statistics.

259 args = {

260 0: (sizeRes, None),

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

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

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

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

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

266 }

267

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

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

270

271 treecorr_config_dict = self.treecorr.toDict()

272

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

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

275 treecorr_config_dict["rng"] = rng

276

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

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

279 for rhoIndex in range(1, 6):

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

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

282 ra,

283 dec,

284 *(args[rhoIndex]),

285 treecorr_config_dict=treecorr_config_dict,

286 )

287

288 _LOG.info("Calculating rho3alt")

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

290 ra,

291 dec,

292 *(args[0]),

293 treecorr_config_dict=treecorr_config_dict,

294 )

295 return cast(KeyedData, rhoStats)

296

297 @classmethod

298 def _corrSpin0(

299 cls,

300 ra: Vector,

301 dec: Vector,

302 k1: Vector,

303 k2: Vector | None = None,

304 raUnits: str = "degrees",

305 decUnits: str = "degrees",

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

307 ) -> KKCorrelation:

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

309

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

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

312

313 Parameters

314 ----------

315 ra : `numpy.array`

316 The right ascension values of entries in the catalog.

317 dec : `numpy.array`

318 The declination values of entries in the catalog.

319 k1 : `numpy.array`

320 The primary scalar field.

321 k2 : `numpy.array`, optional

322 The secondary scalar field.

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

324 raUnits : `str`, optional

325 Unit of the right ascension values. Valid options are

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

327 decUnits : `str`, optional

328 Unit of the declination values. Valid options are

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

330 treecorr_config_dict: `dict`, optional

331 Config dictionary to be passed to `treecorr`

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

333

334 Returns

335 -------

336 xy : `treecorr.KKCorrelation`

337 A `treecorr.KKCorrelation` object containing the correlation

338 function.

339 """

340 _LOG.debug(

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

342 len(ra),

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

344 )

345 xy = treecorr.KKCorrelation(config=treecorr_config_dict)

346 catA = treecorr.Catalog(

347 config=treecorr_config_dict,

348 ra=ra,

349 dec=dec,

350 k=k1,

351 ra_units=raUnits,

352 dec_units=decUnits,

353 logger=_LOG,

354 )

355 if k2 is None:

356 # Calculate the auto-correlation

357 xy.process(catA)

358 else:

359 catB = treecorr.Catalog(

360 config=treecorr_config_dict,

361 ra=ra,

362 dec=dec,

363 k=k2,

364 ra_units=raUnits,

365 dec_units=decUnits,

366 logger=_LOG,

367 patch_centers=catA.patch_centers,

368 )

369 # Calculate the cross-correlation

370 xy.process(catA, catB)

371

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

373 return xy

374

375 @classmethod

376 def _corrSpin2(

377 cls,

378 ra: Vector,

379 dec: Vector,

380 g1a: Vector,

381 g2a: Vector,

382 g1b: Vector | None = None,

383 g2b: Vector | None = None,

384 raUnits: str = "degrees",

385 decUnits: str = "degrees",

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

387 ) -> GGCorrelation:

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

389

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

391 (shear-like) fields.

392

393 Parameters

394 ----------

395 ra : `numpy.array`

396 The right ascension values of entries in the catalog.

397 dec : `numpy.array`

398 The declination values of entries in the catalog.

399 g1a : `numpy.array`

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

401 g2a : `numpy.array`

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

403 g1b : `numpy.array`, optional

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

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

406 g2b : `numpy.array`, optional

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

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

409 raUnits : `str`, optional

410 Unit of the right ascension values. Valid options are

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

412 decUnits : `str`, optional

413 Unit of the declination values. Valid options are

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

415 treecorr_config_dict : `dict`, optional

416 Config dictionary to be passed to `treecorr`

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

418

419 Returns

420 -------

421 xy : `treecorr.GGCorrelation`

422 A `treecorr.GGCorrelation` object containing the correlation

423 function.

424 """

425 _LOG.debug(

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

427 len(ra),

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

429 )

430 xy = treecorr.GGCorrelation(config=treecorr_config_dict)

431 catA = treecorr.Catalog(

432 config=treecorr_config_dict,

433 ra=ra,

434 dec=dec,

435 g1=g1a,

436 g2=g2a,

437 ra_units=raUnits,

438 dec_units=decUnits,

439 logger=_LOG,

440 )

441 if g1b is None or g2b is None:

442 # Calculate the auto-correlation

443 xy.process(catA)

444 else:

445 catB = treecorr.Catalog(

446 config=treecorr_config_dict,

447 ra=ra,

448 dec=dec,

449 g1=g1b,

450 g2=g2b,

451 ra_units=raUnits,

452 dec_units=decUnits,

453 logger=_LOG,

454 patch_centers=catA.patch_centers,

455 )

456 # Calculate the cross-correlation

457 xy.process(catA, catB)

458

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

460 return xy

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

88 statements