Coverage for tests / fgcmcalTestBase.py: 8%

1# See COPYRIGHT file at the top of the source tree. 

2# 

3# This file is part of fgcmcal. 

4# 

5# Developed for the LSST Data Management System. 

6# This product includes software developed by the LSST Project 

7# (https://www.lsst.org). 

8# See the COPYRIGHT file at the top-level directory of this distribution 

9# for details of code ownership. 

10# 

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

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

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

14# (at your option) any later version. 

15# 

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

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

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

19# GNU General Public License for more details. 

20# 

21# You should have received a copy of the GNU General Public License 

22# along with this program. If not, see <https://www.gnu.org/licenses/>. 

23"""General fgcmcal testing class. 

24 

25This class is used as the basis for individual obs package tests using 

26data from testdata_jointcal for Gen3 repos. 

27""" 

28 

29import os 

30import shutil 

31import numpy as np 

32import esutil 

33import warnings 

34 

35import lsst.afw.image 

36import lsst.daf.butler as dafButler 

37import lsst.pipe.base as pipeBase 

38import lsst.geom as geom 

39from lsst.pipe.base import Pipeline, ExecutionResources 

40from lsst.pipe.base.simple_pipeline_executor import SimplePipelineExecutor 

41 

42import lsst.fgcmcal as fgcmcal 

43 

44ROOT = os.path.abspath(os.path.dirname(__file__)) 

45 

46 

47class FgcmcalTestBase(object): 

48 """Base class for gen3 fgcmcal tests, to genericize some test running and setup. 

49 

50 Derive from this first, then from TestCase. 

51 """ 

52 @classmethod 

53 def _importRepository(cls, instrument, exportPath, exportFile): 

54 """Import a test repository into self.testDir 

55 

56 Parameters 

57 ---------- 

58 instrument : `str` 

59 Full string name for the instrument. 

60 exportPath : `str` 

61 Path to location of repository to export. 

62 exportFile : `str` 

63 Filename of export data. 

64 """ 

65 cls.repo = os.path.join(cls.testDir, 'testrepo') 

66 

67 # Make the repo and retrieve a writeable Butler 

68 _ = dafButler.Butler.makeRepo(cls.repo) 

69 with dafButler.Butler.from_config(cls.repo, writeable=True) as butler: 

70 # Register the instrument 

71 instrInstance = pipeBase.Instrument.from_string(instrument) 

72 instrInstance.register(butler.registry) 

73 # Import the exportFile 

74 butler.import_(directory=exportPath, filename=exportFile, 

75 transfer='symlink', 

76 skip_dimensions={'instrument', 'detector', 'physical_filter'}) 

77 

78 def _runPipeline(self, repo, pipelineFile, queryString='', 

79 inputCollections=None, outputCollection=None, 

80 configFiles={}, configOptions={}, 

81 registerDatasetTypes=False): 

82 """Run a pipeline via pipetask. 

83 

84 Parameters 

85 ---------- 

86 repo : `str` 

87 Gen3 repository yaml file. 

88 pipelineFile : `str` 

89 Pipeline definition file. 

90 queryString : `str`, optional 

91 Where query that defines the data to use. 

92 inputCollections : `list` [`str`], optional 

93 Input collections list. 

94 outputCollection : `str`, optional 

95 Output collection name. 

96 configFiles : `dict` [`list`], optional 

97 Dictionary of config files. The key of the ``configFiles`` 

98 dict is the relevant task label. The value of ``configFiles`` 

99 is a list of config files to apply (in order) to that task. 

100 configOptions : `dict` [`dict`], optional 

101 Dictionary of individual config options. The key of the 

102 ``configOptions`` dict is the relevant task label. The value 

103 of ``configOptions`` is another dict that contains config 

104 key/value overrides to apply. 

105 configOptions : `list` [`str`], optional 

106 List of individual config options to use. Each string will 

107 be of the form ``taskName:configField=value``. 

108 registerDatasetTypes : `bool`, optional 

109 Register new dataset types? 

110 

111 Returns 

112 ------- 

113 exit_code : `int` 

114 Exit code for pipetask run. 

115 

116 Raises 

117 ------ 

118 RuntimeError : Raised if the "pipetask" call fails. 

119 """ 

120 butler = SimplePipelineExecutor.prep_butler(repo, 

121 inputs=inputCollections, 

122 output=outputCollection) 

123 self.enterContext(butler) 

124 pipeline = Pipeline.fromFile(pipelineFile) 

125 for taskName, fileList in configFiles.items(): 

126 for fileName in fileList: 

127 pipeline.addConfigFile(taskName, fileName) 

128 for taskName, configDict in configOptions.items(): 

129 for option, value in configDict.items(): 

130 pipeline.addConfigOverride(taskName, option, value) 

131 

132 resources = ExecutionResources(num_cores=1) 

133 

134 executor = SimplePipelineExecutor.from_pipeline(pipeline, 

135 where=queryString, 

136 butler=butler, 

137 resources=resources) 

138 quanta = executor.run(register_dataset_types=registerDatasetTypes) 

139 

140 return len(quanta) 

141 

142 def _testFgcmMakeLut(self, instName, testName, nBand, i0Std, i0Recon, i10Std, i10Recon): 

143 """Test running of FgcmMakeLutTask 

144 

145 Parameters 

146 ---------- 

147 instName : `str` 

148 Short name of the instrument. 

149 testName : `str` 

150 Base name of the test collection. 

151 nBand : `int` 

152 Number of bands tested. 

153 i0Std : `np.ndarray' 

154 Values of i0Std to compare to. 

155 i10Std : `np.ndarray` 

156 Values of i10Std to compare to. 

157 i0Recon : `np.ndarray` 

158 Values of reconstructed i0 to compare to. 

159 i10Recon : `np.ndarray` 

160 Values of reconsntructed i10 to compare to. 

161 """ 

162 instCamel = instName.title() 

163 

164 configFiles = {'fgcmMakeLut': [os.path.join(ROOT, 

165 'config', 

166 f'fgcmMakeLut{instCamel}.py')]} 

167 outputCollection = f'{instName}/{testName}/lut' 

168 

169 self._runPipeline(self.repo, 

170 os.path.join(ROOT, 

171 'pipelines', 

172 f'fgcmMakeLut{instCamel}.yaml'), 

173 configFiles=configFiles, 

174 inputCollections=[f'{instName}/calib', f'{instName}/testdata'], 

175 outputCollection=outputCollection, 

176 registerDatasetTypes=True) 

177 

178 # Check output values 

179 butler = dafButler.Butler(self.repo) 

180 self.enterContext(butler) 

181 lutCat = butler.get('fgcmLookUpTable', 

182 collections=[outputCollection], 

183 instrument=instName) 

184 fgcmLut, lutIndexVals, lutStd = fgcmcal.utilities.translateFgcmLut(lutCat, {}) 

185 

186 self.assertEqual(nBand, len(lutIndexVals[0]['FILTERNAMES'])) 

187 

188 indices = fgcmLut.getIndices(np.arange(nBand, dtype=np.int32), 

189 np.zeros(nBand) + np.log(lutStd[0]['PWVSTD']), 

190 np.zeros(nBand) + lutStd[0]['O3STD'], 

191 np.zeros(nBand) + np.log(lutStd[0]['TAUSTD']), 

192 np.zeros(nBand) + lutStd[0]['ALPHASTD'], 

193 np.zeros(nBand) + 1./np.cos(np.radians(lutStd[0]['ZENITHSTD'])), 

194 np.zeros(nBand, dtype=np.int32), 

195 np.zeros(nBand) + lutStd[0]['PMBSTD']) 

196 i0 = fgcmLut.computeI0(np.zeros(nBand) + np.log(lutStd[0]['PWVSTD']), 

197 np.zeros(nBand) + lutStd[0]['O3STD'], 

198 np.zeros(nBand) + np.log(lutStd[0]['TAUSTD']), 

199 np.zeros(nBand) + lutStd[0]['ALPHASTD'], 

200 np.zeros(nBand) + 1./np.cos(np.radians(lutStd[0]['ZENITHSTD'])), 

201 np.zeros(nBand) + lutStd[0]['PMBSTD'], 

202 indices) 

203 

204 self.assertFloatsAlmostEqual(i0Recon, i0, msg='i0Recon', rtol=1e-5) 

205 

206 i1 = fgcmLut.computeI1(np.zeros(nBand) + np.log(lutStd[0]['PWVSTD']), 

207 np.zeros(nBand) + lutStd[0]['O3STD'], 

208 np.zeros(nBand) + np.log(lutStd[0]['TAUSTD']), 

209 np.zeros(nBand) + lutStd[0]['ALPHASTD'], 

210 np.zeros(nBand) + 1./np.cos(np.radians(lutStd[0]['ZENITHSTD'])), 

211 np.zeros(nBand) + lutStd[0]['PMBSTD'], 

212 indices) 

213 

214 self.assertFloatsAlmostEqual(i10Recon, i1/i0, msg='i10Recon', rtol=1e-5) 

215 

216 # Check that the standard atmosphere was output and non-zero. 

217 atmStd = butler.get('fgcm_standard_atmosphere', 

218 collections=[outputCollection], 

219 instrument=instName) 

220 bounds = atmStd.getWavelengthBounds() 

221 lambdas = np.linspace(bounds[0], bounds[1], 1000) 

222 tputs = atmStd.sampleAt(position=geom.Point2D(0.0, 0.0), wavelengths=lambdas) 

223 self.assertGreater(np.min(tputs), 0.0) 

224 

225 # Check that the standard passbands were output and non-zero. 

226 for physical_filter in fgcmLut.filterNames: 

227 passband = butler.get('fgcm_standard_passband', 

228 collections=[outputCollection], 

229 instrument=instName, 

230 physical_filter=physical_filter) 

231 tputs = passband.sampleAt(position=geom.Point2D(0.0, 0.0), wavelengths=lambdas) 

232 self.assertEqual(np.min(tputs), 0.0) 

233 self.assertGreater(np.max(tputs), 0.0) 

234 

235 # Check that the standard passband tables were output and non-zero. 

236 refs = butler.query_datasets( 

237 "standard_passband", 

238 collections=[outputCollection], 

239 instrument=instName, 

240 ) 

241 self.assertEqual(len(refs), len(fgcmLut.filterNames)) 

242 for ref in refs: 

243 passbandTable = butler.get(ref) 

244 self.assertEqual(np.min(passbandTable["throughput"]), 0.0) 

245 self.assertGreater(np.max(passbandTable["throughput"]), 0.0) 

246 

247 def _testFgcmBuildStarsTable(self, instName, testName, queryString, visits, nStar, nObs, 

248 refcatCollection="refcats/gen2"): 

249 """Test running of FgcmBuildStarsTableTask 

250 

251 Parameters 

252 ---------- 

253 instName : `str` 

254 Short name of the instrument. 

255 testName : `str` 

256 Base name of the test collection. 

257 queryString : `str` 

258 Query to send to the pipetask. 

259 visits : `list` 

260 List of visits to calibrate. 

261 nStar : `int` 

262 Number of stars expected. 

263 nObs : `int` 

264 Number of observations of stars expected. 

265 refcatCollection : `str`, optional 

266 Name of reference catalog collection. 

267 """ 

268 instCamel = instName.title() 

269 

270 configFiles = {'fgcmBuildStarsTable': [os.path.join(ROOT, 

271 'config', 

272 f'fgcmBuildStarsTable{instCamel}.py')]} 

273 outputCollection = f'{instName}/{testName}/buildstars' 

274 

275 self._runPipeline(self.repo, 

276 os.path.join(ROOT, 

277 'pipelines', 

278 'fgcmBuildStarsTable%s.yaml' % (instCamel)), 

279 configFiles=configFiles, 

280 inputCollections=[f'{instName}/{testName}/lut', 

281 refcatCollection], 

282 outputCollection=outputCollection, 

283 queryString=queryString, 

284 registerDatasetTypes=True) 

285 

286 butler = dafButler.Butler(self.repo) 

287 self.enterContext(butler) 

288 

289 visitCat = butler.get('fgcmVisitCatalog', collections=[outputCollection], 

290 instrument=instName) 

291 self.assertEqual(len(visits), len(visitCat)) 

292 

293 starIds = butler.get('fgcmStarIds', collections=[outputCollection], 

294 instrument=instName) 

295 self.assertEqual(len(starIds), nStar) 

296 

297 starObs = butler.get('fgcmStarObservations', collections=[outputCollection], 

298 instrument=instName) 

299 self.assertEqual(len(starObs), nObs) 

300 

301 def _testFgcmBuildFromIsolatedStars(self, instName, testName, queryString, visits, nStar, nObs, 

302 refcatCollection="refcats/gen2"): 

303 """Test running of FgcmBuildFromIsolatedStarsTask. 

304 

305 Parameters 

306 ---------- 

307 instName : `str` 

308 Short name of the instrument. 

309 testName : `str` 

310 Base name of the test collection. 

311 queryString : `str` 

312 Query to send to the pipetask. 

313 visits : `list` 

314 List of visits to calibrate. 

315 nStar : `int` 

316 Number of stars expected. 

317 nObs : `int` 

318 Number of observations of stars expected. 

319 refcatCollection : `str`, optional 

320 Name of reference catalog collection. 

321 """ 

322 instCamel = instName.title() 

323 

324 configFiles = {'fgcmBuildFromIsolatedStars': [ 

325 os.path.join(ROOT, 

326 'config', 

327 f'fgcmBuildFromIsolatedStars{instCamel}.py') 

328 ]} 

329 outputCollection = f'{instName}/{testName}/buildstars' 

330 

331 self._runPipeline(self.repo, 

332 os.path.join(ROOT, 

333 'pipelines', 

334 'fgcmBuildFromIsolatedStars%s.yaml' % (instCamel)), 

335 configFiles=configFiles, 

336 inputCollections=[f'{instName}/{testName}/lut', 

337 refcatCollection], 

338 outputCollection=outputCollection, 

339 queryString=queryString, 

340 registerDatasetTypes=True) 

341 

342 butler = dafButler.Butler(self.repo) 

343 self.enterContext(butler) 

344 

345 visitCat = butler.get('fgcmVisitCatalog', collections=[outputCollection], 

346 instrument=instName) 

347 self.assertEqual(len(visits), len(visitCat)) 

348 

349 starIds = butler.get('fgcm_star_ids', collections=[outputCollection], 

350 instrument=instName) 

351 self.assertEqual(len(starIds), nStar) 

352 

353 starObs = butler.get('fgcm_star_observations', collections=[outputCollection], 

354 instrument=instName) 

355 self.assertEqual(len(starObs), nObs) 

356 

357 def _testFgcmFitCycle(self, instName, testName, cycleNumber, 

358 nZp, nGoodZp, nOkZp, nBadZp, nStdStars, nPlots, 

359 skipChecks=False, extraConfig=None): 

360 """Test running of FgcmFitCycleTask 

361 

362 Parameters 

363 ---------- 

364 instName : `str` 

365 Short name of the instrument. 

366 testName : `str` 

367 Base name of the test collection. 

368 cycleNumber : `int` 

369 Fit cycle number. 

370 nZp : `int` 

371 Number of zeropoints created by the task. 

372 nGoodZp : `int` 

373 Number of good (photometric) zeropoints created. 

374 nOkZp : `int` 

375 Number of constrained zeropoints (photometric or not). 

376 nBadZp : `int` 

377 Number of unconstrained (bad) zeropoints. 

378 nStdStars : `int` 

379 Number of standard stars produced. 

380 nPlots : `int` 

381 Number of plots produced. 

382 skipChecks : `bool`, optional 

383 Skip number checks, when running less-than-final cycle. 

384 extraConfig : `str`, optional 

385 Name of an extra config file to apply. 

386 """ 

387 instCamel = instName.title() 

388 

389 configFiles = {'fgcmFitCycle': [os.path.join(ROOT, 

390 'config', 

391 f'fgcmFitCycle{instCamel}.py')]} 

392 if extraConfig is not None: 

393 configFiles['fgcmFitCycle'].append(extraConfig) 

394 

395 outputCollection = f'{instName}/{testName}/fit' 

396 

397 if cycleNumber == 0: 

398 inputCollections = [f'{instName}/{testName}/buildstars'] 

399 else: 

400 # In these tests we are running the fit cycle task multiple 

401 # times into the same output collection. This code allows 

402 # us to find the correct chained input collections to use 

403 # so that we can both read from previous runs in the output 

404 # collection and write to a new run in the output collection. 

405 # Note that this behavior is handled automatically by the 

406 # pipetask command-line interface, but not by the python 

407 # API. 

408 butler = dafButler.Butler(self.repo) 

409 self.enterContext(butler) 

410 inputCollections = list(butler.registry.getCollectionChain(outputCollection)) 

411 

412 cwd = os.getcwd() 

413 runDir = os.path.join(self.testDir, testName) 

414 os.makedirs(runDir, exist_ok=True) 

415 os.chdir(runDir) 

416 

417 configOptions = {'fgcmFitCycle': 

418 {'cycleNumber': f'{cycleNumber}', 

419 'connections.previousCycleNumber': f'{cycleNumber - 1}', 

420 'connections.cycleNumber': f'{cycleNumber}'}} 

421 self._runPipeline(self.repo, 

422 os.path.join(ROOT, 

423 'pipelines', 

424 f'fgcmFitCycle{instCamel}.yaml'), 

425 configFiles=configFiles, 

426 inputCollections=inputCollections, 

427 outputCollection=outputCollection, 

428 configOptions=configOptions, 

429 registerDatasetTypes=True) 

430 

431 os.chdir(cwd) 

432 

433 if skipChecks: 

434 return 

435 

436 butler = dafButler.Butler(self.repo) 

437 self.enterContext(butler) 

438 

439 # Check that the expected number of plots are there. 

440 plotDatasets = list(butler.registry.queryDatasets( 

441 f"fgcm_Cycle{cycleNumber}_*Plot", 

442 collections=[outputCollection], 

443 )) 

444 self.assertEqual(len(plotDatasets), nPlots) 

445 

446 zps = butler.get(f'fgcm_Cycle{cycleNumber}_Zeropoints', 

447 collections=[outputCollection], 

448 instrument=instName) 

449 self.assertEqual(len(zps), nZp) 

450 

451 gd, = np.where(zps['fgcmFlag'] == 1) 

452 self.assertEqual(len(gd), nGoodZp) 

453 

454 ok, = np.where(zps['fgcmFlag'] < 16) 

455 self.assertEqual(len(ok), nOkZp) 

456 

457 bd, = np.where(zps['fgcmFlag'] >= 16) 

458 self.assertEqual(len(bd), nBadZp) 

459 

460 # Check that there are no illegal values with the ok zeropoints 

461 test, = np.where(zps['fgcmZpt'][gd] < -9000.0) 

462 self.assertEqual(len(test), 0) 

463 

464 stds = butler.get(f'fgcm_Cycle{cycleNumber}_StandardStars', 

465 collections=[outputCollection], 

466 instrument=instName) 

467 

468 self.assertEqual(len(stds), nStdStars) 

469 

470 def _testFgcmOutputProducts(self, instName, testName, 

471 zpOffsets, testVisit, testCcd, testFilter, testBandIndex, 

472 skymapName=None, testSrc=True): 

473 """Test running of FgcmOutputProductsTask. 

474 

475 Parameters 

476 ---------- 

477 instName : `str` 

478 Short name of the instrument. 

479 testName : `str` 

480 Base name of the test collection. 

481 zpOffsets : `np.ndarray` 

482 Zeropoint offsets expected. 

483 testVisit : `int` 

484 Visit id to check for round-trip computations. 

485 testCcd : `int` 

486 Ccd id to check for round-trip computations. 

487 testFilter : `str` 

488 Filtername for testVisit/testCcd. 

489 testBandIndex : `int` 

490 Band index for testVisit/testCcd. 

491 skymapName : `str`, optional 

492 Name of a skymap to use if sharding by tract. 

493 testSrc : `bool`, optional 

494 Test the source catalogs? (Only if available in dataset.) 

495 """ 

496 instCamel = instName.title() 

497 

498 configFiles = {'fgcmOutputProducts': [os.path.join(ROOT, 

499 'config', 

500 f'fgcmOutputProducts{instCamel}.py')]} 

501 inputCollection = f'{instName}/{testName}/fit' 

502 outputCollection = f'{instName}/{testName}/fit/output' 

503 if skymapName is not None: 

504 queryString = f"skymap='{skymapName:s}'" 

505 else: 

506 queryString = "" 

507 

508 self._runPipeline(self.repo, 

509 os.path.join(ROOT, 

510 'pipelines', 

511 'fgcmOutputProducts%s.yaml' % (instCamel)), 

512 queryString=queryString, 

513 configFiles=configFiles, 

514 inputCollections=[inputCollection, 'skymaps'], 

515 outputCollection=outputCollection, 

516 registerDatasetTypes=True) 

517 

518 butler = dafButler.Butler(self.repo, collections=[outputCollection], instrument=instName) 

519 self.enterContext(butler) 

520 offsetCat = butler.get('fgcmReferenceCalibrationOffsets') 

521 offsets = offsetCat['offset'][:] 

522 self.assertFloatsAlmostEqual(offsets, zpOffsets, atol=1e-6) 

523 

524 config = butler.get('fgcmOutputProducts_config', 

525 collections=[outputCollection], instrument=instName) 

526 

527 rawStars = butler.get(f'fgcm_Cycle{config.connections.cycleNumber}_StandardStars', 

528 collections=[inputCollection]) 

529 

530 # Test the fgcm_photoCalib output 

531 zptCat = butler.get(f'fgcm_Cycle{config.connections.cycleNumber}_Zeropoints') 

532 

533 good = (zptCat['fgcmFlag'] < 16) 

534 bad = (zptCat['fgcmFlag'] >= 16) 

535 

536 # Read in all the calibrations, these should all be there 

537 # This test is simply to ensure that all the photoCalib files exist 

538 visits = np.unique(zptCat['visit'][good]) 

539 photoCalibDict = {} 

540 for visit in visits: 

541 expCat = butler.get('fgcmPhotoCalibCatalog', 

542 visit=visit) 

543 for row in expCat: 

544 if row['visit'] == visit: 

545 photoCalibDict[(visit, row['id'])] = row.getPhotoCalib() 

546 

547 # Check that all of the good photocalibs are there. 

548 for rec in zptCat[good]: 

549 self.assertTrue((rec['visit'], rec['detector']) in photoCalibDict) 

550 

551 # Check that none of the bad photocalibs are there. 

552 for rec in zptCat[bad]: 

553 self.assertFalse((rec['visit'], rec['detector']) in photoCalibDict) 

554 

555 # We do round-trip value checking on just the final one (chosen arbitrarily) 

556 testCal = photoCalibDict[(testVisit, testCcd)] 

557 

558 # Check the output standard star catalogs. 

559 if skymapName is not None: 

560 refs = butler.query_datasets("fgcm_standard_star") 

561 self.assertEqual(len(refs), 1) 

562 catalog = butler.get(refs[0]) 

563 self.assertEqual(len(catalog), len(rawStars)) 

564 

565 # Check that the fgcm_id array is what we expect. 

566 np.testing.assert_array_equal(catalog["fgcm_id"], rawStars["id"]) 

567 

568 # Check ids against inputs. 

569 isoConfig = butler.get("fgcmBuildFromIsolatedStars_config") 

570 isolatedCatalog = butler.get( 

571 isoConfig.connections.isolated_star_cats, 

572 tract=refs[0].dataId["tract"], 

573 ) 

574 

575 a, b = esutil.numpy_util.match( 

576 np.asarray(catalog["isolated_star_id"]), 

577 np.asarray(isolatedCatalog["isolated_star_id"]), 

578 ) 

579 

580 # All of the stars should be matched, and we can confirm 

581 # that the RA values are equal as an additional check. 

582 self.assertEqual(len(a), len(catalog)) 

583 np.testing.assert_array_almost_equal(catalog["ra"][a], isolatedCatalog["ra"][b]) 

584 

585 if testSrc: 

586 src = butler.get('src', visit=int(testVisit), detector=int(testCcd), 

587 collections=[outputCollection], instrument=instName) 

588 

589 # Only test sources with positive flux 

590 gdSrc = (src['slot_CalibFlux_instFlux'] > 0.0) 

591 

592 # We need to apply the calibration offset to the fgcmzpt (which is internal 

593 # and doesn't know about that yet) 

594 testZpInd, = np.where((zptCat['visit'] == testVisit) 

595 & (zptCat['detector'] == testCcd)) 

596 fgcmZpt = (zptCat['fgcmZpt'][testZpInd] + offsets[testBandIndex] 

597 + zptCat['fgcmDeltaChrom'][testZpInd]) 

598 fgcmZptGrayErr = np.sqrt(zptCat['fgcmZptVar'][testZpInd]) 

599 

600 if config.doComposeWcsJacobian: 

601 # The raw zeropoint needs to be modified to know about the wcs jacobian 

602 refs = butler.registry.queryDatasets('camera', dimensions=['instrument'], 

603 collections=...) 

604 camera = butler.get(list(refs)[0]) 

605 approxPixelAreaFields = fgcmcal.utilities.computeApproxPixelAreaFields(camera) 

606 center = approxPixelAreaFields[testCcd].getBBox().getCenter() 

607 pixAreaCorr = approxPixelAreaFields[testCcd].evaluate(center) 

608 fgcmZpt += -2.5*np.log10(pixAreaCorr) 

609 

610 # This is the magnitude through the mean calibration 

611 photoCalMeanCalMags = np.zeros(gdSrc.sum()) 

612 # This is the magnitude through the full focal-plane variable mags 

613 photoCalMags = np.zeros_like(photoCalMeanCalMags) 

614 # This is the magnitude with the FGCM (central-ccd) zeropoint 

615 zptMeanCalMags = np.zeros_like(photoCalMeanCalMags) 

616 

617 for i, rec in enumerate(src[gdSrc]): 

618 photoCalMeanCalMags[i] = testCal.instFluxToMagnitude(rec['slot_CalibFlux_instFlux']) 

619 photoCalMags[i] = testCal.instFluxToMagnitude(rec['slot_CalibFlux_instFlux'], 

620 rec.getCentroid()) 

621 zptMeanCalMags[i] = fgcmZpt[0] - 2.5*np.log10(rec['slot_CalibFlux_instFlux']) 

622 

623 # These should be very close but some tiny differences because the fgcm value 

624 # is defined at the center of the bbox, and the photoCal is the mean over the box 

625 self.assertFloatsAlmostEqual(photoCalMeanCalMags, 

626 zptMeanCalMags, rtol=1e-6) 

627 # These should be roughly equal, but not precisely because of the focal-plane 

628 # variation. However, this is a useful sanity check for something going totally 

629 # wrong. 

630 self.assertFloatsAlmostEqual(photoCalMeanCalMags, 

631 photoCalMags, rtol=1e-2) 

632 

633 # The next test compares the "FGCM standard magnitudes" (which are output 

634 # from the fgcm code itself) to the "calibrated magnitudes" that are 

635 # obtained from running photoCalib.calibrateCatalog() on the original 

636 # src catalogs. This summary comparison ensures that using photoCalibs 

637 # yields the same results as what FGCM is computing internally. 

638 # Note that we additionally need to take into account the post-processing 

639 # offsets used in the tests. 

640 

641 # For decent statistics, we are matching all the sources from one visit 

642 # (multiple ccds) 

643 whereClause = f"instrument='{instName:s}' and visit={testVisit:d}" 

644 srcRefs = butler.registry.queryDatasets('src', dimensions=['visit'], 

645 collections='%s/testdata' % (instName), 

646 where=whereClause, 

647 findFirst=True) 

648 photoCals = [] 

649 for srcRef in srcRefs: 

650 photoCals.append(photoCalibDict[(testVisit, srcRef.dataId['detector'])]) 

651 

652 matchMag, matchDelta = self._getMatchedVisitCat(butler, srcRefs, photoCals, 

653 rawStars, testBandIndex, offsets) 

654 

655 st = np.argsort(matchMag) 

656 # Compare the brightest 25% of stars. No matter the setting of 

657 # deltaMagBkgOffsetPercentile, we want to ensure that these stars 

658 # match on average. 

659 brightest, = np.where(matchMag < matchMag[st[int(0.25*st.size)]]) 

660 self.assertFloatsAlmostEqual(np.median(matchDelta[brightest]), 0.0, atol=0.002) 

661 

662 # And the photoCal error is just the zeropoint gray error 

663 self.assertFloatsAlmostEqual(testCal.getCalibrationErr(), 

664 (np.log(10.0)/2.5)*testCal.getCalibrationMean()*fgcmZptGrayErr) 

665 

666 # Test the transmission output 

667 visitCatalog = butler.get('fgcmVisitCatalog', collections=[inputCollection], 

668 instrument=instName) 

669 lutCat = butler.get('fgcmLookUpTable', collections=[inputCollection], 

670 instrument=instName) 

671 

672 testTrans = butler.get('transmission_atmosphere_fgcm', 

673 visit=visitCatalog[0]['visit'], 

674 collections=[outputCollection], instrument=instName) 

675 testResp = testTrans.sampleAt(position=geom.Point2D(0, 0), 

676 wavelengths=lutCat[0]['atmLambda']) 

677 

678 # The test fit is performed with the atmosphere parameters frozen 

679 # (freezeStdAtmosphere = True). Thus the only difference between 

680 # these output atmospheres and the standard is the different 

681 # airmass. Furthermore, this is a very rough comparison because 

682 # the look-up table is computed with very coarse sampling for faster 

683 # testing. 

684 

685 # To account for overall throughput changes, we scale by the median ratio, 

686 # we only care about the shape 

687 ratio = np.median(testResp/lutCat[0]['atmStdTrans']) 

688 self.assertFloatsAlmostEqual(testResp/ratio, lutCat[0]['atmStdTrans'], atol=0.2) 

689 

690 # The second should be close to the first, but there is the airmass 

691 # difference so they aren't identical. 

692 testTrans2 = butler.get('transmission_atmosphere_fgcm', 

693 visit=visitCatalog[1]['visit'], 

694 collections=[outputCollection], instrument=instName) 

695 testResp2 = testTrans2.sampleAt(position=geom.Point2D(0, 0), 

696 wavelengths=lutCat[0]['atmLambda']) 

697 

698 # As above, we scale by the ratio to compare the shape of the curve. 

699 ratio = np.median(testResp/testResp2) 

700 self.assertFloatsAlmostEqual(testResp/ratio, testResp2, atol=0.04) 

701 

702 def _testFgcmOutputIlluminationCorrection(self, instName, testName, detector): 

703 """Test running of FgcmOutputIlluminationCorrectionTask. 

704 

705 Parameters 

706 ---------- 

707 instName : `str` 

708 Short name of the instrument. 

709 testName : `str` 

710 Base name of the test collection. 

711 detector : `int` 

712 Detector ID to test illumination correction output. 

713 """ 

714 instCamel = instName.title() 

715 

716 configFiles = { 

717 "fgcmOutputIlluminationCorrection": [ 

718 os.path.join( 

719 ROOT, 

720 "config", 

721 f"fgcmOutputIlluminationCorrection{instCamel}.py", 

722 ), 

723 ], 

724 } 

725 inputCollection = f"{instName}/{testName}/fit" 

726 outputCollection = f"{instName}/{testName}/fit/illumcorr" 

727 

728 self._runPipeline( 

729 self.repo, 

730 os.path.join( 

731 ROOT, 

732 "pipelines", 

733 f"fgcmOutputIlluminationCorrection{instCamel}.yaml", 

734 ), 

735 configFiles=configFiles, 

736 inputCollections=[inputCollection], 

737 outputCollection=outputCollection, 

738 registerDatasetTypes=True, 

739 queryString=f"detector = {detector}", 

740 ) 

741 

742 butler = dafButler.Butler(self.repo) 

743 self.enterContext(butler) 

744 

745 # Check for illumination correction files. 

746 illumCorrRefs = butler.query_datasets( 

747 "illuminationCorrection", 

748 collections=[outputCollection], 

749 ) 

750 self.assertGreater(len(illumCorrRefs), 0) 

751 

752 with warnings.catch_warnings(record=True) as w: 

753 warnings.simplefilter("always") 

754 testCorr = butler.get(illumCorrRefs[0]) 

755 

756 self.assertEqual(len(w), 0, "Reading an illumination correction had a warning!") 

757 self.assertIsInstance(testCorr, lsst.afw.image.ExposureF) 

758 

759 self.assertIsNotNone(testCorr.metadata.get("CALIB_CREATION_DATETIME")) 

760 self.assertIsNotNone(testCorr.metadata.get("CALIB_CREATION_DATE")) 

761 self.assertIsNotNone(testCorr.metadata.get("CALIB_CREATION_TIME")) 

762 self.assertIsNotNone(testCorr.metadata.get("LSST CALIB UUID FLAT")) 

763 self.assertIsNotNone(testCorr.metadata.get("LSST ISR UNITS")) 

764 self.assertIsNotNone(testCorr.getDetector()) 

765 self.assertEqual(testCorr.getDetector().getId(), illumCorrRefs[0].dataId["detector"]) 

766 

767 def _testFgcmMultiFit(self, instName, testName, queryString, visits, zpOffsets, 

768 nPlotsPenultimateCycle, nPlotsFinalCycle, 

769 refcatCollection="refcats/gen2"): 

770 """Test running the full pipeline with multiple fit cycles. 

771 

772 Parameters 

773 ---------- 

774 instName : `str` 

775 Short name of the instrument. 

776 testName : `str` 

777 Base name of the test collection. 

778 queryString : `str` 

779 Query to send to the pipetask. 

780 visits : `list` 

781 List of visits to calibrate. 

782 zpOffsets : `np.ndarray` 

783 Zeropoint offsets expected. 

784 nPlotsPenultimateCycle : `int` 

785 Number of plots expected in penultimate cycle. 

786 nPlotsFinalCycle : `int` 

787 Number of plots expected in final cycle. 

788 refcatCollection : `str`, optional 

789 Name of reference catalog collection.