Coverage for python / lsst / meas / base / diaCalculationPlugins.py: 42%

1# This file is part of ap_association. 

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/>. 

21 

22"""Plugins for use in DiaSource summary statistics. 

23 

24Output columns must be 

25as defined in the schema of the Apdb both in name and units. 

26""" 

27 

28import functools 

29import warnings 

30 

31from astropy.stats import median_absolute_deviation 

32import numpy as np 

33import pandas as pd 

34from scipy.optimize import lsq_linear 

35 

36import lsst.geom as geom 

37import lsst.pex.config as pexConfig 

38import lsst.pipe.base as pipeBase 

39import lsst.sphgeom as sphgeom 

40from astropy.timeseries import LombScargle 

41from astropy.timeseries import LombScargleMultiband 

42import math 

43import statistics 

44 

45from .diaCalculation import ( 

46 DiaObjectCalculationPluginConfig, 

47 DiaObjectCalculationPlugin) 

48from .pluginRegistry import register 

49 

50 

51__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition", 

52 "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig", 

53 "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig", 

54 "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig", 

55 "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux", 

56 "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig", 

57 "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig", 

58 "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig", 

59 "MadDiaPsfFlux", "MadDiaPsfFluxConfig", 

60 "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig", 

61 "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig", 

62 "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig", 

63 "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig", 

64 "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig", 

65 "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig", 

66 "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig", 

67 "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig", 

68 "LombScarglePeriodogram", "LombScarglePeriodogramConfig", 

69 "LombScarglePeriodogramMulti", "LombScarglePeriodogramMultiConfig", 

70 "UnphysicalDiaSourceSeparation") 

71 

72 

73def catchWarnings(_func=None, *, warns=[]): 

74 """Decorator for generically catching numpy warnings. 

75 """ 

76 def decoratorCatchWarnings(func): 

77 @functools.wraps(func) 

78 def wrapperCatchWarnings(*args, **kwargs): 

79 with warnings.catch_warnings(): 

80 for val in warns: 

81 warnings.filterwarnings("ignore", val) 

82 return func(*args, **kwargs) 

83 return wrapperCatchWarnings 

84 

85 if _func is None: 85 ↛ 88line 85 didn't jump to line 88 because the condition on line 85 was always true

86 return decoratorCatchWarnings 

87 else: 

88 return decoratorCatchWarnings(_func) 

89 

90 

91def typeSafePandasAssignment( 

92 target, 

93 source, 

94 columns, 

95 default_dtype=np.float64, 

96 int_fill_value=0, 

97 # TODO DM-53254: Remove the force_int_to_float hack. 

98 force_int_to_float=False, 

99): 

100 """ 

101 Assign from a source dataframe to a target dataframe in a type safe way. 

102 

103 Parameters 

104 ---------- 

105 target : `pd.DataFrame` 

106 Target pandas dataframe. 

107 source : `pd.DataFrame` or `pd.Series` 

108 Grouped source dataframe. 

109 columns : `list` [`str`] 

110 List of columns to transfer. 

111 default_dtype : `np.dtype`, optional 

112 Default datatype (if not in target). 

113 int_fill_value : `int`, optional 

114 Fill value for integer columns to avoid pandas insisting 

115 that everything should be float-ified as nans. 

116 force_int_to_float : `bool`, optional 

117 Force integer columns to float columns? Use this option 

118 for backwards compatibility for old pandas misfeatures which 

119 are expected by some downstream processes. 

120 """ 

121 is_series = isinstance(source, pd.Series) 

122 for col in columns: 

123 if is_series: 

124 source_col = source 

125 else: 

126 source_col = source[col] 

127 

128 matched_length = False 

129 if col in target.columns: 

130 target_dtype = target[col].dtype 

131 matched_length = len(target) == len(source) 

132 else: 

133 target_dtype = default_dtype 

134 

135 if (matched_length or pd.api.types.is_float_dtype(target_dtype)) and not force_int_to_float: 

136 # If we have a matched length or float, we can do a 

137 # straight assignment here. 

138 target.loc[:, col] = source_col.astype(target_dtype) 

139 else: 

140 # With mis-matched integers, we must do this dance to preserve types. 

141 # Note that this may lose precision with very large numbers. 

142 

143 # Convert to float 

144 target[col] = target[col].astype(np.float64) 

145 # Set the column, casting to float. 

146 target.loc[:, col] = source_col.astype(np.float64) 

147 if not force_int_to_float: 

148 # Convert back to integer 

149 target[col] = target[col].fillna(int_fill_value).astype(target_dtype) 

150 

151 

152def compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100): 

153 """ 

154 Computes an optimized periodogram frequency grid for a given time series. 

155 

156 Parameters 

157 ---------- 

158 x0 : `array` 

159 The input time axis. 

160 oversampling_factor : `int`, optional 

161 The oversampling factor for frequency grid. 

162 nyquist_factor : `int`, optional 

163 The Nyquist factor for frequency grid. 

164 

165 Returns 

166 ------- 

167 frequencies : `array` 

168 The computed optimized periodogram frequency grid. 

169 """ 

170 

171 num_points = len(x0) 

172 baseline = np.max(x0) - np.min(x0) 

173 

174 # Calculate the frequency resolution based on oversampling factor and baseline 

175 frequency_resolution = 1. / baseline / oversampling_factor 

176 

177 num_frequencies = int( 

178 0.5 * oversampling_factor * nyquist_factor * num_points) 

179 frequencies = frequency_resolution + \ 

180 frequency_resolution * np.arange(num_frequencies) 

181 

182 return frequencies 

183 

184 

185class LombScarglePeriodogramConfig(DiaObjectCalculationPluginConfig): 

186 pass 

187 

188 

189@register("ap_lombScarglePeriodogram") 

190class LombScarglePeriodogram(DiaObjectCalculationPlugin): 

191 """Compute the single-band period of a DiaObject given a set of DiaSources. 

192 """ 

193 ConfigClass = LombScarglePeriodogramConfig 

194 

195 plugType = "multi" 

196 outputCols = ["period", "power"] 

197 needsFilter = True 

198 

199 @classmethod 

200 def getExecutionOrder(cls): 

201 return cls.DEFAULT_CATALOGCALCULATION 

202 

203 @catchWarnings(warns=["All-NaN slice encountered"]) 

204 def calculate(self, 

205 diaObjects, 

206 diaSources, 

207 filterDiaSources, 

208 band): 

209 """Compute the periodogram. 

210 

211 Parameters 

212 ---------- 

213 diaObjects : `pandas.DataFrame` 

214 Summary objects to store values in. 

215 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy` 

216 Catalog of DiaSources summarized by this DiaObject. 

217 """ 

218 

219 # Check and initialize output columns in diaObjects. 

220 if (periodCol := f"{band}_period") not in diaObjects.columns: 

221 diaObjects[periodCol] = np.nan 

222 if (powerCol := f"{band}_power") not in diaObjects.columns: 

223 diaObjects[powerCol] = np.nan 

224 

225 def _calculate_period(df, min_detections=5, nterms=1, oversampling_factor=5, nyquist_factor=100): 

226 """Compute the Lomb-Scargle periodogram given a set of DiaSources. 

227 

228 Parameters 

229 ---------- 

230 df : `pandas.DataFrame` 

231 The input DataFrame. 

232 min_detections : `int`, optional 

233 The minimum number of detections. 

234 nterms : `int`, optional 

235 The number of terms in the Lomb-Scargle model. 

236 oversampling_factor : `int`, optional 

237 The oversampling factor for frequency grid. 

238 nyquist_factor : `int`, optional 

239 The Nyquist factor for frequency grid. 

240 

241 Returns 

242 ------- 

243 pd_tab : `pandas.Series` 

244 The output DataFrame with the Lomb-Scargle parameters. 

245 """ 

246 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]), 

247 np.isnan(df["midpointMjdTai"]))] 

248 

249 if len(tmpDf) < min_detections: 

250 return pd.Series({periodCol: np.nan, powerCol: np.nan}) 

251 

252 time = tmpDf["midpointMjdTai"].to_numpy() 

253 flux = tmpDf["psfFlux"].to_numpy() 

254 flux_err = tmpDf["psfFluxErr"].to_numpy() 

255 

256 lsp = LombScargle(time, flux, dy=flux_err, nterms=nterms) 

257 f_grid = compute_optimized_periodogram_grid( 

258 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor) 

259 period = 1/f_grid 

260 power = lsp.power(f_grid) 

261 

262 return pd.Series({periodCol: period[np.argmax(power)], 

263 powerCol: np.max(power)}) 

264 

265 diaObjects.loc[:, [periodCol, powerCol] 

266 ] = filterDiaSources.apply(_calculate_period) 

267 

268 

269class LombScarglePeriodogramMultiConfig(DiaObjectCalculationPluginConfig): 

270 pass 

271 

272 

273@register("ap_lombScarglePeriodogramMulti") 

274class LombScarglePeriodogramMulti(DiaObjectCalculationPlugin): 

275 """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources. 

276 """ 

277 ConfigClass = LombScarglePeriodogramMultiConfig 

278 

279 plugType = "multi" 

280 outputCols = ["multiPeriod", "multiPower", 

281 "multiFap", "multiAmp", "multiPhase"] 

282 needsFilter = True 

283 

284 @classmethod 

285 def getExecutionOrder(cls): 

286 return cls.DEFAULT_CATALOGCALCULATION 

287 

288 @staticmethod 

289 def calculate_baluev_fap(time, n, maxPeriod, zmax): 

290 """Calculate the False-Alarm probability using the Baluev approximation. 

291 

292 Parameters 

293 ---------- 

294 time : `array` 

295 The input time axis. 

296 n : `int` 

297 The number of detections. 

298 maxPeriod : `float` 

299 The maximum period in the grid. 

300 zmax : `float` 

301 The maximum power in the grid. 

302 

303 Returns 

304 ------- 

305 fap_estimate : `float` 

306 The False-Alarm probability Baluev approximation. 

307 

308 Notes 

309 ---------- 

310 .. [1] Baluev, R. V. 2008, MNRAS, 385, 1279 

311 .. [2] Süveges, M., Guy, L.P., Eyer, L., et al. 2015, MNRAS, 450, 2052 

312 """ 

313 if n <= 2: 

314 return np.nan 

315 

316 gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2)) 

317 fu = 1/maxPeriod 

318 return gam_ratio * np.sqrt( 

319 4*np.pi*statistics.variance(time) 

320 ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax) 

321 

322 @staticmethod 

323 def generate_lsp_params(lsp_model, fbest, bands): 

324 """Generate the Lomb-Scargle parameters. 

325 Parameters 

326 ---------- 

327 lsp_model : `astropy.timeseries.LombScargleMultiband` 

328 The Lomb-Scargle model. 

329 fbest : `float` 

330 The best period. 

331 bands : `array` 

332 The bands of the time series. 

333 

334 Returns 

335 ------- 

336 Amp : `array` 

337 The amplitude of the time series. 

338 Ph : `array` 

339 The phase of the time series. 

340 

341 Notes 

342 ---------- 

343 .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18 

344 """ 

345 best_params = lsp_model.model_parameters(fbest, units=True) 

346 

347 name_params = [f"theta_base_{i}" for i in range(3)] 

348 name_params += [f"theta_band_{band}_{i}" for band in np.unique(bands) for i in range(3)] 

349 

350 df_params = pd.DataFrame([best_params], columns=name_params) 

351 

352 unique_bands = np.unique(bands) 

353 

354 amplitude_band = [np.sqrt(df_params[f"theta_band_{band}_1"]**2 

355 + df_params[f"theta_band_{band}_2"]**2) 

356 for band in unique_bands] 

357 phase_bands = [np.arctan2(df_params[f"theta_band_{band}_2"], 

358 df_params[f"theta_band_{band}_1"]) for band in unique_bands] 

359 

360 amp = [a[0] for a in amplitude_band] 

361 ph = [p[0] for p in phase_bands] 

362 

363 return amp, ph 

364 

365 @catchWarnings(warns=["All-NaN slice encountered"]) 

366 def calculate(self, 

367 diaObjects, 

368 diaSources, 

369 **kwargs): 

370 """Compute the multi-band LombScargle periodogram of a DiaObject given 

371 a set of DiaSources. 

372 

373 Parameters 

374 ---------- 

375 diaObjects : `pandas.DataFrame` 

376 Summary objects to store values in. 

377 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy` 

378 Catalog of DiaSources summarized by this DiaObject. 

379 **kwargs : `dict` 

380 Unused kwargs that are always passed to a plugin. 

381 """ 

382 

383 bands_arr = diaSources['band'].unique().values 

384 unique_bands = np.unique(np.concatenate(bands_arr)) 

385 # Check and initialize output columns in diaObjects. 

386 if (periodCol := "multiPeriod") not in diaObjects.columns: 

387 diaObjects[periodCol] = np.nan 

388 if (powerCol := "multiPower") not in diaObjects.columns: 

389 diaObjects[powerCol] = np.nan 

390 if (fapCol := "multiFap") not in diaObjects.columns: 

391 diaObjects[fapCol] = np.nan 

392 ampCol = "multiAmp" 

393 phaseCol = "multiPhase" 

394 for i in range(len(unique_bands)): 

395 ampCol_band = f"{unique_bands[i]}_{ampCol}" 

396 if ampCol_band not in diaObjects.columns: 

397 diaObjects[ampCol_band] = np.nan 

398 phaseCol_band = f"{unique_bands[i]}_{phaseCol}" 

399 if phaseCol_band not in diaObjects.columns: 

400 diaObjects[phaseCol_band] = np.nan 

401 

402 def _calculate_period_multi(df, all_unique_bands, 

403 min_detections=9, oversampling_factor=5, nyquist_factor=100): 

404 """Calculate the multi-band Lomb-Scargle periodogram. 

405 

406 Parameters 

407 ---------- 

408 df : `pandas.DataFrame` 

409 The input DataFrame. 

410 all_unique_bands : `list` of `str` 

411 List of all bands present in the diaSource table that is being worked on. 

412 min_detections : `int`, optional 

413 The minimum number of detections, including all bands. 

414 oversampling_factor : `int`, optional 

415 The oversampling factor for frequency grid. 

416 nyquist_factor : `int`, optional 

417 The Nyquist factor for frequency grid. 

418 

419 Returns 

420 ------- 

421 pd_tab : `pandas.Series` 

422 The output DataFrame with the Lomb-Scargle parameters. 

423 """ 

424 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]), 

425 np.isnan(df["midpointMjdTai"]))] 

426 

427 if (len(tmpDf)) < min_detections: 

428 pd_tab_nodet = pd.Series({periodCol: np.nan, 

429 powerCol: np.nan, 

430 fapCol: np.nan}) 

431 for band in all_unique_bands: 

432 pd_tab_nodet[f"{band}_{ampCol}"] = np.nan 

433 pd_tab_nodet[f"{band}_{phaseCol}"] = np.nan 

434 

435 return pd_tab_nodet 

436 

437 time = tmpDf["midpointMjdTai"].to_numpy() 

438 flux = tmpDf["psfFlux"].to_numpy() 

439 flux_err = tmpDf["psfFluxErr"].to_numpy() 

440 bands = tmpDf["band"].to_numpy() 

441 

442 lsp = LombScargleMultiband(time, flux, bands, dy=flux_err, 

443 nterms_base=1, nterms_band=1) 

444 

445 f_grid = compute_optimized_periodogram_grid( 

446 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor) 

447 period = 1/f_grid 

448 power = lsp.power(f_grid) 

449 

450 fap_estimate = self.calculate_baluev_fap( 

451 time, len(time), period[np.argmax(power)], np.max(power)) 

452 

453 params_table_new = self.generate_lsp_params(lsp, f_grid[np.argmax(power)], bands) 

454 

455 pd_tab = pd.Series({periodCol: period[np.argmax(power)], 

456 powerCol: np.max(power), 

457 fapCol: fap_estimate 

458 }) 

459 

460 # Initialize the per-band amplitude/phase columns as NaNs 

461 for band in all_unique_bands: 

462 pd_tab[f"{band}_{ampCol}"] = np.nan 

463 pd_tab[f"{band}_{phaseCol}"] = np.nan 

464 

465 # Populate the values of only the bands that have data for this diaSource 

466 unique_bands = np.unique(bands) 

467 for i in range(len(unique_bands)): 

468 pd_tab[f"{unique_bands[i]}_{ampCol}"] = params_table_new[0][i] 

469 pd_tab[f"{unique_bands[i]}_{phaseCol}"] = params_table_new[1][i] 

470 

471 return pd_tab 

472 

473 columns_list = [periodCol, powerCol, fapCol] 

474 for i in range(len(unique_bands)): 

475 columns_list.append(f"{unique_bands[i]}_{ampCol}") 

476 columns_list.append(f"{unique_bands[i]}_{phaseCol}") 

477 

478 diaObjects.loc[:, columns_list 

479 ] = diaSources.apply(_calculate_period_multi, unique_bands) 

480 

481 

482class UnphysicalDiaSourceSeparation(pipeBase.AlgorithmError): 

483 """Raised if associated DiaSources are unphysically separated. 

484 

485 Parameters 

486 ---------- 

487 separation : `float` 

488 Observed separation in arseconds. 

489 max_allowed_separation : `float` 

490 Configured maximum separation in arcseconds. 

491 """ 

492 

493 def __init__(self, separation, max_allowed_separation) -> None: 

494 self._separation = separation 

495 self._max_allowed_separation = max_allowed_separation 

496 super().__init__(f"Observed DiaSource separation {separation} exceeds allowed value of " 

497 f"{max_allowed_separation}") 

498 

499 @property 

500 def metadata(self) -> dict: 

501 return { 

502 "separation": self._separation, 

503 "max_allowed_separation": self._max_allowed_separation, 

504 } 

505 

506 

507class MeanDiaPositionConfig(DiaObjectCalculationPluginConfig): 

508 MaxAllowedDiaSourceSeparation = pexConfig.Field( 

509 dtype=float, 

510 default=3.0, 

511 doc="Max allowed separation of associated DiaSources in arcsec. " 

512 "Raises if unphysical separation is found. " 

513 ) 

514 

515 

516@register("ap_meanPosition") 

517class MeanDiaPosition(DiaObjectCalculationPlugin): 

518 """Compute the mean position of a DiaObject given a set of DiaSources. 

519 """ 

520 

521 ConfigClass = MeanDiaPositionConfig 

522 

523 plugType = 'multi' 

524 

525 outputCols = ["ra", "dec"] 

526 needsFilter = False 

527 

528 @classmethod 

529 def getExecutionOrder(cls): 

530 return cls.DEFAULT_CATALOGCALCULATION 

531 

532 def calculate(self, diaObjects, diaSources, **kwargs): 

533 """Compute the mean ra/dec position of the diaObject given the 

534 diaSource locations. 

535 

536 Parameters 

537 ---------- 

538 diaObjects : `pandas.DataFrame` 

539 Summary objects to store values in. 

540 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy` 

541 Catalog of DiaSources summarized by this DiaObject. 

542 **kwargs 

543 Any additional keyword arguments that may be passed to the plugin. 

544 """ 

545 for outCol in self.outputCols: 

546 if outCol not in diaObjects.columns: 

547 diaObjects[outCol] = np.nan 

548 

549 def _computeMeanPos(df): 

550 coords = list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees) 

551 for idx, src in df.iterrows()) 

552 aveCoord = geom.averageSpherePoint(coords) 

553 

554 # We don't want the DIAObject position to move due to misassociated sources 

555 maxSep = max(aveCoord.separation(coord).asArcseconds() for coord in coords) 

556 

557 if maxSep > self.config.MaxAllowedDiaSourceSeparation: 

558 raise UnphysicalDiaSourceSeparation(maxSep, 

559 self.config.MaxAllowedDiaSourceSeparation) 

560 

561 return pd.Series({"ra": aveCoord.getRa().asDegrees(), 

562 "dec": aveCoord.getDec().asDegrees()}) 

563 

564 ans = diaSources.apply(_computeMeanPos) 

565 typeSafePandasAssignment(diaObjects, ans, ["ra", "dec"]) 

566 

567 

568class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig): 

569 

570 htmLevel = pexConfig.Field( 

571 dtype=int, 

572 doc="Level of the HTM pixelization.", 

573 default=20, 

574 ) 

575 

576 

577@register("ap_HTMIndex") 

578class HTMIndexDiaPosition(DiaObjectCalculationPlugin): 

579 """Compute the mean position of a DiaObject given a set of DiaSources. 

580 

581 Notes 

582 ----- 

583 This plugin was implemented to satisfy requirements of old APDB interface 

584 which required ``pixelId`` column in DiaObject with HTM20 index. APDB 

585 interface had migrated to not need that information, but we keep this 

586 plugin in case it may be useful for something else. 

587 """ 

588 ConfigClass = HTMIndexDiaPositionConfig 

589 

590 plugType = 'single' 

591 

592 inputCols = ["ra", "dec"] 

593 outputCols = ["pixelId"] 

594 needsFilter = False 

595 

596 def __init__(self, config, name, metadata): 

597 DiaObjectCalculationPlugin.__init__(self, config, name, metadata) 

598 self.pixelator = sphgeom.HtmPixelization(self.config.htmLevel) 

599 

600 @classmethod 

601 def getExecutionOrder(cls): 

602 return cls.FLUX_MOMENTS_CALCULATED 

603 

604 def calculate(self, diaObjects, diaObjectId, **kwargs): 

605 """Compute the mean position of a DiaObject given a set of DiaSources 

606 

607 Parameters 

608 ---------- 

609 diaObjects : `pandas.dataFrame` 

610 Summary objects to store values in and read ra/dec from. 

611 diaObjectId : `int` 

612 Id of the diaObject to update. 

613 **kwargs 

614 Any additional keyword arguments that may be passed to the plugin. 

615 """ 

616 sphPoint = geom.SpherePoint( 

617 diaObjects.at[diaObjectId, "ra"] * geom.degrees, 

618 diaObjects.at[diaObjectId, "dec"] * geom.degrees) 

619 diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index( 

620 sphPoint.getVector()) 

621 

622 

623class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig): 

624 pass 

625 

626 

627@register("ap_nDiaSources") 

628class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin): 

629 """Compute the total number of DiaSources associated with this DiaObject. 

630 """ 

631 

632 ConfigClass = NumDiaSourcesDiaPluginConfig 

633 outputCols = ["nDiaSources"] 

634 plugType = "multi" 

635 needsFilter = False 

636 

637 @classmethod 

638 def getExecutionOrder(cls): 

639 return cls.DEFAULT_CATALOGCALCULATION 

640 

641 def calculate(self, diaObjects, diaSources, **kwargs): 

642 """Compute the total number of DiaSources associated with this DiaObject. 

643 

644 Parameters 

645 ---------- 

646 diaObject : `dict` 

647 Summary object to store values in and read ra/dec from. 

648 **kwargs 

649 Any additional keyword arguments that may be passed to the plugin. 

650 """ 

651 typeSafePandasAssignment(diaObjects, diaSources.diaObjectId.count(), ["nDiaSources"]) 

652 

653 

654class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig): 

655 pass 

656 

657 

658@register("ap_diaObjectFlag") 

659class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin): 

660 """Find if any DiaSource is flagged. 

661 

662 Set the DiaObject flag if any DiaSource is flagged. 

663 """ 

664 

665 ConfigClass = NumDiaSourcesDiaPluginConfig 

666 outputCols = ["flags"] 

667 plugType = "multi" 

668 needsFilter = False 

669 

670 @classmethod 

671 def getExecutionOrder(cls): 

672 return cls.DEFAULT_CATALOGCALCULATION 

673 

674 def calculate(self, diaObjects, diaSources, **kwargs): 

675 """Find if any DiaSource is flagged. 

676 

677 Set the DiaObject flag if any DiaSource is flagged. 

678 

679 Parameters 

680 ---------- 

681 diaObject : `dict` 

682 Summary object to store values in and read ra/dec from. 

683 **kwargs 

684 Any additional keyword arguments that may be passed to the plugin. 

685 """ 

686 typeSafePandasAssignment(diaObjects, diaSources.flags.any(), ["flags"], default_dtype=np.uint64) 

687 

688 

689class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig): 

690 pass 

691 

692 

693@register("ap_meanFlux") 

694class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin): 

695 """Compute the weighted mean and mean error on the point source fluxes 

696 of the DiaSource measured on the difference image. 

697 

698 Additionally store number of usable data points. 

699 """ 

700 

701 ConfigClass = WeightedMeanDiaPsfFluxConfig 

702 outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"] 

703 plugType = "multi" 

704 needsFilter = True 

705 

706 @classmethod 

707 def getExecutionOrder(cls): 

708 return cls.DEFAULT_CATALOGCALCULATION 

709 

710 @catchWarnings(warns=["invalid value encountered", 

711 "divide by zero"]) 

712 def calculate(self, 

713 diaObjects, 

714 diaSources, 

715 filterDiaSources, 

716 band, 

717 **kwargs): 

718 """Compute the weighted mean and mean error of the point source flux. 

719 

720 Parameters 

721 ---------- 

722 diaObject : `dict` 

723 Summary object to store values in. 

724 diaSources : `pandas.DataFrame` 

725 DataFrame representing all diaSources associated with this 

726 diaObject. 

727 filterDiaSources : `pandas.DataFrame` 

728 DataFrame representing diaSources associated with this 

729 diaObject that are observed in the band pass ``band``. 

730 band : `str` 

731 Simple, string name of the filter for the flux being calculated. 

732 **kwargs 

733 Any additional keyword arguments that may be passed to the plugin. 

734 """ 

735 meanName = "{}_psfFluxMean".format(band) 

736 errName = "{}_psfFluxMeanErr".format(band) 

737 nDataName = "{}_psfFluxNdata".format(band) 

738 if meanName not in diaObjects.columns: 

739 diaObjects[meanName] = np.nan 

740 if errName not in diaObjects.columns: 

741 diaObjects[errName] = np.nan 

742 if nDataName not in diaObjects.columns: 

743 diaObjects[nDataName] = 0 

744 

745 def _weightedMean(df): 

746 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]), 

747 np.isnan(df["psfFluxErr"]))] 

748 tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2) 

749 fluxMean = np.nansum(tmpDf["psfFlux"] 

750 / tmpDf["psfFluxErr"] ** 2) 

751 fluxMean /= tot_weight 

752 if tot_weight > 0: 

753 fluxMeanErr = np.sqrt(1 / tot_weight) 

754 else: 

755 fluxMeanErr = np.nan 

756 nFluxData = len(tmpDf) 

757 

758 return pd.Series({meanName: fluxMean, 

759 errName: fluxMeanErr, 

760 nDataName: nFluxData}, 

761 dtype="object") 

762 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]]) 

763 # TODO DM-53254: Remove the force_int_to_float hack. 

764 typeSafePandasAssignment(diaObjects, df, [meanName, errName, nDataName], force_int_to_float=True) 

765 

766 

767class PercentileDiaPsfFluxConfig(DiaObjectCalculationPluginConfig): 

768 percentiles = pexConfig.ListField( 

769 dtype=int, 

770 default=[5, 25, 50, 75, 95], 

771 doc="Percentiles to calculate to compute values for. Should be " 

772 "integer values." 

773 ) 

774 

775 

776@register("ap_percentileFlux") 

777class PercentileDiaPsfFlux(DiaObjectCalculationPlugin): 

778 """Compute percentiles of diaSource fluxes. 

779 """ 

780 

781 ConfigClass = PercentileDiaPsfFluxConfig 

782 # Output columns are created upon instantiation of the class. 

783 outputCols = [] 

784 plugType = "multi" 

785 needsFilter = True 

786 

787 def __init__(self, config, name, metadata, **kwargs): 

788 DiaObjectCalculationPlugin.__init__(self, 

789 config, 

790 name, 

791 metadata, 

792 **kwargs) 

793 self.outputCols = ["psfFluxPercentile{:02d}".format(percent) 

794 for percent in self.config.percentiles] 

795 

796 @classmethod 

797 def getExecutionOrder(cls): 

798 return cls.DEFAULT_CATALOGCALCULATION 

799 

800 @catchWarnings(warns=["All-NaN slice encountered"]) 

801 def calculate(self, 

802 diaObjects, 

803 diaSources, 

804 filterDiaSources, 

805 band, 

806 **kwargs): 

807 """Compute the percentile fluxes of the point source flux. 

808 

809 Parameters 

810 ---------- 

811 diaObject : `dict` 

812 Summary object to store values in. 

813 diaSources : `pandas.DataFrame` 

814 DataFrame representing all diaSources associated with this 

815 diaObject. 

816 filterDiaSources : `pandas.DataFrame` 

817 DataFrame representing diaSources associated with this 

818 diaObject that are observed in the band pass ``band``. 

819 band : `str` 

820 Simple, string name of the filter for the flux being calculated. 

821 **kwargs 

822 Any additional keyword arguments that may be passed to the plugin. 

823 """ 

824 pTileNames = [] 

825 for tilePercent in self.config.percentiles: 

826 pTileName = "{}_psfFluxPercentile{:02d}".format(band, 

827 tilePercent) 

828 pTileNames.append(pTileName) 

829 if pTileName not in diaObjects.columns: 

830 diaObjects[pTileName] = np.nan 

831 

832 def _fluxPercentiles(df): 

833 pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles) 

834 return pd.Series( 

835 dict((tileName, pTile) 

836 for tileName, pTile in zip(pTileNames, pTiles))) 

837 

838 typeSafePandasAssignment( 

839 diaObjects, 

840 filterDiaSources.apply(_fluxPercentiles), 

841 pTileNames, 

842 ) 

843 

844 

845class SigmaDiaPsfFluxConfig(DiaObjectCalculationPluginConfig): 

846 pass 

847 

848 

849@register("ap_sigmaFlux") 

850class SigmaDiaPsfFlux(DiaObjectCalculationPlugin): 

851 """Compute scatter of diaSource fluxes. 

852 """ 

853 

854 ConfigClass = SigmaDiaPsfFluxConfig 

855 # Output columns are created upon instantiation of the class. 

856 outputCols = ["psfFluxSigma"] 

857 plugType = "multi" 

858 needsFilter = True 

859 

860 @classmethod 

861 def getExecutionOrder(cls): 

862 return cls.DEFAULT_CATALOGCALCULATION 

863 

864 def calculate(self, 

865 diaObjects, 

866 diaSources, 

867 filterDiaSources, 

868 band, 

869 **kwargs): 

870 """Compute the sigma fluxes of the point source flux. 

871 

872 Parameters 

873 ---------- 

874 diaObject : `dict` 

875 Summary object to store values in. 

876 diaSources : `pandas.DataFrame` 

877 DataFrame representing all diaSources associated with this 

878 diaObject. 

879 filterDiaSources : `pandas.DataFrame` 

880 DataFrame representing diaSources associated with this 

881 diaObject that are observed in the band pass ``band``. 

882 band : `str` 

883 Simple, string name of the filter for the flux being calculated. 

884 **kwargs 

885 Any additional keyword arguments that may be passed to the plugin. 

886 """ 

887 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased 

888 # estimator of scatter (i.e. 'N - 1' instead of 'N'). 

889 column = "{}_psfFluxSigma".format(band) 

890 

891 typeSafePandasAssignment( 

892 diaObjects, 

893 filterDiaSources.psfFlux.std(), 

894 [column], 

895 ) 

896 

897 

898class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig): 

899 pass 

900 

901 

902@register("ap_chi2Flux") 

903class Chi2DiaPsfFlux(DiaObjectCalculationPlugin): 

904 """Compute chi2 of diaSource fluxes. 

905 """ 

906 

907 ConfigClass = Chi2DiaPsfFluxConfig 

908 

909 # Required input Cols 

910 inputCols = ["psfFluxMean"] 

911 # Output columns are created upon instantiation of the class. 

912 outputCols = ["psfFluxChi2"]