Coverage for python / lsst / ap / association / diaPipe.py: 19%

1# 

2# LSST Data Management System 

3# Copyright 2008-2016 AURA/LSST. 

4# 

5# This product includes software developed by the 

6# LSST Project (http://www.lsst.org/). 

7# 

8# This program is free software: you can redistribute it and/or modify 

9# it under the terms of the GNU General Public License as published by 

10# the Free Software Foundation, either version 3 of the License, or 

11# (at your option) any later version. 

12# 

13# This program is distributed in the hope that it will be useful, 

14# but WITHOUT ANY WARRANTY; without even the implied warranty of 

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <https://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23"""PipelineTask for associating DiaSources with previous DiaObjects. 

24 

25Additionally performs forced photometry on the calibrated and difference 

26images at the updated locations of DiaObjects. 

27""" 

28 

29__all__ = ("DiaPipelineConfig", 

30 "DiaPipelineTask", 

31 "DiaPipelineConnections") 

32 

33 

34import lsst.dax.apdb as daxApdb 

35import lsst.pex.config as pexConfig 

36import lsst.pipe.base as pipeBase 

37import lsst.pipe.base.connectionTypes as connTypes 

38import lsst.sphgeom 

39 

40from astropy.table import Table 

41import numpy as np 

42import pandas as pd 

43from lsst.ap.association import ( 

44 AssociationTask, 

45 DiaForcedSourceTask, 

46 PackageAlertsTask) 

47 

48from lsst.ap.association.utils import makeEmptyForcedSourceTable, getRegion, paddedRegion, readSchemaFromApdb 

49from lsst.daf.base import DateTime 

50from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

51 DiaObjectCalculationTask 

52from lsst.pipe.tasks.schemaUtils import convertDataFrameToSdmSchema, checkSdmSchemaColumns, \ 

53 dropEmptyColumns, make_empty_catalog 

54from lsst.pipe.tasks.ssoAssociation import SolarSystemAssociationTask 

55from lsst.utils.timer import timeMethod, duration_from_timeMethod 

56 

57 

58class TooManyDiaObjectsError(pipeBase.AlgorithmError): 

59 """Raised if there are an unusually large number of unassociated DiaSources. 

60 This is usually indicative of an image subtraction error, and needs to be 

61 caught before updating the APDB with a large number of spurious entries. 

62 """ 

63 def __init__(self, *, nNewDiaObjects, threshold): 

64 msg = ("Aborting processing before writing to the APDB." 

65 f" {nNewDiaObjects} new DiaObjects would be created, which exceeds the" 

66 f" configured maximum of {threshold}") 

67 super().__init__(msg) 

68 self.nNewDiaObjects = nNewDiaObjects 

69 self.threshold = threshold 

70 

71 @property 

72 def metadata(self): 

73 return {"nNewDiaObjects": self.nNewDiaObjects, 

74 "threshold": self.threshold 

75 } 

76 

77 

78class PostApdbUpdateError(pipeBase.AlgorithmError): 

79 """Raised for any error that occurs after the APDB has been updated. 

80 This allows partial outputs to be written, and signals that processing 

81 can't be retried. 

82 """ 

83 def __init__(self, *, errorMsg): 

84 msg = ("Aborting processing after writing to the APDB. " 

85 "This image cannot be retried! " 

86 "The original error was:\n" 

87 f"{errorMsg}") 

88 super().__init__(msg) 

89 

90 @property 

91 def metadata(self): 

92 return {} 

93 

94 

95class DiaPipelineConnections( 

96 pipeBase.PipelineTaskConnections, 

97 dimensions=("instrument", "visit", "detector"), 

98 defaultTemplates={"coaddName": "deep", "fakesType": ""}): 

99 """Butler connections for DiaPipelineTask. 

100 """ 

101 diaSourceTable = connTypes.Input( 

102 doc="Catalog of calibrated DiaSources.", 

103 name="{fakesType}{coaddName}Diff_diaSrcTable", 

104 storageClass="DataFrame", 

105 dimensions=("instrument", "visit", "detector"), 

106 ) 

107 solarSystemObjectTable = connTypes.Input( 

108 doc="Catalog of SolarSolarSystem objects expected to be observable in " 

109 "this detectorVisit.", 

110 name="preloaded_SsObjects", 

111 storageClass="ArrowAstropy", 

112 dimensions=("instrument", "group", "detector"), 

113 minimum=0, 

114 ) 

115 diffIm = connTypes.Input( 

116 doc="Difference image on which the DiaSources were detected.", 

117 name="{fakesType}{coaddName}Diff_differenceExp", 

118 storageClass="ExposureF", 

119 dimensions=("instrument", "visit", "detector"), 

120 ) 

121 exposure = connTypes.Input( 

122 doc="Calibrated exposure differenced with a template image during " 

123 "image differencing.", 

124 name="{fakesType}calexp", 

125 storageClass="ExposureF", 

126 dimensions=("instrument", "visit", "detector"), 

127 ) 

128 template = connTypes.Input( 

129 doc="Warped template used to create `subtractedExposure`. Not PSF " 

130 "matched.", 

131 dimensions=("instrument", "visit", "detector"), 

132 storageClass="ExposureF", 

133 name="{fakesType}{coaddName}Diff_templateExp", 

134 ) 

135 preloadedDiaObjects = connTypes.Input( 

136 doc="DiaObjects preloaded from the APDB.", 

137 name="preloaded_diaObjects", 

138 storageClass="DataFrame", 

139 dimensions=("instrument", "group", "detector"), 

140 ) 

141 preloadedDiaSources = connTypes.Input( 

142 doc="DiaSources preloaded from the APDB.", 

143 name="preloaded_diaSources", 

144 storageClass="DataFrame", 

145 dimensions=("instrument", "group", "detector"), 

146 ) 

147 preloadedDiaForcedSources = connTypes.Input( 

148 doc="DiaForcedSources preloaded from the APDB.", 

149 name="preloaded_diaForcedSources", 

150 storageClass="DataFrame", 

151 dimensions=("instrument", "group", "detector"), 

152 ) 

153 apdbMarker = connTypes.Output( 

154 doc="Marker dataset storing the configuration of the Apdb for each " 

155 "visit/detector. Used to signal the completion of the pipeline.", 

156 name="apdb_marker", 

157 storageClass="Config", 

158 dimensions=("instrument", "visit", "detector"), 

159 ) 

160 associatedDiaSources = connTypes.Output( 

161 doc="Optional output storing the DiaSource catalog after matching, " 

162 "calibration, and standardization for insertion into the Apdb.", 

163 name="{fakesType}{coaddName}Diff_assocDiaSrc", 

164 storageClass="ArrowAstropy", 

165 dimensions=("instrument", "visit", "detector"), 

166 ) 

167 associatedSsSources = connTypes.Output( 

168 doc="Optional output storing ssSource data computed during association.", 

169 name="{fakesType}{coaddName}Diff_associatedSsSources", 

170 storageClass="ArrowAstropy", 

171 dimensions=("instrument", "visit", "detector"), 

172 ) 

173 unassociatedSsObjects = connTypes.Output( 

174 doc="Expected locations of an ssObject with no source", 

175 name="ssUnassociatedObjects", 

176 storageClass="ArrowAstropy", 

177 dimensions=("instrument", "visit", "detector"), 

178 ) 

179 

180 diaForcedSources = connTypes.Output( 

181 doc="Optional output storing the forced sources computed at the diaObject positions.", 

182 name="{fakesType}{coaddName}Diff_diaForcedSrc", 

183 storageClass="ArrowAstropy", 

184 dimensions=("instrument", "visit", "detector"), 

185 ) 

186 diaObjects = connTypes.Output( 

187 doc="Optional output storing the updated diaObjects associated to these sources.", 

188 name="{fakesType}{coaddName}Diff_diaObject", 

189 storageClass="ArrowAstropy", 

190 dimensions=("instrument", "visit", "detector"), 

191 ) 

192 newDiaSources = connTypes.Output( 

193 doc="New diaSources not associated with an existing diaObject that" 

194 " were used to create a new diaObject", 

195 name="{fakesType}new_dia_source", 

196 storageClass="ArrowAstropy", 

197 dimensions=("instrument", "visit", "detector"), 

198 ) 

199 marginalDiaSources = connTypes.Output( 

200 doc="Low SNR diaSources not associated with an existing diaObject that" 

201 " were rejected instead of creating a new diaObject", 

202 name="{fakesType}marginal_new_dia_source", 

203 storageClass="ArrowAstropy", 

204 dimensions=("instrument", "visit", "detector"), 

205 ) 

206 

207 def __init__(self, *, config=None): 

208 super().__init__(config=config) 

209 

210 if not config.doWriteAssociatedSources: 

211 self.outputs.remove("associatedDiaSources") 

212 self.outputs.remove("diaForcedSources") 

213 self.outputs.remove("diaObjects") 

214 self.outputs.remove("newDiaSources") 

215 self.outputs.remove("marginalDiaSources") 

216 else: 

217 if not config.doRunForcedMeasurement: 

218 self.outputs.remove("diaForcedSources") 

219 if not config.filterUnAssociatedSources: 

220 self.outputs.remove("newDiaSources") 

221 self.outputs.remove("marginalDiaSources") 

222 if not config.doSolarSystemAssociation: 

223 self.inputs.remove("solarSystemObjectTable") 

224 if (not config.doWriteAssociatedSources) or (not config.doSolarSystemAssociation): 

225 self.outputs.remove("associatedSsSources") 

226 self.outputs.remove("unassociatedSsObjects") 

227 

228 def adjustQuantum(self, inputs, outputs, label, dataId): 

229 """Override to make adjustments to `lsst.daf.butler.DatasetRef` objects 

230 in the `lsst.daf.butler.core.Quantum` during the graph generation stage 

231 of the activator. 

232 

233 This implementation checks to make sure that the filters in the dataset 

234 are compatible with AP processing as set by the Apdb/DPDD schema. 

235 

236 Parameters 

237 ---------- 

238 inputs : `dict` 

239 Dictionary whose keys are an input (regular or prerequisite) 

240 connection name and whose values are a tuple of the connection 

241 instance and a collection of associated `DatasetRef` objects. 

242 The exact type of the nested collections is unspecified; it can be 

243 assumed to be multi-pass iterable and support `len` and ``in``, but 

244 it should not be mutated in place. In contrast, the outer 

245 dictionaries are guaranteed to be temporary copies that are true 

246 `dict` instances, and hence may be modified and even returned; this 

247 is especially useful for delegating to `super` (see notes below). 

248 outputs : `dict` 

249 Dict of output datasets, with the same structure as ``inputs``. 

250 label : `str` 

251 Label for this task in the pipeline (should be used in all 

252 diagnostic messages). 

253 data_id : `lsst.daf.butler.DataCoordinate` 

254 Data ID for this quantum in the pipeline (should be used in all 

255 diagnostic messages). 

256 

257 Returns 

258 ------- 

259 adjusted_inputs : `dict` 

260 Dict of the same form as ``inputs`` with updated containers of 

261 input `DatasetRef` objects. Connections that are not changed 

262 should not be returned at all. Datasets may only be removed, not 

263 added. Nested collections may be of any multi-pass iterable type, 

264 and the order of iteration will set the order of iteration within 

265 `PipelineTask.runQuantum`. 

266 adjusted_outputs : `dict` 

267 Dict of updated output datasets, with the same structure and 

268 interpretation as ``adjusted_inputs``. 

269 

270 Raises 

271 ------ 

272 ScalarError 

273 Raised if any `Input` or `PrerequisiteInput` connection has 

274 ``multiple`` set to `False`, but multiple datasets. 

275 NoWorkFound 

276 Raised to indicate that this quantum should not be run; not enough 

277 datasets were found for a regular `Input` connection, and the 

278 quantum should be pruned or skipped. 

279 FileNotFoundError 

280 Raised to cause QuantumGraph generation to fail (with the message 

281 included in this exception); not enough datasets were found for a 

282 `PrerequisiteInput` connection. 

283 """ 

284 _, refs = inputs["diffIm"] 

285 for ref in refs: 

286 if ref.dataId["band"] not in self.config.validBands: 

287 raise ValueError( 

288 f"Requested '{ref.dataId['band']}' not in " 

289 "DiaPipelineConfig.validBands. To process bands not in " 

290 "the standard Rubin set (ugrizy) you must add the band to " 

291 "the validBands list in DiaPipelineConfig and add the " 

292 "appropriate columns to the Apdb schema.") 

293 return super().adjustQuantum(inputs, outputs, label, dataId) 

294 

295 

296class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

297 pipelineConnections=DiaPipelineConnections): 

298 """Config for DiaPipelineTask. 

299 """ 

300 coaddName = pexConfig.Field( 

301 doc="coadd name: typically one of deep, goodSeeing, or dcr", 

302 dtype=str, 

303 default="deep", 

304 ) 

305 apdb_config_url = pexConfig.Field( 

306 dtype=str, 

307 default=None, 

308 optional=False, 

309 doc="A config file specifying the APDB and its connection parameters, " 

310 "typically written by the apdb-cli command-line utility. " 

311 "The database must already be initialized.", 

312 ) 

313 validBands = pexConfig.ListField( 

314 dtype=str, 

315 default=["u", "g", "r", "i", "z", "y"], 

316 doc="List of bands that are valid for AP processing. To process a " 

317 "band not on this list, the appropriate band specific columns " 

318 "must be added to the Apdb schema in dax_apdb.", 

319 ) 

320 associator = pexConfig.ConfigurableField( 

321 target=AssociationTask, 

322 doc="Task used to associate DiaSources with DiaObjects.", 

323 ) 

324 doSolarSystemAssociation = pexConfig.Field( 

325 dtype=bool, 

326 default=True, 

327 doc="Process SolarSystem objects through the pipeline.", 

328 ) 

329 solarSystemAssociator = pexConfig.ConfigurableField( 

330 target=SolarSystemAssociationTask, 

331 doc="Task used to associate DiaSources with SolarSystemObjects.", 

332 ) 

333 diaCalculation = pexConfig.ConfigurableField( 

334 target=DiaObjectCalculationTask, 

335 doc="Task to compute summary statistics for DiaObjects.", 

336 ) 

337 doReloadDiaObjects = pexConfig.Field( 

338 dtype=bool, 

339 default=True, 

340 doc="Drop preloaded DiaObjects and reload them from the APDB?" 

341 "Used in production when the very latest objects from the APDB " 

342 "are needed.", 

343 ) 

344 angleMargin = pexConfig.RangeField( 

345 doc="Padding to add when loading diaObjects if `doReloadDiaObjects=True`", 

346 dtype=float, 

347 default=2, 

348 min=0, 

349 ) 

350 doRunForcedMeasurement = pexConfig.Field( 

351 dtype=bool, 

352 default=True, 

353 deprecated="Added to allow disabling forced sources for performance " 

354 "reasons during the ops rehearsal. " 

355 "It is expected to be removed.", 

356 doc="Run forced measurement on all of the diaObjects? " 

357 "This should only be turned off for debugging purposes.", 

358 ) 

359 diaForcedSource = pexConfig.ConfigurableField( 

360 target=DiaForcedSourceTask, 

361 doc="Task used for force photometer DiaObject locations in direct and " 

362 "difference images.", 

363 ) 

364 forcedReliabilityThreshold = pexConfig.Field( 

365 dtype=float, 

366 default=0.5, 

367 doc="Minimum reliability score of constituent diaSources to use when " 

368 "determining whether to calculate forced photometry for a " 

369 "diaObject.", 

370 ) 

371 forcedTrailLengthThreshold = pexConfig.Field( 

372 dtype=float, 

373 default=1.4, 

374 doc="Maximum trail length (arseconds) of constituent diaSources to use " 

375 "when determining whether to calculate forced photometry for a " 

376 "diaObject.", 

377 ) 

378 forcedBadFlags = pexConfig.ListField( 

379 dtype=str, 

380 default=["glint_trail", "isDipole"], 

381 doc="Flags of constituent diaSources to exclude when determining " 

382 "whether to calculate forced photometry for a diaObject.", 

383 ) 

384 alertPackager = pexConfig.ConfigurableField( 

385 target=PackageAlertsTask, 

386 doc="Subtask for packaging Ap data into alerts.", 

387 ) 

388 doPackageAlerts = pexConfig.Field( 

389 dtype=bool, 

390 default=False, 

391 doc="Package Dia-data into serialized alerts for distribution and " 

392 "write them to disk.", 

393 ) 

394 doWriteAssociatedSources = pexConfig.Field( 

395 dtype=bool, 

396 default=True, 

397 doc="Write out associated DiaSources, DiaForcedSources, and DiaObjects, " 

398 "formatted following the Science Data Model.", 

399 ) 

400 imagePixelMargin = pexConfig.RangeField( 

401 dtype=int, 

402 default=10, 

403 min=0, 

404 doc="Pad the image by this many pixels before removing off-image " 

405 "diaObjects for association.", 

406 ) 

407 filterUnAssociatedSources = pexConfig.Field( 

408 dtype=bool, 

409 default=True, 

410 doc="Check unassociated diaSources for quality before creating new ." 

411 "diaObjects. DiaSources that are associated with existing diaObjects " 

412 "will not be affected." 

413 ) 

414 newObjectBadFlags = pexConfig.ListField( 

415 dtype=str, 

416 default=("centroid_flag", 

417 "pixelFlags_crCenter", 

418 "pixelFlags_nodataCenter", 

419 "pixelFlags_interpolatedCenter", 

420 "pixelFlags_saturatedCenter", 

421 "pixelFlags_suspectCenter", 

422 "pixelFlags_streakCenter", 

423 "glint_trail"), 

424 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

425 "these flags set before creating new diaObjects.", 

426 ) 

427 maxNewDiaObjects = pexConfig.RangeField( 

428 dtype=float, 

429 default=0, 

430 min=0, 

431 doc="Maximum number of new DiaObjects to create before raising an error." 

432 "Set to zero to disable.", 

433 ) 

434 newObjectSnrThreshold = pexConfig.RangeField( 

435 dtype=float, 

436 default=0, 

437 min=0, 

438 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

439 "Abs(flux/fluxErr) less than this threshold before creating new " 

440 "diaObjects." 

441 "Set to zero to disable.", 

442 ) 

443 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

444 dtype=float, 

445 default=0, 

446 min=0, 

447 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

448 "signal-to-noise ratios less than this threshold if they have" 

449 " low reliability scores before creating new diaObjects." 

450 "Set to zero to disable.", 

451 ) 

452 newObjectReliabilityThreshold = pexConfig.RangeField( 

453 dtype=float, 

454 default=0.1, 

455 min=0, 

456 max=1, 

457 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

458 "reliability scores less than this threshold before creating new " 

459 "diaObjects." 

460 "Set to zero to disable.", 

461 ) 

462 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

463 dtype=float, 

464 default=0.1, 

465 min=0, 

466 max=1, 

467 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

468 "low signal-to-noise and reliability scores below this threshold " 

469 "before creating new diaObjects." 

470 "Set to zero to disable.", 

471 ) 

472 newObjectFluxField = pexConfig.Field( 

473 dtype=str, 

474 default="apFlux", 

475 doc="Name of the flux field in the standardized diaSource catalog to " 

476 "use for checking the signal-to-noise before creating new diaObjects.", 

477 ) 

478 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

479 

480 def setDefaults(self): 

481 self.diaCalculation.plugins = ["ap_meanPosition", 

482 "ap_nDiaSources", 

483 "ap_meanFlux", 

484 "ap_sigmaFlux", 

485 "ap_minMaxFlux", 

486 "ap_maxSlopeFlux", 

487 "ap_meanErrFlux", 

488 "ap_meanTotFlux"] 

489 

490 

491class DiaPipelineTask(pipeBase.PipelineTask): 

492 """Task for loading, associating and storing Difference Image Analysis 

493 (DIA) Objects and Sources. 

494 """ 

495 ConfigClass = DiaPipelineConfig 

496 _DefaultName = "diaPipe" 

497 

498 def __init__(self, initInputs=None, **kwargs): 

499 super().__init__(**kwargs) 

500 self.apdb = daxApdb.Apdb.from_uri(self.config.apdb_config_url) 

501 self.schema = readSchemaFromApdb(self.apdb) 

502 self.makeSubtask("associator") 

503 self.makeSubtask("diaCalculation") 

504 if self.config.doRunForcedMeasurement: 

505 self.makeSubtask("diaForcedSource") 

506 if self.config.doPackageAlerts: 

507 self.makeSubtask("alertPackager") 

508 if self.config.doSolarSystemAssociation: 

509 self.makeSubtask("solarSystemAssociator") 

510 if self.config.filterUnAssociatedSources: 

511 columns = [self.config.newObjectFluxField, 

512 self.config.newObjectFluxField + "Err", 

513 "reliability"] 

514 columns += self.config.newObjectBadFlags 

515 

516 missing = checkSdmSchemaColumns(self.schema, columns, "DiaSource") 

517 if missing: 

518 raise pipeBase.InvalidQuantumError("Field %s not in the DiaSource schema" % missing) 

519 

520 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

521 inputs = butlerQC.get(inputRefs) 

522 inputs["idGenerator"] = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

523 inputs["band"] = butlerQC.quantum.dataId["band"] 

524 inputs["legacySolarSystemTable"] = None 

525 if not self.config.doSolarSystemAssociation: 

526 inputs["solarSystemObjectTable"] = None 

527 

528 associationResults = pipeBase.Struct( 

529 apdbMarker=None, 

530 associatedDiaSources=None, 

531 diaForcedSources=None, 

532 diaObjects=None, 

533 associatedSsSources=None, 

534 unassociatedSsObjects=None, 

535 newDiaSources=None, 

536 marginalDiaSources=None, 

537 ) 

538 try: 

539 self.run(**inputs, associationResults=associationResults) 

540 except pipeBase.AlgorithmError as e: 

541 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

542 e, 

543 self, 

544 log=self.log 

545 ) 

546 butlerQC.put(associationResults, outputRefs) 

547 raise error from e 

548 butlerQC.put(associationResults, outputRefs) 

549 

550 @timeMethod 

551 def run(self, 

552 diaSourceTable, 

553 legacySolarSystemTable, 

554 diffIm, 

555 exposure, 

556 template, 

557 preloadedDiaObjects, 

558 preloadedDiaSources, 

559 preloadedDiaForcedSources, 

560 band, 

561 idGenerator, 

562 solarSystemObjectTable=None, 

563 associationResults=None): 

564 """Process DiaSources and DiaObjects. 

565 

566 Load previous DiaObjects and their DiaSource history. Calibrate the 

567 values in the diaSourceCat. Associate new DiaSources with previous 

568 DiaObjects. Run forced photometry at the updated DiaObject locations. 

569 Store the results in the Alert Production Database (Apdb). 

570 

571 Parameters 

572 ---------- 

573 diaSourceTable : `pandas.DataFrame` 

574 Newly detected DiaSources. 

575 legacySolarSystemTable : `pandas.DataFrame` 

576 Not used 

577 diffIm : `lsst.afw.image.ExposureF` 

578 Difference image exposure in which the sources in ``diaSourceCat`` 

579 were detected. 

580 exposure : `lsst.afw.image.ExposureF` 

581 Calibrated exposure differenced with a template to create 

582 ``diffIm``. 

583 template : `lsst.afw.image.ExposureF` 

584 Template exposure used to create diffIm. 

585 preloadedDiaObjects : `pandas.DataFrame` 

586 Previously detected DiaObjects, loaded from the APDB. 

587 preloadedDiaSources : `pandas.DataFrame` 

588 Previously detected DiaSources, loaded from the APDB. 

589 preloadedDiaForcedSources : `pandas.DataFrame` 

590 Catalog of previously detected forced DiaSources, from the APDB 

591 band : `str` 

592 The band in which the new DiaSources were detected. 

593 idGenerator : `lsst.meas.base.IdGenerator` 

594 Object that generates source IDs and random number generator seeds. 

595 solarSystemObjectTable : `astropy.table.Table` 

596 Preloaded Solar System objects expected to be visible in the image. 

597 associationResults : `lsst.pipe.base.Struct`, optional 

598 Result struct that is modified to allow saving of partial outputs 

599 for some failure conditions. If the task completes successfully, 

600 this is also returned. 

601 

602 Returns 

603 ------- 

604 associationResults : `lsst.pipe.base.Struct` 

605 Results struct with components. 

606 

607 - ``apdbMarker`` : Marker dataset to store in the Butler indicating 

608 that this ccdVisit has completed successfully. 

609 (`lsst.dax.apdb.ApdbConfig`) 

610 - ``associatedDiaSources`` : Catalog of newly associated 

611 DiaSources. (`pandas.DataFrame`) 

612 - ``diaForcedSources`` : Catalog of new and previously detected 

613 forced DiaSources. (`pandas.DataFrame`) 

614 - ``diaObjects`` : Updated table of DiaObjects. (`pandas.DataFrame`) 

615 - ``associatedSsSources`` : Catalog of ssSource records. 

616 (`pandas.DataFrame`) 

617 

618 Raises 

619 ------ 

620 RuntimeError 

621 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

622 """ 

623 

624 if associationResults is None: 

625 associationResults = pipeBase.Struct() 

626 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

627 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 

628 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

629 if self.config.doReloadDiaObjects: 

630 try: 

631 preloadedDiaObjects = self.loadRefreshedDiaObjects(getRegion(exposure), preloadedDiaObjects) 

632 except Exception as e: 

633 self.log.warning("Error encountered while attempting to load " 

634 "the latest diaObjects from the APDB. Processing " 

635 "will continue with only the diaObjects from " 

636 "preload.", exc_info=e) 

637 finally: 

638 self.metadata["loadDiaObjectsDuration"] = duration_from_timeMethod( 

639 self.metadata, "loadRefreshedDiaObjects", clock="Utc" 

640 ) 

641 self.log.verbose("Re-loading DiaObjects: Took %.4f seconds", 

642 self.metadata["loadDiaObjectsDuration"]) 

643 

644 else: 

645 self.metadata["loadDiaObjectsDuration"] = -1 

646 

647 self.checkTableIndex(preloadedDiaSources, index=["diaObjectId", "band", "diaSourceId"]) 

648 self.checkTableIndex(preloadedDiaObjects, index="diaObjectId") 

649 self.checkTableIndex(preloadedDiaForcedSources, index=["diaObjectId", "diaForcedSourceId"]) 

650 

651 if not preloadedDiaObjects.empty: 

652 # Include a small buffer outside the image so that we can associate sources near the edge 

653 diaObjects, _ = self.purgeDiaObjects(diffIm.getBBox(), diffIm.getWcs(), preloadedDiaObjects, 

654 buffer=self.config.imagePixelMargin) 

655 else: 

656 self.log.info("Preloaded DiaObject table is empty.") 

657 diaObjects = preloadedDiaObjects 

658 

659 # Associate DiaSources with DiaObjects 

660 assocResults = self.associateDiaSources(diaSourceTable, solarSystemObjectTable, diffIm, diaObjects) 

661 

662 # Set unassociated diaObjectIds and ssObjectIds to NULL, and convert to SDM schema format 

663 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

664 assocResults.associatedDiaSources, 

665 "DiaSource", 

666 nullColumns=["diaObjectId", "ssObjectId"] 

667 ) 

668 

669 # Merge associated diaSources 

670 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

671 preloadedDiaSources, 

672 assocResults.associatedDiaSources, 

673 diaObjects, 

674 assocResults.newDiaObjects, 

675 diffIm 

676 ) 

677 

678 # Compute DiaObject Summary statistics from their full DiaSource 

679 # history. 

680 diaCalResult = self.diaCalculation.run( 

681 mergedDiaObjects, 

682 mergedDiaSourceHistory, 

683 updatedDiaObjectIds, 

684 [band]) 

685 updatedDiaObjects = convertDataFrameToSdmSchema(self.schema, diaCalResult.updatedDiaObjects, 

686 tableName="DiaObject", skipIndex=True) 

687 

688 # Test for duplication in the updated DiaObjects. 

689 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 

690 raise RuntimeError( 

691 "Duplicate DiaObjects (loaded + updated) created after " 

692 "DiaCalculation. This is unexpected behavior and should be " 

693 "reported. Exiting.") 

694 if self.testDataFrameIndex(updatedDiaObjects): 

695 raise RuntimeError( 

696 "Duplicate DiaObjects (updated) created after " 

697 "DiaCalculation. This is unexpected behavior and should be " 

698 "reported. Exiting.") 

699 

700 # Forced source measurement 

701 if self.config.doRunForcedMeasurement: 

702 diaObjectsForced = self._selectGoodDiaObjects(diaCalResult.diaObjectCat, mergedDiaSourceHistory) 

703 diaForcedSources = self.runForcedMeasurement( 

704 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

705 ) 

706 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

707 else: 

708 # alertPackager needs correct columns 

709 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

710 forcedSourceHistoryThreshold = 0 

711 

712 # Write results to Alert Production Database (APDB) 

713 try: 

714 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

715 diaForcedSources) 

716 finally: 

717 self.metadata["writeToApdbDuration"] = duration_from_timeMethod(self.metadata, "writeToApdb", 

718 clock="Utc") 

719 # A single log message is easier for Loki to parse than timeMethod's start+end pairs. 

720 self.log.verbose("writeToApdb: Took %.4f seconds", self.metadata["writeToApdbDuration"]) 

721 

722 # For historical reasons, apdbMarker is a Config even if it's not meant to be read. 

723 # A default Config is the cheapest way to satisfy the storage class. 

724 associationResults.apdbMarker = pexConfig.Config() 

725 associationResults.associatedDiaSources = assocResults.associatedDiaSources 

726 associationResults.diaForcedSources = diaForcedSources 

727 associationResults.diaObjects = diaCalResult.diaObjectCat 

728 associationResults.unassociatedSsObjects = assocResults.unassociatedSsObjects 

729 associationResults.newDiaSources = assocResults.newDiaSources 

730 associationResults.marginalDiaSources = assocResults.marginalDiaSources 

731 

732 # Catch *any* error after we have updated the APDB 

733 try: 

734 associatedSsSources = assocResults.associatedSsSources 

735 associatedSsSourcesPlusMpcorb = None 

736 if associatedSsSources is not None: 

737 mpcorbColumns = [col for col in associatedSsSources.columns if col[:7] == 'MPCORB_'] 

738 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

739 associatedSsSources = self.standardizeTable(associatedSsSources, "SSSource", nullColumns=[]) 

740 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

741 for mpcorbColumn in mpcorbColumns: 

742 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

743 associationResults.associatedSsSources = associatedSsSources 

744 

745 # patch the otherwise-empty validityStart field for the alerts 

746 updatedDiaObjects['validityStartMjdTai'] = validityStart.get(system=DateTime.MJD, 

747 scale=DateTime.TAI) 

748 

749 # Package alerts 

750 if self.config.doPackageAlerts: 

751 # Append new forced sources to the full history 

752 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

753 tableName="DiaForcedSource") 

754 if self.testDataFrameIndex(diaForcedSources): 

755 self.log.warning( 

756 "Duplicate DiaForcedSources created after merge with " 

757 "history and new sources. This may cause downstream " 

758 "problems. Dropping duplicates.") 

759 # Drop duplicates via index and keep the first appearance. 

760 # Reset due to the index shape being slight different than 

761 # expected. 

762 diaForcedSourcesFull = diaForcedSourcesFull.groupby( 

763 diaForcedSourcesFull.index).first() 

764 diaForcedSourcesFull.reset_index(drop=True, inplace=True) 

765 diaForcedSourcesFull.set_index( 

766 ["diaObjectId", "diaForcedSourceId"], 

767 drop=False, 

768 inplace=True) 

769 

770 self.alertPackager.run(assocResults.associatedDiaSources, 

771 updatedDiaObjects, 

772 preloadedDiaSources, 

773 diaForcedSourcesFull, 

774 diffIm, 

775 exposure, 

776 template, 

777 ssSrc=associatedSsSourcesPlusMpcorb, 

778 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

779 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

780 ) 

781 except Exception as e: 

782 # Catch *any* error after we have updated the APDB 

783 raise PostApdbUpdateError(errorMsg=repr(e)) 

784 

785 return associationResults 

786 

787 def _selectGoodDiaObjects(self, diaObjects, diaSources): 

788 """ 

789 Extract a subset of diaObjects with multiple associated diaSources that 

790 pass selection requirements. 

791 """ 

792 n_history = self.diaForcedSource.config.historyThreshold 

793 required_columns = ["reliability", "trailLength"] 

794 for col in required_columns: 

795 if col not in diaSources.columns: 

796 raise RuntimeError("Required column '%s' missing from diaSource table." % col) 

797 

798 rel_ok = diaSources["reliability"] >= self.config.forcedReliabilityThreshold 

799 trailed_ok = diaSources["trailLength"] <= self.config.forcedTrailLengthThreshold 

800 badFlags = [f for f in self.config.forcedBadFlags if f in diaSources.columns] 

801 if badFlags: 

802 self.log.info("Excluding diaSources with %s flags set when calculating" 

803 " diaObject history for forced photometry." % badFlags) 

804 # Exclude rows where ANY bad flag is True. 

805 # If bad flag columns can contain NA, treat NA as True (i.e., bad) via fillna(True). 

806 bad_any = diaSources[badFlags].fillna(True).any(axis=1) 

807 mask = rel_ok & trailed_ok & ~bad_any 

808 ids = diaSources.index.get_level_values("diaObjectId")[mask] 

809 

810 # Count occurrences in diaSource 

811 counts = pd.Series(ids, copy=False).value_counts() 

812 

813 # IDs that exceed threshold 

814 keep_ids = counts[counts >= n_history].index 

815 

816 # Since diaObject index is unique, direct index filtering is safe and efficient 

817 return diaObjects.loc[diaObjects.index.intersection(keep_ids)] 

818 

819 def createNewDiaObjects(self, unassociatedDiaSources): 

820 """Loop through the set of DiaSources and create new DiaObjects 

821 for unassociated DiaSources. 

822 

823 Parameters 

824 ---------- 

825 unassociatedDiaSources : `pandas.DataFrame` 

826 Set of DiaSources to create new DiaObjects from. 

827 

828 Returns 

829 ------- 

830 results : `lsst.pipe.base.Struct` 

831 Results struct containing: 

832 

833 - diaSources : `pandas.DataFrame` 

834 DiaSource catalog with updated DiaObject ids. 

835 - newDiaObjects : `pandas.DataFrame` 

836 Newly created DiaObjects from the unassociated DiaSources. 

837 - nNewDiaObjects : `int` 

838 Number of newly created diaObjects. 

839 - marginalDiaSources : `pandas.DataFrame` 

840 Unassociated diaSources with low signal-to-noise and/or 

841 reliability, which are excluded from the new DiaObjects. 

842 """ 

843 marginalDiaSources = make_empty_catalog(self.schema, tableName="DiaSource") 

844 if len(unassociatedDiaSources) == 0: 

845 newDiaObjects = make_empty_catalog(self.schema, tableName="DiaObject") 

846 else: 

847 if self.config.filterUnAssociatedSources: 

848 results = self.filterSources(unassociatedDiaSources) 

849 unassociatedDiaSources = results.goodSources 

850 marginalDiaSources = results.badSources 

851 unassociatedDiaSources["diaObjectId"] = unassociatedDiaSources["diaSourceId"] 

852 newDiaObjects = convertDataFrameToSdmSchema(self.schema, unassociatedDiaSources, 

853 tableName="DiaObject", skipIndex=True) 

854 self.metadata["nRejectedNewDiaObjects"] = len(marginalDiaSources) 

855 return pipeBase.Struct(diaSources=unassociatedDiaSources, 

856 newDiaObjects=newDiaObjects, 

857 nNewDiaObjects=len(newDiaObjects), 

858 marginalDiaSources=marginalDiaSources) 

859 

860 def filterSources(self, sources, snrThreshold=None, lowReliabilitySnrThreshold=None, 

861 reliabilityThreshold=None, lowSnrReliabilityThreshold=None, badFlags=None): 

862 """Select good sources out of a catalog. 

863 

864 Parameters 

865 ---------- 

866 sources : `pandas.DataFrame` 

867 Set of DiaSources to check. 

868 snrThreshold : `float`, optional 

869 The minimum signal to noise diaSource to make a new diaObject. 

870 Included for unit tests. Uses the task config value if not set. 

871 lowReliabilitySnrThreshold : `float`, optional 

872 Use ``lowSnrReliabilityThreshold`` as the reliability threshold for 

873 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

874 Included for unit tests. Uses the task config value if not set. 

875 reliabilityThreshold : `float`, optional 

876 The minimum reliability score diaSource to make a new diaObject 

877 Included for unit tests. Uses the task config value if not set. 

878 lowSnrReliabilityThreshold : `float`, optional 

879 Use ``lowSnrReliabilityThreshold`` as the reliability threshold for 

880 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

881 Included for unit tests. Uses the task config value if not set. 

882 badFlags : `list` of `str`, optional 

883 Do not create new diaObjects for any diaSource with any of these 

884 flags set. 

885 Included for unit tests. Uses the task config value if not set. 

886 

887 Returns 

888 ------- 

889 results : `lsst.pipe.base.Struct` 

890 Results struct containing: 

891 

892 - goodSources : `pandas.DataFrame` 

893 Subset of the input sources that pass all checks. 

894 - badSources : `pandas.DataFrame` 

895 Subset of the input sources that fail any check. 

896 """ 

897 if snrThreshold is None: 

898 snrThreshold = self.config.newObjectSnrThreshold 

899 if lowReliabilitySnrThreshold is None: 

900 lowReliabilitySnrThreshold = self.config.newObjectLowReliabilitySnrThreshold 

901 if reliabilityThreshold is None: 

902 reliabilityThreshold = self.config.newObjectReliabilityThreshold 

903 if lowSnrReliabilityThreshold is None: 

904 lowSnrReliabilityThreshold = self.config.newObjectLowSnrReliabilityThreshold 

905 if badFlags is None: 

906 badFlags = self.config.newObjectBadFlags 

907 flagged = sources[badFlags].fillna(False).any(axis=1) 

908 fluxField = self.config.newObjectFluxField 

909 fluxErrField = fluxField + "Err" 

910 signalToNoise = np.abs(np.array(sources[fluxField]/sources[fluxErrField])) 

911 reliability = np.array(sources['reliability']) 

912 nFlagged = np.count_nonzero(flagged)