Coverage for tests / test_linearity.py: 7%
657 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 08:44 +0000
1#!/usr/bin/env python
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 linearity code."""
27from astropy.table import Table
28import logging
29import unittest
30import numpy as np
31import copy
32from scipy.interpolate import Akima1DInterpolator
34import lsst.utils
35import lsst.utils.tests
37from lsst.ip.isr import PhotonTransferCurveDataset, IsrMockLSST
39import lsst.afw.image
40import lsst.afw.math
41from lsst.cp.pipe import LinearitySolveTask, LinearityNormalizeTask, LinearityDoubleSplineSolveTask
42from lsst.cp.pipe.cpLinearitySolve import _computeTurnoffAndMax, _noderator, _determineInputGroups
43from lsst.cp.pipe.ptc import PhotonTransferCurveSolveTask, PhotonTransferCurveExtractPairTask
44from lsst.cp.pipe.utils import funcPolynomial, funcAstier
45from lsst.daf.base import DateTime
46from lsst.ip.isr.isrMock import FlatMock, IsrMock
47from lsst.pipe.base import InMemoryDatasetHandle
50class FakeCamera(list):
51 def getName(self):
52 return "FakeCam"
55class LinearityTaskTestCase(lsst.utils.tests.TestCase):
56 """Test case for the linearity tasks."""
58 def setUp(self):
59 mock_image_config = IsrMock.ConfigClass()
60 mock_image_config.flatDrop = 0.99999
61 mock_image_config.isTrimmed = True
63 self.dummy_exposure = FlatMock(config=mock_image_config).run()
64 self.detector = self.dummy_exposure.getDetector()
65 self.input_dims = {"detector": 0}
67 self.camera = FakeCamera([self.detector])
69 self.amp_names = []
70 for amp in self.detector:
71 self.amp_names.append(amp.getName())
73 self.sequencerMetadata = {
74 "SEQNAME": "a_sequencer",
75 "SEQFILE": "a_sequencer_file",
76 "SEQCKSUM": "deadbeef",
77 }
79 def _create_ptc(
80 self,
81 amp_names,
82 exp_times,
83 means,
84 ccobcurr=None,
85 photo_charges=None,
86 temperatures=None,
87 ptc_turnoff=None,
88 ):
89 """
90 Create a PTC with values for linearity tests.
92 Parameters
93 ----------
94 amp_names : `list` [`str`]
95 Names of amps.
96 exp_times : `np.ndarray`
97 Array of exposure times.
98 means : `np.ndarray`
99 Array of means.
100 ccobcurr : `np.ndarray`, optional
101 Array of CCOBCURR to put into auxiliary values.
102 photo_charges : `np.ndarray`, optional
103 Array of photoCharges to put into ptc.
104 temperatures : `np.ndarray`, optional
105 Array of temperatures (TEMP6) to put into ptc.
106 ptc_turnoff : `float`, optional
107 Turnoff value to set (by hand) for testing.
109 Returns
110 -------
111 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
112 PTC filled with relevant values.
113 """
114 exp_id_pairs = np.arange(len(exp_times)*2).reshape((len(exp_times), 2)).tolist()
116 if photo_charges is None:
117 photo_charges = np.full(len(exp_times), np.nan)
119 datasets = []
120 for i in range(len(exp_times)):
121 partial = PhotonTransferCurveDataset(amp_names, ptcFitType="PARTIAL", covMatrixSide=1)
122 for amp_name in amp_names:
123 # For the first amp, we add a few bad points.
124 if amp_name == amp_names[0] and i >= 5 and i < 7:
125 exp_id_mask = False
126 raw_mean = np.nan
127 else:
128 exp_id_mask = True
129 raw_mean = means[i]
131 partial.setAmpValuesPartialDataset(
132 amp_name,
133 inputExpIdPair=exp_id_pairs[i],
134 rawExpTime=exp_times[i],
135 rawMean=raw_mean,
136 rawVar=raw_mean,
137 kspValue=1.0,
138 expIdMask=exp_id_mask,
139 photoCharge=photo_charges[i],
140 overscanMedianLevel=100.0,
141 )
143 aux_dict = {}
144 if ccobcurr is not None:
145 aux_dict["CCOBCURR"] = ccobcurr[i]
146 if temperatures is not None:
147 aux_dict["TEMP6"] = temperatures[i]
149 if aux_dict:
150 partial.setAuxValuesPartialDataset(aux_dict)
152 datasets.append(partial)
154 datasets.append(PhotonTransferCurveDataset(amp_names, ptcFitType="DUMMY"))
156 config = PhotonTransferCurveSolveTask.ConfigClass()
157 config.maximumRangeCovariancesAstier = 1
158 config.maxDeltaInitialPtcOutlierFit = 100_000.0
159 solve_task = PhotonTransferCurveSolveTask(config=config)
160 # Suppress logging here.
161 with self.assertNoLogs(level=logging.CRITICAL):
162 ptc = solve_task.run(datasets).outputPtcDataset
164 # Make the last amp a bad amp.
165 ptc.badAmps = [amp_names[-1]]
167 if ptc_turnoff is not None:
168 for amp_name in amp_names:
169 if amp_name in ptc.badAmps:
170 ptc.ptcTurnoff[amp_name] = np.nan
171 ptc.finalMeans[amp_name][:] = np.nan
172 else:
173 ptc.ptcTurnoff[amp_name] = ptc_turnoff
174 high = (ptc.rawMeans[amp_name] > ptc_turnoff)
175 ptc.expIdMask[amp_name][high] = False
176 ptc.finalMeans[amp_name][high] = np.nan
178 ptc.updateMetadata(**self.sequencerMetadata, setCalibInfo=True)
179 ptc.updateMetadata(camera=self.camera, detector=self.detector)
181 return ptc
183 def _check_linearity(self, linearity_type, min_adu=0.0, max_adu=100000.0):
184 """Run and check linearity.
186 Parameters
187 ----------
188 linearity_type : `str`
189 Must be ``Polynomial``, ``Squared``, or ``LookupTable``.
190 min_adu : `float`, optional
191 Minimum cut on ADU for fit.
192 max_adu : `float`, optional
193 Maximum cut on ADU for fit.
194 """
195 flux = 1000.
196 time_vec = np.arange(1., 101., 5)
197 k2_non_linearity = -5e-6
198 coeff = k2_non_linearity/(flux**2.)
200 mu_vec = flux * time_vec + k2_non_linearity * time_vec**2.
202 ptc = self._create_ptc(self.amp_names, time_vec, mu_vec)
204 config = LinearitySolveTask.ConfigClass()
205 config.linearityType = linearity_type
206 config.minLinearAdu = min_adu
207 config.maxLinearAdu = max_adu
209 task = LinearitySolveTask(config=config)
210 linearizer = task.run(ptc, [self.dummy_exposure], self.camera, self.input_dims).outputLinearizer
212 if linearity_type == "LookupTable":
213 t_max = config.maxLookupTableAdu / flux
214 time_range = np.linspace(0.0, t_max, config.maxLookupTableAdu)
215 signal_ideal = time_range * flux
216 signal_uncorrected = funcPolynomial(np.array([0.0, flux, k2_non_linearity]), time_range)
217 linearizer_table_row = signal_ideal - signal_uncorrected
219 # Skip the last amp which is marked bad.
220 for i, amp_name in enumerate(ptc.ampNames[:-1]):
221 if linearity_type in ["Squared", "Polynomial"]:
222 self.assertFloatsAlmostEqual(linearizer.fitParams[amp_name][0], 0.0, atol=1e-2)
223 self.assertFloatsAlmostEqual(linearizer.fitParams[amp_name][1], 1.0, rtol=1e-5)
224 self.assertFloatsAlmostEqual(linearizer.fitParams[amp_name][2], coeff, rtol=1e-6)
226 if linearity_type == "Polynomial":
227 self.assertFloatsAlmostEqual(linearizer.fitParams[amp_name][3], 0.0)
229 if linearity_type == "Squared":
230 self.assertEqual(len(linearizer.linearityCoeffs[amp_name]), 1)
231 self.assertFloatsAlmostEqual(linearizer.linearityCoeffs[amp_name][0], -coeff, rtol=1e-6)
232 else:
233 self.assertEqual(len(linearizer.linearityCoeffs[amp_name]), 2)
234 self.assertFloatsAlmostEqual(linearizer.linearityCoeffs[amp_name][0], -coeff, rtol=1e-6)
235 self.assertFloatsAlmostEqual(linearizer.linearityCoeffs[amp_name][1], 0.0)
237 else:
238 index = linearizer.linearityCoeffs[amp_name][0]
239 self.assertEqual(index, i)
240 self.assertEqual(len(linearizer.tableData[index, :]), len(linearizer_table_row))
241 self.assertFloatsAlmostEqual(linearizer.tableData[index, :], linearizer_table_row, rtol=1e-4)
243 lin_mask = np.isfinite(linearizer.fitResiduals[amp_name])
244 lin_mask_expected = (mu_vec > min_adu) & (mu_vec < max_adu) & ptc.expIdMask[amp_name]
246 self.assertListEqual(lin_mask.tolist(), lin_mask_expected.tolist())
247 self.assertFloatsAlmostEqual(linearizer.fitResiduals[amp_name][lin_mask], 0.0, atol=1e-2)
249 # If we apply the linearity correction, we should get the true
250 # linear values out.
251 image = lsst.afw.image.ImageF(len(mu_vec), 1)
252 image.array[:, :] = mu_vec
253 lin_func = linearizer.getLinearityTypeByName(linearizer.linearityType[amp_name])
254 lin_func()(
255 image,
256 coeffs=linearizer.linearityCoeffs[amp_name],
257 table=linearizer.tableData,
258 log=None,
259 )
261 linear_signal = flux * time_vec
262 self.assertFloatsAlmostEqual(image.array[0, :] / linear_signal, 1.0, rtol=1e-6)
264 self._check_linearizer_lengths(linearizer)
266 def _check_linearizer_lengths(self, linearizer):
267 # Check that the lengths of all the fields match.
268 lenCoeffs = -1
269 lenParams = -1
270 lenParamsErr = -1
271 lenResiduals = -1
272 lenFit = -1
273 for ampName in linearizer.ampNames:
274 if lenCoeffs < 0:
275 lenCoeffs = len(linearizer.linearityCoeffs[ampName])
276 lenParams = len(linearizer.fitParams[ampName])
277 lenParamsErr = len(linearizer.fitParamsErr[ampName])
278 lenResiduals = len(linearizer.fitResiduals[ampName])
279 lenFit = len(linearizer.linearFit[ampName])
280 else:
281 self.assertEqual(
282 len(linearizer.linearityCoeffs[ampName]),
283 lenCoeffs,
284 msg=f"amp {ampName} linearityCoeffs length mismatch",
285 )
286 self.assertEqual(
287 len(linearizer.fitParams[ampName]),
288 lenParams,
289 msg=f"amp {ampName} fitParams length mismatch",
290 )
291 self.assertEqual(
292 len(linearizer.fitParamsErr[ampName]),
293 lenParamsErr,
294 msg=f"amp {ampName} fitParamsErr length mismatch",
295 )
296 self.assertEqual(
297 len(linearizer.fitResiduals[ampName]),
298 lenResiduals,
299 msg=f"amp {ampName} fitResiduals length mismatch",
300 )
301 self.assertEqual(
302 len(linearizer.linearFit[ampName]),
303 lenFit,
304 msg=f"amp {ampName} linearFit length mismatch",
305 )
307 def test_linearity_polynomial(self):
308 """Test linearity with polynomial fit."""
309 self._check_linearity("Polynomial")
311 def test_linearity_squared(self):
312 """Test linearity with a single order squared solution."""
313 self._check_linearity("Squared")
315 def test_linearity_table(self):
316 """Test linearity with a lookup table solution."""
317 self._check_linearity("LookupTable")
319 def test_linearity_polynomial_aducuts(self):
320 """Test linearity with polynomial and ADU cuts."""
321 self._check_linearity("Polynomial", min_adu=10000.0, max_adu=90000.0)
323 def _check_linearity_spline(
324 self,
325 do_pd_offsets=False,
326 n_points=200,
327 do_mu_offset=False,
328 do_weight_fit=False,
329 do_temperature_fit=False,
330 ):
331 """Check linearity with a spline solution.
333 Parameters
334 ----------
335 do_pd_offsets : `bool`, optional
336 Apply offsets to the photodiode data.
337 do_mu_offset : `bool`, optional
338 Apply constant offset to mu data.
339 do_weight_fit : `bool`, optional
340 Fit the weight parameters?
341 do_temperature_fit : `bool`, optional
342 Apply a temperature dependence and fit it?
343 """
344 np.random.seed(12345)
346 # Create a test dataset representative of real data.
347 pd_values = np.linspace(1e-8, 2e-5, n_points)
348 time_values = pd_values * 1000000.
349 linear_ratio = 5e9
350 mu_linear = linear_ratio * pd_values
352 # Test spline parameters are taken from a test fit to LSSTCam
353 # data, run 7193D, detector 22, amp C00. The exact fit is not
354 # important, but this is only meant to be representative of
355 # the shape of the non-linearity that we see.
357 n_nodes = 10
359 non_lin_spline_nodes = np.linspace(0, mu_linear.max(), n_nodes)
360 non_lin_spline_values = np.array(
361 [0.0, -8.87, 1.46, 1.69, -6.92, -68.23, -78.01, -11.56, 80.26, 185.01]
362 )
364 spl = lsst.afw.math.makeInterpolate(
365 non_lin_spline_nodes,
366 non_lin_spline_values,
367 lsst.afw.math.stringToInterpStyle("AKIMA_SPLINE"),
368 )
370 mu_values = mu_linear + spl.interpolate(mu_linear)
372 # Add a temperature dependence if necessary.
373 if do_temperature_fit:
374 temp_coeff = 0.0006
375 temperatures = np.random.normal(scale=0.5, size=len(mu_values)) - 100.0
377 # We use a negative sign here because we are doing the
378 # opposite of the correction.
379 mu_values *= (1 - temp_coeff*(temperatures - (-100.0)))
380 else:
381 temperatures = None
383 # Add a constant offset if necessary.
384 if do_mu_offset:
385 offset_value = 2.0
386 mu_values += offset_value
387 else:
388 offset_value = 0.0
390 # Add some noise.
391 mu_values += np.random.normal(scale=mu_values, size=len(mu_values)) / 10000.
393 # Add some outlier values.
394 if n_points >= 200:
395 outlier_indices = np.arange(5) + 170
396 else:
397 outlier_indices = []
398 mu_values[outlier_indices] += 200.0
400 # Add some small offsets to the pd_values if requested.
401 pd_values_offset = pd_values.copy()
402 ccobcurr = None
403 if do_pd_offsets:
404 ccobcurr = np.zeros(pd_values.size)
405 n_points_group = n_points//4
406 group0 = np.arange(n_points_group)
407 group1 = np.arange(n_points_group) + n_points_group
408 group2 = np.arange(n_points_group) + 2*n_points_group
409 group3 = np.arange(n_points_group) + 3*n_points_group
410 ccobcurr[group0] = 0.01
411 ccobcurr[group1] = 0.02
412 ccobcurr[group2] = 0.03
413 ccobcurr[group3] = 0.04
415 pd_offset_factors = [0.995, 1.0, 1.005, 1.002]
416 pd_values_offset[group0] *= pd_offset_factors[0]
417 pd_values_offset[group2] *= pd_offset_factors[2]
418 pd_values_offset[group3] *= pd_offset_factors[3]
420 # Add one bad photodiode value, but don't put it at the very
421 # end because that would change the spline node positions
422 # and make comparisons to the "truth" here in the tests
423 # more difficult.
424 pd_values_offset[-2] = np.nan
426 ptc = self._create_ptc(
427 self.amp_names,
428 time_values,
429 mu_values,
430 ccobcurr=ccobcurr,
431 photo_charges=pd_values_offset,
432 temperatures=temperatures,
433 )
435 config = LinearitySolveTask.ConfigClass()
436 config.linearityType = "Spline"
437 config.usePhotodiode = True
438 config.minLinearAdu = 0.0
439 config.splineKnots = n_nodes
440 config.splineGroupingMinPoints = 101
441 config.doSplineFitOffset = do_mu_offset
442 config.doSplineFitWeights = do_weight_fit
443 config.splineFitWeightParsStart = [7.2e-5, 1e-4]
444 config.doSplineFitTemperature = do_temperature_fit
445 config.maxFracLinearityDeviation = 0.05
447 if do_pd_offsets:
448 config.splineGroupingColumn = "CCOBCURR"
450 if do_temperature_fit:
451 config.splineFitTemperatureColumn = "TEMP6"
453 task = LinearitySolveTask(config=config)
454 linearizer = task.run(
455 ptc,
456 [self.dummy_exposure],
457 self.camera,
458 self.input_dims,
459 ).outputLinearizer
461 for key, value in self.sequencerMetadata.items():
462 self.assertEqual(linearizer.metadata[key], value)
464 for key in ["INSTRUME", "DETECTOR", "DET_NAME", "DET_SER"]:
465 self.assertEqual(linearizer.metadata[key], ptc.metadata[key])
467 if do_weight_fit:
468 # These checks currently fail, and weight fitting is not
469 # recommended.
470 return
472 # Skip the last amp which is marked bad.
473 for amp_name in ptc.ampNames[:-1]:
474 # This test data doesn't have a real turnoff.
475 self.assertEqual(linearizer.linearityTurnoff[amp_name], np.nanmax(ptc.rawMeans[amp_name]))
476 self.assertEqual(linearizer.linearityMaxSignal[amp_name], np.nanmax(ptc.rawMeans[amp_name]))
478 lin_mask = np.isfinite(linearizer.fitResiduals[amp_name])
480 # Make sure that non-finite initial values in range are also
481 # masked.
482 check, = np.where(~np.isfinite(ptc.rawMeans[amp_name]))
483 if len(check) > 0:
484 np.testing.assert_array_equal(lin_mask[check], False)
486 # Make sure the outliers are masked.
487 np.testing.assert_array_equal(lin_mask[outlier_indices], False)
489 # Check the turnoff and max values.
490 resid_atol = 4e-4
491 self.assertFloatsAlmostEqual(
492 linearizer.fitResiduals[amp_name][lin_mask] / mu_linear[lin_mask],
493 0.0,
494 atol=resid_atol,
495 )
497 # Loose check on the chi-squared.
498 self.assertLess(linearizer.fitChiSq[amp_name], 2.0)
500 # Check the residual sigma_mad.
501 self.assertLess(linearizer.fitResidualsSigmaMad[amp_name], 1.2e-4)
503 # If we apply the linearity correction, we should get the true
504 # linear values out.
505 image = lsst.afw.image.ImageF(len(mu_values), 1)
506 image.array[:, :] = mu_values
507 lin_func = linearizer.getLinearityTypeByName(linearizer.linearityType[amp_name])
508 lin_func()(
509 image,
510 coeffs=linearizer.linearityCoeffs[amp_name],
511 log=None,
512 )
514 # We scale by the median because of ambiguity in the overall
515 # gain parameter which is not part of the non-linearity.
516 ratio = image.array[0, lin_mask]/mu_linear[lin_mask]
517 # When we have an offset, this test gets a bit confused
518 # mixing truth and offset values.
519 ratio_rtol = 5e-2 if do_mu_offset else 5e-4
520 self.assertFloatsAlmostEqual(
521 ratio / np.median(ratio),
522 1.0,
523 rtol=ratio_rtol,
524 )
526 # Check that the spline parameters recovered are consistent,
527 # with input to some low-grade precision.
528 # The first element should be identically zero.
529 self.assertFloatsEqual(linearizer.linearityCoeffs[amp_name][0], 0.0)
531 # We have two different comparisons here; for the terms that are
532 # |value| < 20 (offset) or |value| > 20 (ratio), to avoid
533 # divide-by-small-number problems. In all cases these are
534 # approximate, and the real test is in the residuals.
535 small = (np.abs(non_lin_spline_values) < 20)
537 spline_atol = 6.0 if do_pd_offsets else 2.0
538 spline_rtol = 0.14 if do_pd_offsets else 0.05
540 self.assertFloatsAlmostEqual(
541 linearizer.linearityCoeffs[amp_name][n_nodes:][small],
542 non_lin_spline_values[small],
543 atol=spline_atol,
544 )
545 self.assertFloatsAlmostEqual(
546 linearizer.linearityCoeffs[amp_name][n_nodes:][~small],
547 non_lin_spline_values[~small],
548 rtol=spline_rtol,
549 )
551 # And check the offsets if they were included.
552 if do_pd_offsets:
553 # The relative scaling is to group 1.
554 fit_offset_factors = linearizer.fitParams[amp_name][1] / linearizer.fitParams[amp_name]
555 extra_pars = 0
556 if do_mu_offset:
557 extra_pars += 1
558 if do_temperature_fit:
559 extra_pars += 1
561 if extra_pars > 0:
562 fit_offset_factors = fit_offset_factors[:-extra_pars]
564 self.assertFloatsAlmostEqual(fit_offset_factors, np.array(pd_offset_factors), rtol=6e-4)
566 # And check if the offset is fit well.
567 fit_offset = None
568 fit_temp_coeff = None
569 if do_mu_offset and do_temperature_fit:
570 fit_offset = linearizer.fitParams[amp_name][-2]
571 fit_temp_coeff = linearizer.fitParams[amp_name][-1]
572 elif do_mu_offset:
573 fit_offset = linearizer.fitParams[amp_name][-1]
574 elif do_temperature_fit:
575 fit_temp_coeff = linearizer.fitParams[amp_name][-1]
577 if fit_offset is not None:
578 self.assertFloatsAlmostEqual(fit_offset, offset_value, rtol=6e-3)
580 if fit_temp_coeff is not None:
581 self.assertFloatsAlmostEqual(fit_temp_coeff, temp_coeff, rtol=2e-2)
583 self._check_linearizer_lengths(linearizer)
585 def test_linearity_spline(self):
586 self._check_linearity_spline(do_pd_offsets=False, do_mu_offset=False)
588 def test_linearity_spline_offsets(self):
589 self._check_linearity_spline(do_pd_offsets=True, do_mu_offset=False)
591 def test_linearity_spline_mu_offset(self):
592 self._check_linearity_spline(do_pd_offsets=True, do_mu_offset=True)
594 def test_linearity_spline_fit_weights(self):
595 self._check_linearity_spline(do_pd_offsets=True, do_mu_offset=True, do_weight_fit=True)
597 def test_linearity_spline_fit_temperature(self):
598 self._check_linearity_spline(do_pd_offsets=True, do_mu_offset=True, do_temperature_fit=True)
600 def test_linearity_spline_offsets_too_few_points(self):
601 with self.assertRaisesRegex(RuntimeError, "too few points"):
602 self._check_linearity_spline(do_pd_offsets=True, n_points=100)
604 def test_linearity_turnoff(self):
605 # Use some real LSSTComCam linearity data to measure the turnoff.
606 abscissa, ordinate, ptc_mask = self._comcam_raw_linearity_data()
608 config = LinearitySolveTask.ConfigClass()
609 task = LinearitySolveTask(config=config)
611 with self.assertNoLogs(level=logging.WARNING):
612 turnoff_index, turnoff, max_signal, _ = _computeTurnoffAndMax(
613 abscissa,
614 ordinate,
615 ptc_mask,
616 np.zeros(len(abscissa)),
617 )
619 # This was visually inspected such that these are reasonable.
620 self.assertEqual(turnoff_index, 90)
621 np.testing.assert_almost_equal(turnoff, 99756.30512572)
622 np.testing.assert_almost_equal(max_signal, 108730.32842316)
624 # Do the linearity fit with these data.
625 ptc = self._create_ptc(
626 self.amp_names,
627 abscissa * 1000000000,
628 ordinate,
629 photo_charges=abscissa,
630 ptc_turnoff=np.max(ordinate[ptc_mask]),
631 )
632 config = LinearitySolveTask.ConfigClass()
633 config.linearityType = "Spline"
634 config.usePhotodiode = True
635 config.minLinearAdu = 30.0
636 config.splineKnots = 10
637 config.doSplineFitOffset = False
638 config.doSplineFitWeights = False
639 config.splineFitWeightParsStart = [7.2e-5, 1e-4]
640 config.doSplineFitTemperature = False
642 task = LinearitySolveTask(config=config)
643 linearizer = task.run(
644 ptc,
645 [self.dummy_exposure],
646 self.camera,
647 self.input_dims,
648 ).outputLinearizer
650 # Confirm that the turnoff for the good amps is the same.
651 for amp_name in self.amp_names:
652 if amp_name in ptc.badAmps:
653 self.assertTrue(np.isnan(linearizer.linearityTurnoff[amp_name]))
654 self.assertTrue(np.isnan(linearizer.linearityMaxSignal[amp_name]))
655 else:
656 self.assertEqual(linearizer.linearityTurnoff[amp_name], turnoff)
657 self.assertEqual(linearizer.linearityMaxSignal[amp_name], max_signal)
659 # Check that the linearizer gives reasonable values over the
660 # range up to the ptc turnoff.
661 nodes, values = np.split(linearizer.linearityCoeffs[amp_name], 2)
662 self.assertEqual(values[0], 0.0)
663 to_test = (nodes > 0.0) & (nodes < ptc.ptcTurnoff[amp_name])
664 np.testing.assert_array_less(np.abs(values[to_test]/nodes[to_test]), 0.002)
666 # Check the residuals are reasonable up to the linearity
667 # turnoff.
668 to_test = ((ptc.rawMeans[amp_name] <= linearizer.linearityTurnoff[amp_name])
669 & np.isfinite(ptc.rawMeans[amp_name]))
670 residuals_scaled = linearizer.fitResiduals[amp_name][to_test]/ptc.rawMeans[amp_name][to_test]
671 np.testing.assert_array_less(np.abs(residuals_scaled), 0.0015)
673 # Try again after cutting it off, make sure it warns.
674 cutoff = (ordinate < turnoff)
676 with self.assertLogs(level=logging.INFO) as cm:
677 turnoff_index2, turnoff2, max_signal2, _ = _computeTurnoffAndMax(
678 abscissa[cutoff],
679 ordinate[cutoff],
680 ptc_mask[cutoff],
681 np.zeros(len(abscissa))[cutoff],
682 )
683 self.assertIn("No linearity turnoff", cm.output[0])
684 self.assertEqual(turnoff_index2, len(ptc_mask[cutoff]) - 1)
686 def test_linearity_turnoff_lsstcam(self):
687 # Use some real LSSTCam linearity data to measure the turnoff.
688 exp_times, photo_charges, raw_means, ptc_mask = self._lsstcam_raw_linearity_data()
690 ptc = PhotonTransferCurveDataset()
691 ptc.ampNames = ["Amp1", "Amp2"]
692 for amp_name in ptc.ampNames:
693 ptc.inputExpIdPairs[amp_name] = np.zeros((len(exp_times), 2), dtype=np.int64).tolist()
694 ptc.rawExpTimes[amp_name] = exp_times
695 ptc.photoCharges[amp_name] = photo_charges
696 ptc.rawMeans[amp_name] = raw_means
697 ptc.expIdMask[amp_name] = ptc_mask
698 ptc.gain[amp_name] = 1.0
700 # Construct an auxiliary array by hand since we know these data.
701 aux_arr = np.zeros(len(exp_times))
702 grouping_values_truth = np.zeros(len(exp_times), dtype=np.int64)
703 ratio = raw_means / exp_times
704 low, = np.where(ratio < 100)
705 grouping_values_truth[low] = 0
706 aux_arr[low] = 0.15
707 med, = np.where((ratio > 800) & (800 < 1000))
708 grouping_values_truth[med] = 1
709 aux_arr[med] = 0.5
710 high, = np.where((ratio > 1200) | (exp_times > 70.0))
711 grouping_values_truth[high] = 2
712 aux_arr[high] = 0.8
714 ptc.auxValues["AUX"] = aux_arr
716 def _compare_grouping_values(arr_a, arr_b):
717 self.assertEqual(len(arr_a), len(arr_b))
718 u_a, ind_a = np.unique(arr_a, return_index=True)
719 u_a = arr_a[np.sort(ind_a)]
720 u_b, ind_b = np.unique(arr_b, return_index=True)
721 u_b = arr_b[np.sort(ind_b)]
722 self.assertEqual(len(u_a), len(u_b))
723 for value_a, value_b in zip(u_a, u_b):
724 np.testing.assert_array_equal(arr_a == value_a, arr_b == value_b)
726 # First test: do grouping by aux value.
727 config = LinearitySolveTask.ConfigClass()
728 config.splineGroupingColumn = "AUX"
730 grouping_values = _determineInputGroups(
731 ptc,
732 config.doAutoGrouping,
733 config.autoGroupingUseExptime,
734 config.autoGroupingMaxSignalFraction,
735 config.autoGroupingThreshold,
736 config.splineGroupingColumn,
737 config.minPhotodiodeCurrent,
738 )
739 _compare_grouping_values(grouping_values, grouping_values_truth)
741 # Second test: do grouping automatically, with exp time.
742 config = LinearitySolveTask.ConfigClass()
743 config.doAutoGrouping = True
745 grouping_values = _determineInputGroups(
746 ptc,
747 config.doAutoGrouping,
748 config.autoGroupingUseExptime,
749 config.autoGroupingMaxSignalFraction,
750 config.autoGroupingThreshold,
751 config.splineGroupingColumn,
752 config.minPhotodiodeCurrent,
753 )
754 _compare_grouping_values(grouping_values, grouping_values_truth)
756 # Third test: do grouping automatically, with photodiode.
757 config = LinearitySolveTask.ConfigClass()
758 config.usePhotodiode = True
759 config.doAutoGrouping = True
760 config.autoGroupingUseExptime = False
761 config.autoGroupingThreshold = 0.008
762 config.minPhotodiodeCurrent = 1e-10
764 # We have to modify here to allow the lowest point to be "ungrouped"
765 # because it has an illegally low photocharge.
766 grouping_values = _determineInputGroups(
767 ptc,
768 config.doAutoGrouping,
769 config.autoGroupingUseExptime,
770 config.autoGroupingMaxSignalFraction,
771 config.autoGroupingThreshold,
772 config.splineGroupingColumn,
773 config.minPhotodiodeCurrent,
774 )
775 grouping_values_truth_mod = grouping_values_truth.copy()
776 grouping_values_truth_mod[0] = -1
777 _compare_grouping_values(grouping_values, grouping_values_truth_mod)
779 def test_linearity_fit_lsstcam(self):
780 """
781 Check that fitting the linearity works reasonably well for
782 both the photodiode and exposure time.
783 """
784 exp_times, photo_charges, raw_means, ptc_mask = self._lsstcam_raw_linearity_data()
786 # Cut off the fake data at the start and end.
787 exp_times = exp_times[1: -2]
788 photo_charges = photo_charges[1: -2]
789 raw_means = raw_means[1: -2]
790 ptc_mask = ptc_mask[1: -2]
792 ptc = PhotonTransferCurveDataset()
793 ptc.ampNames = self.amp_names
794 for amp_name in ptc.ampNames:
795 ptc.inputExpIdPairs[amp_name] = np.zeros((len(exp_times), 2), dtype=np.int64).tolist()
796 ptc.rawExpTimes[amp_name] = exp_times
797 ptc.photoCharges[amp_name] = photo_charges
798 ptc.rawMeans[amp_name] = raw_means
799 ptc.expIdMask[amp_name] = ptc_mask
800 ptc.gain[amp_name] = 1.0
801 ptc.ptcTurnoff[amp_name] = 75565.8
802 ptc.inputExpPairMjdStartList[amp_name] = np.cumsum(exp_times)
804 # Use the photodiode to solve for linearity.
805 config = LinearitySolveTask.ConfigClass()
806 config.usePhotodiode = True
807 config.doAutoGrouping = True
808 config.splineKnots = 5
809 config.trimmedState = False
810 config.linearityType = "Spline"
811 config.doSplineFitWeights = False
812 config.doSplineFitTemperature = False
813 config.doSplineFitOffset = True
814 config.splineFitMaxIter = 40
815 config.splineGroupingMinPoints = 50
816 task = LinearitySolveTask(config=config)
818 linearizer_pd = task.run(ptc, [self.dummy_exposure], self.camera, self.input_dims).outputLinearizer
820 # Use the exposure time to solve for linearity.
821 config = LinearitySolveTask.ConfigClass()
822 config.usePhotodiode = False
823 config.doAutoGrouping = True
824 config.splineKnots = 5
825 config.trimmedState = False
826 config.linearityType = "Spline"
827 config.doSplineFitWeights = False
828 config.doSplineFitTemperature = False
829 config.doSplineFitOffset = True
830 config.splineFitMaxIter = 40
831 config.splineGroupingMinPoints = 50
832 task = LinearitySolveTask(config=config)
834 ptc2 = copy.deepcopy(ptc)
835 for amp_name in self.amp_names:
836 ptc2.photoCharges[amp_name][:] = 0.0
838 linearizer_et = task.run(ptc2, [self.dummy_exposure], self.camera, self.input_dims).outputLinearizer
840 # Check that the linearizer with the photodiode seems reasonable.
841 amp_name = self.amp_names[0]
842 nodes_pd, values_pd = np.split(linearizer_pd.linearityCoeffs[amp_name], 2)
843 np.testing.assert_array_less(np.abs(values_pd), 100.0)
844 in_range_pd = np.isfinite(linearizer_pd.fitResiduals[amp_name])
845 np.testing.assert_array_less(
846 np.abs(linearizer_pd.fitResiduals[amp_name][in_range_pd] / ptc.rawMeans[amp_name][in_range_pd]),
847 0.004,
848 )
850 # Check that the linearizer with the exposure time seems reasonable.
851 nodes_et, values_et = np.split(linearizer_et.linearityCoeffs[amp_name], 2)
852 np.testing.assert_array_less(np.abs(values_et), 100.0)
853 in_range_et = np.isfinite(linearizer_et.fitResiduals[amp_name])
854 np.testing.assert_array_less(
855 np.abs(linearizer_et.fitResiduals[amp_name][in_range_et] / ptc.rawMeans[amp_name][in_range_et]),
856 0.004,
857 )
859 # Check consistency. This is very rough
860 self.assertTrue(np.allclose(nodes_et, nodes_pd))
861 self.assertTrue(np.allclose(values_et[1:] / nodes_et[1:], values_pd[1:] / nodes_pd[1:], atol=5e-4))
863 def test_linearity_renormalization_lsstcam(self):
864 # In this test we take the sample data, and set a bunch
865 # of amps with the same values. After renormalization and
866 # refitting there should be no residuals. (This is not a
867 # recommended way of using the renormalization code).
868 exp_times, photo_charges, raw_means, ptc_mask = self._lsstcam_raw_linearity_data()
870 # Cut off the fake data at the start and end.
871 exp_times = exp_times[1: -2]
872 photo_charges = photo_charges[1: -2]
873 raw_means = raw_means[1: -2]
874 ptc_mask = ptc_mask[1: -2]
876 ptc = PhotonTransferCurveDataset()
877 ptc.ampNames = self.amp_names
878 exp_ids = np.arange(len(exp_times)*2).reshape((len(exp_times), 2))
879 for amp_name in ptc.ampNames:
880 ptc.inputExpIdPairs[amp_name] = exp_ids.tolist()
881 ptc.rawExpTimes[amp_name] = exp_times
882 ptc.photoCharges[amp_name] = photo_charges
883 ptc.rawMeans[amp_name] = raw_means
884 ptc.expIdMask[amp_name] = ptc_mask
885 ptc.gain[amp_name] = 1.0
886 ptc.ptcTurnoff[amp_name] = 75565.8
887 ptc.inputExpPairMjdStartList[amp_name] = np.cumsum(exp_times)
889 # Use the exposure time to solve for linearity.
890 config = LinearitySolveTask.ConfigClass()
891 config.usePhotodiode = False
892 config.doAutoGrouping = True
893 config.splineKnots = 5
894 config.trimmedState = False
895 config.linearityType = "Spline"
896 config.doSplineFitWeights = False
897 config.doSplineFitTemperature = False
898 config.doSplineFitOffset = True
899 config.splineFitMaxIter = 40
900 config.splineGroupingMinPoints = 50
901 task = LinearitySolveTask(config=config)
903 linearizer = task.run(ptc, [self.dummy_exposure], self.camera, self.input_dims).outputLinearizer
905 # Now run the normalization task.
906 nconfig = LinearityNormalizeTask.ConfigClass()
907 nconfig.normalizeDetectors = [0]
908 nconfig.referenceDetector = 0
909 nconfig.doNormalizeAbsoluteLinearizer = False
910 ntask = LinearityNormalizeTask(config=nconfig)
912 ptc_handles = [
913 InMemoryDatasetHandle(
914 ptc,
915 dataId={"detector": 0},
916 ),
917 ]
918 linearizer_handles = [
919 InMemoryDatasetHandle(
920 linearizer,
921 dataId={"detector": 0},
922 ),
923 ]
925 normalization = ntask.run(
926 camera=self.camera,
927 inputPtcHandles=ptc_handles,
928 inputLinearizerHandles=linearizer_handles,
929 ).outputNormalization
931 # And re-solve the linearity.
932 config.useFocalPlaneNormalization = True
933 task2 = LinearitySolveTask(config=config)
935 linearizer2 = task2.run(
936 ptc,
937 [self.dummy_exposure],
938 self.camera,
939 self.input_dims,
940 inputNormalization=normalization,
941 ).outputLinearizer
943 # Check for zero residuals.
944 resids = linearizer2.fitResiduals[self.amp_names[0]]/linearizer2.fitResidualsModel[self.amp_names[0]]
945 self.assertFloatsAlmostEqual(np.nan_to_num(resids), 0.0, atol=5e-6)
947 def _comcam_raw_linearity_data(self):
948 # These are LSSTComCam measurements taken from a calibration
949 # run as part of DM-46357.
950 photo_charges = np.array(
951 [
952 3.22636000e-10, 3.38780400e-10, 3.54841300e-10, 3.87126000e-10,
953 4.03360000e-10, 4.35709800e-10, 4.51693200e-10, 4.83922500e-10,
954 5.16297600e-10, 5.48525400e-10, 5.80748400e-10, 6.13082500e-10,
955 6.45148000e-10, 6.77922000e-10, 7.25721750e-10, 7.74712800e-10,
956 8.22586650e-10, 8.70876900e-10, 9.19831800e-10, 9.68568000e-10,
957 1.03278400e-09, 1.08122255e-09, 1.14522645e-09, 1.22590280e-09,
958 1.29087600e-09, 1.37134750e-09, 1.45214100e-09, 1.53323350e-09,
959 1.62903910e-09, 1.72676065e-09, 1.83921900e-09, 1.93578600e-09,
960 2.06456960e-09, 2.17876500e-09, 2.30696180e-09, 2.45230720e-09,
961 2.59714735e-09, 2.74327300e-09, 2.91997345e-09, 3.08153670e-09,
962 3.27340599e-09, 3.46962700e-09, 3.67778820e-09, 3.88937850e-09,
963 4.12976640e-09, 4.37225980e-09, 4.62984095e-09, 4.90212160e-09,
964 5.19730540e-09, 5.50027885e-09, 5.83860750e-09, 6.17903475e-09,
965 6.54889207e-09, 6.93710400e-09, 7.35386640e-09, 7.79136960e-09,
966 8.25815040e-09, 8.74769030e-09, 9.27673375e-09, 9.82420530e-09,
967 1.04087197e-08, 1.10388366e-08, 1.17030225e-08, 1.23862272e-08,
968 1.31366576e-08, 1.39122921e-08, 1.47334462e-08, 1.56154856e-08,
969 1.65528171e-08, 1.75371688e-08, 1.85675291e-08, 1.96824430e-08,
970 2.08605509e-08, 2.21063200e-08, 2.34186546e-08, 2.48260115e-08,
971 2.62975235e-08, 2.78649723e-08, 2.95408665e-08, 3.12853772e-08,
972 3.31612268e-08, 3.51149011e-08, 3.72152551e-08, 3.94439604e-08,
973 4.17678940e-08, 4.42723820e-08, 4.69153456e-08, 4.96932949e-08,
974 5.27096702e-08, 5.58084160e-08, 5.91300138e-08, 6.26865948e-08,
975 6.64351212e-08, 7.03429300e-08, 7.45695391e-08, 7.89900791e-08,
976 8.37413792e-08, 9.40025100e-08, 9.95942560e-08, 8.87289232e-08,
977 ],
978 )
979 raw_means = np.array(
980 [
981 545.61773431, 572.92351724, 599.73600889, 654.80522687,
982 681.89236561, 736.36468948, 763.35233159, 818.3558025,
983 872.28757051, 926.97891772, 981.93722274, 1036.0641088,
984 1091.15536336, 1145.34900099, 1227.09504139, 1308.5247877,
985 1390.45828326, 1472.05043362, 1553.69250543, 1636.30579999,
986 1744.80644944, 1826.53524631, 1936.27610568, 2072.63450198,
987 2181.56004372, 2316.47487095, 2453.76647885, 2589.54660646,
988 2754.19117479, 2916.91523806, 3107.95392897, 3272.65202763,
989 3490.14569333, 3680.430583, 3896.17664859, 4146.03684447,
990 4391.418256, 4633.96031196, 4935.13618128, 5208.9852533,
991 5533.58336929, 5860.25295412, 6218.8609866, 6570.97901761,
992 6977.63082612, 7386.3804372, 7823.50205095, 8291.10153869,
993 8780.51405056, 9300.4324222, 9870.51463262, 10446.11494902,
994 11069.51440385, 11733.55548663, 12437.83099445, 13179.9215084,
995 13972.28308964, 14789.63243003, 15692.2466057, 16622.22691148,
996 17592.14597643, 18670.98290396, 19778.68937623, 20960.06668402,
997 22216.68121994, 23518.23906809, 24907.61374144, 26424.82837902,
998 28002.38769651, 29681.8198974, 31422.27519355, 33285.31574829,
999 35303.14661517, 37381.69611718, 39633.10759127, 41958.90250067,
1000 44458.19627488, 47100.85288959, 49902.33562771, 52837.68322699,
1001 55976.05887806, 59291.58601354, 62822.08965674, 66519.7625026,
1002 70480.78728045, 74721.82649122, 79220.00783343, 83928.39988665,
1003 88959.09479122, 94289.83758648, 99756.30512572, 104244.9446884,
1004 106938.47971957, 108134.21868844, 108602.08685706, 108702.48845118,
1005 108730.32842316, 108799.56402709, 108803.45364906, 108851.90458027,
1006 ],
1007 )
1009 ptc_mask = np.array(
1010 [
1011 True, True, True, True, True, True, True, True, True,
1012 True, True, True, True, True, True, True, True, True,
1013 True, True, True, True, True, True, True, True, True,
1014 True, True, True, True, True, True, True, True, True,
1015 True, True, True, True, True, True, True, True, True,
1016 True, True, True, True, True, True, True, True, True,
1017 True, True, True, True, True, True, True, True, True,
1018 True, True, True, True, True, True, True, True, True,
1019 True, True, True, True, True, True, True, True, True,
1020 True, True, True, True, True, True, True, False, False,
1021 False, False, False, False, False, False, False, False, False,
1022 False,
1023 ],
1024 )
1026 return photo_charges, raw_means, ptc_mask
1028 def _lsstcam_raw_linearity_data(self):
1029 # These are LSSTCam in-dome measurements taken from a calibration
1030 # run as part of DM-51470. Additional fake measurements have been
1031 # added that are below the photodiode limit or above the linearity
1032 # turnoff to test those aspects of the code.
1034 exp_times = np.array(
1035 [
1036 1.0,
1037 1.30007792, 2.60008502, 2.90007734, 3.20008755, 3.50007701,
1038 3.80058908, 4.10007262, 4.70008993, 6.00008655, 6.30008602,
1039 7.19907308, 7.50007725, 8.79909039, 9.10008216, 9.40107679,
1040 9.70009685, 10.9000566, 11.30008197, 11.90007639, 11.6015799,
1041 12.20007157, 12.80010366, 13.10108137, 13.70007062, 1.16707611,
1042 15.00008154, 15.30057025, 16.20007229, 1.3335712, 17.20008039,
1043 1.50008011, 18.4000895, 19.00007415, 21.20008397, 22.10008192,
1044 1.00008583, 1.00108552, 1.00058579, 24.00007939, 2.16806793,
1045 2.33307958, 1.32606268, 2.50008583, 2.66807985, 1.76707816,
1046 3.33358383, 3.50056648, 1.98808122, 3.83307695, 2.21007681,
1047 4.00059676, 4.16754985, 2.42907929, 4.66706753, 2.65107918,
1048 4.83408904, 5.00008345, 5.3330698, 5.50058341, 5.66706872,
1049 5.83259177, 3.31405592, 6.16708112, 3.53410554, 6.66809487,
1050 6.83407307, 3.97509313, 7.16708589, 4.19707441, 7.66709065,
1051 4.41806412, 8.50006843, 5.30106926, 9.83357477, 5.52107573,
1052 5.74307823, 10.33309817, 10.83407974, 11.50107312, 6.40508437,
1053 11.66708422, 11.83308768, 12.16760278, 6.84805465, 12.667593,
1054 7.06858993, 12.83307147, 7.50908542, 13.50007677, 13.83408427,
1055 14.00057459, 14.16807103, 7.95108128, 15.00008368, 15.00008249,
1056 15.00109196, 15.00109196, 15.00007319, 15.00058889, 15.00008273,
1057 15.00009346, 15.00007486, 15.00108337, 15.00109553, 15.00105524,
1058 15.00056362, 15.16708612, 15.33409524, 8.83507133, 9.05507565,
1059 16.33209324, 16.50058484, 16.66708279, 16.8340795, 17.33309412,
1060 9.71907902, 17.50007749, 17.6670754, 18.00008154, 18.16807771,
1061 10.38105488, 10.60207963, 19.16808653, 19.66707206, 20.1670773,
1062 20.66808057, 20.8330586, 11.7070713, 21.16708922, 21.33307409,
1063 11.92708254, 21.83359718, 22.00107551, 22.33308864, 12.59008217,
1064 22.83308482, 23.00007963, 12.8110795, 23.16708589, 23.33407331,
1065 13.03108048, 23.50059509, 13.47310114, 24.33407569, 13.69507217,
1066 24.66806412, 13.91508937, 25.00106311, 25.16805458, 25.33307934,
1067 25.50108433, 14.35706139, 26.167588, 14.57808304, 26.33309245,
1068 26.50007701, 26.8335743, 27.00106406, 15.24107933, 27.50059843,
1069 27.83407807, 28.001091, 28.16757345, 28.33307838, 28.50007868,
1070 15.90408254, 29.00007534, 16.124089, 29.16807365, 29.33310032,
1071 16.34508967, 29.50060129, 16.56608367, 30.00009513, 16.78757954,
1072 17.00807858, 17.45009923, 17.89208889, 18.11209965, 18.55409193,
1073 18.77606678, 19.21658659, 19.66007686, 20.54309058, 20.76208639,
1074 21.86608076, 22.30808902, 23.41357875, 23.63407302, 24.07609272,
1075 24.29708791, 24.51807809, 24.95957422, 25.18107915, 25.40207815,
1076 25.84308672, 26.0640893, 26.50607157, 26.72708106, 26.94707298,
1077 27.39009333, 27.83209324, 29.37858009, 30.04057312, 30.48257208,
1078 30.9245832, 31.14458227, 31.36508775, 31.58659363, 31.80709291,
1079 32.47107959, 32.91309094, 33.35208154, 34.01608729, 34.68008041,
1080 35.56258774, 36.4460969, 36.66707397, 36.88811016, 37.10908842,
1081 37.99308038, 38.65557122, 38.87708688, 39.75906563, 40.42207551,
1082 40.8640964, 41.30507493, 41.52558637, 42.18909407, 42.41057515,
1083 42.63108134, 43.07208395, 43.29357743, 43.51308179, 43.73508239,
1084 45.2810781, 45.50158954, 45.72306991, 46.1655736, 46.3850894,
1085 46.60510993, 47.04907584, 47.27006769, 47.48958015, 47.93157744,
1086 48.37308311, 48.59308457, 48.81608295, 49.47805619, 49.69907236,
1087 50.1410749, 50.5805645, 51.02406359, 51.68808389, 52.56959724,
1088 53.01307774, 53.67507911, 54.77957559, 55.00056076, 80.0, 90.0])
1090 photo_charges = np.array(
1091 [
1092 1e-11,
1093 1.82817260e-09, 3.53148649e-09, 4.06277898e-09, 4.35869544e-09,
1094 4.78330872e-09, 5.22915865e-09, 5.78226268e-09, 6.48815270e-09,
1095 8.25712380e-09, 8.51615745e-09, 1.02197526e-08, 1.05099569e-08,
1096 1.22837309e-08, 1.28114263e-08, 1.30927807e-08, 1.33760842e-08,
1097 1.53687893e-08, 1.56255057e-08, 1.61675487e-08, 1.64177802e-08,
1098 1.68976483e-08, 1.75646793e-08, 1.81330449e-08, 1.87132138e-08,
1099 2.91637877e-08, 2.07739826e-08, 2.18182545e-08, 2.19715436e-08,
1100 3.32051817e-08, 2.36684166e-08, 3.72810404e-08, 2.53565939e-08,
1101 2.59972840e-08, 2.96734901e-08, 3.02007963e-08, 4.55297856e-08,
1102 4.55746585e-08, 4.56246682e-08, 3.37802654e-08, 5.42455205e-08,
1103 5.81404431e-08, 6.04291412e-08, 6.24474643e-08, 6.64574931e-08,
1104 8.04680210e-08, 8.31232983e-08, 8.72382377e-08, 9.05033796e-08,
1105 9.55119486e-08, 1.00627749e-07, 9.97876553e-08, 1.04031649e-07,
1106 1.10573123e-07, 1.16371623e-07, 1.20965846e-07, 1.20569281e-07,
1107 1.24911511e-07, 1.32978879e-07, 1.37649727e-07, 1.41153092e-07,
1108 1.45213041e-07, 1.50986279e-07, 1.54133358e-07, 1.60948534e-07,
1109 1.66481769e-07, 1.70432685e-07, 1.81069871e-07, 1.78186815e-07,
1110 1.91063717e-07, 1.91086417e-07, 2.01389150e-07, 2.11345064e-07,
1111 2.41348459e-07, 2.45181656e-07, 2.51154663e-07, 2.61397147e-07,
1112 2.56946642e-07, 2.69402724e-07, 2.86605066e-07, 2.91289314e-07,
1113 2.90809466e-07, 2.95092192e-07, 3.03471554e-07, 3.11797168e-07,
1114 3.16873891e-07, 3.22205728e-07, 3.19714411e-07, 3.41500913e-07,
1115 3.36233190e-07, 3.44508625e-07, 3.49027053e-07, 3.53546367e-07,
1116 3.62263658e-07, 3.73349576e-07, 3.73605323e-07, 3.73355323e-07,
1117 3.73424493e-07, 3.73690267e-07, 3.74063761e-07, 3.73328827e-07,
1118 3.73513232e-07, 3.74292450e-07, 3.74058846e-07, 3.74199516e-07,
1119 3.73805125e-07, 3.74801019e-07, 3.77549794e-07, 3.82109671e-07,
1120 4.02813124e-07, 4.12015519e-07, 4.06590291e-07, 4.11174216e-07,
1121 4.15620699e-07, 4.20346548e-07, 4.31620648e-07, 4.42191673e-07,
1122 4.36834459e-07, 4.40553709e-07, 4.48188921e-07, 4.52683265e-07,
1123 4.72329423e-07, 4.82356891e-07, 4.77541263e-07, 4.90793892e-07,
1124 5.04259245e-07, 5.15785320e-07, 5.20004090e-07, 5.32583692e-07,
1125 5.27353277e-07, 5.31060399e-07, 5.43193778e-07, 5.44255280e-07,
1126 5.47986412e-07, 5.56257861e-07, 5.73223115e-07, 5.69229732e-07,
1127 5.73483192e-07, 5.84122304e-07, 5.78339239e-07, 5.80922371e-07,
1128 5.93536614e-07, 5.85343924e-07, 6.13633727e-07, 6.07300226e-07,
1129 6.23672682e-07, 6.14832531e-07, 6.33686470e-07, 6.24144708e-07,
1130 6.26834243e-07, 6.31470139e-07, 6.35002684e-07, 6.54419975e-07,
1131 6.52114981e-07, 6.63381949e-07, 6.56440208e-07, 6.63513673e-07,
1132 6.69092777e-07, 6.73678035e-07, 6.95448734e-07, 6.86136092e-07,
1133 6.94088711e-07, 6.97076014e-07, 7.02235128e-07, 7.06743404e-07,
1134 7.12541252e-07, 7.23700469e-07, 7.22175554e-07, 7.34407902e-07,
1135 7.26219718e-07, 7.30123718e-07, 7.43971514e-07, 7.36034974e-07,
1136 7.54263790e-07, 7.48980267e-07, 7.65759441e-07, 7.75726865e-07,
1137 7.94317280e-07, 8.14849394e-07, 8.24539049e-07, 8.44726672e-07,
1138 8.54993598e-07, 8.74385137e-07, 8.95097268e-07, 9.34897394e-07,
1139 9.44723938e-07, 9.96059623e-07, 1.01658226e-06, 1.06624643e-06,
1140 1.07704428e-06, 1.09541320e-06, 1.10813448e-06, 1.11542569e-06,
1141 1.13651858e-06, 1.14503830e-06, 1.15640469e-06, 1.17643071e-06,
1142 1.18651568e-06, 1.20725029e-06, 1.21732323e-06, 1.22769579e-06,
1143 1.24728449e-06, 1.26572090e-06, 1.33844128e-06, 1.36705783e-06,
1144 1.38812903e-06, 1.40766463e-06, 1.41726953e-06, 1.42920091e-06,
1145 1.43961190e-06, 1.44785313e-06, 1.47873846e-06, 1.49792958e-06,
1146 1.51710029e-06, 1.55136912e-06, 1.57902098e-06, 1.62106226e-06,
1147 1.65860885e-06, 1.66904934e-06, 1.68223145e-06, 1.68977232e-06,
1148 1.72980750e-06, 1.75960946e-06, 1.76905651e-06, 1.80953552e-06,
1149 1.83924183e-06, 1.86149564e-06, 1.88067160e-06, 1.88919825e-06,
1150 1.92234685e-06, 1.93096053e-06, 1.94121685e-06, 1.95972805e-06,
1151 1.97182165e-06, 1.98003698e-06, 1.99059056e-06, 2.06132937e-06,
1152 2.07290637e-06, 2.08140970e-06, 2.10187884e-06, 2.11122988e-06,
1153 2.12235208e-06, 2.14279142e-06, 2.15065476e-06, 2.16372892e-06,
1154 2.18410909e-06, 2.20113980e-06, 2.21079153e-06, 2.22210751e-06,
1155 2.25089289e-06, 2.26316028e-06, 2.28389811e-06, 2.30118685e-06,
1156 2.32201405e-06, 2.35181662e-06, 2.39407598e-06, 2.41469194e-06,
1157 2.44405239e-06, 2.49260152e-06, 2.50544044e-06, 3.64427700e-06,
1158 4.09981162e-06,
1159 ],
1160 )
1162 raw_means = np.array(
1163 [
1164 70.0,
1165 92.86846201, 182.20115781, 210.17843191, 224.75373652,
1166 249.41990014, 270.3750848, 298.13356004, 335.2491762,
1167 419.85418654, 440.20870837, 523.35053806, 533.95229561,
1168 624.33520762, 654.19741753, 663.41872989, 691.00929737,
1169 781.82840105, 797.53228322, 829.48856791, 829.77401764,
1170 869.06926893, 910.33500907, 924.56998291, 961.33369162,
1171 1028.74632532, 1080.23191006, 1115.72725114, 1145.29298889,
1172 1171.88485425, 1215.27614848, 1316.63175163, 1319.90794786,
1173 1329.20039219, 1513.53524709, 1544.58715074, 1577.72509245,
1174 1580.42422803, 1582.08035824, 1730.68368453, 1909.06875316,
1175 2051.60147106, 2095.65486231, 2200.50081135, 2346.93082746,
1176 2790.61021016, 2933.16186159, 3077.82799132, 3136.72492514,
1177 3371.89652643, 3491.51091289, 3520.01335135, 3670.82553208,
1178 3838.56300847, 4102.6067899, 4195.75368619, 4247.1381654,
1179 4403.42219284, 4690.74436914, 4847.29352829, 4976.48877486,
1180 5120.97366735, 5240.94946909, 5433.09472747, 5590.18821641,
1181 5873.91979779, 6013.31482841, 6283.71507483, 6297.6863956,
1182 6631.55638028, 6747.12360371, 6992.05043094, 7473.70409509,
1183 8383.09997091, 8651.84396759, 8724.92350749, 9081.20376784,
1184 9088.49929752, 9526.30681728, 10113.24561528, 10125.09736262,
1185 10275.02009773, 10413.30226, 10716.08347878, 10831.62417085,
1186 11174.97345771, 11187.10433017, 11297.94913483, 11874.65293352,
1187 11877.8316926, 12174.65186236, 12321.831863, 12476.93862506,
1188 12595.30191395, 13183.70211442, 13186.68422581, 13187.2122269,
1189 13191.52699678, 13196.79290796, 13204.41871041, 13205.11971342,
1190 13210.91337662, 13213.64928761, 13215.36094497, 13215.79616817,
1191 13217.20001814, 13228.6272711, 13341.71183545, 13492.87662206,
1192 13983.88641078, 14320.43230657, 14365.05001829, 14521.47003941,
1193 14676.59292711, 14836.45065264, 15256.72039757, 15383.54413924,
1194 15429.46336546, 15552.75261379, 15829.4470338, 15980.36094398,
1195 16422.11284161, 16773.36075251, 16873.76417919, 17347.93050569,
1196 17780.51617957, 18213.68136006, 18370.11017442, 18520.97641651,
1197 18626.94414553, 18774.53387747, 18872.50742297, 19220.11365054,
1198 19371.40480846, 19669.37342574, 19918.48776394, 20119.78165825,
1199 20260.25387413, 20297.01090622, 20420.94501903, 20540.65397603,
1200 20647.36826825, 20695.08978438, 21326.45152597, 21450.48862807,
1201 21693.06777509, 21746.54361219, 22045.21577418, 22057.42176734,
1202 22169.88708478, 22322.21832224, 22461.40321731, 22752.13471875,
1203 23049.46907456, 23089.19566699, 23176.810312, 23435.02471521,
1204 23659.65784493, 23804.24440484, 24177.15459574, 24239.82590366,
1205 24538.08527224, 24653.42322741, 24834.63478483, 24968.23972901,
1206 25166.10351596, 25177.0618337, 25535.14004216, 25536.46094104,
1207 25672.96722423, 25819.18373405, 25881.30357988, 26032.17886968,
1208 26236.72558439, 26474.25109961, 26624.66379018, 26976.051577,
1209 27641.68356156, 28354.41545864, 28684.89234527, 29401.04773183,
1210 29749.42931579, 30424.52582892, 31124.41787207, 32531.15694883,
1211 32876.38317561, 34657.47503525, 35358.55517438, 37065.5660748,
1212 37453.44032424, 38121.45997741, 38516.02714121, 38806.56101841,
1213 39499.85472998, 39834.64346431, 40243.25907892, 40942.40514966,
1214 41279.92463902, 41977.41767446, 42332.91963243, 42700.34604663,
1215 43372.00824829, 44037.57319312, 46538.15651649, 47533.00991698,
1216 48252.54686708, 48926.33644579, 49278.174449, 49662.66006463,
1217 49986.82486074, 50323.45507686, 51397.80697509, 52070.1823656,
1218 52704.70512386, 53869.17336706, 54765.71986416, 56274.31527709,
1219 57552.6631826, 57973.25092965, 58357.98782968, 58660.00904752,
1220 60045.31679453, 61092.97983081, 61429.86712469, 62836.09931594,
1221 63824.97166359, 64603.41083149, 65266.46358084, 65563.45761106,
1222 66707.66379327, 66974.76300745, 67353.24759166, 68038.21337638,
1223 68400.62102102, 68744.27455369, 69097.6239373, 71540.01858064,
1224 71916.28448469, 72279.98815642, 72939.45375775, 73289.46074577,
1225 73692.17736103, 74419.50733063, 74659.92625471, 75082.30090212,
1226 75797.60997538, 76437.89738215, 76797.10336102, 77171.44260012,
1227 78130.37592659, 78568.4843529, 79279.27081724, 79917.26494941,
1228 80621.12449606, 81643.09849749, 83131.19278591, 83798.73997235,
1229 84860.43286242, 86483.42983764, 86958.02382474, 100000.0,
1230 100000.0,
1231 ],
1232 )
1234 ptc_mask = np.array(
1235 [
1236 False,
1237 False, True, True, True, True, True, True, True, True,
1238 True, True, True, True, True, True, True, True, True,
1239 True, True, True, True, True, True, True, True, True,
1240 True, True, True, True, True, True, True, True, True,
1241 True, True, True, True, True, True, True, True, True,
1242 True, True, True, True, True, True, True, True, True,
1243 True, True, True, True, True, True, True, True, True,
1244 True, True, True, True, True, True, True, True, True,
1245 True, True, True, True, True, True, True, True, True,
1246 True, True, True, True, True, True, True, True, True,
1247 True, True, True, True, True, True, True, True, True,
1248 True, True, True, True, True, True, True, True, True,
1249 True, True, True, True, True, True, True, True, True,
1250 True, True, True, True, True, True, True, True, True,
1251 True, True, True, True, True, True, True, True, True,
1252 True, True, True, True, True, True, True, True, True,
1253 True, True, True, True, True, True, True, True, True,
1254 True, True, True, True, True, True, True, True, True,
1255 True, True, True, True, True, True, True, True, True,
1256 True, True, True, True, True, True, True, True, True,
1257 True, True, True, True, True, True, True, True, True,
1258 True, True, True, True, True, True, True, True, True,
1259 True, True, True, True, True, True, True, True, True,
1260 True, True, True, True, True, True, True, True, True,
1261 True, True, True, True, True, True, True, True, True,
1262 True, True, True, True, True, True, True, True, True,
1263 True, True, False, False, False, False, False, False, False,
1264 False, False, False, False, False, False, False, False, False,
1265 False, False, False, False, False, False, False, False, False,
1266 ],
1267 )
1269 return exp_times, photo_charges, raw_means, ptc_mask
1272class DoubleSplineLinearityTestCase(lsst.utils.tests.TestCase):
1273 def setUp(self):
1274 mock = IsrMockLSST()
1275 self.camera = mock.getCamera()
1276 self.detector = self.camera[20]
1277 self.detector_id = self.detector.getId()
1279 def _compute_logistic_nonlinearity(self, xvals, midpoint, amplitude, transition=3000.0):
1280 """Compute a simple non-linearity with a logistic curve.
1282 Parameters
1283 ----------
1284 xvals : `np.ndarray`
1285 Input count values.
1286 midpoint : `float`
1287 Transition point for logistic curve.
1288 amplitude : `float`
1289 Fractional amplitude of logistic curve.
1290 transition : `float`, optional
1291 Transition value (sharpness).
1293 Returns
1294 -------
1295 offsets : `np.ndarray`
1296 Offset count values.
1297 """
1298 frac_offset = amplitude / 2. - (amplitude / (1. + np.exp(-(1./transition)*(xvals - midpoint))))
1300 return xvals * frac_offset
1302 def test_linearity_doublespline(self):
1303 n_pair = 100
1304 pair_sigma = 0.005 # Fractional variation.
1306 rng = np.random.RandomState(seed=12345)
1308 amp_names = [amp.getName() for amp in self.detector]
1309 n_amps = len(amp_names)
1311 rel_amplitudes = rng.uniform(low=0.001, high=0.003, size=n_amps)
1312 rel_midpoints = rng.uniform(low=20000.0, high=50000.0, size=n_amps)
1313 linearity_turnoffs = rng.uniform(low=90000.0, high=100000.0, size=n_amps)
1315 noises = rng.uniform(low=5.0, high=10.0, size=n_amps)
1316 gains = rng.uniform(low=1.4, high=1.6, size=n_amps)
1317 a00s = rng.uniform(low=-2e-6, high=-4e-6, size=n_amps)
1318 ptc_turnoffs = rng.uniform(low=70000.0, high=80000.0, size=n_amps)
1320 ref_amp_index = np.argmax(linearity_turnoffs)
1321 ref_amp_name = amp_names[ref_amp_index]
1323 # Make sure this one is significantly higher for consistency.
1324 linearity_turnoffs[ref_amp_index] *= 1.1
1326 # Reset the relative linearity for the reference amp.
1327 rel_amplitudes[ref_amp_index] = 0.0
1329 # Create an absolute linearizer.
1330 abs_amplitude = -0.001
1331 abs_midpoint = 10000.0
1333 range_e = np.asarray([50.0, linearity_turnoffs[ref_amp_index]]) * gains[ref_amp_index]
1335 pair_levels_e = np.linspace(range_e[0], range_e[1], n_pair)
1337 ptc_extract_config = PhotonTransferCurveExtractPairTask.ConfigClass()
1338 ptc_extract_config.doOutputBinnedImages = True
1339 ptc_extract_config.maximumRangeCovariancesAstier = 1
1340 ptc_extract_config.minNumberGoodPixelsForCovariance = 100
1341 ptc_extract_config.numEdgeSuspect = 4
1342 ptc_extract_config.edgeMaskLevel = "AMP"
1343 ptc_extract_config.doGain = False
1345 ptc_extract_task = PhotonTransferCurveExtractPairTask(config=ptc_extract_config)
1347 pair_handles = []
1348 binned_handles = []
1350 normalization_exposures = []
1351 normalization_values = []
1353 for i in range(n_pair):
1354 # We generate two images with the same level for
1355 # simpler testing.
1356 # The level is done in electrons then we apply gain,
1357 # and compute the variance from the Astier function.
1359 levels_e = pair_levels_e[i] + rng.normal(loc=0.0, scale=pair_sigma * pair_levels_e[i], size=2)
1360 # levels_e = pair_levels_e[i] + np.array([0.0, 0.0])
1361 exptime = np.mean(levels_e) / 10.0
1363 # normalization_exposures.extend([i * 2
1365 flat_pair = []
1366 for j in range(2):
1367 flat = lsst.afw.image.ExposureF(self.detector.getBBox())
1368 flat.setDetector(self.detector)
1369 visit_info = lsst.afw.image.VisitInfo(
1370 exposureTime=exptime,
1371 date=DateTime("2025-09-25T00:00:00", DateTime.TAI),
1372 )
1373 flat.info.id = i * 2 + j
1374 flat.info.setVisitInfo(visit_info)
1375 flat_pair.append(flat)
1377 normalization_exposures.append(i * 2 + j)
1378 normalization_values.append((levels_e[j] / 10.) / exptime)
1380 for j in range(2):
1381 flat = flat_pair[j]
1382 # Compute the variance for all the amps at this level.
1383 var_adu = funcAstier([a00s, gains, noises**2.], levels_e[j] / gains)
1385 # Adjust the variance for those above the ptc turnoff
1386 var_adu[(levels_e[j] / gains) > ptc_turnoffs] *= 0.5
1388 # Build the flat.
1389 for k, amp in enumerate(self.detector):
1391 # Offset things above the linearizer turnoff.
1392 # I don't know if this will actually work ...
1393 level_e = levels_e[j]
1394 if level_e > linearity_turnoffs[k] * gains[k]:
1395 level_e *= 0.8
1396 var_adu[k] *= 0.8
1398 noise_key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {amp.getName()}"
1399 flat.metadata[noise_key] = noises[k] / gains[k]
1400 median_key = f"LSST ISR OVERSCAN SERIAL MEDIAN {amp.getName()}"
1401 flat.metadata[median_key] = 0.0
1402 pedestal_key = f"LSST ISR AMPOFFSET PEDESTAL {amp.getName()}"
1403 flat.metadata[pedestal_key] = 0.0
1405 bbox = amp.getBBox()
1406 flat[bbox].image.array[:, :] = rng.normal(
1407 loc=level_e,
1408 scale=np.sqrt(var_adu[k]) * gains[k],
1409 size=flat[bbox].image.array.shape,
1410 ) / gains[k]
1412 # Apply the non-linearity... first the absolute.
1413 abs_offset = self._compute_logistic_nonlinearity(
1414 flat[bbox].image.array,
1415 abs_midpoint,
1416 abs_amplitude,
1417 )
1418 flat[bbox].image.array[:, :] += abs_offset
1420 # And then the relative.
1421 rel_offset = self._compute_logistic_nonlinearity(
1422 flat[bbox].image.array,
1423 rel_midpoints[k],
1424 rel_amplitudes[k],
1425 )
1426 flat[bbox].image.array[:, :] += rel_offset
1428 # Run the PTC extraction on the pair and store the rebinned images.
1429 # How long does this take? Memory?
1431 handles = [
1432 InMemoryDatasetHandle(
1433 flat,
1434 dataId={"exposure": flat.info.id, "detector": flat.getDetector().getId()},
1435 )
1436 for flat in flat_pair
1437 ]
1438 inputDims = [flat.info.id for flat in flat_pair]
1439 results = ptc_extract_task.run(inputExp=handles, inputDims=inputDims)
1441 data_id = {"exposure": flat_pair[0].info.id, "detector": flat_pair[0].getDetector().getId()}
1442 pair_handles.append(
1443 InMemoryDatasetHandle(
1444 results.outputCovariance,
1445 dataId=data_id,
1446 )
1447 )
1448 binned_handles.append(
1449 InMemoryDatasetHandle(
1450 results.outputBinnedImages,
1451 dataId=data_id,
1452 )
1453 )
1455 normalization = Table({"exposure": normalization_exposures, "normalization": normalization_values})
1457 # Build the linearizer PTC.
1458 ptc_solve_config = PhotonTransferCurveSolveTask.ConfigClass()
1459 ptc_solve_config.ptcFitType = "EXPAPPROXIMATION"
1460 ptc_solve_config.maximumRangeCovariancesAstier = 1
1461 ptc_solve_config.maximumRangeCovariancesAstierFullCovFit = 1
1462 ptc_solve_config.ksTestMinPvalue = 0.0
1463 # This is a large number because these small amplifiers are noisy.
1464 ptc_solve_config.sigmaCutPtcOutliers = 20.0
1466 ptc_solve_task = PhotonTransferCurveSolveTask(config=ptc_solve_config)
1467 input_covariances = [handle.get() for handle in pair_handles]
1468 results = ptc_solve_task.run(input_covariances, camera=self.camera, detId=self.detector_id)
1469 linearizer_ptc = results.outputPtcDataset
1471 # Now we run the linearize solving code.
1473 linearity_solve_config = LinearityDoubleSplineSolveTask.ConfigClass()
1474 linearity_solve_config.relativeSplineMinimumSignalNode = 0.0
1475 linearity_solve_config.relativeSplineLowThreshold = 0.0
1476 linearity_solve_config.relativeSplineMidNodeSize = 10000.0
1477 linearity_solve_config.relativeSplineHighNodeSize = 10000.0
1478 linearity_solve_config.absoluteSplineNodeSize = 10000.0
1479 linearity_solve_config.useFocalPlaneNormalization = True
1480 linearity_solve_config.maxLinearityTurnoffRelativeToPtcTurnoff = 1.5
1481 linearity_solve_task = LinearityDoubleSplineSolveTask(config=linearity_solve_config)
1483 results = linearity_solve_task.run(
1484 inputPtc=linearizer_ptc,
1485 camera=self.camera,
1486 inputBinnedImagesHandles=binned_handles,
1487 inputNormalization=normalization,
1488 )
1490 linearizer = results.outputLinearizer
1492 # Check that we chose the correct reference amplifier.
1493 self.assertEqual(linearizer.absoluteReferenceAmplifier, ref_amp_name)
1495 # Check the relative residuals.
1496 for amp_name in linearizer.ampNames:
1497 if amp_name == ref_amp_name:
1498 continue
1500 frac_resid = linearizer.fitResiduals[amp_name] / linearizer.inputOrdinate[amp_name]
1502 self.assertFloatsAlmostEqual(np.nan_to_num(frac_resid), 0.0, atol=1e-3)
1504 # Check the absolute residuals.
1505 abs_frac_resid = linearizer.fitResiduals[ref_amp_name] / linearizer.inputOrdinate[ref_amp_name]
1506 self.assertFloatsAlmostEqual(np.nan_to_num(abs_frac_resid), 0.0, atol=5e-4)
1508 def _check_spline_midpoint(coeffs, expected_midpoint, turnoff, tolerance=6000.0):
1509 centers, values = np.split(coeffs, 2)
1511 xvals = np.linspace(1.0, centers[-1], 1000)
1512 spl = Akima1DInterpolator(centers, values, method="akima")
1513 yvals = spl(xvals) / xvals
1515 grad = np.gradient(yvals, xvals)
1516 # This is the point of maximum gradient, which should match
1517 # the input midpoints.
1518 max_grad = xvals[np.argmax(np.abs(grad))]
1520 self.assertFloatsAlmostEqual(max_grad, expected_midpoint, atol=tolerance)
1522 self.assertLessEqual(centers[-1], turnoff)
1524 # Check that the linearizers are consistent with expectations.
1525 coeffs = linearizer.linearityCoeffs[ref_amp_name]
1526 n_nodes1 = int(coeffs[0])
1527 n_nodes2 = int(coeffs[1])
1528 self.assertEqual(n_nodes1, 0)
1529 abs_coeff = coeffs[2 + 2 * n_nodes1: 2 + 2 * n_nodes1 + 2 * n_nodes2]
1531 _check_spline_midpoint(abs_coeff, abs_midpoint, linearity_turnoffs[ref_amp_index], tolerance=4000.0)
1533 for i, amp_name in enumerate(linearizer.ampNames):
1534 if amp_name == ref_amp_name:
1535 continue
1537 coeffs = linearizer.linearityCoeffs[amp_name]
1539 n_nodes1 = int(coeffs[0])
1540 n_nodes2 = int(coeffs[1])
1542 spline_coeff1 = coeffs[2: 2 + 2 * n_nodes1]
1544 _check_spline_midpoint(spline_coeff1, rel_midpoints[i], linearity_turnoffs[i])
1546 # Confirm that the absolute parameters are matched.
1547 spline_coeff2 = coeffs[2 + 2 * n_nodes1: 2 + 2 * n_nodes1 + 2 * n_nodes2]
1548 self.assertFloatsAlmostEqual(spline_coeff2, abs_coeff)
1550 def test_noderator(self):
1551 # Test "regular" usage.
1552 low = 5000.0
1553 mid = 70000.0
1554 high = 100000.0
1555 nodes = _noderator(low, mid, high, 50.0, 750.0, 5000.0, 2000.0)
1557 self.assertEqual(nodes[0], 0.0)
1558 self.assertEqual(nodes[1], 50.0)
1559 self.assertEqual(nodes[-1], high)
1560 self.assertTrue(np.any(nodes == low))
1561 self.assertTrue(np.any(nodes == mid))
1562 self.assertTrue(np.all(np.arange(len(nodes)) == np.argsort(nodes)))
1563 self.assertGreater(nodes[2] - nodes[1], 750.0)
1564 self.assertGreater(nodes[10] - nodes[9], 5000.0)
1565 self.assertGreater(nodes[-1] - nodes[-2], 2000.0)
1567 # Test no "low turnoff"
1568 low = 0.0
1569 nodes = _noderator(low, mid, high, 0.0, 750.0, 5000.0, 2000.0)
1570 self.assertEqual(nodes[0], 0.0)
1571 self.assertEqual(nodes[-1], high)
1572 self.assertTrue(np.any(nodes == mid))
1573 self.assertGreater(nodes[1] - nodes[0], 5000.0)
1574 self.assertGreater(nodes[-1] - nodes[-2], 2000.0)
1576 # Test no "high turnoff"
1577 # This includes when the ptc turnoff is above the linearity
1578 # turnoff.
1579 low = 5000.0
1580 high = mid - 1.0
1581 nodes = _noderator(low, mid, high, 50.0, 750.0, 5000.0, 2000.0)
1582 self.assertEqual(nodes[0], 0.0)
1583 self.assertEqual(nodes[1], 50.0)
1584 self.assertEqual(nodes[-1], mid)
1585 self.assertTrue(np.any(nodes == low))
1586 self.assertTrue(np.any(nodes == mid))
1587 self.assertTrue(np.all(np.arange(len(nodes)) == np.argsort(nodes)))
1588 self.assertGreater(nodes[2] - nodes[1], 750.0)
1589 self.assertGreater(nodes[10] - nodes[9], 5000.0)
1590 self.assertGreater(nodes[-1] - nodes[-2], 5000.0)
1593class TestMemory(lsst.utils.tests.MemoryTestCase):
1594 pass
1597def setup_module(module):
1598 lsst.utils.tests.init()
1601if __name__ == "__main__": 1601 ↛ 1602line 1601 didn't jump to line 1602 because the condition on line 1601 was never true
1602 lsst.utils.tests.init()
1603 unittest.main()