Inheritance diagram for lsst.cp.pipe.linearity.LinearitySolveTask:

Public Member Functions
def	runQuantum (self, butlerQC, inputRefs, outputRefs)

def	run (self, inputPtc, dummy, camera, inputDims)

def	debugFit (self, stepname, xVector, yVector, yModel, mask, ampName)

Static Public Attributes
	ConfigClass = LinearitySolveConfig

Detailed Description

Fit the linearity from the PTC dataset.

Definition at line 127 of file linearity.py.

Member Function Documentation

◆ debugFit()

def lsst.cp.pipe.linearity.LinearitySolveTask.debugFit	(	self,
		stepname,
		xVector,
		yVector,
		yModel,
		mask,
		ampName
	)

Debug method for linearity fitting.

Parameters
----------
stepname : `str`
    A label to use to check if we care to debug at a given
    line of code.
xVector : `numpy.array`
    The values to use as the independent variable in the
    linearity fit.
yVector : `numpy.array`
    The values to use as the dependent variable in the
    linearity fit.
yModel : `numpy.array`
    The values to use as the linearized result.
mask : `numpy.array` [ `bool` ], optional
    A mask to indicate which entries of ``xVector`` and
    ``yVector`` to keep.
ampName : `str`
    Amplifier name to lookup linearity correction values.

Definition at line 377 of file linearity.py.

     def debugFit(self, stepname, xVector, yVector, yModel, mask, ampName):
         """Debug method for linearity fitting.
  
         Parameters
         ----------
         stepname : `str`
             A label to use to check if we care to debug at a given
             line of code.
         xVector : `numpy.array`
             The values to use as the independent variable in the
             linearity fit.
         yVector : `numpy.array`
             The values to use as the dependent variable in the
             linearity fit.
         yModel : `numpy.array`
             The values to use as the linearized result.
         mask : `numpy.array` [ `bool` ], optional
             A mask to indicate which entries of ``xVector`` and
             ``yVector`` to keep.
         ampName : `str`
             Amplifier name to lookup linearity correction values.
  
         """
         frame = getDebugFrame(self._display, stepname)
         if frame:
             import matplotlib.pyplot as plt
             fig, axs = plt.subplots(2)
  
             if mask is None:
                 mask = np.ones_like(xVector, dtype=bool)
  
             fig.suptitle(f"{stepname} {ampName} {self.config.linearityType}")
             if stepname == 'linearFit':
                 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
                 axs[0].set_ylabel("Input Ordinate (flux)")
                 axs[1].set_xlabel("Linear Ordinate (linear flux)")
                 axs[1].set_ylabel("Flux Difference: (input - linear)")
             elif stepname in ('polyFit', 'splineFit'):
                 axs[0].set_xlabel("Linear Abscissa (linear flux)")
                 axs[0].set_ylabel("Input Ordinate (flux)")
                 axs[1].set_xlabel("Linear Ordinate (linear flux)")
                 axs[1].set_ylabel("Flux Difference: (input - full model fit)")
             elif stepname == 'solution':
                 axs[0].set_xlabel("Input Abscissa (time or mondiode)")
                 axs[0].set_ylabel("Linear Ordinate (linear flux)")
                 axs[1].set_xlabel("Model flux (linear flux)")
                 axs[1].set_ylabel("Flux Difference: (linear - model)")
  
             axs[0].set_yscale('log')
             axs[0].set_xscale('log')
             axs[0].scatter(xVector, yVector)
             axs[0].scatter(xVector[~mask], yVector[~mask], c='red', marker='x')
             axs[1].set_xscale('log')
  
             axs[1].scatter(yModel, yVector[mask] - yModel)
             fig.show()
  
             prompt = "Press Enter or c to continue [chpx]..."
             while True:
                 ans = input(prompt).lower()
                 if ans in ("", " ", "c",):
                     break
                 elif ans in ("p", ):
                     import pdb
                     pdb.set_trace()
                 elif ans in ("h", ):
                     print("[h]elp [c]ontinue [p]db")
                 elif ans in ('x', ):
                     exit()
             plt.close()
  
  

◆ run()

def lsst.cp.pipe.linearity.LinearitySolveTask.run	(	self,
		inputPtc,
		dummy,
		camera,
		inputDims
	)

Fit non-linearity to PTC data, returning the correct Linearizer
object.

Parameters
----------
inputPtc : `lsst.cp.pipe.PtcDataset`
    Pre-measured PTC dataset.
dummy : `lsst.afw.image.Exposure`
    The exposure used to select the appropriate PTC dataset.
    In almost all circumstances, one of the input exposures
    used to generate the PTC dataset is the best option.
camera : `lsst.afw.cameraGeom.Camera`
    Camera geometry.
inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
    DataIds to use to populate the output calibration.

Returns
-------
results : `lsst.pipe.base.Struct`
    The results struct containing:

    ``outputLinearizer`` : `lsst.ip.isr.Linearizer`
        Final linearizer calibration.
    ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
        Provenance data for the new calibration.

Notes
-----
This task currently fits only polynomial-defined corrections,
where the correction coefficients are defined such that:
    corrImage = uncorrImage + sum_i c_i uncorrImage^(2 + i)
These `c_i` are defined in terms of the direct polynomial fit:
    meanVector ~ P(x=timeVector) = sum_j k_j x^j
such that c_(j-2) = -k_j/(k_1^j) in units of DN^(1-j) (c.f.,
Eq. 37 of 2003.05978). The `config.polynomialOrder` or
`config.splineKnots` define the maximum order of x^j to fit.
As k_0 and k_1 are degenerate with bias level and gain, they
are not included in the non-linearity correction.

Definition at line 153 of file linearity.py.

     def run(self, inputPtc, dummy, camera, inputDims):
         """Fit non-linearity to PTC data, returning the correct Linearizer
         object.
  
         Parameters
         ----------
         inputPtc : `lsst.cp.pipe.PtcDataset`
             Pre-measured PTC dataset.
         dummy : `lsst.afw.image.Exposure`
             The exposure used to select the appropriate PTC dataset.
             In almost all circumstances, one of the input exposures
             used to generate the PTC dataset is the best option.
         camera : `lsst.afw.cameraGeom.Camera`
             Camera geometry.
         inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
             DataIds to use to populate the output calibration.
  
         Returns
         -------
         results : `lsst.pipe.base.Struct`
             The results struct containing:
  
             ``outputLinearizer`` : `lsst.ip.isr.Linearizer`
                 Final linearizer calibration.
             ``outputProvenance`` : `lsst.ip.isr.IsrProvenance`
                 Provenance data for the new calibration.
  
         Notes
         -----
         This task currently fits only polynomial-defined corrections,
         where the correction coefficients are defined such that:
             corrImage = uncorrImage + sum_i c_i uncorrImage^(2 + i)
         These `c_i` are defined in terms of the direct polynomial fit:
             meanVector ~ P(x=timeVector) = sum_j k_j x^j
         such that c_(j-2) = -k_j/(k_1^j) in units of DN^(1-j) (c.f.,
         Eq. 37 of 2003.05978). The `config.polynomialOrder` or
         `config.splineKnots` define the maximum order of x^j to fit.
         As k_0 and k_1 are degenerate with bias level and gain, they
         are not included in the non-linearity correction.
         """
         if len(dummy) == 0:
             self.log.warn("No dummy exposure found.")
  
         detector = camera[inputDims['detector']]
         if self.config.linearityType == 'LookupTable':
             table = np.zeros((len(detector), self.config.maxLookupTableAdu), dtype=np.float32)
             tableIndex = 0
         else:
             table = None
             tableIndex = None  # This will fail if we increment it.
  
         if self.config.linearityType == 'Spline':
             fitOrder = self.config.splineKnots
         else:
             fitOrder = self.config.polynomialOrder
  
         # Initialize the linearizer.
         linearizer = Linearizer(detector=detector, table=table, log=self.log)
  
         for i, amp in enumerate(detector):
             ampName = amp.getName()
             if ampName in inputPtc.badAmps:
                 nEntries = 1
                 pEntries = 1
                 if self.config.linearityType in ['Polynomial']:
                     nEntries = fitOrder + 1
                     pEntries = fitOrder + 1
                 elif self.config.linearityType in ['Spline']:
                     nEntries = fitOrder * 2
                 elif self.config.linearityType in ['Squared', 'None']:
                     nEntries = 1
                     pEntries = fitOrder + 1
                 elif self.config.linearityType in ['LookupTable']:
                     nEntries = 2
                     pEntries = fitOrder + 1
  
                 linearizer.linearityType[ampName] = "None"
                 linearizer.linearityCoeffs[ampName] = np.zeros(nEntries)
                 linearizer.linearityBBox[ampName] = amp.getBBox()
                 linearizer.fitParams[ampName] = np.zeros(pEntries)
                 linearizer.fitParamsErr[ampName] = np.zeros(pEntries)
                 linearizer.fitChiSq[ampName] = np.nan
                 self.log.warn("Amp %s has no usable PTC information.  Skipping!", ampName)
                 continue
  
             if (len(inputPtc.expIdMask[ampName]) == 0) or self.config.ignorePtcMask:
                 self.log.warn(f"Mask not found for {ampName} in non-linearity fit. Using all points.")
                 mask = np.repeat(True, len(inputPtc.expIdMask[ampName]))
             else:
                 mask = np.array(inputPtc.expIdMask[ampName], dtype=bool)
  
             inputAbscissa = np.array(inputPtc.rawExpTimes[ampName])[mask]
             inputOrdinate = np.array(inputPtc.rawMeans[ampName])[mask]
  
             # Determine proxy-to-linear-flux transformation
             fluxMask = inputOrdinate < self.config.maxLinearAdu
             lowMask = inputOrdinate > self.config.minLinearAdu
             fluxMask = fluxMask & lowMask
             linearAbscissa = inputAbscissa[fluxMask]
             linearOrdinate = inputOrdinate[fluxMask]
  
             linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 100.0], linearAbscissa,
                                                               linearOrdinate, funcPolynomial)
             # Convert this proxy-to-flux fit into an expected linear flux
             linearOrdinate = linearFit[0] + linearFit[1] * inputAbscissa
  
             # Exclude low end outliers
             threshold = self.config.nSigmaClipLinear * np.sqrt(linearOrdinate)
             fluxMask = np.abs(inputOrdinate - linearOrdinate) < threshold
             linearOrdinate = linearOrdinate[fluxMask]
             fitOrdinate = inputOrdinate[fluxMask]
             self.debugFit('linearFit', inputAbscissa, inputOrdinate, linearOrdinate, fluxMask, ampName)
             # Do fits
             if self.config.linearityType in ['Polynomial', 'Squared', 'LookupTable']:
                 polyFit = np.zeros(fitOrder + 1)
                 polyFit[1] = 1.0
                 polyFit, polyFitErr, chiSq, weights = irlsFit(polyFit, linearOrdinate,
                                                               fitOrdinate, funcPolynomial)
  
                 # Truncate the polynomial fit
                 k1 = polyFit[1]
                 linearityFit = [-coeff/(k1**order) for order, coeff in enumerate(polyFit)]
                 significant = np.where(np.abs(linearityFit) > 1e-10, True, False)
                 self.log.info(f"Significant polynomial fits: {significant}")
  
                 modelOrdinate = funcPolynomial(polyFit, linearAbscissa)
                 self.debugFit('polyFit', linearAbscissa, fitOrdinate, modelOrdinate, None, ampName)
  
                 if self.config.linearityType == 'Squared':
                     linearityFit = [linearityFit[2]]
                 elif self.config.linearityType == 'LookupTable':
                     # Use linear part to get time at wich signal is maxAduForLookupTableLinearizer DN
                     tMax = (self.config.maxLookupTableAdu - polyFit[0])/polyFit[1]
                     timeRange = np.linspace(0, tMax, self.config.maxLookupTableAdu)
                     signalIdeal = polyFit[0] + polyFit[1]*timeRange
                     signalUncorrected = funcPolynomial(polyFit, timeRange)
                     lookupTableRow = signalIdeal - signalUncorrected  # LinearizerLookupTable has correction
  
                     linearizer.tableData[tableIndex, :] = lookupTableRow
                     linearityFit = [tableIndex, 0]
                     tableIndex += 1
             elif self.config.linearityType in ['Spline']:
                 # See discussion in `lsst.ip.isr.linearize.py` before modifying.
                 numPerBin, binEdges = np.histogram(linearOrdinate, bins=fitOrder)
                 with np.errstate(invalid="ignore"):
                     # Algorithm note: With the counts of points per
                     # bin above, the next histogram calculates the
                     # values to put in each bin by weighting each
                     # point by the correction value.
                     values = np.histogram(linearOrdinate, bins=fitOrder,
                                           weights=(inputOrdinate[fluxMask] - linearOrdinate))[0]/numPerBin
  
                     # After this is done, the binCenters are
                     # calculated by weighting by the value we're
                     # binning over.  This ensures that widely
                     # spaced/poorly sampled data aren't assigned to
                     # the midpoint of the bin (as could be done using
                     # the binEdges above), but to the weighted mean of
                     # the inputs.  Note that both histograms are
                     # scaled by the count per bin to normalize what
                     # the histogram returns (a sum of the points
                     # inside) into an average.
                     binCenters = np.histogram(linearOrdinate, bins=fitOrder,
                                               weights=linearOrdinate)[0]/numPerBin
                     values = values[numPerBin > 0]
                     binCenters = binCenters[numPerBin > 0]
  
                 self.debugFit('splineFit', binCenters, np.abs(values), values, None, ampName)
                 interp = afwMath.makeInterpolate(binCenters.tolist(), values.tolist(),
                                                  afwMath.stringToInterpStyle("AKIMA_SPLINE"))
                 modelOrdinate = linearOrdinate + interp.interpolate(linearOrdinate)
                 self.debugFit('splineFit', linearOrdinate, fitOrdinate, modelOrdinate, None, ampName)
  
                 # If we exclude a lot of points, we may end up with
                 # less than fitOrder points.  Pad out the low-flux end
                 # to ensure equal lengths.
                 if len(binCenters) != fitOrder:
                     padN = fitOrder - len(binCenters)
                     binCenters = np.pad(binCenters, (padN, 0), 'linear_ramp',
                                         end_values=(binCenters.min() - 1.0, ))
                     # This stores the correction, which is zero at low values.
                     values = np.pad(values, (padN, 0))
  
                 # Pack the spline into a single array.
                 linearityFit = np.concatenate((binCenters.tolist(), values.tolist())).tolist()
                 polyFit = [0.0]
                 polyFitErr = [0.0]
                 chiSq = np.nan
             else:
                 polyFit = [0.0]
                 polyFitErr = [0.0]
                 chiSq = np.nan
                 linearityFit = [0.0]
  
             linearizer.linearityType[ampName] = self.config.linearityType
             linearizer.linearityCoeffs[ampName] = np.array(linearityFit)
             linearizer.linearityBBox[ampName] = amp.getBBox()
             linearizer.fitParams[ampName] = np.array(polyFit)
             linearizer.fitParamsErr[ampName] = np.array(polyFitErr)
             linearizer.fitChiSq[ampName] = chiSq
  
             image = afwImage.ImageF(len(inputOrdinate), 1)
             image.getArray()[:, :] = inputOrdinate
             linearizeFunction = linearizer.getLinearityTypeByName(linearizer.linearityType[ampName])
             linearizeFunction()(image,
                                 **{'coeffs': linearizer.linearityCoeffs[ampName],
                                    'table': linearizer.tableData,
                                    'log': linearizer.log})
             linearizeModel = image.getArray()[0, :]
  
             self.debugFit('solution', inputOrdinate[fluxMask], linearOrdinate,
                           linearizeModel[fluxMask], None, ampName)
  
         linearizer.hasLinearity = True
         linearizer.validate()
         linearizer.updateMetadata(camera=camera, detector=detector, filterName='NONE')
         linearizer.updateMetadata(setDate=True, setCalibId=True)
         provenance = IsrProvenance(calibType='linearizer')
  
         return pipeBase.Struct(
             outputLinearizer=linearizer,
             outputProvenance=provenance,
         )
  

◆ runQuantum()

def lsst.cp.pipe.linearity.LinearitySolveTask.runQuantum	(	self,
		butlerQC,
		inputRefs,
		outputRefs
	)

Ensure that the input and output dimensions are passed along.

Parameters
----------
butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
    Butler to operate on.
inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
    Input data refs to load.
ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
    Output data refs to persist.

Definition at line 133 of file linearity.py.

     def runQuantum(self, butlerQC, inputRefs, outputRefs):
         """Ensure that the input and output dimensions are passed along.
  
         Parameters
         ----------
         butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
             Butler to operate on.
         inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
             Input data refs to load.
         ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
             Output data refs to persist.
         """
         inputs = butlerQC.get(inputRefs)
  
         # Use the dimensions to set calib/provenance information.
         inputs['inputDims'] = inputRefs.inputPtc.dataId.byName()
  
         outputs = self.run(**inputs)
         butlerQC.put(outputs, outputRefs)
  

Member Data Documentation

◆ ConfigClass

lsst.cp.pipe.linearity.LinearitySolveTask.ConfigClass = LinearitySolveConfig

static

Definition at line 130 of file linearity.py.

The documentation for this class was generated from the following file:

/j/snowflake/release/lsstsw/stack/lsst-scipipe-0.6.0/Linux64/cp_pipe/22.0.1-18-gb17765a+2264247a6b/python/lsst/cp/pipe/linearity.py

Public Member Functions