Coverage for tests/test_ptc.py: 6%

1#!/usr/bin/env python 

2 

3# 

4# LSST Data Management System 

5# 

6# Copyright 2008-2017 AURA/LSST. 

7# 

8# This product includes software developed by the 

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

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 LSST License Statement and 

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

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

24# 

25"""Test cases for cp_pipe.""" 

26 

27import unittest 

28import numpy as np 

29import copy 

30import tempfile 

31import logging 

32 

33import lsst.log 

34import lsst.utils 

35import lsst.utils.tests 

36 

37import lsst.cp.pipe as cpPipe 

38import lsst.ip.isr.isrMock as isrMock 

39from lsst.ip.isr import PhotonTransferCurveDataset, PhotodiodeCalib, AmpOffsetTask 

40from lsst.cp.pipe.utils import makeMockFlats, ampOffsetGainRatioFixup 

41 

42from lsst.pipe.base import InMemoryDatasetHandle 

43 

44 

45class FakeCamera(list): 

46 def getName(self): 

47 return "FakeCam" 

48 

49 

50class PretendRef: 

51 "A class to act as a mock exposure reference" 

52 

53 def __init__(self, exposure): 

54 self.exp = exposure 

55 

56 def get(self, component=None): 

57 if component == "visitInfo": 

58 return self.exp.info.getVisitInfo() 

59 elif component == "detector": 

60 return self.exp.getDetector() 

61 elif component == "metadata": 

62 return self.exp.getMetadata() 

63 else: 

64 return self.exp 

65 

66 

67class MeasurePhotonTransferCurveTaskTestCase(lsst.utils.tests.TestCase): 

68 """A test case for the PTC tasks.""" 

69 

70 def setUp(self): 

71 self.defaultConfigExtract = ( 

72 cpPipe.ptc.PhotonTransferCurveExtractPairTask.ConfigClass() 

73 ) 

74 self.defaultTaskExtract = cpPipe.ptc.PhotonTransferCurveExtractTask( 

75 config=self.defaultConfigExtract 

76 ) 

77 

78 self.defaultConfigSolve = cpPipe.ptc.PhotonTransferCurveSolveTask.ConfigClass() 

79 self.defaultTaskSolve = cpPipe.ptc.PhotonTransferCurveSolveTask( 

80 config=self.defaultConfigSolve 

81 ) 

82 

83 self.flatMean = 2000 

84 self.readNoiseAdu = 10 

85 mockImageConfig = isrMock.IsrMock.ConfigClass() 

86 

87 # flatDrop is not really relevant as we replace the data 

88 # but good to note it in case we change how this image is made 

89 mockImageConfig.flatDrop = 0.99999 

90 mockImageConfig.isTrimmed = True 

91 

92 self.flatExp1 = isrMock.FlatMock(config=mockImageConfig).run() 

93 self.flatExp2 = self.flatExp1.clone() 

94 (shapeY, shapeX) = self.flatExp1.getDimensions() 

95 

96 self.flatWidth = np.sqrt(self.flatMean) + self.readNoiseAdu 

97 

98 self.rng1 = np.random.RandomState(1984) 

99 flatData1 = self.rng1.normal(self.flatMean, self.flatWidth, (shapeX, shapeY)) 

100 self.rng2 = np.random.RandomState(666) 

101 flatData2 = self.rng2.normal(self.flatMean, self.flatWidth, (shapeX, shapeY)) 

102 

103 self.flatExp1.image.array[:] = flatData1 

104 self.flatExp2.image.array[:] = flatData2 

105 

106 # create fake PTC data to see if fit works, for one amp ('amp') 

107 self.flux = 1000.0 # ADU/sec 

108 self.timeVec = np.arange(1.0, 101.0, 5) 

109 self.k2NonLinearity = -5e-6 

110 # quadratic signal-chain non-linearity 

111 muVec = self.flux * self.timeVec + self.k2NonLinearity * self.timeVec**2 

112 self.gain = 0.75 # e-/ADU 

113 self.c1 = 1.0 / self.gain 

114 self.noiseSq = 2 * self.gain # 7.5 (e-)^2 

115 self.a00 = -1.2e-6 

116 self.c2 = -1.5e-6 

117 self.c3 = -4.7e-12 # tuned so that it turns over for 200k mean 

118 self.photoCharges = np.linspace(1e-8, 1e-5, len(self.timeVec)) 

119 

120 self.ampNames = [ 

121 amp.getName() for amp in self.flatExp1.getDetector().getAmplifiers() 

122 ] 

123 self.dataset = PhotonTransferCurveDataset(self.ampNames, ptcFitType="PARTIAL") 

124 self.covariancesSqrtWeights = {} 

125 for ( 

126 ampName 

127 ) in self.ampNames: # just the expTimes and means here - vars vary per function 

128 self.dataset.rawExpTimes[ampName] = self.timeVec 

129 self.dataset.rawMeans[ampName] = muVec 

130 self.dataset.covariancesSqrtWeights[ampName] = np.zeros( 

131 (1, self.dataset.covMatrixSide, self.dataset.covMatrixSide) 

132 ) 

133 

134 def test_covAstier(self): 

135 """Test to check getCovariancesAstier 

136 

137 We check that the gain is the same as the input gain from the 

138 mock data, that the covariances via FFT (as it is in 

139 MeasurePhotonTransferCurveTask when doCovariancesAstier=True) 

140 are the same as calculated in real space, and that Cov[0, 0] 

141 (i.e., the variances) are similar to the variances calculated 

142 with the standard method (when doCovariancesAstier=false), 

143 

144 """ 

145 extractConfig = self.defaultConfigExtract 

146 extractConfig.minNumberGoodPixelsForCovariance = 5000 

147 extractConfig.detectorMeasurementRegion = "FULL" 

148 extractConfig.doExtractPhotodiodeData = True 

149 extractConfig.doKsHistMeasurement = True 

150 extractConfig.auxiliaryHeaderKeys = ["CCOBCURR", "CCDTEMP"] 

151 extractPairTask = cpPipe.ptc.PhotonTransferCurveExtractPairTask(config=extractConfig) 

152 

153 # Create solve task config 

154 solveConfig = self.defaultConfigSolve 

155 # Cut off the low-flux point which is a bad fit, and this 

156 # also exercises this functionality and makes the tests 

157 # run a lot faster. 

158 solveConfig.minMeanSignal["ALL_AMPS"] = 2000.0 

159 # Set the outlier fit threshold higher than the default appropriate 

160 # for this test dataset. 

161 solveConfig.maxSignalInitialPtcOutlierFit = 90000.0 

162 # Given the limited nature of the test data, we can only 

163 # reliably fit out to a 3x3 covariance matrix. 

164 # Improvements will be investigated on DM-46131. 

165 solveConfig.maximumRangeCovariancesAstierFullCovFit = 3 

166 

167 inputGain = self.gain 

168 

169 muStandard, varStandard = {}, {} 

170 nPixelCovarianceStandard = {} 

171 expHandles = [] 

172 expIds = [] 

173 pdHandles = [] 

174 idCounter = 0 

175 for i, expTime in enumerate(self.timeVec): 

176 mockExp1, mockExp2 = makeMockFlats( 

177 expTime, 

178 gain=inputGain, 

179 readNoiseElectrons=3, 

180 expId1=idCounter, 

181 expId2=idCounter + 1, 

182 ampNames=self.ampNames, 

183 ) 

184 for mockExp in [mockExp1, mockExp2]: 

185 md = mockExp.getMetadata() 

186 # These values are chosen to be easily compared after 

187 # processing for correct ordering. 

188 md['CCOBCURR'] = float(idCounter) 

189 md['CCDTEMP'] = float(idCounter + 1) 

190 mockExp.setMetadata(md) 

191 md['SEQFILE'] = 'a_seqfile' 

192 md['SEQNAME'] = 'a_seqfile' 

193 md['SEQCKSUM'] = 'deadbeef' 

194 

195 mockExpRef1 = PretendRef(mockExp1) 

196 mockExpRef2 = PretendRef(mockExp2) 

197 expHandles.extend((mockExpRef1, mockExpRef2)) 

198 expIds.append(idCounter) 

199 expIds.append(idCounter + 1) 

200 maskValue = mockExp1.mask.getPlaneBitMask(extractPairTask.config.maskNameList) 

201 for ampNumber, ampName in enumerate(self.ampNames): 

202 # cov has (i, j, var, cov, npix) 

203 ( 

204 im1Area, 

205 im2Area, 

206 imStatsCtrl, 

207 mu1, 

208 mu2, 

209 ) = extractPairTask.getImageAreasMasksStats(mockExp1, mockExp2) 

210 nPixelCovariance = ((im1Area.mask.array & maskValue) == 0).sum() 

211 nPixelCovarianceStandard.setdefault(ampName, []).append(nPixelCovariance) 

212 muDiff, varDiff, covAstier, rowMeanVariance = extractPairTask.measureMeanVarCov( 

213 im1Area, im2Area, imStatsCtrl, mu1, mu2 

214 ) 

215 muStandard.setdefault(ampName, []).append(muDiff) 

216 varStandard.setdefault(ampName, []).append(varDiff) 

217 

218 # Make a photodiode dataset to integrate. 

219 timeSamples = np.linspace(0, 20.0, 100) 

220 currentSamples1 = np.zeros(100) 

221 currentSamples1[50] = -1.0*self.photoCharges[i] 

222 currentSamples2 = np.zeros(100) 

223 currentSamples2[50] = -1.0*(self.photoCharges[i] + 1e-12) 

224 

225 pdCalib1 = PhotodiodeCalib(timeSamples=timeSamples, currentSamples=currentSamples1) 

226 pdCalib1.currentScale = -1.0 

227 pdCalib1.integrationMethod = "CHARGE_SUM" 

228 

229 pdCalib2 = PhotodiodeCalib(timeSamples=timeSamples, currentSamples=currentSamples2) 

230 pdCalib2.currentScale = -1.0 

231 pdCalib2.integrationMethod = "CHARGE_SUM" 

232 

233 pdHandles.append( 

234 InMemoryDatasetHandle( 

235 pdCalib1, 

236 dataId={"exposure": idCounter}, 

237 ) 

238 ) 

239 pdHandles.append( 

240 InMemoryDatasetHandle( 

241 pdCalib2, 

242 dataId={"exposure": idCounter + 1}, 

243 ) 

244 ) 

245 idCounter += 2 

246 

247 # Divide out the pairs. 

248 outputCovariances = [] 

249 for i in range(len(expHandles) // 2): 

250 inputExp = [expHandles[2*i], expHandles[2*i + 1]] 

251 inputDims = [expIds[2*i], expIds[2*i + 1]] 

252 inputPhotodiodeData = [pdHandles[2*i], pdHandles[2*i + 1]] 

253 results = extractPairTask.run( 

254 inputExp=inputExp, 

255 inputDims=inputDims, 

256 inputPhotodiodeData=inputPhotodiodeData, 

257 ) 

258 outputCovariances.append(results.outputCovariance) 

259 

260 # Force the last PTC dataset to have a NaN, and ensure that the 

261 # task runs (DM-38029). This is a minor perturbation and does not 

262 # affect the output comparison. 

263 outputCovariances[-1].rawMeans["C:0,0"] = np.array([np.nan]) 

264 outputCovariances[-1].rawVars["C:0,0"] = np.array([np.nan]) 

265 

266 # Force the next-to-last PTC dataset to have a decreased variance to 

267 # ensure that the outlier fit rejection works. 

268 rawVar = outputCovariances[-2].rawVars["C:0,0"] 

269 outputCovariances[-2].rawVars["C:0,0"] = rawVar * 0.9 

270 

271 # Reorganize the outputCovariances so we can confirm they come 

272 # out sorted afterwards. 

273 outputCovariancesRev = outputCovariances[::-1] 

274 

275 # Some expected values for noise matrix, just to check that 

276 # it was calculated. 

277 expectedNoiseMatrix = { 

278 "FULLCOVARIANCE": np.array( 

279 [[8.99474598, 9.94916264, -27.90587299], 

280 [-2.95079527, -17.11827641, -47.88156244], 

281 [5.24915021, -3.25786165, 26.09634067]], 

282 ), 

283 "FULLCOVARIANCE_NO_B": np.array( 

284 [[8.71049338, 12.48584043, -37.06585088], 

285 [-4.80523971, -23.29102809, -66.37815343], 

286 [7.48654766, -4.10168337, 35.64469824]], 

287 ), 

288 } 

289 

290 for fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]: 

291 solveConfig.ptcFitType = fitType 

292 solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig) 

293 

294 # Ensure no warnings re: read noise mismatches are logged. 

295 with self.assertNoLogs(level=logging.WARNING): 

296 resultsSolve = solveTask.run( 

297 outputCovariancesRev, camera=FakeCamera([self.flatExp1.getDetector()]) 

298 ) 

299 

300 ptc = resultsSolve.outputPtcDataset 

301 

302 # Confirm that metadata keywords are set. 

303 for key in ["SEQFILE", "SEQNAME", "SEQCKSUM"]: 

304 self.assertEqual(ptc.metadata[key], expHandles[0].get().metadata[key]) 

305 self.assertEqual(ptc.metadata["INSTRUME"], FakeCamera().getName()) 

306 det = expHandles[0].get().getDetector() 

307 self.assertEqual(ptc.metadata["DETECTOR"], det.getId()) 

308 self.assertEqual(ptc.metadata["DET_NAME"], det.getName()) 

309 self.assertEqual(ptc.metadata["DET_SER"], det.getSerial()) 

310 

311 for amp in self.ampNames: 

312 self.assertAlmostEqual(ptc.gain[amp], inputGain, places=2) 

313 self.assertFloatsAlmostEqual( 

314 np.asarray(varStandard[amp])[ptc.expIdMask[amp]] / ptc.finalVars[amp][ptc.expIdMask[amp]], 

315 1.0, 

316 rtol=1e-4, 

317 ) 

318 np.testing.assert_array_equal(ptc.nPixelCovariances[amp], nPixelCovarianceStandard[amp][0]) 

319 

320 # Check that the PTC turnoff is correctly computed. 

321 # This will be different for the C:0,0 amp. 

322 if amp == "C:0,0": 

323 self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[amp][-3]) 

324 else: 

325 self.assertAlmostEqual(ptc.ptcTurnoff[amp], ptc.rawMeans[amp][-1]) 

326 

327 # Test that all the quantities are correctly ordered and 

328 # have not accidentally been masked. 

329 for i, extractPtc in enumerate(outputCovariances): 

330 self.assertFloatsAlmostEqual( 

331 extractPtc.rawExpTimes[ampName][0], 

332 ptc.rawExpTimes[ampName][i], 

333 ) 

334 self.assertFloatsAlmostEqual( 

335 extractPtc.rawMeans[ampName][0], 

336 ptc.rawMeans[ampName][i], 

337 ) 

338 self.assertFloatsAlmostEqual( 

339 extractPtc.rawVars[ampName][0], 

340 ptc.rawVars[ampName][i], 

341 ) 

342 self.assertFloatsAlmostEqual( 

343 extractPtc.rawDeltas[ampName][0], 

344 ptc.rawDeltas[ampName][i], 

345 ) 

346 # Ensure these are set to a non-zero value. 

347 self.assertFloatsNotEqual( 

348 extractPtc.rawDeltas[ampName][0], 

349 0.0, 

350 ) 

351 self.assertFloatsAlmostEqual( 

352 extractPtc.photoCharges[ampName][0], 

353 ptc.photoCharges[ampName][i], 

354 ) 

355 self.assertFloatsAlmostEqual( 

356 extractPtc.photoChargeDeltas[ampName][0], 

357 ptc.photoChargeDeltas[ampName][i], 

358 ) 

359 # Ensure these are set to a non-zero value. 

360 self.assertFloatsNotEqual( 

361 extractPtc.photoChargeDeltas[ampName][0], 

362 0.0, 

363 ) 

364 self.assertFloatsAlmostEqual( 

365 extractPtc.histVars[ampName][0], 

366 ptc.histVars[ampName][i], 

367 ) 

368 self.assertFloatsAlmostEqual( 

369 extractPtc.histChi2Dofs[ampName][0], 

370 ptc.histChi2Dofs[ampName][i], 

371 ) 

372 self.assertFloatsAlmostEqual( 

373 extractPtc.kspValues[ampName][0], 

374 ptc.kspValues[ampName][i], 

375 ) 

376 self.assertFloatsAlmostEqual( 

377 extractPtc.covariances[ampName][0], 

378 ptc.covariances[ampName][i], 

379 ) 

380 self.assertFloatsAlmostEqual( 

381 extractPtc.covariancesSqrtWeights[ampName][0], 

382 ptc.covariancesSqrtWeights[ampName][i], 

383 ) 

384 self.assertFloatsAlmostEqual( 

385 ptc.noiseMatrix[ampName], 

386 expectedNoiseMatrix[fitType], 

387 atol=1e-8, 

388 rtol=None, 

389 ) 

390 self.assertFloatsAlmostEqual( 

391 ptc.ampOffsets[ampName], 

392 0.0, 

393 ) 

394 self.assertFloatsAlmostEqual( 

395 ptc.noise[ampName], 

396 np.nanmedian(ptc.noiseList[ampName]) * ptc.gain[ampName], 

397 rtol=0.05, 

398 ) 

399 # If the noise error is greater than the noise, 

400 # something is seriously wrong. Possibly some 

401 # kind of gain application mismatch. 

402 self.assertLess( 

403 ptc.noiseErr[ampName], 

404 ptc.noise[ampName], 

405 ) 

406 

407 mask = ptc.getGoodPoints(amp) 

408 

409 values = ( 

410 ptc.covariancesModel[amp][mask, 0, 0] - ptc.covariances[amp][mask, 0, 0] 

411 ) / ptc.covariancesModel[amp][mask, 0, 0] 

412 np.testing.assert_array_less(np.abs(values), 2e-3) 

413 

414 if ptc.ptcFitType == "FULLCOVARIANCE": 

415 values = ( 

416 ptc.covariancesModel[amp][mask, 0, 1] - ptc.covariances[amp][mask, 0, 1] 

417 ) / ptc.covariancesModel[amp][mask, 0, 1] 

418 np.testing.assert_array_less(np.abs(values), 0.3) 

419 

420 values = ( 

421 ptc.covariancesModel[amp][mask, 1, 0] - ptc.covariances[amp][mask, 1, 0] 

422 ) / ptc.covariancesModel[amp][mask, 1, 0] 

423 np.testing.assert_array_less(np.abs(values), 0.3) 

424 

425 # And test that the auxiliary values are there and 

426 # correctly ordered. 

427 self.assertIn('CCOBCURR', ptc.auxValues) 

428 self.assertIn('CCDTEMP', ptc.auxValues) 

429 firstExpIds = np.array([i for i, _ in ptc.inputExpIdPairs['C:0,0']], dtype=np.float64) 

430 self.assertFloatsAlmostEqual(ptc.auxValues['CCOBCURR'], firstExpIds) 

431 self.assertFloatsAlmostEqual(ptc.auxValues['CCDTEMP'], firstExpIds + 1) 

432 

433 expIdsUsed = ptc.getExpIdsUsed("C:0,0") 

434 # Check that these are the same as the inputs, paired up, with the 

435 # first two (low flux) and final four (outliers, nans) removed. 

436 self.assertTrue( 

437 np.all(expIdsUsed == np.array(expIds).reshape(len(expIds) // 2, 2)[1:-2]) 

438 ) 

439 

440 goodAmps = ptc.getGoodAmps() 

441 self.assertEqual(goodAmps, self.ampNames) 

442 

443 # Check that every possibly modified field has the same length. 

444 covShape = None 

445 covSqrtShape = None 

446 covModelShape = None 

447 

448 for ampName in self.ampNames: 

449 if covShape is None: 

450 covShape = ptc.covariances[ampName].shape 

451 covSqrtShape = ptc.covariancesSqrtWeights[ampName].shape 

452 covModelShape = ptc.covariancesModel[ampName].shape 

453 else: 

454 self.assertEqual(ptc.covariances[ampName].shape, covShape) 

455 self.assertEqual( 

456 ptc.covariancesSqrtWeights[ampName].shape, covSqrtShape 

457 ) 

458 self.assertEqual(ptc.covariancesModel[ampName].shape, covModelShape) 

459 

460 # Check if evalPtcModel produces expected values 

461 nanMask = ~np.isnan(ptc.finalMeans[ampName]) 

462 means = ptc.finalMeans[ampName][nanMask] 

463 covModel = ptc.covariancesModel[ampName][nanMask] 

464 covariancesModel = ptc.evalPtcModel(means)[ampName] 

465 self.assertFloatsAlmostEqual(covariancesModel, covModel, atol=1e-12) 

466 

467 # And check that this is serializable 

468 with tempfile.NamedTemporaryFile(suffix=".fits") as f: 

469 usedFilename = ptc.writeFits(f.name) 

470 fromFits = PhotonTransferCurveDataset.readFits(usedFilename) 

471 self.assertEqual(fromFits, ptc) 

472 

473 def ptcFitAndCheckPtc( 

474 self, 

475 fitType=None, 

476 doFitBootstrap=False, 

477 ): 

478 localDataset = copy.deepcopy(self.dataset) 

479 localDataset.ptcFitType = fitType 

480 configSolve = copy.copy(self.defaultConfigSolve) 

481 if doFitBootstrap: 

482 configSolve.doFitBootstrap = True 

483 

484 if fitType == "EXPAPPROXIMATION": 

485 g = self.gain 

486 for ampName in self.ampNames: 

487 localDataset.rawVars[ampName] = [ 

488 ( 

489 0.5 / (self.a00 * g**2) * (np.exp(2 * self.a00 * mu * g) - 1) 

490 + self.noiseSq / (g * g) 

491 ) 

492 for mu in localDataset.rawMeans[ampName] 

493 ] 

494 else: 

495 raise RuntimeError( 

496 "Fit type must be 'EXPAPPROXIMATION'" 

497 ) 

498 

499 # Initialize mask and covariance weights that will be used in fits. 

500 # Covariance weights values empirically determined from one of 

501 # the cases in test_covAstier. 

502 matrixSize = localDataset.covMatrixSide 

503 maskLength = len(localDataset.rawMeans[ampName]) 

504 for ampName in self.ampNames: 

505 localDataset.expIdMask[ampName] = np.repeat(True, maskLength) 

506 localDataset.covariancesSqrtWeights[ampName] = np.repeat( 

507 np.ones((matrixSize, matrixSize)), maskLength 

508 ).reshape((maskLength, matrixSize, matrixSize)) 

509 localDataset.covariancesSqrtWeights[ampName][:, 0, 0] = [ 

510 0.07980188, 

511 0.01339653, 

512 0.0073118, 

513 0.00502802, 

514 0.00383132, 

515 0.00309475, 

516 0.00259572, 

517 0.00223528, 

518 0.00196273, 

519 0.00174943, 

520 0.00157794, 

521 0.00143707, 

522 0.00131929, 

523 0.00121935, 

524 0.0011334, 

525 0.00105893, 

526 0.00099357, 

527 0.0009358, 

528 0.00088439, 

529 0.00083833, 

530 ] 

531 

532 solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=configSolve) 

533 

534 if fitType == "EXPAPPROXIMATION": 

535 localDataset = solveTask.fitPtc(localDataset) 

536 else: 

537 localDataset = solveTask.fitDataFullCovariance(localDataset) 

538 

539 # Check entries in localDataset, which was modified by the function 

540 for ampName in self.ampNames: 

541 self.assertEqual(fitType, localDataset.ptcFitType) 

542 self.assertAlmostEqual(self.gain, localDataset.gain[ampName]) 

543 if fitType == "EXPAPPROXIMATION": 

544 self.assertAlmostEqual( 

545 self.a00, localDataset.ptcFitPars[ampName][0] 

546 ) 

547 # Noise already in electrons 

548 self.assertAlmostEqual( 

549 np.sqrt(self.noiseSq), localDataset.noise[ampName] 

550 ) 

551 # If the noise error is greater than the noise, something 

552 # is seriously wrong. Possibly some kind of gain application 

553 # mismatch. 

554 self.assertLess( 

555 localDataset.noiseErr[ampName], 

556 np.sqrt(self.noiseSq) 

557 ) 

558 

559 # Check if evalColModel produces expected values 

560 if not doFitBootstrap: 

561 nanMask = ~np.isnan(localDataset.finalMeans[ampName]) 

562 means = localDataset.finalMeans[ampName][nanMask] 

563 model = localDataset.finalModelVars[ampName][nanMask] 

564 evalVarModel = localDataset.evalPtcModel(means)[ampName] 

565 self.assertFloatsEqual(evalVarModel, model) 

566 

567 def test_lsstcam_samples(self): 

568 for dense in [False, True]: 

569 for mode in ["normal", "upturn", "dip"]: 

570 for doModelPtcRolloff in [False, True]: 

571 rawMeans, rawVars, ptcTurnoff = self._getSampleMeanAndVar(dense=dense, mode=mode) 

572 

573 # We only need a single amp. 

574 ampName = self.ampNames[0] 

575 dataset = PhotonTransferCurveDataset([ampName], ptcFitType="EXPAPPROXIMATION") 

576 dataset.rawExpTimes[ampName] = np.arange(len(rawMeans), dtype=np.float64) + 1.0 

577 dataset.rawMeans[ampName] = rawMeans 

578 dataset.rawVars[ampName] = rawVars 

579 dataset.covariancesSqrtWeights[ampName] = np.repeat( 

580 np.ones((dataset.covMatrixSide, dataset.covMatrixSide)), len(rawMeans) 

581 ).reshape((len(rawMeans), dataset.covMatrixSide, dataset.covMatrixSide)) 

582 dataset.covariancesSqrtWeights[ampName][:, 0, 0] = np.sqrt(rawVars) 

583 mask = np.ones(len(rawMeans), dtype=np.bool_) 

584 mask[~np.isfinite(rawMeans) | ~np.isfinite(rawVars)] = False 

585 dataset.expIdMask[ampName] = mask 

586 

587 configSolve = copy.copy(self.defaultConfigSolve) 

588 configSolve.doModelPtcRolloff = doModelPtcRolloff 

589 configSolve.doFitBootstrap = False 

590 

591 solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=configSolve) 

592 

593 initialFitDataset = solveTask.fitPtc(dataset, computePtcTurnoff=True) 

594 

595 # Model the PTC rolloff 

596 if doModelPtcRolloff: 

597 initialFitDataset = solveTask.fitPtcRolloff(initialFitDataset) 

598 

599 solvedDataset = solveTask.fitPtc(initialFitDataset, computePtcTurnoff=False) 

600 # Check that the ptcTurnoff is what is expected. 

601 self.assertFloatsAlmostEqual( 

602 solvedDataset.ptcTurnoff[ampName], 

603 ptcTurnoff, 

604 msg=( 

605 f"Dense: {dense}; Mode: {mode}, " 

606 f"doModelPtcRolloff: {doModelPtcRolloff}; Amp: {ampName}" 

607 ), 

608 ) 

609 

610 # Check that no values above the 

611 # turnoff/rolloff are "good". 

612 maxMean = solvedDataset.ptcTurnoff[ampName] 

613 if doModelPtcRolloff: 

614 # Check that the PTC rolloff is less 

615 # than the PTC turnoff and within 10% 

616 # of the turnoff. 

617 if dense: 

618 # Check that 

619 self.assertLess( 

620 solvedDataset.ptcRolloff[ampName], 

621 solvedDataset.ptcTurnoff[ampName], 

622 msg=( 

623 f"Dense: {dense}; Mode: {mode}, " 

624 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

625 f"Amp: {ampName}" 

626 ), 

627 ) 

628 self.assertFloatsAlmostEqual( 

629 solvedDataset.ptcRolloff[ampName], 

630 solvedDataset.ptcTurnoff[ampName], 

631 atol=0.1 * solvedDataset.ptcTurnoff[ampName], 

632 msg=( 

633 f"Dense: {dense}; Mode: {mode}, " 

634 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

635 f"Amp: {ampName}" 

636 ), 

637 ) 

638 else: 

639 # In all other cases, it should not be able to 

640 # find points beyond the turnoff, and the 

641 # rolloff should be set to the PTC turnoff. 

642 self.assertEqual( 

643 solvedDataset.ptcRolloff[ampName], 

644 solvedDataset.ptcTurnoff[ampName], 

645 msg=( 

646 f"Dense: {dense}; Mode: {mode}, " 

647 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

648 f"Amp: {ampName}" 

649 ), 

650 ) 

651 

652 # If that passes set the maxMean to the PTC Rolloff 

653 maxMean = solvedDataset.ptcRolloff[ampName] 

654 

655 above = (solvedDataset.finalMeans[ampName] > maxMean) 

656 self.assertEqual( 

657 np.sum(above), 

658 0, 

659 msg=( 

660 f"Dense: {dense}; Mode: {mode}, " 

661 f"doModelPtcRolloff: {doModelPtcRolloff}; Amp: {ampName}" 

662 ), 

663 ) 

664 

665 # Check the sampling error on the turnoff. 

666 turnoffIdx = np.argwhere( 

667 solvedDataset.rawMeans[ampName] == solvedDataset.ptcTurnoff[ampName] 

668 ) 

669 samplingError = (rawMeans[turnoffIdx + 1] - rawMeans[turnoffIdx - 1])/2. 

670 self.assertFloatsAlmostEqual( 

671 solvedDataset.ptcTurnoffSamplingError[ampName], 

672 samplingError, 

673 msg=f"Dense: {dense}; Mode: {mode}", 

674 ) 

675 

676 def test_ptcFit(self): 

677 for fitType in ["EXPAPPROXIMATION"]: 

678 self.ptcFitAndCheckPtc( 

679 fitType=fitType 

680 ) 

681 

682 def test_meanVarMeasurement(self): 

683 task = self.defaultTaskExtract 

684 im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats( 

685 self.flatExp1, self.flatExp2 

686 ) 

687 mu, varDiff, _, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2) 

688 

689 self.assertLess(self.flatWidth - np.sqrt(varDiff), 1) 

690 self.assertLess(self.flatMean - mu, 1) 

691 

692 def test_meanVarMeasurementWithNans(self): 

693 task = self.defaultTaskExtract 

694 

695 flatExp1 = self.flatExp1.clone() 

696 flatExp2 = self.flatExp2.clone() 

697 

698 flatExp1.image.array[20:30, :] = np.nan 

699 flatExp2.image.array[20:30, :] = np.nan 

700 

701 im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats( 

702 flatExp1, flatExp2 

703 ) 

704 mu, varDiff, _, _ = task.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2) 

705 

706 expectedMu1 = np.nanmean(flatExp1.image.array) 

707 expectedMu2 = np.nanmean(flatExp2.image.array) 

708 expectedMu = 0.5 * (expectedMu1 + expectedMu2) 

709 

710 # Now the variance of the difference. First, create the diff image. 

711 im1 = flatExp1.maskedImage 

712 im2 = flatExp2.maskedImage 

713 

714 temp = im2.clone() 

715 temp *= expectedMu1 

716 diffIm = im1.clone() 

717 diffIm *= expectedMu2 

718 diffIm -= temp 

719 diffIm /= expectedMu 

720 

721 # Divide by two as it is what measureMeanVarCov returns 

722 # (variance of difference) 

723 expectedVar = 0.5 * np.nanvar(diffIm.image.array) 

724 

725 # Check that the standard deviations and the emans agree to 

726 # less than 1 ADU 

727 self.assertLess(np.sqrt(expectedVar) - np.sqrt(varDiff), 1) 

728 self.assertLess(expectedMu - mu, 1) 

729 

730 def test_meanVarMeasurementAllNan(self): 

731 task = self.defaultTaskExtract 

732 flatExp1 = self.flatExp1.clone() 

733 flatExp2 = self.flatExp2.clone() 

734 

735 flatExp1.image.array[:, :] = np.nan 

736 flatExp2.image.array[:, :] = np.nan 

737 

738 im1Area, im2Area, imStatsCtrl, mu1, mu2 = task.getImageAreasMasksStats( 

739 flatExp1, flatExp2 

740 ) 

741 mu, varDiff, covDiff, rowMeanVariance = task.measureMeanVarCov( 

742 im1Area, im2Area, imStatsCtrl, mu1, mu2 

743 )