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 butler = dafButler.Butler(cls.repo, writeable=True) 

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 pipeline = Pipeline.fromFile(pipelineFile) 

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

125 for fileName in fileList: 

126 pipeline.addConfigFile(taskName, fileName) 

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

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

129 pipeline.addConfigOverride(taskName, option, value) 

130 

131 resources = ExecutionResources(num_cores=1) 

132 

133 executor = SimplePipelineExecutor.from_pipeline(pipeline, 

134 where=queryString, 

135 butler=butler, 

136 resources=resources) 

137 quanta = executor.run(register_dataset_types=registerDatasetTypes) 

138 

139 return len(quanta) 

140 

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

142 """Test running of FgcmMakeLutTask 

143 

144 Parameters 

145 ---------- 

146 instName : `str` 

147 Short name of the instrument. 

148 testName : `str` 

149 Base name of the test collection. 

150 nBand : `int` 

151 Number of bands tested. 

152 i0Std : `np.ndarray' 

153 Values of i0Std to compare to. 

154 i10Std : `np.ndarray` 

155 Values of i10Std to compare to. 

156 i0Recon : `np.ndarray` 

157 Values of reconstructed i0 to compare to. 

158 i10Recon : `np.ndarray` 

159 Values of reconsntructed i10 to compare to. 

160 """ 

161 instCamel = instName.title() 

162 

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

164 'config', 

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

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

167 

168 self._runPipeline(self.repo, 

169 os.path.join(ROOT, 

170 'pipelines', 

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

172 configFiles=configFiles, 

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

174 outputCollection=outputCollection, 

175 registerDatasetTypes=True) 

176 

177 # Check output values 

178 butler = dafButler.Butler(self.repo) 

179 lutCat = butler.get('fgcmLookUpTable', 

180 collections=[outputCollection], 

181 instrument=instName) 

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

183 

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

185 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

200 indices) 

201 

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

203 

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

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

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

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

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

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

210 indices) 

211 

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

213 

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

215 atmStd = butler.get('fgcm_standard_atmosphere', 

216 collections=[outputCollection], 

217 instrument=instName) 

218 bounds = atmStd.getWavelengthBounds() 

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

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

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

222 

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

224 for physical_filter in fgcmLut.filterNames: 

225 passband = butler.get('fgcm_standard_passband', 

226 collections=[outputCollection], 

227 instrument=instName, 

228 physical_filter=physical_filter) 

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

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

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

232 

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

234 refs = butler.query_datasets( 

235 "standard_passband", 

236 collections=[outputCollection], 

237 instrument=instName, 

238 ) 

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

240 for ref in refs: 

241 passbandTable = butler.get(ref) 

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

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

244 

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

246 refcatCollection="refcats/gen2"): 

247 """Test running of FgcmBuildStarsTableTask 

248 

249 Parameters 

250 ---------- 

251 instName : `str` 

252 Short name of the instrument. 

253 testName : `str` 

254 Base name of the test collection. 

255 queryString : `str` 

256 Query to send to the pipetask. 

257 visits : `list` 

258 List of visits to calibrate. 

259 nStar : `int` 

260 Number of stars expected. 

261 nObs : `int` 

262 Number of observations of stars expected. 

263 refcatCollection : `str`, optional 

264 Name of reference catalog collection. 

265 """ 

266 instCamel = instName.title() 

267 

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

269 'config', 

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

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

272 

273 self._runPipeline(self.repo, 

274 os.path.join(ROOT, 

275 'pipelines', 

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

277 configFiles=configFiles, 

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

279 refcatCollection], 

280 outputCollection=outputCollection, 

281 queryString=queryString, 

282 registerDatasetTypes=True) 

283 

284 butler = dafButler.Butler(self.repo) 

285 

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

287 instrument=instName) 

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

289 

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

291 instrument=instName) 

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

293 

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

295 instrument=instName) 

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

297 

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

299 refcatCollection="refcats/gen2"): 

300 """Test running of FgcmBuildFromIsolatedStarsTask. 

301 

302 Parameters 

303 ---------- 

304 instName : `str` 

305 Short name of the instrument. 

306 testName : `str` 

307 Base name of the test collection. 

308 queryString : `str` 

309 Query to send to the pipetask. 

310 visits : `list` 

311 List of visits to calibrate. 

312 nStar : `int` 

313 Number of stars expected. 

314 nObs : `int` 

315 Number of observations of stars expected. 

316 refcatCollection : `str`, optional 

317 Name of reference catalog collection. 

318 """ 

319 instCamel = instName.title() 

320 

321 configFiles = {'fgcmBuildFromIsolatedStars': [ 

322 os.path.join(ROOT, 

323 'config', 

324 f'fgcmBuildFromIsolatedStars{instCamel}.py') 

325 ]} 

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

327 

328 self._runPipeline(self.repo, 

329 os.path.join(ROOT, 

330 'pipelines', 

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

332 configFiles=configFiles, 

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

334 refcatCollection], 

335 outputCollection=outputCollection, 

336 queryString=queryString, 

337 registerDatasetTypes=True) 

338 

339 butler = dafButler.Butler(self.repo) 

340 

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

342 instrument=instName) 

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

344 

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

346 instrument=instName) 

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

348 

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

350 instrument=instName) 

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

352 

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

354 nZp, nGoodZp, nOkZp, nBadZp, nStdStars, nPlots, 

355 skipChecks=False, extraConfig=None): 

356 """Test running of FgcmFitCycleTask 

357 

358 Parameters 

359 ---------- 

360 instName : `str` 

361 Short name of the instrument. 

362 testName : `str` 

363 Base name of the test collection. 

364 cycleNumber : `int` 

365 Fit cycle number. 

366 nZp : `int` 

367 Number of zeropoints created by the task. 

368 nGoodZp : `int` 

369 Number of good (photometric) zeropoints created. 

370 nOkZp : `int` 

371 Number of constrained zeropoints (photometric or not). 

372 nBadZp : `int` 

373 Number of unconstrained (bad) zeropoints. 

374 nStdStars : `int` 

375 Number of standard stars produced. 

376 nPlots : `int` 

377 Number of plots produced. 

378 skipChecks : `bool`, optional 

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

380 extraConfig : `str`, optional 

381 Name of an extra config file to apply. 

382 """ 

383 instCamel = instName.title() 

384 

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

386 'config', 

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

388 if extraConfig is not None: 

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

390 

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

392 

393 if cycleNumber == 0: 

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

395 else: 

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

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

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

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

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

401 # Note that this behavior is handled automatically by the 

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

403 # API. 

404 butler = dafButler.Butler(self.repo) 

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

406 

407 cwd = os.getcwd() 

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

409 os.makedirs(runDir, exist_ok=True) 

410 os.chdir(runDir) 

411 

412 configOptions = {'fgcmFitCycle': 

413 {'cycleNumber': f'{cycleNumber}', 

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

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

416 self._runPipeline(self.repo, 

417 os.path.join(ROOT, 

418 'pipelines', 

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

420 configFiles=configFiles, 

421 inputCollections=inputCollections, 

422 outputCollection=outputCollection, 

423 configOptions=configOptions, 

424 registerDatasetTypes=True) 

425 

426 os.chdir(cwd) 

427 

428 if skipChecks: 

429 return 

430 

431 butler = dafButler.Butler(self.repo) 

432 

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

434 plotDatasets = list(butler.registry.queryDatasets( 

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

436 collections=[outputCollection], 

437 )) 

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

439 

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

441 collections=[outputCollection], 

442 instrument=instName) 

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

444 

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

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

447 

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

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

450 

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

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

453 

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

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

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

457 

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

459 collections=[outputCollection], 

460 instrument=instName) 

461 

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

463 

464 def _testFgcmOutputProducts(self, instName, testName, 

465 zpOffsets, testVisit, testCcd, testFilter, testBandIndex, 

466 skymapName=None, testSrc=True): 

467 """Test running of FgcmOutputProductsTask. 

468 

469 Parameters 

470 ---------- 

471 instName : `str` 

472 Short name of the instrument. 

473 testName : `str` 

474 Base name of the test collection. 

475 zpOffsets : `np.ndarray` 

476 Zeropoint offsets expected. 

477 testVisit : `int` 

478 Visit id to check for round-trip computations. 

479 testCcd : `int` 

480 Ccd id to check for round-trip computations. 

481 testFilter : `str` 

482 Filtername for testVisit/testCcd. 

483 testBandIndex : `int` 

484 Band index for testVisit/testCcd. 

485 skymapName : `str`, optional 

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

487 testSrc : `bool`, optional 

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

489 """ 

490 instCamel = instName.title() 

491 

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

493 'config', 

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

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

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

497 if skymapName is not None: 

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

499 else: 

500 queryString = "" 

501 

502 self._runPipeline(self.repo, 

503 os.path.join(ROOT, 

504 'pipelines', 

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

506 queryString=queryString, 

507 configFiles=configFiles, 

508 inputCollections=[inputCollection, 'skymaps'], 

509 outputCollection=outputCollection, 

510 registerDatasetTypes=True) 

511 

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

513 offsetCat = butler.get('fgcmReferenceCalibrationOffsets') 

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

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

516 

517 config = butler.get('fgcmOutputProducts_config', 

518 collections=[outputCollection], instrument=instName) 

519 

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

521 collections=[inputCollection]) 

522 

523 # Test the fgcm_photoCalib output 

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

525 

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

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

528 

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

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

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

532 photoCalibDict = {} 

533 for visit in visits: 

534 expCat = butler.get('fgcmPhotoCalibCatalog', 

535 visit=visit) 

536 for row in expCat: 

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

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

539 

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

541 for rec in zptCat[good]: 

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

543 

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

545 for rec in zptCat[bad]: 

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

547 

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

549 testCal = photoCalibDict[(testVisit, testCcd)] 

550 

551 # Check the output standard star catalogs. 

552 if skymapName is not None: 

553 refs = butler.query_datasets("fgcm_standard_star") 

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

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

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

557 

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

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

560 

561 # Check ids against inputs. 

562 isoConfig = butler.get("fgcmBuildFromIsolatedStars_config") 

563 isolatedCatalog = butler.get( 

564 isoConfig.connections.isolated_star_cats, 

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

566 ) 

567 

568 a, b = esutil.numpy_util.match( 

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

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

571 ) 

572 

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

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

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

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

577 

578 if testSrc: 

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

580 collections=[outputCollection], instrument=instName) 

581 

582 # Only test sources with positive flux 

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

584 

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

586 # and doesn't know about that yet) 

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

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

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

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

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

592 

593 if config.doComposeWcsJacobian: 

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

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

596 collections=...) 

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

598 approxPixelAreaFields = fgcmcal.utilities.computeApproxPixelAreaFields(camera) 

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

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

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

602 

603 # This is the magnitude through the mean calibration 

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

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

606 photoCalMags = np.zeros_like(photoCalMeanCalMags) 

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

608 zptMeanCalMags = np.zeros_like(photoCalMeanCalMags) 

609 

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

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

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

613 rec.getCentroid()) 

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

615 

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

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

618 self.assertFloatsAlmostEqual(photoCalMeanCalMags, 

619 zptMeanCalMags, rtol=1e-6) 

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

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

622 # wrong. 

623 self.assertFloatsAlmostEqual(photoCalMeanCalMags, 

624 photoCalMags, rtol=1e-2) 

625 

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

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

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

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

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

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

632 # offsets used in the tests. 

633 

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

635 # (multiple ccds) 

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

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

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

639 where=whereClause, 

640 findFirst=True) 

641 photoCals = [] 

642 for srcRef in srcRefs: 

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

644 

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

646 rawStars, testBandIndex, offsets) 

647 

648 st = np.argsort(matchMag) 

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

650 # deltaMagBkgOffsetPercentile, we want to ensure that these stars 

651 # match on average. 

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

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

654 

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

656 self.assertFloatsAlmostEqual(testCal.getCalibrationErr(), 

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

658 

659 # Test the transmission output 

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

661 instrument=instName) 

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

663 instrument=instName) 

664 

665 testTrans = butler.get('transmission_atmosphere_fgcm', 

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

667 collections=[outputCollection], instrument=instName) 

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

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

670 

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

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

673 # these output atmospheres and the standard is the different 

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

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

676 # testing. 

677 

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

679 # we only care about the shape 

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

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

682 

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

684 # difference so they aren't identical. 

685 testTrans2 = butler.get('transmission_atmosphere_fgcm', 

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

687 collections=[outputCollection], instrument=instName) 

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

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

690 

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

692 ratio = np.median(testResp/testResp2) 

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

694 

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

696 """Test running of FgcmOutputIlluminationCorrectionTask. 

697 

698 Parameters 

699 ---------- 

700 instName : `str` 

701 Short name of the instrument. 

702 testName : `str` 

703 Base name of the test collection. 

704 detector : `int` 

705 Detector ID to test illumination correction output. 

706 """ 

707 instCamel = instName.title() 

708 

709 configFiles = { 

710 "fgcmOutputIlluminationCorrection": [ 

711 os.path.join( 

712 ROOT, 

713 "config", 

714 f"fgcmOutputIlluminationCorrection{instCamel}.py", 

715 ), 

716 ], 

717 } 

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

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

720 

721 self._runPipeline( 

722 self.repo, 

723 os.path.join( 

724 ROOT, 

725 "pipelines", 

726 f"fgcmOutputIlluminationCorrection{instCamel}.yaml", 

727 ), 

728 configFiles=configFiles, 

729 inputCollections=[inputCollection], 

730 outputCollection=outputCollection, 

731 registerDatasetTypes=True, 

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

733 ) 

734 

735 butler = dafButler.Butler(self.repo) 

736 

737 # Check for illumination correction files. 

738 illumCorrRefs = butler.query_datasets( 

739 "illuminationCorrection", 

740 collections=[outputCollection], 

741 ) 

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

743 

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

745 warnings.simplefilter("always") 

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

747 

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

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

750 

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

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

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

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

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

756 self.assertIsNotNone(testCorr.getDetector()) 

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

758 

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

760 refcatCollection="refcats/gen2"): 

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

762 

763 Parameters 

764 ---------- 

765 instName : `str` 

766 Short name of the instrument. 

767 testName : `str` 

768 Base name of the test collection. 

769 queryString : `str` 

770 Query to send to the pipetask. 

771 visits : `list` 

772 List of visits to calibrate. 

773 zpOffsets : `np.ndarray` 

774 Zeropoint offsets expected. 

775 refcatCollection : `str`, optional 

776 Name of reference catalog collection. 

777 """ 

778 instCamel = instName.title() 

779 

780 configFiles = { 

781 'fgcmBuildFromIsolatedStars': [ 

782 os.path.join( 

783 ROOT, 

784 'config', 

785 f'fgcmBuildFromIsolatedStars{instCamel}.py', 

786 ), 

787 ], 

788 'fgcmFitCycle': [ 

789 os.path.join( 

790 ROOT, 

791 'config',