LSST Applications  21.0.0-131-g8cabc107+528f53ee53,22.0.0+00495a2688,22.0.0+0ef2527977,22.0.0+11a2aa21cd,22.0.0+269b7e55e3,22.0.0+2c6b6677a3,22.0.0+64c1bc5aa5,22.0.0+7b3a3f865e,22.0.0+e1b6d2281c,22.0.0+ff3c34362c,22.0.1-1-g1b65d06+c95cbdf3df,22.0.1-1-g7058be7+1cf78af69b,22.0.1-1-g7dab645+2a65e40b06,22.0.1-1-g8760c09+64c1bc5aa5,22.0.1-1-g949febb+64c1bc5aa5,22.0.1-1-ga324b9c+269b7e55e3,22.0.1-1-gf9d8b05+ff3c34362c,22.0.1-10-g781e53d+9b51d1cd24,22.0.1-10-gba590ab+b9624b875d,22.0.1-13-g76f9b8d+2c6b6677a3,22.0.1-14-g22236948+57af756299,22.0.1-18-g3db9cf4b+9b7092c56c,22.0.1-18-gb17765a+2264247a6b,22.0.1-2-g8ef0a89+2c6b6677a3,22.0.1-2-gcb770ba+c99495d3c6,22.0.1-24-g2e899d296+4206820b0d,22.0.1-3-g7aa11f2+2c6b6677a3,22.0.1-3-g8c1d971+f253ffa91f,22.0.1-3-g997b569+ff3b2f8649,22.0.1-4-g1930a60+6871d0c7f6,22.0.1-4-g5b7b756+6b209d634c,22.0.1-6-ga02864e+6871d0c7f6,22.0.1-7-g3402376+a1a2182ac4,22.0.1-7-g65f59fa+54b92689ce,master-gcc5351303a+e1b6d2281c,w.2021.32
LSST Data Management Base Package
linearity.py
Go to the documentation of this file.
1 # This file is part of cp_pipe.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
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 GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 #
22 import numpy as np
23 
24 import lsst.afw.image as afwImage
25 import lsst.afw.math as afwMath
26 import lsst.pipe.base as pipeBase
27 import lsst.pipe.base.connectionTypes as cT
28 import lsst.pex.config as pexConfig
29 
30 from lsstDebug import getDebugFrame
31 from lsst.ip.isr import (Linearizer, IsrProvenance)
32 
33 from .utils import (funcPolynomial, irlsFit)
34 from ._lookupStaticCalibration import lookupStaticCalibration
35 
36 __all__ = ["LinearitySolveTask", "LinearitySolveConfig", "MeasureLinearityTask"]
37 
38 
39 class LinearitySolveConnections(pipeBase.PipelineTaskConnections,
40  dimensions=("instrument", "exposure", "detector")):
41  dummy = cT.Input(
42  name="raw",
43  doc="Dummy exposure.",
44  storageClass='Exposure',
45  dimensions=("instrument", "exposure", "detector"),
46  multiple=True,
47  deferLoad=True,
48  )
49  camera = cT.PrerequisiteInput(
50  name="camera",
51  doc="Camera Geometry definition.",
52  storageClass="Camera",
53  dimensions=("instrument", ),
54  isCalibration=True,
55  lookupFunction=lookupStaticCalibration,
56  )
57  inputPtc = cT.PrerequisiteInput(
58  name="ptc",
59  doc="Input PTC dataset.",
60  storageClass="PhotonTransferCurveDataset",
61  dimensions=("instrument", "detector"),
62  isCalibration=True,
63  )
64 
65  outputLinearizer = cT.Output(
66  name="linearity",
67  doc="Output linearity measurements.",
68  storageClass="Linearizer",
69  dimensions=("instrument", "detector"),
70  isCalibration=True,
71  )
72 
73 
74 class LinearitySolveConfig(pipeBase.PipelineTaskConfig,
75  pipelineConnections=LinearitySolveConnections):
76  """Configuration for solving the linearity from PTC dataset.
77  """
78  linearityType = pexConfig.ChoiceField(
79  dtype=str,
80  doc="Type of linearizer to construct.",
81  default="Squared",
82  allowed={
83  "LookupTable": "Create a lookup table solution.",
84  "Polynomial": "Create an arbitrary polynomial solution.",
85  "Squared": "Create a single order squared solution.",
86  "Spline": "Create a spline based solution.",
87  "None": "Create a dummy solution.",
88  }
89  )
90  polynomialOrder = pexConfig.Field(
91  dtype=int,
92  doc="Degree of polynomial to fit.",
93  default=3,
94  )
95  splineKnots = pexConfig.Field(
96  dtype=int,
97  doc="Number of spline knots to use in fit.",
98  default=10,
99  )
100  maxLookupTableAdu = pexConfig.Field(
101  dtype=int,
102  doc="Maximum DN value for a LookupTable linearizer.",
103  default=2**18,
104  )
105  maxLinearAdu = pexConfig.Field(
106  dtype=float,
107  doc="Maximum DN value to use to estimate linear term.",
108  default=20000.0,
109  )
110  minLinearAdu = pexConfig.Field(
111  dtype=float,
112  doc="Minimum DN value to use to estimate linear term.",
113  default=2000.0,
114  )
115  nSigmaClipLinear = pexConfig.Field(
116  dtype=float,
117  doc="Maximum deviation from linear solution for Poissonian noise.",
118  default=5.0,
119  )
120  ignorePtcMask = pexConfig.Field(
121  dtype=bool,
122  doc="Ignore the expIdMask set by the PTC solver?",
123  default=False,
124  )
125 
126 
127 class LinearitySolveTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
128  """Fit the linearity from the PTC dataset.
129  """
130  ConfigClass = LinearitySolveConfig
131  _DefaultName = 'cpLinearitySolve'
132 
133  def runQuantum(self, butlerQC, inputRefs, outputRefs):
134  """Ensure that the input and output dimensions are passed along.
135 
136  Parameters
137  ----------
138  butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
139  Butler to operate on.
140  inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
141  Input data refs to load.
142  ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
143  Output data refs to persist.
144  """
145  inputs = butlerQC.get(inputRefs)
146 
147  # Use the dimensions to set calib/provenance information.
148  inputs['inputDims'] = inputRefs.inputPtc.dataId.byName()
149 
150  outputs = self.runrun(**inputs)
151  butlerQC.put(outputs, outputRefs)
152 
153  def run(self, inputPtc, dummy, camera, inputDims):
154  """Fit non-linearity to PTC data, returning the correct Linearizer
155  object.
156 
157  Parameters
158  ----------
159  inputPtc : `lsst.cp.pipe.PtcDataset`
160  Pre-measured PTC dataset.
161  dummy : `lsst.afw.image.Exposure`
162  The exposure used to select the appropriate PTC dataset.
163  In almost all circumstances, one of the input exposures
164  used to generate the PTC dataset is the best option.
165  camera : `lsst.afw.cameraGeom.Camera`
166  Camera geometry.
167  inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
168  DataIds to use to populate the output calibration.
169 
170  Returns
171  -------
172  results : `lsst.pipe.base.Struct`
173  The results struct containing:
174 
175  ``outputLinearizer`` : `lsst.ip.isr.Linearizer`
176  Final linearizer calibration.
177  ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
178  Provenance data for the new calibration.
179 
180  Notes
181  -----
182  This task currently fits only polynomial-defined corrections,
183  where the correction coefficients are defined such that:
184  corrImage = uncorrImage + sum_i c_i uncorrImage^(2 + i)
185  These `c_i` are defined in terms of the direct polynomial fit:
186  meanVector ~ P(x=timeVector) = sum_j k_j x^j
187  such that c_(j-2) = -k_j/(k_1^j) in units of DN^(1-j) (c.f.,
188  Eq. 37 of 2003.05978). The `config.polynomialOrder` or
189  `config.splineKnots` define the maximum order of x^j to fit.
190  As k_0 and k_1 are degenerate with bias level and gain, they
191  are not included in the non-linearity correction.
192  """
193  if len(dummy) == 0:
194  self.log.warn("No dummy exposure found.")
195 
196  detector = camera[inputDims['detector']]
197  if self.config.linearityType == 'LookupTable':
198  table = np.zeros((len(detector), self.config.maxLookupTableAdu), dtype=np.float32)
199  tableIndex = 0
200  else:
201  table = None
202  tableIndex = None # This will fail if we increment it.
203 
204  if self.config.linearityType == 'Spline':
205  fitOrder = self.config.splineKnots
206  else:
207  fitOrder = self.config.polynomialOrder
208 
209  # Initialize the linearizer.
210  linearizer = Linearizer(detector=detector, table=table, log=self.log)
211 
212  for i, amp in enumerate(detector):
213  ampName = amp.getName()
214  if ampName in inputPtc.badAmps:
215  nEntries = 1
216  pEntries = 1
217  if self.config.linearityType in ['Polynomial']:
218  nEntries = fitOrder + 1
219  pEntries = fitOrder + 1
220  elif self.config.linearityType in ['Spline']:
221  nEntries = fitOrder * 2
222  elif self.config.linearityType in ['Squared', 'None']:
223  nEntries = 1
224  pEntries = fitOrder + 1
225  elif self.config.linearityType in ['LookupTable']:
226  nEntries = 2
227  pEntries = fitOrder + 1
228 
229  linearizer.linearityType[ampName] = "None"
230  linearizer.linearityCoeffs[ampName] = np.zeros(nEntries)
231  linearizer.linearityBBox[ampName] = amp.getBBox()
232  linearizer.fitParams[ampName] = np.zeros(pEntries)
233  linearizer.fitParamsErr[ampName] = np.zeros(pEntries)
234  linearizer.fitChiSq[ampName] = np.nan
235  self.log.warn("Amp %s has no usable PTC information. Skipping!", ampName)
236  continue
237 
238  if (len(inputPtc.expIdMask[ampName]) == 0) or self.config.ignorePtcMask:
239  self.log.warn(f"Mask not found for {ampName} in non-linearity fit. Using all points.")
240  mask = np.repeat(True, len(inputPtc.expIdMask[ampName]))
241  else:
242  mask = np.array(inputPtc.expIdMask[ampName], dtype=bool)
243 
244  inputAbscissa = np.array(inputPtc.rawExpTimes[ampName])[mask]
245  inputOrdinate = np.array(inputPtc.rawMeans[ampName])[mask]
246 
247  # Determine proxy-to-linear-flux transformation
248  fluxMask = inputOrdinate < self.config.maxLinearAdu
249  lowMask = inputOrdinate > self.config.minLinearAdu
250  fluxMask = fluxMask & lowMask
251  linearAbscissa = inputAbscissa[fluxMask]
252  linearOrdinate = inputOrdinate[fluxMask]
253 
254  linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 100.0], linearAbscissa,
255  linearOrdinate, funcPolynomial)
256  # Convert this proxy-to-flux fit into an expected linear flux
257  linearOrdinate = linearFit[0] + linearFit[1] * inputAbscissa
258 
259  # Exclude low end outliers
260  threshold = self.config.nSigmaClipLinear * np.sqrt(linearOrdinate)
261  fluxMask = np.abs(inputOrdinate - linearOrdinate) < threshold
262  linearOrdinate = linearOrdinate[fluxMask]
263  fitOrdinate = inputOrdinate[fluxMask]
264  self.debugFitdebugFit('linearFit', inputAbscissa, inputOrdinate, linearOrdinate, fluxMask, ampName)
265  # Do fits
266  if self.config.linearityType in ['Polynomial', 'Squared', 'LookupTable']:
267  polyFit = np.zeros(fitOrder + 1)
268  polyFit[1] = 1.0
269  polyFit, polyFitErr, chiSq, weights = irlsFit(polyFit, linearOrdinate,
270  fitOrdinate, funcPolynomial)
271 
272  # Truncate the polynomial fit
273  k1 = polyFit[1]
274  linearityFit = [-coeff/(k1**order) for order, coeff in enumerate(polyFit)]
275  significant = np.where(np.abs(linearityFit) > 1e-10, True, False)
276  self.log.info(f"Significant polynomial fits: {significant}")
277 
278  modelOrdinate = funcPolynomial(polyFit, linearAbscissa)
279  self.debugFitdebugFit('polyFit', linearAbscissa, fitOrdinate, modelOrdinate, None, ampName)
280 
281  if self.config.linearityType == 'Squared':
282  linearityFit = [linearityFit[2]]
283  elif self.config.linearityType == 'LookupTable':
284  # Use linear part to get time at wich signal is maxAduForLookupTableLinearizer DN
285  tMax = (self.config.maxLookupTableAdu - polyFit[0])/polyFit[1]
286  timeRange = np.linspace(0, tMax, self.config.maxLookupTableAdu)
287  signalIdeal = polyFit[0] + polyFit[1]*timeRange
288  signalUncorrected = funcPolynomial(polyFit, timeRange)
289  lookupTableRow = signalIdeal - signalUncorrected # LinearizerLookupTable has correction
290 
291  linearizer.tableData[tableIndex, :] = lookupTableRow
292  linearityFit = [tableIndex, 0]
293  tableIndex += 1
294  elif self.config.linearityType in ['Spline']:
295  # See discussion in `lsst.ip.isr.linearize.py` before modifying.
296  numPerBin, binEdges = np.histogram(linearOrdinate, bins=fitOrder)
297  with np.errstate(invalid="ignore"):
298  # Algorithm note: With the counts of points per
299  # bin above, the next histogram calculates the
300  # values to put in each bin by weighting each
301  # point by the correction value.
302  values = np.histogram(linearOrdinate, bins=fitOrder,
303  weights=(inputOrdinate[fluxMask] - linearOrdinate))[0]/numPerBin
304 
305  # After this is done, the binCenters are
306  # calculated by weighting by the value we're
307  # binning over. This ensures that widely
308  # spaced/poorly sampled data aren't assigned to
309  # the midpoint of the bin (as could be done using
310  # the binEdges above), but to the weighted mean of
311  # the inputs. Note that both histograms are
312  # scaled by the count per bin to normalize what
313  # the histogram returns (a sum of the points
314  # inside) into an average.
315  binCenters = np.histogram(linearOrdinate, bins=fitOrder,
316  weights=linearOrdinate)[0]/numPerBin
317  values = values[numPerBin > 0]
318  binCenters = binCenters[numPerBin > 0]
319 
320  self.debugFitdebugFit('splineFit', binCenters, np.abs(values), values, None, ampName)
321  interp = afwMath.makeInterpolate(binCenters.tolist(), values.tolist(),
322  afwMath.stringToInterpStyle("AKIMA_SPLINE"))
323  modelOrdinate = linearOrdinate + interp.interpolate(linearOrdinate)
324  self.debugFitdebugFit('splineFit', linearOrdinate, fitOrdinate, modelOrdinate, None, ampName)
325 
326  # If we exclude a lot of points, we may end up with
327  # less than fitOrder points. Pad out the low-flux end
328  # to ensure equal lengths.
329  if len(binCenters) != fitOrder:
330  padN = fitOrder - len(binCenters)
331  binCenters = np.pad(binCenters, (padN, 0), 'linear_ramp',
332  end_values=(binCenters.min() - 1.0, ))
333  # This stores the correction, which is zero at low values.
334  values = np.pad(values, (padN, 0))
335 
336  # Pack the spline into a single array.
337  linearityFit = np.concatenate((binCenters.tolist(), values.tolist())).tolist()
338  polyFit = [0.0]
339  polyFitErr = [0.0]
340  chiSq = np.nan
341  else:
342  polyFit = [0.0]
343  polyFitErr = [0.0]
344  chiSq = np.nan
345  linearityFit = [0.0]
346 
347  linearizer.linearityType[ampName] = self.config.linearityType
348  linearizer.linearityCoeffs[ampName] = np.array(linearityFit)
349  linearizer.linearityBBox[ampName] = amp.getBBox()
350  linearizer.fitParams[ampName] = np.array(polyFit)
351  linearizer.fitParamsErr[ampName] = np.array(polyFitErr)
352  linearizer.fitChiSq[ampName] = chiSq
353 
354  image = afwImage.ImageF(len(inputOrdinate), 1)
355  image.getArray()[:, :] = inputOrdinate
356  linearizeFunction = linearizer.getLinearityTypeByName(linearizer.linearityType[ampName])
357  linearizeFunction()(image,
358  **{'coeffs': linearizer.linearityCoeffs[ampName],
359  'table': linearizer.tableData,
360  'log': linearizer.log})
361  linearizeModel = image.getArray()[0, :]
362 
363  self.debugFitdebugFit('solution', inputOrdinate[fluxMask], linearOrdinate,
364  linearizeModel[fluxMask], None, ampName)
365 
366  linearizer.hasLinearity = True
367  linearizer.validate()
368  linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
369  linearizer.updateMetadata(setDate=True, setCalibId=True)
370  provenance = IsrProvenance(calibType='linearizer')
371 
372  return pipeBase.Struct(
373  outputLinearizer=linearizer,
374  outputProvenance=provenance,
375  )
376 
377  def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
378  """Debug method for linearity fitting.
379 
380  Parameters
381  ----------
382  stepname : `str`
383  A label to use to check if we care to debug at a given
384  line of code.
385  xVector : `numpy.array`
386  The values to use as the independent variable in the
387  linearity fit.
388  yVector : `numpy.array`
389  The values to use as the dependent variable in the
390  linearity fit.
391  yModel : `numpy.array`
392  The values to use as the linearized result.
393  mask : `numpy.array` [ `bool` ], optional
394  A mask to indicate which entries of ``xVector`` and
395  ``yVector`` to keep.
396  ampName : `str`
397  Amplifier name to lookup linearity correction values.
398 
399  """
400  frame = getDebugFrame(self._display, stepname)
401  if frame:
402  import matplotlib.pyplot as plt
403  fig, axs = plt.subplots(2)
404 
405  if mask is None:
406  mask = np.ones_like(xVector, dtype=bool)
407 
408  fig.suptitle(f"{stepname} {ampName} {self.config.linearityType}")
409  if stepname == 'linearFit':
410  axs[0].set_xlabel("Input Abscissa (time or mondiode)")
411  axs[0].set_ylabel("Input Ordinate (flux)")
412  axs[1].set_xlabel("Linear Ordinate (linear flux)")
413  axs[1].set_ylabel("Flux Difference: (input - linear)")
414  elif stepname in ('polyFit', 'splineFit'):
415  axs[0].set_xlabel("Linear Abscissa (linear flux)")
416  axs[0].set_ylabel("Input Ordinate (flux)")
417  axs[1].set_xlabel("Linear Ordinate (linear flux)")
418  axs[1].set_ylabel("Flux Difference: (input - full model fit)")
419  elif stepname == 'solution':
420  axs[0].set_xlabel("Input Abscissa (time or mondiode)")
421  axs[0].set_ylabel("Linear Ordinate (linear flux)")
422  axs[1].set_xlabel("Model flux (linear flux)")
423  axs[1].set_ylabel("Flux Difference: (linear - model)")
424 
425  axs[0].set_yscale('log')
426  axs[0].set_xscale('log')
427  axs[0].scatter(xVector, yVector)
428  axs[0].scatter(xVector[~mask], yVector[~mask], c='red', marker='x')
429  axs[1].set_xscale('log')
430 
431  axs[1].scatter(yModel, yVector[mask] - yModel)
432  fig.show()
433 
434  prompt = "Press Enter or c to continue [chpx]..."
435  while True:
436  ans = input(prompt).lower()
437  if ans in ("", " ", "c",):
438  break
439  elif ans in ("p", ):
440  import pdb
441  pdb.set_trace()
442  elif ans in ("h", ):
443  print("[h]elp [c]ontinue [p]db")
444  elif ans in ('x', ):
445  exit()
446  plt.close()
447 
448 
449 class MeasureLinearityConfig(pexConfig.Config):
450  solver = pexConfig.ConfigurableField(
451  target=LinearitySolveTask,
452  doc="Task to convert PTC data to linearity solutions.",
453  )
454 
455 
456 class MeasureLinearityTask(pipeBase.CmdLineTask):
457  """Stand alone Gen2 linearity measurement.
458 
459  This class wraps the Gen3 linearity task to allow it to be run as
460  a Gen2 CmdLineTask.
461  """
462  ConfigClass = MeasureLinearityConfig
463  _DefaultName = "measureLinearity"
464 
465  def __init__(self, **kwargs):
466  super().__init__(**kwargs)
467  self.makeSubtask("solver")
468 
469  def runDataRef(self, dataRef):
470  """Run new linearity code for gen2.
471 
472  Parameters
473  ----------
474  dataRef : `lsst.daf.persistence.ButlerDataRef`
475  Input dataref for the photon transfer curve data.
476 
477  Returns
478  -------
479  results : `lsst.pipe.base.Struct`
480  The results struct containing:
481 
482  ``outputLinearizer`` : `lsst.ip.isr.Linearizer`
483  Final linearizer calibration.
484  ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
485  Provenance data for the new calibration.
486  """
487  ptc = dataRef.get('photonTransferCurveDataset')
488  camera = dataRef.get('camera')
489  inputDims = dataRef.dataId # This is the closest gen2 has.
490  linearityResults = self.solver.run(ptc, camera=camera, inputDims=inputDims)
491 
492  inputDims['calibDate'] = linearityResults.outputLinearizer.getMetadata().get('CALIBDATE')
493  butler = dataRef.getButler()
494  butler.put(linearityResults.outputLinearizer, "linearizer", inputDims)
495  return linearityResults
def runQuantum(self, butlerQC, inputRefs, outputRefs)
Definition: linearity.py:133
def run(self, inputPtc, dummy, camera, inputDims)
Definition: linearity.py:153
def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName)
Definition: linearity.py:377
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:257
std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
A factory function to make Interpolate objects.
Definition: Interpolate.cc:343
def irlsFit(initialParams, dataX, dataY, function, weightsY=None, weightType='Cauchy')
Definition: utils.py:327
def funcPolynomial(pars, x)
Definition: utils.py:558
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95