Coverage for tests/test_ptcAdjustGainRatios.py: 8%

2# LSST Data Management System

6# This product includes software developed by the

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

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

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

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

12# (at your option) any later version.

13#

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

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

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

17# GNU General Public License for more details.

18#

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

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

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

22#

23"""Test cases for cp_pipe gain adjustment code."""

25import copy

26import logging

27import unittest

28import numpy as np

29from scipy.interpolate import Akima1DInterpolator

31import lsst.utils.tests

33import lsst.afw.cameraGeom

34from lsst.afw.image import ExposureF

35from lsst.cp.pipe.ptc import PhotonTransferCurveAdjustGainRatiosTask

36from lsst.ip.isr import IsrMockLSST, PhotonTransferCurveDataset

37from lsst.pipe.base import InMemoryDatasetHandle

40class PtcAdjustGainRatiosTestCase(lsst.utils.tests.TestCase):

41 def setUp(self):

42 mock = IsrMockLSST()

43 self.camera = mock.getCamera()

44 self.detector = self.camera[20]

45 self.amp_names = [amp.getName() for amp in self.detector]

47 self.n_radial_nodes = 10

49 rng = np.random.RandomState(seed=12345)

51 self.gain_true = {}

52 for amp_name in self.amp_names:

53 self.gain_true[amp_name] = float(rng.normal(loc=1.5, scale=0.05, size=1)[0])

55 def _get_adjusted_flat(

56 self,

57 normalization,

58 radial_gradient_outer_value,

59 planar_gradient_x,

60 planar_gradient_y,

61 bad_amps=[],

62 weird_amps=[],

63 ):

64 """Get an adjusted flat.

66 Parameters

67 ----------

68 normalization : `float`

69 Normalization value.

70 radial_gradient_outer_value : `float`

71 Value at max radius (should be < 1.0).

72 planar_gradient_x : `float`

73 Planar gradient x slope.

74 planar_gradient_y : `float`

75 Planar gradient y slope.

76 bad_amps : `list` [`str`], optional

77 List of bad amp names.

79 Returns

80 -------

81 flat : `lsst.afw.image.Exposure`

82 Flat image generated.

83 radial_gradient : `np.ndarray`

84 Radial gradient image applied.

85 planar_gradient : `np.ndarray`

86 Planar gradient image applied.

87 """

88 flat = ExposureF(self.detector.getBBox())

89 flat.setDetector(self.detector)

91 xx = np.arange(flat.image.array.shape[1], dtype=np.float64)

92 yy = np.arange(flat.image.array.shape[0], dtype=np.float64)

93 x, y = np.meshgrid(xx, yy)

94 x = x.ravel()

95 y = y.ravel()

97 transform = self.detector.getTransform(lsst.afw.cameraGeom.PIXELS, lsst.afw.cameraGeom.FOCAL_PLANE)

98 xy = np.vstack((x, y))

99 xf, yf = np.vsplit(transform.getMapping().applyForward(xy), 2)

100 xf = xf.ravel()

101 yf = yf.ravel()

102

103 radius = np.sqrt(xf**2. + yf**2.)

104

105 nodes = np.linspace(radius.min(), radius.max(), self.n_radial_nodes)

106 spline_values = np.linspace(1.0, 0.98, self.n_radial_nodes)

107

108 spl = Akima1DInterpolator(nodes, spline_values, method="akima")

109 radial_gradient = spl(radius).reshape(flat.image.array.shape)

110 flat.image.array[:, :] = radial_gradient.copy()

111

112 # This will be a parameter

113 normalization = 10000.0

114

115 flat.image.array[:, :] *= normalization

116

117 # Add a bit of a planar gradient; approx +/- 1%.

118 planar_gradient = (

119 1

120 - 0.00005 * (x - self.detector.getBBox().getCenter().getX())

121 + 0.00001 * (y - self.detector.getBBox().getCenter().getY())

122 ).reshape(flat.image.array.shape)

123

124 flat.image.array[:, :] *= planar_gradient

125

126 for amp_name in self.amp_names:

127 flat.image[self.detector[amp_name].getBBox()].array[:, :] /= self.gain_true[amp_name]

128

129 # Make the bad amps goofy.

130 for bad_amp in bad_amps:

131 flat.image[self.detector[bad_amp].getBBox()].array[:, :] *= 0.1

132

133 return flat, radial_gradient, planar_gradient

134

135 def _gain_correct_flat(self, flat, ptc, adjust_ratios):

136 """Gain correct a flat to a biased and debiased version.

137

138 Parameters

139 ----------

140 flat : `lsst.afw.image.Exposure`

141 Flat to correct.

142 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`

143 PTC dataset with unadjusted gains.

144 adjust_ratios : `np.ndarray`

145 Array of adjustment ratios.

146

147 Returns

148 -------

149 flat_biased : `lsst.afw.image.Exposure`

150 Flat with the original unadjusted gains.

151 flat_debiased : `lsst.afw.image.Exposure`

152 Flat with adjusted gains.

153 """

154 flat_biased = flat.clone()

155 flat_debiased = flat.clone()

156 for i, amp_name in enumerate(self.amp_names):

157 bbox = self.detector[amp_name].getBBox()

158 if amp_name in ptc.badAmps:

159 flat_biased[bbox].image.array[:, :] = np.nan

160 flat_debiased[bbox].image.array[:, :] = np.nan

161 continue

162

163 flat_biased[bbox].image.array[:, :] *= ptc.gainUnadjusted[amp_name]

164 flat_debiased[bbox].image.array[:, :] *= (ptc.gainUnadjusted[amp_name] / adjust_ratios[i])

165

166 return flat_biased, flat_debiased

167

168 def test_compute_gain_ratios_small_gradient(self):

169 fixed_amp_num = 2

170

171 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

172 10000.0,

173 0.98,

174 -0.00005,

175 0.00001,

176 )

177

178 ptc = PhotonTransferCurveDataset()

179 ptc.ampNames = self.amp_names

180 ptc.badAmps = []

181

182 gain_biased = copy.copy(self.gain_true)

183 gain_biased[self.amp_names[0]] *= 1.002

184 gain_biased[self.amp_names[4]] *= 0.998

185

186 # These are the gains recorded in the PTC.

187 for amp_name in self.amp_names:

188 ptc.gain[amp_name] = gain_biased[amp_name]

189 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

190

191 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

192 config.bin_factor = 2

193 config.n_flat = 10

194 config.amp_boundary = 0

195 config.radial_gradient_n_spline_nodes = 3

196

197 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

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

199 adjust_ratios = task._compute_gain_ratios(ptc, flat.clone(), fixed_amp_num)

200

201 # Make two flats with the old adjustment and the new adjustment.

202 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

203

204 # Remove the gradients...

205 flat_biased.image.array[:, :] /= radial_gradient

206 flat_biased.image.array[:, :] /= planar_gradient

207 flat_debiased.image.array[:, :] /= radial_gradient

208 flat_debiased.image.array[:, :] /= planar_gradient

209

210 # Insist that the standard deviation decreases after

211 # removing gradients.

212 self.assertLess(

213 np.std(flat_debiased.image.array.ravel()),

214 np.std(flat_biased.image.array.ravel()),

215 )

216

217 def test_compute_gain_ratios_large_gradient(self):

218 fixed_amp_num = 2

219

220 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

221 10000.0,

222 0.6,

223 -0.00005,

224 0.00001,

225 )

226

227 ptc = PhotonTransferCurveDataset()

228 ptc.ampNames = self.amp_names

229 ptc.badAmps = []

230

231 gain_biased = copy.copy(self.gain_true)

232 gain_biased[self.amp_names[0]] *= 1.002

233 gain_biased[self.amp_names[4]] *= 0.998

234

235 # These are the gains recorded in the PTC.

236 for amp_name in self.amp_names:

237 ptc.gain[amp_name] = gain_biased[amp_name]

238 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

239

240 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

241 config.bin_factor = 2

242 config.n_flat = 10

243 config.amp_boundary = 0

244 config.radial_gradient_n_spline_nodes = 3

245

246 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

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

248 adjust_ratios = task._compute_gain_ratios(ptc, flat.clone(), fixed_amp_num)

249

250 # Make two flats with the old adjustment and the new adjustment.

251 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

252

253 # Remove the gradients...

254 flat_biased.image.array[:, :] /= radial_gradient

255 flat_biased.image.array[:, :] /= planar_gradient

256 flat_debiased.image.array[:, :] /= radial_gradient

257 flat_debiased.image.array[:, :] /= planar_gradient

258

259 # Insist that the standard deviation decreases after

260 # removing gradients.

261 self.assertLess(

262 np.std(flat_debiased.image.array.ravel()),

263 np.std(flat_biased.image.array.ravel()),

264 )

265

266 def test_compute_gain_ratios_bad_amp(self):

267 fixed_amp_num = 2

268

269 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

270 10000.0,

271 0.6,

272 -0.00005,

273 0.00001,

274 bad_amps=[self.amp_names[1]],

275 )

276

277 ptc = PhotonTransferCurveDataset()

278 ptc.ampNames = self.amp_names

279 ptc.badAmps = [self.amp_names[1]]

280

281 gain_biased = copy.copy(self.gain_true)

282 gain_biased[self.amp_names[0]] *= 1.002

283 gain_biased[self.amp_names[4]] *= 0.998

284

285 # These are the gains recorded in the PTC.

286 for amp_name in self.amp_names:

287 ptc.gain[amp_name] = gain_biased[amp_name]

288 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

289

290 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

291 config.bin_factor = 2

292 config.n_flat = 10

293 config.amp_boundary = 0

294 config.radial_gradient_n_spline_nodes = 3

295

296 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

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

298 adjust_ratios = task._compute_gain_ratios(ptc, flat.clone(), fixed_amp_num)

299

300 # Make two flats with the old adjustment and the new adjustment.

301 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

302

303 # Remove the gradients...

304 flat_biased.image.array[:, :] /= radial_gradient

305 flat_biased.image.array[:, :] /= planar_gradient

306 flat_debiased.image.array[:, :] /= radial_gradient

307 flat_debiased.image.array[:, :] /= planar_gradient

308

309 # Insist that the standard deviation decreases after

310 # removing gradients.

311 # Note that this contains nans which must be filtered.

312 self.assertLess(

313 np.nanstd(flat_debiased.image.array.ravel()),

314 np.nanstd(flat_biased.image.array.ravel()),

315 )

316

317 def test_compute_gain_ratios_weird_amp(self):

318 fixed_amp_num = 2

319

320 # The weird amps have a large offset but are not known

321 # as bad.

322 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

323 10000.0,

324 0.6,

325 -0.00005,

326 0.00001,

327 bad_amps=[self.amp_names[1]],

328 )

329

330 ptc = PhotonTransferCurveDataset()

331 ptc.ampNames = self.amp_names

332

333 gain_biased = copy.copy(self.gain_true)

334 gain_biased[self.amp_names[0]] *= 1.002

335 gain_biased[self.amp_names[4]] *= 0.998

336

337 # These are the gains recorded in the PTC.

338 for amp_name in self.amp_names:

339 ptc.gain[amp_name] = gain_biased[amp_name]

340 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

341

342 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

343 config.bin_factor = 2

344 config.n_flat = 10

345 config.amp_boundary = 0

346 config.radial_gradient_n_spline_nodes = 3

347

348 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

349 with self.assertLogs(level=logging.WARNING) as cm:

350 adjust_ratios = task._compute_gain_ratios(ptc, flat.clone(), fixed_amp_num)

351 self.assertIn("Found bad amp", cm.output[0])

352

353 np.testing.assert_array_less(np.abs(adjust_ratios - 1.0), config.max_fractional_gain_ratio)

354

355 # Make two flats with the old adjustment and the new adjustment.

356 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

357

358 # Remove the gradients...

359 flat_biased.image.array[:, :] /= radial_gradient

360 flat_biased.image.array[:, :] /= planar_gradient

361 flat_debiased.image.array[:, :] /= radial_gradient

362 flat_debiased.image.array[:, :] /= planar_gradient

363

364 # Insist that the standard deviation decreases after

365 # removing gradients.

366 # Note that this contains nans which must be filtered.

367 self.assertLess(

368 np.nanstd(flat_debiased.image.array.ravel()),

369 np.nanstd(flat_biased.image.array.ravel()),

370 )

371

372 def test_task_run(self):

373 levels = np.linspace(500.0, 25000.0, 40)

374

375 flat_handles = []

376 exp_id_pairs = []

377 for i, level in enumerate(levels):

378 exp_id_1 = 100 + 2*i

379 exp_id_2 = exp_id_1 + 1

380

381 data_id_1 = {

382 "detector": self.detector.getId(),

383 "exposure": exp_id_1,

384 }

385 data_id_2 = {

386 "detector": self.detector.getId(),

387 "exposure": exp_id_2,

388 }

389

390 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

391 levels[i],

392 0.6,

393 -0.00005,

394 0.00001,

395 )

396 flat_handles.append(InMemoryDatasetHandle(flat.clone(), dataId=data_id_1))

397 flat_handles.append(InMemoryDatasetHandle(flat.clone(), dataId=data_id_2))

398

399 exp_id_pairs.append([exp_id_1, exp_id_2])

400

401 ptc = PhotonTransferCurveDataset()

402 ptc.ampNames = self.amp_names

403 ptc.badAmps = []

404

405 gain_biased = copy.copy(self.gain_true)

406 gain_biased[self.amp_names[0]] *= 1.002

407 gain_biased[self.amp_names[4]] *= 0.998

408

409 # These are the gains recorded in the PTC.

410 for amp_name in self.amp_names:

411 ptc.gain[amp_name] = gain_biased[amp_name]

412 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

413

414 ptc.rawMeans[amp_name] = levels

415 ptc.finalMeans[amp_name] = levels

416

417 ptc.inputExpIdPairs[amp_name] = exp_id_pairs

418

419 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

420 config.bin_factor = 2

421 config.n_flat = 10

422 config.amp_boundary = 0

423 config.radial_gradient_n_spline_nodes = 3

424 config.n_flat = 10

425

426 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

427 struct = task.run(exposures=flat_handles, input_ptc=ptc)

428

429 output_ptc = struct.output_ptc

430 summary = struct.gain_adjust_summary

431

432 self.assertGreater(summary["mean_adu"].min(), config.min_adu)

433 self.assertLess(summary["mean_adu"].max(), config.max_adu)

434

435 self.assertIn("fixed_amp_index", summary.meta.keys())

436 self.assertIn("fixed_amp_name", summary.meta.keys())

437 self.assertIn("median_correction", summary.meta.keys())

438

439 # In this idealized case, the summary values should be the same

440 # for every flux level.

441 adjust_ratios = np.zeros(len(self.amp_names))

442 for i, amp_name in enumerate(self.amp_names):

443 np.testing.assert_array_equal(

444 summary[f"{amp_name}_gain_ratio"][0],

445 summary[f"{amp_name}_gain_ratio"],

446 )

447 adjust_ratios[i] = output_ptc.gainUnadjusted[amp_name] / output_ptc.gain[amp_name]

448 self.assertFloatsAlmostEqual(

449 adjust_ratios[i],

450 summary[f"{amp_name}_gain_ratio"] / summary.meta["median_correction"],

451 )

452

453 # Confirm this actually works.

454 flat = flat_handles[0].get()

455 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

456

457 # Remove the gradients...

458 flat_biased.image.array[:, :] /= radial_gradient

459 flat_biased.image.array[:, :] /= planar_gradient

460 flat_debiased.image.array[:, :] /= radial_gradient

461 flat_debiased.image.array[:, :] /= planar_gradient

462

463 # Insist that the standard deviation decreases after

464 # removing gradients.

465 self.assertLess(

466 np.std(flat_debiased.image.array.ravel()),

467 np.std(flat_biased.image.array.ravel()),

468 )

469

470 def test_task_run_bad_amp(self):

471 levels = np.linspace(500.0, 25000.0, 40)

472

473 bad_amps = [self.amp_names[1]]

474

475 flat_handles = []

476 exp_id_pairs = []

477 for i, level in enumerate(levels):

478 exp_id_1 = 100 + 2*i

479 exp_id_2 = exp_id_1 + 1

480

481 data_id_1 = {

482 "detector": self.detector.getId(),

483 "exposure": exp_id_1,

484 }

485 data_id_2 = {

486 "detector": self.detector.getId(),

487 "exposure": exp_id_2,

488 }

489

490 flat, radial_gradient, planar_gradient = self._get_adjusted_flat(

491 levels[i],

492 0.6,

493 -0.00005,

494 0.00001,

495 bad_amps=bad_amps,

496 )

497 flat_handles.append(InMemoryDatasetHandle(flat.clone(), dataId=data_id_1))

498 flat_handles.append(InMemoryDatasetHandle(flat.clone(), dataId=data_id_2))

499

500 exp_id_pairs.append([exp_id_1, exp_id_2])

501

502 ptc = PhotonTransferCurveDataset()

503 ptc.ampNames = self.amp_names

504 ptc.badAmps = bad_amps

505

506 gain_biased = copy.copy(self.gain_true)

507 gain_biased[self.amp_names[0]] *= 1.002

508 gain_biased[self.amp_names[4]] *= 0.998

509

510 # These are the gains recorded in the PTC.

511 for amp_name in self.amp_names:

512 ptc.gain[amp_name] = gain_biased[amp_name]

513 ptc.gainUnadjusted[amp_name] = gain_biased[amp_name]

514

515 ptc.rawMeans[amp_name] = levels

516 ptc.finalMeans[amp_name] = levels

517

518 ptc.inputExpIdPairs[amp_name] = exp_id_pairs

519

520 config = PhotonTransferCurveAdjustGainRatiosTask.ConfigClass()

521 config.bin_factor = 2

522 config.n_flat = 10

523 config.amp_boundary = 0

524 config.radial_gradient_n_spline_nodes = 3

525 # Make sure the code works when we have fewer flats.

526 config.n_flat = 100

527

528 task = PhotonTransferCurveAdjustGainRatiosTask(config=config)

529 struct = task.run(exposures=flat_handles, input_ptc=ptc)

530

531 output_ptc = struct.output_ptc

532 summary = struct.gain_adjust_summary

533

534 self.assertGreater(summary["mean_adu"].min(), config.min_adu)

535 self.assertLess(summary["mean_adu"].max(), config.max_adu)

536

537 self.assertIn("fixed_amp_index", summary.meta.keys())

538 self.assertIn("fixed_amp_name", summary.meta.keys())

539 self.assertIn("median_correction", summary.meta.keys())

540

541 # In this idealized case, the summary values should be the same

542 # for every flux level.

543 adjust_ratios = np.zeros(len(self.amp_names))

544 for i, amp_name in enumerate(self.amp_names):

545 np.testing.assert_array_equal(

546 summary[f"{amp_name}_gain_ratio"][0],

547 summary[f"{amp_name}_gain_ratio"],

548 )

549 adjust_ratios[i] = output_ptc.gainUnadjusted[amp_name] / output_ptc.gain[amp_name]

550 self.assertFloatsAlmostEqual(

551 adjust_ratios[i],

552 summary[f"{amp_name}_gain_ratio"] / summary.meta["median_correction"],

553 )

554

555 # Confirm this actually works.

556 flat = flat_handles[0].get()

557 flat_biased, flat_debiased = self._gain_correct_flat(flat, ptc, adjust_ratios)

558

559 # Remove the gradients...

560 flat_biased.image.array[:, :] /= radial_gradient

561 flat_biased.image.array[:, :] /= planar_gradient

562 flat_debiased.image.array[:, :] /= radial_gradient

563 flat_debiased.image.array[:, :] /= planar_gradient

564

565 # Insist that the standard deviation decreases after

566 # removing gradients.

567 self.assertLess(

568 np.nanstd(flat_debiased.image.array.ravel()),

569 np.nanstd(flat_biased.image.array.ravel()),

570 )

571

572

573class TestMemory(lsst.utils.tests.MemoryTestCase):

574 pass

575

576

577def setup_module(module):

578 lsst.utils.tests.init()

579

580

581if __name__ == "__main__": 581 ↛ 582line 581 didn't jump to line 582 because the condition on line 581 was never true

582 lsst.utils.tests.init()

583 unittest.main()