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 

32import warnings 

33 

34import lsst.log 

35import lsst.utils 

36import lsst.utils.tests 

37 

38import lsst.cp.pipe as cpPipe 

39import lsst.ip.isr.isrMock as isrMock 

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

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

42 

43from lsst.pipe.base import InMemoryDatasetHandle 

44 

45 

46class FakeCamera(list): 

47 def getName(self): 

48 return "FakeCam" 

49 

50 

51class PretendRef: 

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

53 

54 def __init__(self, exposure): 

55 self.exp = exposure 

56 

57 def get(self, component=None): 

58 if component == "visitInfo": 

59 return self.exp.info.getVisitInfo() 

60 elif component == "detector": 

61 return self.exp.getDetector() 

62 elif component == "metadata": 

63 return self.exp.getMetadata() 

64 else: 

65 return self.exp 

66 

67 

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

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

70 

71 def setUp(self): 

72 self.defaultConfigExtract = ( 

73 cpPipe.ptc.PhotonTransferCurveExtractPairTask.ConfigClass() 

74 ) 

75 self.defaultTaskExtract = cpPipe.ptc.PhotonTransferCurveExtractTask( 

76 config=self.defaultConfigExtract 

77 ) 

78 

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

80 self.defaultTaskSolve = cpPipe.ptc.PhotonTransferCurveSolveTask( 

81 config=self.defaultConfigSolve 

82 ) 

83 

84 self.flatMean = 2000 

85 self.readNoiseAdu = 10 

86 mockImageConfig = isrMock.IsrMock.ConfigClass() 

87 

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

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

90 mockImageConfig.flatDrop = 0.99999 

91 mockImageConfig.isTrimmed = True 

92 

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

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

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

96 

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

98 

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

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

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

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

103 

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

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

106 

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

108 self.flux = 1000.0 # ADU/sec 

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

110 self.k2NonLinearity = -5e-6 

111 # quadratic signal-chain non-linearity 

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

113 self.gain = 0.75 # e-/ADU 

114 self.c1 = 1.0 / self.gain 

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

116 self.a00 = -1.2e-6 

117 self.c2 = -1.5e-6 

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

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

120 

121 self.ampNames = [ 

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

123 ] 

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

125 self.covariancesSqrtWeights = {} 

126 for ( 

127 ampName 

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

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

130 self.dataset.rawMeans[ampName] = muVec 

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

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

133 ) 

134 

135 def test_covAstier(self): 

136 """Test to check getCovariancesAstier 

137 

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

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

140 MeasurePhotonTransferCurveTask when doCovariancesAstier=True) 

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

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

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

144 

145 """ 

146 extractConfig = self.defaultConfigExtract 

147 extractConfig.minNumberGoodPixelsForCovariance = 5000 

148 extractConfig.detectorMeasurementRegion = "FULL" 

149 extractConfig.doExtractPhotodiodeData = True 

150 extractConfig.doKsHistMeasurement = True 

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

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

153 

154 # Create solve task config 

155 solveConfig = self.defaultConfigSolve 

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

157 # also exercises this functionality and makes the tests 

158 # run a lot faster. 

159 solveConfig.minMeanSignal["ALL_AMPS"] = 2000.0 

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

161 # for this test dataset. 

162 solveConfig.maxSignalInitialPtcOutlierFit = 90000.0 

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

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

165 # Improvements will be investigated on DM-46131. 

166 solveConfig.maximumRangeCovariancesAstierFullCovFit = 3 

167 

168 inputGain = self.gain 

169 

170 muStandard, varStandard = {}, {} 

171 nPixelCovarianceStandard = {} 

172 expHandles = [] 

173 expIds = [] 

174 pdHandles = [] 

175 idCounter = 0 

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

177 mockExp1, mockExp2 = makeMockFlats( 

178 expTime, 

179 gain=inputGain, 

180 readNoiseElectrons=3, 

181 expId1=idCounter, 

182 expId2=idCounter + 1, 

183 ampNames=self.ampNames, 

184 ) 

185 for mockExp in [mockExp1, mockExp2]: 

186 md = mockExp.getMetadata() 

187 # These values are chosen to be easily compared after 

188 # processing for correct ordering. 

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

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

191 mockExp.setMetadata(md) 

192 md['SEQFILE'] = 'a_seqfile' 

193 md['SEQNAME'] = 'a_seqfile' 

194 md['SEQCKSUM'] = 'deadbeef' 

195 

196 mockExpRef1 = PretendRef(mockExp1) 

197 mockExpRef2 = PretendRef(mockExp2) 

198 expHandles.extend((mockExpRef1, mockExpRef2)) 

199 expIds.append(idCounter) 

200 expIds.append(idCounter + 1) 

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

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

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

204 ( 

205 im1Area, 

206 im2Area, 

207 imStatsCtrl, 

208 mu1, 

209 mu2, 

210 ) = extractPairTask.getImageAreasMasksStats(mockExp1, mockExp2) 

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

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

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

214 im1Area, im2Area, imStatsCtrl, mu1, mu2 

215 ) 

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

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

218 

219 # Make a photodiode dataset to integrate. 

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

221 currentSamples1 = np.zeros(100) 

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

223 currentSamples2 = np.zeros(100) 

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

225 

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

227 pdCalib1.currentScale = -1.0 

228 pdCalib1.integrationMethod = "CHARGE_SUM" 

229 

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

231 pdCalib2.currentScale = -1.0 

232 pdCalib2.integrationMethod = "CHARGE_SUM" 

233 

234 pdHandles.append( 

235 InMemoryDatasetHandle( 

236 pdCalib1, 

237 dataId={"exposure": idCounter}, 

238 ) 

239 ) 

240 pdHandles.append( 

241 InMemoryDatasetHandle( 

242 pdCalib2, 

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

244 ) 

245 ) 

246 idCounter += 2 

247 

248 # Divide out the pairs. 

249 outputCovariances = [] 

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

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

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

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

254 results = extractPairTask.run( 

255 inputExp=inputExp, 

256 inputDims=inputDims, 

257 inputPhotodiodeData=inputPhotodiodeData, 

258 ) 

259 outputCovariances.append(results.outputCovariance) 

260 

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

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

263 # affect the output comparison. 

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

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

266 

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

268 # ensure that the outlier fit rejection works. 

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

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

271 

272 # Reorganize the outputCovariances so we can confirm they come 

273 # out sorted afterwards. 

274 outputCovariancesRev = outputCovariances[::-1] 

275 

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

277 # it was calculated. 

278 expectedNoiseMatrix = { 

279 "FULLCOVARIANCE": np.array( 

280 [[8.99474598, 9.94916264, -27.90587299], 

281 [-2.95079527, -17.11827641, -47.88156244], 

282 [5.24915021, -3.25786165, 26.09634067]], 

283 ), 

284 "FULLCOVARIANCE_NO_B": np.array( 

285 [[8.71049338, 12.48584043, -37.06585088], 

286 [-4.80523971, -23.29102809, -66.37815343], 

287 [7.48654766, -4.10168337, 35.64469824]], 

288 ), 

289 } 

290 

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

292 solveConfig.ptcFitType = fitType 

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

294 

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

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

297 resultsSolve = solveTask.run( 

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

299 ) 

300 

301 ptc = resultsSolve.outputPtcDataset 

302 

303 # Confirm that metadata keywords are set. 

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

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

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

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

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

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

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

311 

312 for amp in self.ampNames: 

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

314 self.assertFloatsAlmostEqual( 

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

316 1.0, 

317 rtol=1e-4, 

318 ) 

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

320 

321 # Check that the PTC turnoff is correctly computed. 

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

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

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

325 else: 

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

327 

328 # Test that all the quantities are correctly ordered and 

329 # have not accidentally been masked. 

330 for i, extractPtc in enumerate(outputCovariances): 

331 self.assertFloatsAlmostEqual( 

332 extractPtc.rawExpTimes[ampName][0], 

333 ptc.rawExpTimes[ampName][i], 

334 ) 

335 self.assertFloatsAlmostEqual( 

336 extractPtc.rawMeans[ampName][0], 

337 ptc.rawMeans[ampName][i], 

338 ) 

339 self.assertFloatsAlmostEqual( 

340 extractPtc.rawVars[ampName][0], 

341 ptc.rawVars[ampName][i], 

342 ) 

343 self.assertFloatsAlmostEqual( 

344 extractPtc.rawDeltas[ampName][0], 

345 ptc.rawDeltas[ampName][i], 

346 ) 

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

348 self.assertFloatsNotEqual( 

349 extractPtc.rawDeltas[ampName][0], 

350 0.0, 

351 ) 

352 self.assertFloatsAlmostEqual( 

353 extractPtc.photoCharges[ampName][0], 

354 ptc.photoCharges[ampName][i], 

355 ) 

356 self.assertFloatsAlmostEqual( 

357 extractPtc.photoChargeDeltas[ampName][0], 

358 ptc.photoChargeDeltas[ampName][i], 

359 ) 

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

361 self.assertFloatsNotEqual( 

362 extractPtc.photoChargeDeltas[ampName][0], 

363 0.0, 

364 ) 

365 self.assertFloatsAlmostEqual( 

366 extractPtc.histVars[ampName][0], 

367 ptc.histVars[ampName][i], 

368 ) 

369 self.assertFloatsAlmostEqual( 

370 extractPtc.histChi2Dofs[ampName][0], 

371 ptc.histChi2Dofs[ampName][i], 

372 ) 

373 self.assertFloatsAlmostEqual( 

374 extractPtc.kspValues[ampName][0], 

375 ptc.kspValues[ampName][i], 

376 ) 

377 self.assertFloatsAlmostEqual( 

378 extractPtc.covariances[ampName][0], 

379 ptc.covariances[ampName][i], 

380 ) 

381 self.assertFloatsAlmostEqual( 

382 extractPtc.covariancesSqrtWeights[ampName][0], 

383 ptc.covariancesSqrtWeights[ampName][i], 

384 ) 

385 self.assertFloatsAlmostEqual( 

386 ptc.noiseMatrix[ampName], 

387 expectedNoiseMatrix[fitType], 

388 atol=1e-8, 

389 rtol=None, 

390 ) 

391 self.assertFloatsAlmostEqual( 

392 ptc.ampOffsets[ampName], 

393 0.0, 

394 ) 

395 self.assertFloatsAlmostEqual( 

396 ptc.noise[ampName], 

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

398 rtol=0.05, 

399 ) 

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

401 # something is seriously wrong. Possibly some 

402 # kind of gain application mismatch. 

403 self.assertLess( 

404 ptc.noiseErr[ampName], 

405 ptc.noise[ampName], 

406 ) 

407 

408 mask = ptc.getGoodPoints(amp) 

409 

410 values = ( 

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

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

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

414 

415 if ptc.ptcFitType == "FULLCOVARIANCE": 

416 values = ( 

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

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

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

420 

421 values = ( 

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

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

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

425 

426 # And test that the auxiliary values are there and 

427 # correctly ordered. 

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

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

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

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

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

433 

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

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

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

437 self.assertTrue( 

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

439 ) 

440 

441 goodAmps = ptc.getGoodAmps() 

442 self.assertEqual(goodAmps, self.ampNames) 

443 

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

445 covShape = None 

446 covSqrtShape = None 

447 covModelShape = None 

448 

449 for ampName in self.ampNames: 

450 if covShape is None: 

451 covShape = ptc.covariances[ampName].shape 

452 covSqrtShape = ptc.covariancesSqrtWeights[ampName].shape 

453 covModelShape = ptc.covariancesModel[ampName].shape 

454 else: 

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

456 self.assertEqual( 

457 ptc.covariancesSqrtWeights[ampName].shape, covSqrtShape 

458 ) 

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

460 

461 # Check if evalPtcModel produces expected values 

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

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

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

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

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

467 

468 # And check that this is serializable 

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

470 usedFilename = ptc.writeFits(f.name) 

471 fromFits = PhotonTransferCurveDataset.readFits(usedFilename) 

472 self.assertEqual(fromFits, ptc) 

473 

474 def ptcFitAndCheckPtc( 

475 self, 

476 fitType=None, 

477 doFitBootstrap=False, 

478 ): 

479 localDataset = copy.deepcopy(self.dataset) 

480 localDataset.ptcFitType = fitType 

481 configSolve = copy.copy(self.defaultConfigSolve) 

482 if doFitBootstrap: 

483 configSolve.doFitBootstrap = True 

484 

485 if fitType == "EXPAPPROXIMATION": 

486 g = self.gain 

487 for ampName in self.ampNames: 

488 localDataset.rawVars[ampName] = [ 

489 ( 

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

491 + self.noiseSq / (g * g) 

492 ) 

493 for mu in localDataset.rawMeans[ampName] 

494 ] 

495 else: 

496 raise RuntimeError( 

497 "Fit type must be 'EXPAPPROXIMATION'" 

498 ) 

499 

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

501 # Covariance weights values empirically determined from one of 

502 # the cases in test_covAstier. 

503 matrixSize = localDataset.covMatrixSide 

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

505 for ampName in self.ampNames: 

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

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

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

509 ).reshape((maskLength, matrixSize, matrixSize)) 

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

511 0.07980188, 

512 0.01339653, 

513 0.0073118, 

514 0.00502802, 

515 0.00383132, 

516 0.00309475, 

517 0.00259572, 

518 0.00223528, 

519 0.00196273, 

520 0.00174943, 

521 0.00157794, 

522 0.00143707, 

523 0.00131929, 

524 0.00121935, 

525 0.0011334, 

526 0.00105893, 

527 0.00099357, 

528 0.0009358, 

529 0.00088439, 

530 0.00083833, 

531 ] 

532 

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

534 

535 if fitType == "EXPAPPROXIMATION": 

536 localDataset = solveTask.fitPtc(localDataset) 

537 else: 

538 localDataset = solveTask.fitDataFullCovariance(localDataset) 

539 

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

541 for ampName in self.ampNames: 

542 self.assertEqual(fitType, localDataset.ptcFitType) 

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

544 if fitType == "EXPAPPROXIMATION": 

545 self.assertAlmostEqual( 

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

547 ) 

548 # Noise already in electrons 

549 self.assertAlmostEqual( 

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

551 ) 

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

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

554 # mismatch. 

555 self.assertLess( 

556 localDataset.noiseErr[ampName], 

557 np.sqrt(self.noiseSq) 

558 ) 

559 

560 # Check if evalColModel produces expected values 

561 if not doFitBootstrap: 

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

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

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

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

566 self.assertFloatsEqual(evalVarModel, model) 

567 

568 def test_lsstcam_samples(self): 

569 for dense in [False, True]: 

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

571 for doModelPtcRolloff in [False, True]: 

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

573 

574 # We only need a single amp. 

575 ampName = self.ampNames[0] 

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

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

578 dataset.rawMeans[ampName] = rawMeans 

579 dataset.rawVars[ampName] = rawVars 

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

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

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

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

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

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

586 dataset.expIdMask[ampName] = mask 

587 

588 configSolve = copy.copy(self.defaultConfigSolve) 

589 configSolve.doModelPtcRolloff = doModelPtcRolloff 

590 configSolve.doFitBootstrap = False 

591 

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

593 

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

595 

596 # Model the PTC rolloff 

597 if doModelPtcRolloff: 

598 initialFitDataset = solveTask.fitPtcRolloff(initialFitDataset) 

599 

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

601 # Check that the ptcTurnoff is what is expected. 

602 self.assertFloatsAlmostEqual( 

603 solvedDataset.ptcTurnoff[ampName], 

604 ptcTurnoff, 

605 msg=( 

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

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

608 ), 

609 ) 

610 

611 # Check that no values above the 

612 # turnoff/rolloff are "good". 

613 maxMean = solvedDataset.ptcTurnoff[ampName] 

614 if doModelPtcRolloff: 

615 # Check that the PTC rolloff is less 

616 # than the PTC turnoff and within 10% 

617 # of the turnoff. 

618 if dense: 

619 # Check that 

620 self.assertLess( 

621 solvedDataset.ptcRolloff[ampName], 

622 solvedDataset.ptcTurnoff[ampName], 

623 msg=( 

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

625 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

626 f"Amp: {ampName}" 

627 ), 

628 ) 

629 self.assertFloatsAlmostEqual( 

630 solvedDataset.ptcRolloff[ampName], 

631 solvedDataset.ptcTurnoff[ampName], 

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

633 msg=( 

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

635 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

636 f"Amp: {ampName}" 

637 ), 

638 ) 

639 else: 

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

641 # find points beyond the turnoff, and the 

642 # rolloff should be set to the PTC turnoff. 

643 self.assertEqual( 

644 solvedDataset.ptcRolloff[ampName], 

645 solvedDataset.ptcTurnoff[ampName], 

646 msg=( 

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

648 f"doModelPtcRolloff: {doModelPtcRolloff}; " 

649 f"Amp: {ampName}" 

650 ), 

651 ) 

652 

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

654 maxMean = solvedDataset.ptcRolloff[ampName] 

655 

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

657 self.assertEqual( 

658 np.sum(above), 

659 0, 

660 msg=( 

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

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

663 ), 

664 ) 

665 

666 # Check the sampling error on the turnoff. 

667 turnoffIdx = np.argwhere( 

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

669 ) 

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

671 self.assertFloatsAlmostEqual( 

672 solvedDataset.ptcTurnoffSamplingError[ampName], 

673 samplingError, 

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

675 ) 

676 

677 def test_ptcFit(self): 

678 for fitType in ["EXPAPPROXIMATION"]: 

679 self.ptcFitAndCheckPtc( 

680 fitType=fitType 

681 ) 

682 

683 def test_meanVarMeasurement(self): 

684 task = self.defaultTaskExtract 

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

686 self.flatExp1, self.flatExp2 

687 ) 

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

689 

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

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

692 

693 def test_meanVarMeasurementWithNans(self): 

694 task = self.defaultTaskExtract 

695 

696 flatExp1 = self.flatExp1.clone() 

697 flatExp2 = self.flatExp2.clone() 

698 

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

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

701 

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

703 flatExp1, flatExp2 

704 ) 

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

706 

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

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

709 expectedMu = 0.5 * (expectedMu1 + expectedMu2) 

710 

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

712 im1 = flatExp1.maskedImage 

713 im2 = flatExp2.maskedImage 

714 

715 temp = im2.clone() 

716 temp *= expectedMu1 

717 diffIm = im1.clone() 

718 diffIm *= expectedMu2 

719 diffIm -= temp 

720 diffIm /= expectedMu 

721 

722 # Divide by two as it is what measureMeanVarCov returns 

723 # (variance of difference) 

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

725 

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

727 # less than 1 ADU 

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

729 self.assertLess(expectedMu - mu, 1) 

730 

731 def test_meanVarMeasurementAllNan(self): 

732 task = self.defaultTaskExtract 

733 flatExp1 = self.flatExp1.clone() 

734 flatExp2 = self.flatExp2.clone() 

735 

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

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

738 

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

740 flatExp1, flatExp2 

741 ) 

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

743 im1Area, im2Area, imStatsCtrl, mu1, mu2