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

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: 220 ↛ 222line 220 didn't jump to line 222 because the condition on line 220 was always true

221 diaObjects[periodCol] = np.nan 

222 if (powerCol := f"{band}_power") not in diaObjects.columns: 222 ↛ 225line 222 didn't jump to line 225 because the condition on line 222 was always true

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: 249 ↛ 250line 249 didn't jump to line 250 because the condition on line 249 was never true

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: 313 ↛ 314line 313 didn't jump to line 314 because the condition on line 313 was never true

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: 386 ↛ 388line 386 didn't jump to line 388 because the condition on line 386 was always true

387 diaObjects[periodCol] = np.nan 

388 if (powerCol := "multiPower") not in diaObjects.columns: 388 ↛ 390line 388 didn't jump to line 390 because the condition on line 388 was always true

389 diaObjects[powerCol] = np.nan 

390 if (fapCol := "multiFap") not in diaObjects.columns: 390 ↛ 392line 390 didn't jump to line 392 because the condition on line 390 was always true

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: 396 ↛ 398line 396 didn't jump to line 398 because the condition on line 396 was always true

397 diaObjects[ampCol_band] = np.nan 

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

399 if phaseCol_band not in diaObjects.columns: 399 ↛ 394line 399 didn't jump to line 394 because the condition on line 399 was always true

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 and position uncertainty of a DiaObject 

519 from its associated DiaSources. 

520 

521 The reported (ra, dec) is the inverse-variance weighted spherical 

522 mean of the per-source positions, using only the DiaSources whose 

523 per-source coordinate errors are finite. The reported covariance 

524 is the covariance of that weighted-mean estimator, 

525 

526 C_formal = (sum_i inv(C_i))^-1, 

527 

528 where ``C_i`` is the per-source 2x2 covariance, multiplied by a 

529 *chi-squared scale factor* that inflates the result when the 

530 empirical scatter is larger than the per-source errors predict: 

531 

532 chi2 = sum_i r_i^T inv(C_i) r_i (r_i = tangent-plane 

533 residual from the mean) 

534 dof = 2 * (N - 1) 

535 C_obj = max(1, chi2 / dof) * C_formal 

536 

537 When the data agree with the per-source errors (chi2 ~ dof) the 

538 scale factor is unity and ``C_obj == C_formal``. When the residuals 

539 exceed what the errors predict, the covariance is inflated. 

540 """ 

541 

542 ConfigClass = MeanDiaPositionConfig 

543 

544 plugType = 'multi' 

545 

546 outputCols = ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"] 

547 needsFilter = False 

548 

549 @classmethod 

550 def getExecutionOrder(cls): 

551 return cls.DEFAULT_CATALOGCALCULATION 

552 

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

554 """Compute the mean ra/dec position and its uncertainty for each 

555 DiaObject. 

556 

557 Parameters 

558 ---------- 

559 diaObjects : `pandas.DataFrame` 

560 Summary objects to store values in. 

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

562 Catalog of DiaSources summarized by this DiaObject. 

563 **kwargs 

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

565 """ 

566 for outCol in self.outputCols: 

567 if outCol not in diaObjects.columns: 567 ↛ 566line 567 didn't jump to line 566 because the condition on line 567 was always true

568 diaObjects[outCol] = np.nan 

569 

570 maxAllowedSep = self.config.MaxAllowedDiaSourceSeparation 

571 

572 def _column(df, name): 

573 """Read a column as float64, or return all-NaN if absent.""" 

574 if name in df.columns: 

575 return df[name].to_numpy(dtype=np.float64) 

576 return np.full(len(df), np.nan) 

577 

578 def _computeMeanPos(df): 

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

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

581 # Quick check of the unweighted mean to catch unphysical objects 

582 # and as a fallback. 

583 unweightedAvg = geom.averageSpherePoint(coords) 

584 maxSep = max(unweightedAvg.separation(coord).asArcseconds() for coord in coords) 

585 if maxSep > maxAllowedSep: 

586 raise UnphysicalDiaSourceSeparation(maxSep, maxAllowedSep) 

587 

588 nSrc = len(coords) 

589 raErr = _column(df, "raErr") 

590 decErr = _column(df, "decErr") 

591 raDecCov = _column(df, "ra_dec_Cov") 

592 finiteDiag = (np.isfinite(raErr) & (raErr > 0)) & (np.isfinite(decErr) & (decErr > 0)) 

593 nDiag = int(finiteDiag.sum()) 

594 

595 # Tangent-plane offsets from the unweighted mean reference, 

596 # in degrees (arc-length east, north). 

597 offsets = np.array( 

598 [[off.asDegrees() for off in unweightedAvg.getTangentPlaneOffset(coord)] for coord in coords] 

599 ) if nSrc >= 1 else np.zeros((0, 2)) 

600 

601 aveCoord = unweightedAvg 

602 varRa = np.nan 

603 varDec = np.nan 

604 raDecCovObj = np.nan 

605 

606 if nDiag >= 1: 

607 # Weighted-mean path: include only sources with finite 

608 # per-source coordinate errors. 

609 idx = np.where(finiteDiag)[0] 

610 sigmaRaSq = raErr[idx]**2 

611 sigmaDecSq = decErr[idx]**2 

612 # Use full 2x2 weights only when every included source 

613 # has a finite ra_dec_Cov; otherwise use diagonal weights 

614 # for all of them and emit NaN for the output covariance. 

615 haveFullCov = bool(np.isfinite(raDecCov[idx]).all()) 

616 

617 C = np.zeros((nDiag, 2, 2), dtype=np.float64) 

618 C[:, 0, 0] = sigmaRaSq 

619 C[:, 1, 1] = sigmaDecSq 

620 if haveFullCov: 

621 C[:, 0, 1] = raDecCov[idx] 

622 C[:, 1, 0] = raDecCov[idx] 

623 try: 

624 W = np.linalg.inv(C) 

625 except np.linalg.LinAlgError: 

626 # Singular per-source covariance: drop off-diagonal 

627 # and retry with diagonal weights only. 

628 haveFullCov = False 

629 C[:, 0, 1] = 0.0 

630 C[:, 1, 0] = 0.0 

631 W = np.linalg.inv(C) 

632 

633 W_sum = W.sum(axis=0) 

634 C_formal = np.linalg.inv(W_sum) 

635 

636 # Weighted-mean tangent-plane offset from the unweighted 

637 # reference: mu_offset = (sum W_i)^-1 (sum W_i r_i). Then 

638 # apply that offset to the unweighted mean to get the 

639 # weighted spherical mean. SpherePoint.offset takes a 

640 # bearing measured counter-clockwise from east (bearing 0 

641 # is due east, bearing pi/2 is due north). 

642 mu_offset = C_formal @ np.einsum('nij,nj->i', W, offsets[idx]) 

643 east, north = float(mu_offset[0]), float(mu_offset[1]) 

644 sep_deg = float(np.hypot(east, north)) 

645 if sep_deg > 0.0: 

646 bearing = geom.Angle(float(np.arctan2(north, east)), geom.radians) 

647 aveCoord = unweightedAvg.offset(bearing, sep_deg*geom.degrees) 

648 

649 varRaFormal = float(C_formal[0, 0]) 

650 varDecFormal = float(C_formal[1, 1]) 

651 covFormal = float(C_formal[0, 1]) if haveFullCov else np.nan 

652 

653 if nDiag == 1: 

654 # Single included source: pass its own 2x2 covariance through. 

655 varRa = varRaFormal 

656 varDec = varDecFormal 

657 raDecCovObj = covFormal 

658 else: 

659 # Inflate the variance if the reduced chi**2 fit of the 

660 # diaSource coordinates is greater than 1, to account for 

661 # scatter. 

662 residuals = offsets[idx] - mu_offset 

663 if haveFullCov: 

664 # Vectorized Σₙ rₙᵀ Wₙ rₙ (r = residuals). 

665 chi2 = float(np.einsum('ni,nij,nj->', residuals, W, residuals)) 

666 else: 

667 chi2 = float(np.sum(residuals[:, 0]**2/sigmaRaSq + residuals[:, 1]**2/sigmaDecSq)) 

668 dof = 2*(nDiag - 1) 

669 scale = max(1.0, chi2/dof) 

670 varRa = varRaFormal*scale 

671 varDec = varDecFormal*scale 

672 raDecCovObj = (covFormal*scale if np.isfinite(covFormal) else np.nan) 

673 else: 

674 # No source has finite per-source coordinate errors. Fall 

675 # back to the unweighted mean position plus the empirical 

676 # SEM (when N >= 2), and warn that the per-source errors 

677 # were not usable. 

678 warnings.warn( 

679 "No DiaSources with finite coordinate errors; falling " 

680 "back to the unweighted mean position with scatter-only " 

681 "uncertainty.", 

682 stacklevel=2, 

683 ) 

684 if nSrc >= 2: 

685 sampleCov = np.cov(offsets, rowvar=False, ddof=1) 

686 semCov = sampleCov/nSrc 

687 varRa = float(semCov[0, 0]) 

688 varDec = float(semCov[1, 1]) 

689 raDecCovObj = float(semCov[0, 1]) 

690 # else: nSrc == 1 with no error -- outputs stay NaN. 

691 

692 raErrObj = float(np.sqrt(varRa)) if np.isfinite(varRa) else np.nan 

693 decErrObj = float(np.sqrt(varDec)) if np.isfinite(varDec) else np.nan 

694 

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

696 "dec": aveCoord.getDec().asDegrees(), 

697 "raErr": raErrObj, 

698 "decErr": decErrObj, 

699 "ra_dec_Cov": raDecCovObj}) 

700 

701 ans = diaSources.apply(_computeMeanPos) 

702 typeSafePandasAssignment(diaObjects, ans, ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"]) 

703 

704 

705class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig): 

706 

707 htmLevel = pexConfig.Field( 

708 dtype=int, 

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

710 default=20, 

711 ) 

712 

713 

714@register("ap_HTMIndex") 

715class HTMIndexDiaPosition(DiaObjectCalculationPlugin): 

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

717 

718 Notes 

719 ----- 

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

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

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

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

724 """ 

725 ConfigClass = HTMIndexDiaPositionConfig 

726 

727 plugType = 'single' 

728 

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

730 outputCols = ["pixelId"] 

731 needsFilter = False 

732 

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

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

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

736 

737 @classmethod 

738 def getExecutionOrder(cls): 

739 return cls.FLUX_MOMENTS_CALCULATED 

740 

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

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

743 

744 Parameters 

745 ---------- 

746 diaObjects : `pandas.dataFrame` 

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

748 diaObjectId : `int` 

749 Id of the diaObject to update. 

750 **kwargs 

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

752 """ 

753 sphPoint = geom.SpherePoint( 

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

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

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

757 sphPoint.getVector()) 

758 

759 

760class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig): 

761 pass 

762 

763 

764@register("ap_nDiaSources") 

765class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin): 

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

767 """ 

768 

769 ConfigClass = NumDiaSourcesDiaPluginConfig 

770 outputCols = ["nDiaSources"] 

771 plugType = "multi" 

772 needsFilter = False 

773 

774 @classmethod 

775 def getExecutionOrder(cls): 

776 return cls.DEFAULT_CATALOGCALCULATION 

777 

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

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

780 

781 Parameters 

782 ---------- 

783 diaObject : `dict` 

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

785 **kwargs 

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

787 """ 

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

789 

790 

791class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig): 

792 pass 

793 

794 

795@register("ap_diaObjectFlag") 

796class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin): 

797 """Find if any DiaSource is flagged. 

798 

799 Set the DiaObject flag if any DiaSource is flagged. 

800 """ 

801 

802 ConfigClass = NumDiaSourcesDiaPluginConfig 

803 outputCols = ["flags"] 

804 plugType = "multi" 

805 needsFilter = False 

806 

807 @classmethod 

808 def getExecutionOrder(cls): 

809 return cls.DEFAULT_CATALOGCALCULATION 

810 

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

812 """Find if any DiaSource is flagged. 

813 

814 Set the DiaObject flag if any DiaSource is flagged. 

815 

816 Parameters 

817 ---------- 

818 diaObject : `dict` 

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

820 **kwargs 

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

822 """ 

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

824 

825 

826class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig): 

827 pass 

828 

829 

830@register("ap_meanFlux") 

831class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin): 

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

833 of the DiaSource measured on the difference image. 

834 

835 Additionally store number of usable data points. 

836 """ 

837 

838 ConfigClass = WeightedMeanDiaPsfFluxConfig 

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

840 plugType = "multi" 

841 needsFilter = True 

842 

843 @classmethod 

844 def getExecutionOrder(cls): 

845 return cls.DEFAULT_CATALOGCALCULATION 

846 

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

848 "divide by zero"]) 

849 def calculate(self, 

850 diaObjects, 

851 diaSources, 

852 filterDiaSources, 

853 band, 

854 **kwargs): 

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

856 

857 Parameters 

858 ---------- 

859 diaObject : `dict` 

860 Summary object to store values in. 

861 diaSources : `pandas.DataFrame` 

862 DataFrame representing all diaSources associated with this 

863 diaObject. 

864 filterDiaSources : `pandas.DataFrame` 

865 DataFrame representing diaSources associated with this 

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

867 band : `str` 

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

869 **kwargs 

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

871 """ 

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

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

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