Coverage for tests/test_diaCalculationPlugins.py: 99%

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 

22import warnings 

23 

24from astropy.stats import median_absolute_deviation 

25import numpy as np 

26import pandas as pd 

27from scipy.stats import skew 

28import unittest 

29 

30from lsst.meas.base import ( 

31 MeanDiaPosition, MeanDiaPositionConfig, 

32 HTMIndexDiaPosition, HTMIndexDiaPositionConfig, 

33 NumDiaSourcesDiaPlugin, NumDiaSourcesDiaPluginConfig, 

34 SimpleSourceFlagDiaPlugin, SimpleSourceFlagDiaPluginConfig, 

35 WeightedMeanDiaPsfFlux, WeightedMeanDiaPsfFluxConfig, 

36 PercentileDiaPsfFlux, PercentileDiaPsfFluxConfig, 

37 SigmaDiaPsfFlux, SigmaDiaPsfFluxConfig, 

38 Chi2DiaPsfFlux, Chi2DiaPsfFluxConfig, 

39 MadDiaPsfFlux, MadDiaPsfFluxConfig, 

40 SkewDiaPsfFlux, SkewDiaPsfFluxConfig, 

41 MinMaxDiaPsfFlux, MinMaxDiaPsfFluxConfig, 

42 MaxSlopeDiaPsfFlux, MaxSlopeDiaPsfFluxConfig, 

43 ErrMeanDiaPsfFlux, ErrMeanDiaPsfFluxConfig, 

44 LinearFitDiaPsfFlux, LinearFitDiaPsfFluxConfig, 

45 StetsonJDiaPsfFlux, StetsonJDiaPsfFluxConfig, 

46 WeightedMeanDiaTotFlux, WeightedMeanDiaTotFluxConfig, 

47 SigmaDiaTotFlux, SigmaDiaTotFluxConfig, 

48 LombScarglePeriodogram, LombScarglePeriodogramConfig, 

49 LombScarglePeriodogramMulti, LombScarglePeriodogramMultiConfig, 

50 UnphysicalDiaSourceSeparation) 

51import lsst.utils.tests 

52 

53 

54def run_single_plugin(diaObjectCat, 

55 diaObjectId, 

56 diaSourceCat, 

57 band, 

58 plugin): 

59 """Wrapper for running single plugins. 

60 

61 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run` 

62 

63 Parameters 

64 ---------- 

65 diaObjectCat : `pandas.DataFrame` 

66 Input object catalog to store data into and read from. 

67 diaSourcesCat : `pandas.DataFrame` 

68 DiaSource catalog to read data from and groupby on. 

69 fitlerName : `str` 

70 String name of the filter to process. 

71 plugin : `lsst.ap.association.DiaCalculationPlugin` 

72 Plugin to run. 

73 """ 

74 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False) 

75 diaSourceCat.set_index( 

76 ["diaObjectId", "band", "diaSourceId"], 

77 inplace=True, 

78 drop=False) 

79 

80 objDiaSources = diaSourceCat.loc[diaObjectId] 

81 updatingFilterDiaSources = diaSourceCat.loc[ 

82 (diaObjectId, band), : 

83 ] 

84 

85 plugin.calculate(diaObjects=diaObjectCat, 

86 diaObjectId=diaObjectId, 

87 diaSources=objDiaSources, 

88 filterDiaSources=updatingFilterDiaSources, 

89 band=band) 

90 

91 

92def run_multi_plugin(diaObjectCat, diaSourceCat, band, plugin): 

93 """Wrapper for running multi plugins. 

94 

95 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run` 

96 

97 Parameters 

98 ---------- 

99 diaObjectCat : `pandas.DataFrame` 

100 Input object catalog to store data into and read from. 

101 diaSourcesCat : `pandas.DataFrame` 

102 DiaSource catalog to read data from and groupby on. 

103 filterName : `str` 

104 String name of the filter to process. 

105 plugin : `lsst.ap.association.DiaCalculationPlugin` 

106 Plugin to run. 

107 """ 

108 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False) 

109 diaSourceCat.set_index( 

110 ["diaObjectId", "band", "diaSourceId"], 

111 inplace=True, 

112 drop=False) 

113 

114 updatingFilterDiaSources = diaSourceCat.loc[ 

115 (slice(None), band), : 

116 ] 

117 

118 diaSourcesGB = diaSourceCat.groupby(level=0) 

119 filterDiaSourcesGB = updatingFilterDiaSources.groupby(level=0) 

120 

121 plugin.calculate(diaObjects=diaObjectCat, 

122 diaSources=diaSourcesGB, 

123 filterDiaSources=filterDiaSourcesGB, 

124 band=band) 

125 

126 

127def run_multiband_plugin(diaObjectCat, diaSourceCat, plugin): 

128 """Wrapper for running multi plugins. 

129 

130 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run` 

131 

132 Parameters 

133 ---------- 

134 diaObjectCat : `pandas.DataFrame` 

135 Input object catalog to store data into and read from. 

136 diaSourcesCat : `pandas.DataFrame` 

137 DiaSource catalog to read data from and groupby on. 

138 plugin : `lsst.ap.association.DiaCalculationPlugin` 

139 Plugin to run. 

140 """ 

141 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False) 

142 diaSourceCat.set_index( 

143 ["diaObjectId", "band", "diaSourceId"], 

144 inplace=True, 

145 drop=False) 

146 

147 diaSourcesGB = diaSourceCat.groupby(level=0) 

148 

149 plugin.calculate(diaObjects=diaObjectCat, 

150 diaSources=diaSourcesGB, 

151 ) 

152 

153 

154def make_diaObject_table(objId, plugin, default_value=np.nan, band=None): 

155 """Create a minimal diaObject table with columns required for the plugin 

156 

157 Parameters 

158 ---------- 

159 objId : `int` 

160 The diaObjectId 

161 plugin : `lsst.ap.association.DiaCalculationPlugin` 

162 The plugin that will be run. 

163 default_value : `float` or `int`, optional 

164 Value to set new columns to. 

165 band : `str`, optional 

166 Band designation to append to the plugin columns. 

167 

168 Returns 

169 ------- 

170 diaObjects : `pandas.DataFrame` 

171 Output catalog with the required columns for the plugin. 

172 """ 

173 # Add an extra empty diaObject here. This ensures that 

174 # we properly test the source/object matching implicit 

175 # in the plugin calculations. 

176 diaObjects = {"diaObjectId": [objId, objId + 1]} 

177 for col in plugin.outputCols: 

178 if band is not None: 

179 diaObjects[f"{band}_{col}"] = default_value 

180 else: 

181 diaObjects[col] = default_value 

182 return pd.DataFrame(diaObjects) 

183 

184 

185class TestMeanPosition(unittest.TestCase): 

186 

187 def testCalculate(self): 

188 """Test mean position calculation. 

189 

190 DiaSources here are constructed without per-source coordinate 

191 errors, so each ``run_multi_plugin`` call legitimately triggers 

192 the no-errors warning from the plugin. Suppress it here so the 

193 signal of an unexpected warning elsewhere is not drowned out. 

194 """ 

195 n_sources = 10 

196 objId = 0 

197 

198 # configure a 2 degree max separation 

199 plug = MeanDiaPosition(MeanDiaPositionConfig(MaxAllowedDiaSourceSeparation=7200.0), 

200 "ap_meanPosition", 

201 None) 

202 

203 warnings.filterwarnings("ignore", 

204 message="No DiaSources with finite coordinate errors", 

205 category=UserWarning) 

206 self.addCleanup(warnings.resetwarnings) 

207 

208 # Test expected means in RA. 

209 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

210 diaSources = pd.DataFrame(data={"ra": np.linspace(-1, 1, n_sources), 

211 "dec": np.zeros(n_sources), 

212 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

213 "diaObjectId": n_sources * [objId], 

214 "band": n_sources * ["g"], 

215 "diaSourceId": np.arange(n_sources, 

216 dtype=int)}) 

217 run_multi_plugin(diaObjects, diaSources, "g", plug) 

218 

219 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0) 

220 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0) 

221 

222 # Test expected means in DEC. 

223 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

224 diaSources = pd.DataFrame(data={"ra": np.zeros(n_sources), 

225 "dec": np.linspace(-1, 1, n_sources), 

226 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

227 "diaObjectId": n_sources * [objId], 

228 "band": n_sources * ["g"], 

229 "diaSourceId": np.arange(n_sources, 

230 dtype=int)}) 

231 run_multi_plugin(diaObjects, diaSources, "g", plug) 

232 

233 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0) 

234 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0) 

235 

236 # Test failure mode RA is nan. 

237 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

238 diaSources = pd.DataFrame(data={"ra": np.full(n_sources, np.nan), 

239 "dec": np.zeros(n_sources), 

240 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

241 "diaObjectId": n_sources * [objId], 

242 "band": n_sources * ["g"], 

243 "diaSourceId": np.arange(n_sources, 

244 dtype=int)}) 

245 run_multi_plugin(diaObjects, diaSources, "g", plug) 

246 

247 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra"])) 

248 self.assertTrue(np.isnan(diaObjects.loc[objId, "dec"])) 

249 

250 # Test failure mode DEC is nan. 

251 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

252 diaSources = pd.DataFrame(data={"ra": np.zeros(n_sources), 

253 "dec": np.full(n_sources, np.nan), 

254 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

255 "diaObjectId": n_sources * [objId], 

256 "band": n_sources * ["g"], 

257 "diaSourceId": np.arange(n_sources, 

258 dtype=int)}) 

259 run_multi_plugin(diaObjects, diaSources, "g", plug) 

260 

261 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra"])) 

262 self.assertTrue(np.isnan(diaObjects.loc[objId, "dec"])) 

263 

264 # configure the default 3 arcsecond separation 

265 plug = MeanDiaPosition(MeanDiaPositionConfig(MaxAllowedDiaSourceSeparation=3.0), 

266 "ap_meanPosition", 

267 None) 

268 

269 # These 1 degree separations should raise 

270 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

271 diaSources = pd.DataFrame(data={"ra": np.linspace(-1, 1, n_sources), 

272 "dec": np.zeros(n_sources), 

273 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

274 "diaObjectId": n_sources * [objId], 

275 "band": n_sources * ["g"], 

276 "diaSourceId": np.arange(n_sources, 

277 dtype=int)}) 

278 with self.assertRaises(UnphysicalDiaSourceSeparation): 

279 run_multi_plugin(diaObjects, diaSources, "g", plug) 

280 

281 # 1 arcsecond separations should not raise 

282 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

283 diaSources = pd.DataFrame(data={"ra": np.linspace(-1/3600., 1/3600., n_sources), 

284 "dec": np.zeros(n_sources), 

285 "midpointMjdTai": np.linspace(0, n_sources, n_sources), 

286 "diaObjectId": n_sources * [objId], 

287 "band": n_sources * ["g"], 

288 "diaSourceId": np.arange(n_sources, 

289 dtype=int)}) 

290 run_multi_plugin(diaObjects, diaSources, "g", plug) 

291 

292 def _makeDiaSourcesWithUncertainties(self, raErr, decErr, ra_dec_Cov, ra=None, dec=None, objId=0): 

293 """Build a tiny DiaSource DataFrame with per-source uncertainties. 

294 """ 

295 n = len(raErr) 

296 if ra is None: 

297 ra = np.zeros(n) 

298 if dec is None: 

299 dec = np.zeros(n) 

300 return pd.DataFrame(data={ 

301 "ra": ra, 

302 "dec": dec, 

303 "raErr": raErr, 

304 "decErr": decErr, 

305 "ra_dec_Cov": ra_dec_Cov, 

306 "midpointMjdTai": np.arange(n, dtype=float), 

307 "diaObjectId": n * [objId], 

308 "band": n * ["g"], 

309 "diaSourceId": np.arange(n, dtype=int), 

310 }) 

311 

312 def testUncertaintyDiagonalOnly(self): 

313 """Two coincident sources, no off-diagonal: chi^2 = 0, scale 

314 factor is 1, so the output is the diagonal inverse-variance 

315 weighted-mean covariance. 

316 """ 

317 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

318 objId = 0 

319 raErr = np.array([1e-6, 2e-6]) 

320 decErr = np.array([1e-6, 2e-6]) 

321 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

322 diaSources = self._makeDiaSourcesWithUncertainties( 

323 raErr=raErr, 

324 decErr=decErr, 

325 ra_dec_Cov=np.array([np.nan, np.nan]), 

326 ) 

327 run_multi_plugin(diaObjects, diaSources, "g", plug) 

328 

329 # Var(weighted mean) = 1 / sum(1/sigma^2) per axis. 

330 expectedRaErr = 1.0/np.sqrt(np.sum(1.0/raErr**2)) 

331 expectedDecErr = 1.0/np.sqrt(np.sum(1.0/decErr**2)) 

332 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], expectedRaErr) 

333 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], expectedDecErr) 

334 # No per-source ra_dec_Cov, so output ra_dec_Cov is NaN. 

335 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra_dec_Cov"])) 

336 

337 def testUncertaintyFullCovariance(self): 

338 """Two coincident sources with off-diagonal covariance: chi^2 = 0, 

339 so the output is C_formal = (sum_i inv(C_i))^-1. 

340 """ 

341 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

342 objId = 0 

343 raErr = np.array([1e-6, 2e-6]) 

344 decErr = np.array([1.5e-6, 1.0e-6]) 

345 rho = np.array([0.3, -0.2]) # correlation coefficient 

346 raDecCov = rho * raErr * decErr 

347 

348 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

349 diaSources = self._makeDiaSourcesWithUncertainties( 

350 raErr=raErr, decErr=decErr, ra_dec_Cov=raDecCov) 

351 run_multi_plugin(diaObjects, diaSources, "g", plug) 

352 

353 n = len(raErr) 

354 cov = np.zeros((n, 2, 2)) 

355 cov[:, 0, 0] = raErr**2 

356 cov[:, 1, 1] = decErr**2 

357 cov[:, 0, 1] = raDecCov 

358 cov[:, 1, 0] = raDecCov 

359 covObj = np.linalg.inv(np.linalg.inv(cov).sum(axis=0)) 

360 

361 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], np.sqrt(covObj[0, 0])) 

362 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], np.sqrt(covObj[1, 1])) 

363 self.assertAlmostEqual(diaObjects.loc[objId, "ra_dec_Cov"], covObj[0, 1]) 

364 

365 def testUncertaintySingleSourceCopiesThrough(self): 

366 """A single-source group: outputs equal the source's uncertainties. 

367 """ 

368 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

369 objId = 0 

370 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

371 diaSources = self._makeDiaSourcesWithUncertainties( 

372 raErr=np.array([1.5e-6]), 

373 decErr=np.array([0.8e-6]), 

374 ra_dec_Cov=np.array([2e-13]), 

375 ) 

376 run_multi_plugin(diaObjects, diaSources, "g", plug) 

377 

378 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], 1.5e-6) 

379 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], 0.8e-6) 

380 self.assertAlmostEqual(diaObjects.loc[objId, "ra_dec_Cov"], 2e-13) 

381 

382 def testUncertaintyMissingColumnsFallsBackToScatter(self): 

383 """No raErr/decErr columns + N>=2 spread-out sources -> the 

384 weighted-mean path is unusable, so the DiaObject falls back to 

385 the unweighted mean position with the scatter-only SEM, and a 

386 warning is emitted. 

387 """ 

388 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

389 objId = 0 

390 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

391 n = 3 

392 ra = np.linspace(-1e-4, 1e-4, n) 

393 dec = np.zeros(n) 

394 diaSources = pd.DataFrame(data={ 

395 "ra": ra, 

396 "dec": dec, 

397 "midpointMjdTai": np.arange(n, dtype=float), 

398 "diaObjectId": n * [objId], 

399 "band": n * ["g"], 

400 "diaSourceId": np.arange(n, dtype=int), 

401 }) 

402 with self.assertWarnsRegex(UserWarning, "No DiaSources with finite coordinate errors"): 

403 run_multi_plugin(diaObjects, diaSources, "g", plug) 

404 

405 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0) 

406 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0) 

407 

408 # Expected scatter term: at dec=0, the tangent-plane east offsets 

409 # equal the RA deltas (in degrees) to high precision, so the 

410 # standard error of the mean in RA is sample_std(ra, ddof=1) / 

411 # sqrt(N). Dec is identically zero so its scatter is zero. 

412 expectedRaErr = np.std(ra, ddof=1)/np.sqrt(n) 

413 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], expectedRaErr) 

414 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], 0.0) 

415 self.assertAlmostEqual(diaObjects.loc[objId, "ra_dec_Cov"], 0.0) 

416 

417 def testUncertaintyMissingColumnsSingleSourceEmitsNaN(self): 

418 """With a single source and no raErr/decErr the uncertainty is 

419 undefined but the mean position is still computed; a warning is 

420 emitted that the per-source errors are unusable. 

421 """ 

422 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

423 objId = 0 

424 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

425 diaSources = pd.DataFrame(data={ 

426 "ra": [0.0], 

427 "dec": [0.0], 

428 "midpointMjdTai": [0.0], 

429 "diaObjectId": [objId], 

430 "band": ["g"], 

431 "diaSourceId": [0], 

432 }) 

433 with self.assertWarnsRegex(UserWarning, "No DiaSources with finite coordinate errors"): 

434 run_multi_plugin(diaObjects, diaSources, "g", plug) 

435 

436 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0) 

437 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0) 

438 self.assertTrue(np.isnan(diaObjects.loc[objId, "raErr"])) 

439 self.assertTrue(np.isnan(diaObjects.loc[objId, "decErr"])) 

440 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra_dec_Cov"])) 

441 

442 def testUncertaintyPartialCovarianceFallsBackToDiagonal(self): 

443 """Some sources have NaN ra_dec_Cov but valid raErr/decErr. 

444 

445 When at least one included source lacks a finite ra_dec_Cov, 

446 diagonal weights are used for all included sources, the off- 

447 diagonal output is NaN, and the chi-squared falls back to its 

448 diagonal form. Sources are coincident here, so chi^2 = 0 and 

449 there is no rescaling. 

450 """ 

451 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

452 objId = 0 

453 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

454 raErr = np.array([1e-6, 2e-6, 3e-6]) 

455 decErr = np.array([1e-6, 2e-6, 3e-6]) 

456 # Only one source has a finite covariance 

457 raDecCov = np.array([1e-13, np.nan, np.nan]) 

458 diaSources = self._makeDiaSourcesWithUncertainties(raErr=raErr, decErr=decErr, ra_dec_Cov=raDecCov) 

459 run_multi_plugin(diaObjects, diaSources, "g", plug) 

460 

461 # Diagonal inverse-variance weighted mean covariance. 

462 expectedRaErr = 1.0/np.sqrt(np.sum(1.0/raErr**2)) 

463 expectedDecErr = 1.0/np.sqrt(np.sum(1.0/decErr**2)) 

464 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], expectedRaErr) 

465 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], expectedDecErr) 

466 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra_dec_Cov"])) 

467 

468 def testUncertaintyScaleFactorInflatesWhenScatterExceedsErrors(self): 

469 """Two sources whose positional separation is much larger than 

470 their per-source errors: chi^2 >> dof, so the chi-squared scale 

471 factor inflates the formal covariance by chi^2 / dof. 

472 """ 

473 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

474 objId = 0 

475 # Two sources symmetric around RA = 0, Dec = 0. Spread of 

476 # 0.36 arcsec on each axis (within MaxAllowedDiaSourceSeparation 

477 # of 3 arcsec), still 100x the per-source sigma of 1e-6 deg = 3.6 mas, 

478 # so chi^2 / dof ~ 4e4 / 2 ~ 2e4. 

479 delta = 1e-4 # degrees 

480 ra = np.array([-delta, delta]) 

481 dec = np.array([-delta, delta]) 

482 raErr = np.array([1e-6, 1e-6]) 

483 decErr = np.array([1e-6, 1e-6]) 

484 raDecCov = np.array([np.nan, np.nan]) 

485 

486 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

487 diaSources = self._makeDiaSourcesWithUncertainties( 

488 raErr=raErr, decErr=decErr, ra_dec_Cov=raDecCov, ra=ra, dec=dec) 

489 run_multi_plugin(diaObjects, diaSources, "g", plug) 

490 

491 # With equal per-source errors the weighted mean equals the 

492 # unweighted mean = (0, 0), so the tangent-plane residuals are 

493 # essentially (ra, dec). 

494 # chi^2 = sum_i (r_east^2/sigma_ra^2 + r_north^2/sigma_dec^2) 

495 # dof = 2 * (N - 1) = 2. 

496 chi2 = np.sum(ra**2/raErr**2 + dec**2/decErr**2) 

497 dof = 2*(len(raErr) - 1) 

498 scale = max(1.0, chi2/dof) 

499 # Diagonal weighted-mean variance: 1/sum(1/sigma^2) per axis. 

500 varRaFormal = 1.0/np.sum(1.0/raErr**2) 

501 varDecFormal = 1.0/np.sum(1.0/decErr**2) 

502 expectedRaErr = np.sqrt(scale*varRaFormal) 

503 expectedDecErr = np.sqrt(scale*varDecFormal) 

504 

505 # The inflation should be large (scale ~ 1e6): 

506 self.assertGreater(scale, 1e3) 

507 

508 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"] / expectedRaErr, 1.0, places=5) 

509 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"] / expectedDecErr, 1.0, places=5) 

510 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra_dec_Cov"])) 

511 

512 def testUncertaintyScaleFactorNoInflationWhenConsistent(self): 

513 """When per-source positions are coincident, chi^2 = 0 and the 

514 chi-squared scale factor is exactly 1: the output equals 

515 C_formal = (sum_i inv(C_i))^-1. 

516 """ 

517 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

518 objId = 0 

519 raErr = np.array([1e-6, 2e-6, 1.5e-6]) 

520 decErr = np.array([1e-6, 1.5e-6, 2e-6]) 

521 raDecCov = np.array([2e-13, -1e-13, 3e-14]) 

522 

523 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

524 diaSources = self._makeDiaSourcesWithUncertainties(raErr=raErr, decErr=decErr, ra_dec_Cov=raDecCov) 

525 run_multi_plugin(diaObjects, diaSources, "g", plug) 

526 

527 n = len(raErr) 

528 cov = np.zeros((n, 2, 2)) 

529 cov[:, 0, 0] = raErr**2 

530 cov[:, 1, 1] = decErr**2 

531 cov[:, 0, 1] = raDecCov 

532 cov[:, 1, 0] = raDecCov 

533 covObj = np.linalg.inv(np.linalg.inv(cov).sum(axis=0)) 

534 

535 self.assertAlmostEqual(diaObjects.loc[objId, "raErr"], np.sqrt(covObj[0, 0])) 

536 self.assertAlmostEqual(diaObjects.loc[objId, "decErr"], np.sqrt(covObj[1, 1])) 

537 self.assertAlmostEqual(diaObjects.loc[objId, "ra_dec_Cov"], covObj[0, 1]) 

538 

539 def testWeightedMeanPositionDiffersFromUnweighted(self): 

540 """Two sources at asymmetric (ra) positions with very unequal 

541 per-source errors: the reported position is the inverse-variance 

542 weighted mean, which is much closer to the more-precise source 

543 than to the unweighted midpoint. 

544 """ 

545 plug = MeanDiaPosition(MeanDiaPositionConfig(), "ap_meanPosition", None) 

546 objId = 0 

547 d = 1e-5 # degrees; well inside MaxAllowedDiaSourceSeparation. 

548 ra = np.array([-d, d]) 

549 dec = np.array([0.0, 0.0]) 

550 # Source 0 is 100x more precise than source 1. 

551 raErr = np.array([1e-6, 1e-4]) 

552 decErr = np.array([1e-6, 1e-4]) 

553 raDecCov = np.array([np.nan, np.nan]) 

554 

555 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

556 diaSources = self._makeDiaSourcesWithUncertainties( 

557 raErr=raErr, decErr=decErr, ra_dec_Cov=raDecCov, ra=ra, dec=dec) 

558 run_multi_plugin(diaObjects, diaSources, "g", plug) 

559 

560 # Diagonal weighted mean of (ra_0, ra_1) with weights w_i = 1/raErr_i^2. 

561 w = 1.0/raErr**2 

562 expectedRa = np.sum(w*ra)/np.sum(w) 

563 # Output RA is wrapped to [0, 360); normalize back to a signed 

564 # offset near zero before comparing. 

565 outRa = ((diaObjects.loc[objId, "ra"] + 180.0) % 360.0) - 180.0 

566 # The unweighted mean would be 0; the weighted mean should be 

567 # very close to source 0 at ra = -d. 

568 self.assertAlmostEqual(outRa, expectedRa) 

569 self.assertLess(abs(outRa - (-d)), 0.01*d) 

570 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0) 

571 

572 

573class TestHTMIndexPosition(unittest.TestCase): 

574 

575 def testCalculate(self): 

576 """Test HTMPixel assignment calculation. 

577 """ 

578 # Test expected pixelId at RA, DEC = 0 

579 objId = 0 

580 n_sources = 10 

581 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

582 diaObjects.loc[objId, "ra"] = 0. 

583 diaObjects.loc[objId, "dec"] = 0. 

584 diaSources = pd.DataFrame( 

585 data={"diaObjectId": n_sources * [objId], 

586 "band": n_sources * ["g"], 

587 "diaSourceId": np.arange(n_sources, dtype=int)}) 

588 plug = HTMIndexDiaPosition(HTMIndexDiaPositionConfig(), 

589 "ap_HTMIndex", 

590 None) 

591 

592 run_single_plugin(diaObjectCat=diaObjects, 

593 diaObjectId=objId, 

594 diaSourceCat=diaSources, 

595 band="g", 

596 plugin=plug) 

597 self.assertEqual(diaObjects.at[objId, "pixelId"], 

598 17042430230528) 

599 

600 # Test expected pixelId at some value of RA and DEC. 

601 diaObjects = pd.DataFrame({"diaObjectId": [objId]}) 

602 diaObjects.loc[objId, "ra"] = 45.37 

603 diaObjects.loc[objId, "dec"] = 13.67 

604 diaSources = pd.DataFrame( 

605 data={"diaObjectId": n_sources * [objId], 

606 "band": n_sources * ["g"], 

607 "diaSourceId": np.arange(n_sources, dtype=int)}) 

608 run_single_plugin(diaObjectCat=diaObjects, 

609 diaObjectId=objId, 

610 diaSourceCat=diaSources, 

611 band="g", 

612 plugin=plug) 

613 self.assertEqual(diaObjects.at[objId, "pixelId"], 

614 17450571968473) 

615 

616 

617class TestNDiaSourcesDiaPlugin(unittest.TestCase): 

618 

619 def testCalculate(self): 

620 """Test that the number of DiaSources is correct. 

621 """ 

622 

623 for n_sources in [1, 8, 10]: 

624 # Test expected number of sources per object. 

625 objId = 0 

626 diaSources = pd.DataFrame( 

627 data={"diaObjectId": n_sources * [objId], 

628 "band": n_sources * ["g"], 

629 "diaSourceId": np.arange(n_sources, dtype=int)}) 

630 plug = NumDiaSourcesDiaPlugin(NumDiaSourcesDiaPluginConfig(), 

631 "ap_nDiaSources", 

632 None) 

633 diaObjects = make_diaObject_table(objId, plug, default_value=int(-1)) 

634 run_multi_plugin(diaObjects, diaSources, "g", plug) 

635 

636 self.assertEqual(n_sources, diaObjects.at[objId, "nDiaSources"]) 

637 self.assertEqual(diaObjects["nDiaSources"].dtype, np.int64) 

638 

639 

640class TestSimpleSourceFlagDiaPlugin(unittest.TestCase): 

641 

642 def testCalculate(self): 

643 """Test that DiaObject flags are set. 

644 """ 

645 objId = 0 

646 n_sources = 10 

647 

648 # Test expected flags, no flags set. 

649 diaSources = pd.DataFrame( 

650 data={"diaObjectId": n_sources * [objId], 

651 "band": n_sources * ["g"], 

652 "diaSourceId": np.arange(n_sources, dtype=int), 

653 "flags": np.zeros(n_sources, dtype=np.uint64)}) 

654 plug = SimpleSourceFlagDiaPlugin(SimpleSourceFlagDiaPluginConfig(), 

655 "ap_diaObjectFlag", 

656 None) 

657 

658 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64(0)) 

659 run_multi_plugin(diaObjects, diaSources, "g", plug) 

660 self.assertEqual(diaObjects.at[objId, "flags"], 0) 

661 self.assertEqual(diaObjects["flags"].dtype, np.uint64) 

662 

663 # Test expected flags, all flags set. 

664 diaSources = pd.DataFrame( 

665 data={"diaObjectId": n_sources * [objId], 

666 "band": n_sources * ["g"], 

667 "diaSourceId": np.arange(n_sources, dtype=int), 

668 "flags": np.ones(n_sources, dtype=np.uint64)}) 

669 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64(0)) 

670 run_multi_plugin(diaObjects, diaSources, "g", plug) 

671 self.assertEqual(diaObjects.at[objId, "flags"], 1) 

672 self.assertEqual(diaObjects["flags"].dtype, np.uint64) 

673 

674 # Test expected flags, random flags. 

675 diaSources = pd.DataFrame( 

676 data={"diaObjectId": n_sources * [objId], 

677 "band": n_sources * ["g"], 

678 "diaSourceId": np.arange(n_sources, dtype=int), 

679 "flags": np.random.randint(0, 2 ** 16, size=n_sources)}) 

680 

681 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64(0)) 

682 run_multi_plugin(diaObjects, diaSources, "g", plug) 

683 self.assertEqual(diaObjects.at[objId, "flags"], 1) 

684 self.assertEqual(diaObjects["flags"].dtype, np.uint64) 

685 

686 # Test expected flags, one flag set. 

687 flag_array = np.zeros(n_sources, dtype=np.uint64) 

688 flag_array[4] = 256 

689 diaSources = pd.DataFrame( 

690 data={"diaObjectId": n_sources * [objId], 

691 "band": n_sources * ["g"], 

692 "diaSourceId": np.arange(n_sources, dtype=int), 

693 "flags": flag_array}) 

694 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64(0)) 

695 run_multi_plugin(diaObjects, diaSources, "g", plug) 

696 self.assertEqual(diaObjects.at[objId, "flags"], 1) 

697 self.assertEqual(diaObjects["flags"].dtype, np.uint64) 

698 

699 

700class TestWeightedMeanDiaPsfFlux(unittest.TestCase): 

701 

702 def testCalculate(self): 

703 """Test mean value calculation. 

704 """ 

705 n_sources = 10 

706 objId = 0 

707 

708 # Test expected mean. 

709 # In the first test, we have only one object. 

710 diaObjects = pd.DataFrame({"diaObjectId": [objId]})