LSST Applications  22.0.1,22.0.1+01bcf6a671,22.0.1+046ee49490,22.0.1+05c7de27da,22.0.1+0c6914dbf6,22.0.1+1220d50b50,22.0.1+12fd109e95,22.0.1+1a1dd69893,22.0.1+1c910dc348,22.0.1+1ef34551f5,22.0.1+30170c3d08,22.0.1+39153823fd,22.0.1+611137eacc,22.0.1+771eb1e3e8,22.0.1+94e66cc9ed,22.0.1+9a075d06e2,22.0.1+a5ff6e246e,22.0.1+a7db719c1a,22.0.1+ba0d97e778,22.0.1+bfe1ee9056,22.0.1+c4e1e0358a,22.0.1+cc34b8281e,22.0.1+d640e2c0fa,22.0.1+d72a2e677a,22.0.1+d9a6b571bd,22.0.1+e485e9761b,22.0.1+ebe8d3385e
LSST Data Management Base Package
Public Member Functions | Public Attributes | Static Public Attributes | List of all members
lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask Class Reference
Inheritance diagram for lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask:

Public Member Functions

def __init__ (self, *args, **kwargs)
 
def validateIsrConfig (self)
 
def runDataRef (self, dataRef, visitPairs)
 
def estimateGains (self, dataRef, visitPairs)
 
def generateKernel (self, corrs, means, objId, rejectLevel=None)
 
def successiveOverRelax (self, source, maxIter=None, eLevel=None)
 

Public Attributes

 debug
 
 disp1
 
 disp2
 

Static Public Attributes

 RunnerClass = PairedVisitListTaskRunner
 
 ConfigClass = MakeBrighterFatterKernelTaskConfig
 

Detailed Description

Brighter-fatter effect correction-kernel calculation task.

A command line task for calculating the brighter-fatter correction
kernel from pairs of flat-field images (with the same exposure length).

The following operations are performed:

- The configurable isr task is called, which unpersists and assembles the
  raw images, and performs the selected instrument signature removal tasks.
  For the purpose of brighter-fatter coefficient calculation is it
  essential that certain components of isr are *not* performed, and
  recommended that certain others are. The task checks the selected isr
  configuration before it is run, and if forbidden components have been
  selected task will raise, and if recommended ones have not been selected,
  warnings are logged.

- The gain of the each amplifier in the detector is calculated using
  the photon transfer curve (PTC) method and used to correct the images
  so that all calculations are done in units of electrons, and so that the
  level across amplifier boundaries is continuous.
  Outliers in the PTC are iteratively rejected
  before fitting, with the nSigma rejection level set by
  config.nSigmaClipRegression. Individual pixels are ignored in the input
  images the image based on config.nSigmaClipGainCalc.

- Each image is then cross-correlated with the one it's paired with
  (with the pairing defined by the --visit-pairs command line argument),
  which is done either the whole-image to whole-image,
  or amplifier-by-amplifier, depending on config.level.

- Once the cross-correlations have been calculated for each visit pair,
  these are used to generate the correction kernel.
  The maximum lag used, in pixels, and hence the size of the half-size
  of the kernel generated, is given by config.maxLag,
  i.e. a value of 10 will result in a kernel of size 2n-1 = 19x19 pixels.
  Outlier values in these cross-correlations are rejected by using a
  pixel-wise sigma-clipped thresholding to each cross-correlation in
  the visit-pairs-length stack of cross-correlations.
  The number of sigma clipped to is set by config.nSigmaClipKernelGen.

- Once DM-15277 has been completed, a method will exist to calculate the
  empirical correction factor, config.biasCorr.
  TODO: DM-15277 update this part of the docstring once the ticket is done.

Definition at line 330 of file makeBrighterFatterKernel.py.

Constructor & Destructor Documentation

◆ __init__()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.__init__ (   self,
args,
**  kwargs 
)

Definition at line 380 of file makeBrighterFatterKernel.py.

380  def __init__(self, *args, **kwargs):
381  pipeBase.CmdLineTask.__init__(self, *args, **kwargs)
382  self.makeSubtask("isr")
383 
384  self.debug = lsstDebug.Info(__name__)
385  if self.debug.enabled:
386  self.log.info("Running with debug enabled...")
387  # If we're displaying, test it works and save displays for later.
388  # It's worth testing here as displays are flaky and sometimes
389  # can't be contacted, and given processing takes a while,
390  # it's a shame to fail late due to display issues.
391  if self.debug.display:
392  try:
393  afwDisp.setDefaultBackend(self.debug.displayBackend)
394  afwDisp.Display.delAllDisplays()
395  self.disp1 = afwDisp.Display(0, open=True)
396  self.disp2 = afwDisp.Display(1, open=True)
397 
398  im = afwImage.ImageF(1, 1)
399  im.array[:] = [[1]]
400  self.disp1.mtv(im)
401  self.disp1.erase()
402  except NameError:
403  self.debug.display = False
404  self.log.warn('Failed to setup/connect to display! Debug display has been disabled')
405 
406  plt.interactive(False) # stop windows popping up when plotting. When headless, use 'agg' backend too
407  self.validateIsrConfig()
408  self.config.validate()
409  self.config.freeze()
410 
def erase(frame=None)
Definition: ds9.py:96
def mtv(data, frame=None, title="", wcs=None, *args, **kwargs)
Definition: ds9.py:92

Member Function Documentation

◆ estimateGains()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.estimateGains (   self,
  dataRef,
  visitPairs 
)
Estimate the amplifier gains using the specified visits.

Given a dataRef and list of flats of varying intensity,
calculate the gain for each amplifier in the detector
using the photon transfer curve (PTC) method.

The config.fixPtcThroughOrigin option determines whether the iterative
fitting is forced to go through the origin or not.
This defaults to True, fitting var=1/gain * mean.
If set to False then var=1/g * mean + const is fitted.

This is really a photo transfer curve (PTC) gain measurement task.
See DM-14063 for results from of a comparison between
this task's numbers and the gain values in the HSC camera model,
and those measured by the PTC task in eotest.

Parameters
----------
dataRef : `lsst.daf.persistence.butler.Butler.dataRef`
    dataRef for the detector for the flats to be used
visitPairs : `list` of `tuple`
    List of visit-pairs to use, as [(v1,v2), (v3,v4)...]

Returns
-------
gains : `lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterGain`
    Object holding the per-amplifier gains, essentially a
    dict of the as-calculated amplifier gain values, keyed by amp name
nominalGains : `dict` [`str`, `float`]
    Dict of the amplifier gains, as reported by the `detector` object,
    keyed by amplifier name

Definition at line 899 of file makeBrighterFatterKernel.py.

899  def estimateGains(self, dataRef, visitPairs):
900  """Estimate the amplifier gains using the specified visits.
901 
902  Given a dataRef and list of flats of varying intensity,
903  calculate the gain for each amplifier in the detector
904  using the photon transfer curve (PTC) method.
905 
906  The config.fixPtcThroughOrigin option determines whether the iterative
907  fitting is forced to go through the origin or not.
908  This defaults to True, fitting var=1/gain * mean.
909  If set to False then var=1/g * mean + const is fitted.
910 
911  This is really a photo transfer curve (PTC) gain measurement task.
912  See DM-14063 for results from of a comparison between
913  this task's numbers and the gain values in the HSC camera model,
914  and those measured by the PTC task in eotest.
915 
916  Parameters
917  ----------
918  dataRef : `lsst.daf.persistence.butler.Butler.dataRef`
919  dataRef for the detector for the flats to be used
920  visitPairs : `list` of `tuple`
921  List of visit-pairs to use, as [(v1,v2), (v3,v4)...]
922 
923  Returns
924  -------
925  gains : `lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterGain`
926  Object holding the per-amplifier gains, essentially a
927  dict of the as-calculated amplifier gain values, keyed by amp name
928  nominalGains : `dict` [`str`, `float`]
929  Dict of the amplifier gains, as reported by the `detector` object,
930  keyed by amplifier name
931  """
932  # NB: don't use dataRef.get('raw_detector') due to composites
933  detector = dataRef.get('camera')[dataRef.dataId[self.config.ccdKey]]
934  amps = detector.getAmplifiers()
935  ampNames = [amp.getName() for amp in amps]
936 
937  ampMeans = {key: [] for key in ampNames} # these get turned into np.arrays later
938  ampCoVariances = {key: [] for key in ampNames}
939  ampVariances = {key: [] for key in ampNames}
940 
941  # Loop over the amps in the detector,
942  # calculating a PTC for each amplifier.
943  # The amplifier iteration is performed in _calcMeansAndVars()
944  # NB: no gain correction is applied
945  for visPairNum, visPair in enumerate(visitPairs):
946  _means, _vars, _covars = self._calcMeansAndVars(dataRef, visPair[0], visPair[1])
947 
948  # Do sanity checks; if these are failed more investigation is needed
949  breaker = 0
950  for amp in detector:
951  ampName = amp.getName()
952  if _means[ampName]*10 < _vars[ampName] or _means[ampName]*10 < _covars[ampName]:
953  msg = 'Sanity check failed; check visit pair %s amp %s' % (visPair, ampName)
954  self.log.warn(msg)
955  breaker += 1
956  if breaker:
957  continue
958 
959  # having made sanity checks
960  # pull the values out into the respective dicts
961  for k in _means.keys(): # keys are necessarily the same
962  if _vars[k]*1.3 < _covars[k] or _vars[k]*0.7 > _covars[k]:
963  self.log.warn('Dropped a value')
964  continue
965  ampMeans[k].append(_means[k])
966  ampVariances[k].append(_vars[k])
967  ampCoVariances[k].append(_covars[k])
968 
969  gains = {}
970  nomGains = {}
971  ptcResults = {}
972  for amp in detector:
973  ampName = amp.getName()
974  if ampMeans[ampName] == []: # all the data was dropped, amp is presumed bad
975  gains[ampName] = 1.0
976  ptcResults[ampName] = (0, 0, 1, 0)
977  continue
978 
979  nomGains[ampName] = amp.getGain()
980  slopeRaw, interceptRaw, rVal, pVal, stdErr = \
981  stats.linregress(np.asarray(ampMeans[ampName]), np.asarray(ampCoVariances[ampName]))
982  slopeFix, _ = self._iterativeRegression(np.asarray(ampMeans[ampName]),
983  np.asarray(ampCoVariances[ampName]),
984  fixThroughOrigin=True)
985  slopeUnfix, intercept = self._iterativeRegression(np.asarray(ampMeans[ampName]),
986  np.asarray(ampCoVariances[ampName]),
987  fixThroughOrigin=False)
988  self.log.info("Slope of raw fit: %s, intercept: %s p value: %s" % (slopeRaw,
989  interceptRaw, pVal))
990  self.log.info("slope of fixed fit: %s, difference vs raw:%s" % (slopeFix,
991  slopeFix - slopeRaw))
992  self.log.info("slope of unfixed fit: %s, difference vs fix:%s" % (slopeUnfix,
993  slopeFix - slopeUnfix))
994  if self.config.fixPtcThroughOrigin:
995  slopeToUse = slopeFix
996  else:
997  slopeToUse = slopeUnfix
998 
999  if self.debug.enabled:
1000  fig = plt.figure()
1001  ax = fig.add_subplot(111)
1002  ax.plot(np.asarray(ampMeans[ampName]),
1003  np.asarray(ampCoVariances[ampName]), linestyle='None', marker='x', label='data')
1004  if self.config.fixPtcThroughOrigin:
1005  ax.plot(np.asarray(ampMeans[ampName]),
1006  np.asarray(ampMeans[ampName])*slopeToUse, label='Fit through origin')
1007  else:
1008  ax.plot(np.asarray(ampMeans[ampName]),
1009  np.asarray(ampMeans[ampName])*slopeToUse + intercept,
1010  label='Fit (intercept unconstrained')
1011 
1012  dataRef.put(fig, "plotBrighterFatterPtc", amp=ampName)
1013  self.log.info('Saved PTC for detector %s amp %s' % (detector.getId(), ampName))
1014  gains[ampName] = 1.0/slopeToUse
1015  # change the fit to use a cubic and match parameters with Lage method
1016  # or better, use the PTC task here too
1017  ptcResults[ampName] = (0, 0, 1, 0)
1018 
1019  return BrighterFatterGain(gains, ptcResults), nomGains
1020 
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33

◆ generateKernel()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.generateKernel (   self,
  corrs,
  means,
  objId,
  rejectLevel = None 
)
Generate the full kernel from a list of cross-correlations and means.

Taking a list of quarter-image, gain-corrected cross-correlations,
do a pixel-wise sigma-clipped mean of each,
and tile into the full-sized kernel image.

Each corr in corrs is one quarter of the full cross-correlation,
and has been gain-corrected. Each mean in means is a tuple of the means
of the two individual images, corresponding to that corr.

Parameters:
-----------
corrs : `list` of `numpy.ndarray`, (Ny, Nx)
    A list of the quarter-image cross-correlations
means : `list` of `tuples` of `floats`
    The means of the input images for each corr in corrs
rejectLevel : `float`, optional
    This is essentially is a sanity check parameter.
    If this condition is violated there is something unexpected
    going on in the image, and it is discarded from the stack before
    the clipped-mean is calculated.
    If not provided then config.xcorrCheckRejectLevel is used

Returns:
--------
kernel : `numpy.ndarray`, (Ny, Nx)
    The output kernel

Definition at line 1270 of file makeBrighterFatterKernel.py.

1270  def generateKernel(self, corrs, means, objId, rejectLevel=None):
1271  """Generate the full kernel from a list of cross-correlations and means.
1272 
1273  Taking a list of quarter-image, gain-corrected cross-correlations,
1274  do a pixel-wise sigma-clipped mean of each,
1275  and tile into the full-sized kernel image.
1276 
1277  Each corr in corrs is one quarter of the full cross-correlation,
1278  and has been gain-corrected. Each mean in means is a tuple of the means
1279  of the two individual images, corresponding to that corr.
1280 
1281  Parameters:
1282  -----------
1283  corrs : `list` of `numpy.ndarray`, (Ny, Nx)
1284  A list of the quarter-image cross-correlations
1285  means : `list` of `tuples` of `floats`
1286  The means of the input images for each corr in corrs
1287  rejectLevel : `float`, optional
1288  This is essentially is a sanity check parameter.
1289  If this condition is violated there is something unexpected
1290  going on in the image, and it is discarded from the stack before
1291  the clipped-mean is calculated.
1292  If not provided then config.xcorrCheckRejectLevel is used
1293 
1294  Returns:
1295  --------
1296  kernel : `numpy.ndarray`, (Ny, Nx)
1297  The output kernel
1298  """
1299  self.log.info('Calculating kernel for %s'%objId)
1300 
1301  if not rejectLevel:
1302  rejectLevel = self.config.xcorrCheckRejectLevel
1303 
1304  if self.config.correlationQuadraticFit:
1305  xcorrList = []
1306  fluxList = []
1307 
1308  for corrNum, ((mean1, mean2), corr) in enumerate(zip(means, corrs)):
1309  msg = 'For item %s, initial corr[0,0] = %g, corr[1,0] = %g'%(corrNum, corr[0, 0], corr[1, 0])
1310  self.log.info(msg)
1311  if self.config.level == 'DETECTOR':
1312  # This is now done in _applyGains() but only if level is not DETECTOR
1313  corr[0, 0] -= (mean1 + mean2)
1314  fullCorr = self._tileArray(corr)
1315 
1316  # Craig Lage says he doesn't understand the negative sign, but it needs to be there
1317  xcorrList.append(-fullCorr / 2.0)
1318  flux = (mean1 + mean2) / 2.0
1319  fluxList.append(flux * flux)
1320  # We're using the linear fit algorithm to find a quadratic fit,
1321  # so we square the x-axis.
1322  # The step below does not need to be done, but is included
1323  # so that correlations can be compared
1324  # directly to existing code. We might want to take it out.
1325  corr /= -1.0*(mean1**2 + mean2**2)
1326 
1327  if not xcorrList:
1328  raise RuntimeError("Cannot generate kernel because all inputs were discarded. "
1329  "Either the data is bad, or config.xcorrCheckRejectLevel is too low")
1330 
1331  # This method fits a quadratic vs flux and keeps only the quadratic term.
1332  meanXcorr = np.zeros_like(fullCorr)
1333  xcorrList = np.asarray(xcorrList)
1334 
1335  for i in range(np.shape(meanXcorr)[0]):
1336  for j in range(np.shape(meanXcorr)[1]):
1337  # Note the i,j inversion. This serves the same function as the transpose step in
1338  # the base code. I don't understand why it is there, but I put it in to be consistent.
1339  slopeRaw, interceptRaw, rVal, pVal, stdErr = stats.linregress(np.asarray(fluxList),
1340  xcorrList[:, j, i])
1341  try:
1342  slope, intercept = self._iterativeRegression(np.asarray(fluxList),
1343  xcorrList[:, j, i],
1344  fixThroughOrigin=True)
1345  msg = "(%s,%s):Slope of raw fit: %s, intercept: %s p value: %s" % (i, j, slopeRaw,
1346  interceptRaw, pVal)
1347  self.log.debug(msg)
1348  self.log.debug("(%s,%s):Slope of fixed fit: %s" % (i, j, slope))
1349 
1350  meanXcorr[i, j] = slope
1351  except ValueError:
1352  meanXcorr[i, j] = slopeRaw
1353 
1354  msg = f"i={i}, j={j}, slope = {slope:.6g}, slopeRaw = {slopeRaw:.6g}"
1355  self.log.debug(msg)
1356  self.log.info('Quad Fit meanXcorr[0,0] = %g, meanXcorr[1,0] = %g'%(meanXcorr[8, 8],
1357  meanXcorr[9, 8]))
1358 
1359  else:
1360  # Try to average over a set of possible inputs.
1361  # This generates a simple function of the kernel that
1362  # should be constant across the images, and averages that.
1363  xcorrList = []
1364  sctrl = afwMath.StatisticsControl()
1365  sctrl.setNumSigmaClip(self.config.nSigmaClipKernelGen)
1366 
1367  for corrNum, ((mean1, mean2), corr) in enumerate(zip(means, corrs)):
1368  corr[0, 0] -= (mean1 + mean2)
1369  if corr[0, 0] > 0:
1370  self.log.warn('Skipped item %s due to unexpected value of (variance-mean)' % corrNum)
1371  continue
1372  corr /= -1.0*(mean1**2 + mean2**2)
1373 
1374  fullCorr = self._tileArray(corr)
1375 
1376  xcorrCheck = np.abs(np.sum(fullCorr))/np.sum(np.abs(fullCorr))
1377  if xcorrCheck > rejectLevel:
1378  self.log.warn("Sum of the xcorr is unexpectedly high. Investigate item num %s for %s. \n"
1379  "value = %s" % (corrNum, objId, xcorrCheck))
1380  continue
1381  xcorrList.append(fullCorr)
1382 
1383  if not xcorrList:
1384  raise RuntimeError("Cannot generate kernel because all inputs were discarded. "
1385  "Either the data is bad, or config.xcorrCheckRejectLevel is too low")
1386 
1387  # stack the individual xcorrs and apply a per-pixel clipped-mean
1388  meanXcorr = np.zeros_like(fullCorr)
1389  xcorrList = np.transpose(xcorrList)
1390  for i in range(np.shape(meanXcorr)[0]):
1391  for j in range(np.shape(meanXcorr)[1]):
1392  meanXcorr[i, j] = afwMath.makeStatistics(xcorrList[i, j],
1393  afwMath.MEANCLIP, sctrl).getValue()
1394 
1395  if self.config.correlationModelRadius < (meanXcorr.shape[0] - 1) / 2:
1396  sumToInfinity = self._buildCorrelationModel(meanXcorr, self.config.correlationModelRadius,
1397  self.config.correlationModelSlope)
1398  self.log.info("SumToInfinity = %s" % sumToInfinity)
1399  else:
1400  sumToInfinity = 0.0
1401  if self.config.forceZeroSum:
1402  self.log.info("Forcing sum of correlation matrix to zero")
1403  meanXcorr = self._forceZeroSum(meanXcorr, sumToInfinity)
1404 
1405  return meanXcorr, self.successiveOverRelax(meanXcorr)
1406 
Pass parameters to a Statistics object.
Definition: Statistics.h:93
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:354

◆ runDataRef()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.runDataRef (   self,
  dataRef,
  visitPairs 
)
Run the brighter-fatter measurement task.

For a dataRef (which is each detector here),
and given a list of visit pairs, calculate the
brighter-fatter kernel for the detector.

Parameters
----------
dataRef : `list` of `lsst.daf.persistence.ButlerDataRef`
    dataRef for the detector for the visits to be fit.
visitPairs : `iterable` of `tuple` of `int`
    Pairs of visit numbers to be processed together

Definition at line 460 of file makeBrighterFatterKernel.py.

460  def runDataRef(self, dataRef, visitPairs):
461  """Run the brighter-fatter measurement task.
462 
463  For a dataRef (which is each detector here),
464  and given a list of visit pairs, calculate the
465  brighter-fatter kernel for the detector.
466 
467  Parameters
468  ----------
469  dataRef : `list` of `lsst.daf.persistence.ButlerDataRef`
470  dataRef for the detector for the visits to be fit.
471  visitPairs : `iterable` of `tuple` of `int`
472  Pairs of visit numbers to be processed together
473  """
474  np.random.seed(0) # used in the PTC fit bootstrap
475 
476  # setup necessary objects
477  # NB: don't use dataRef.get('raw_detector')
478  # this currently doesn't work for composites because of the way
479  # composite objects (i.e. LSST images) are handled/constructed
480  # these need to be retrieved from the camera and dereferenced
481  # rather than accessed directly
482  detNum = dataRef.dataId[self.config.ccdKey]
483  detector = dataRef.get('camera')[dataRef.dataId[self.config.ccdKey]]
484  amps = detector.getAmplifiers()
485  ampNames = [amp.getName() for amp in amps]
486 
487  if self.config.level == 'DETECTOR':
488  kernels = {detNum: []}
489  means = {detNum: []}
490  xcorrs = {detNum: []}
491  meanXcorrs = {detNum: []}
492  elif self.config.level == 'AMP':
493  kernels = {key: [] for key in ampNames}
494  means = {key: [] for key in ampNames}
495  xcorrs = {key: [] for key in ampNames}
496  meanXcorrs = {key: [] for key in ampNames}
497  else:
498  raise RuntimeError("Unsupported level: {}".format(self.config.level))
499 
500  # we must be able to get the gains one way or the other, so check early
501  if not self.config.doCalcGains:
502  deleteMe = None
503  try:
504  deleteMe = dataRef.get('photonTransferCurveDataset')
505  except butlerExceptions.NoResults:
506  try:
507  deleteMe = dataRef.get('brighterFatterGain')
508  except butlerExceptions.NoResults:
509  pass
510  if not deleteMe:
511  raise RuntimeError("doCalcGains == False and gains could not be got from butler") from None
512  else:
513  del deleteMe
514 
515  # if the level is DETECTOR we need to have the gains first so that each
516  # amp can be gain corrected in order to treat the detector as a single
517  # imaging area. However, if the level is AMP we can wait, calculate
518  # the correlations and correct for the gains afterwards
519  if self.config.level == 'DETECTOR':
520  if self.config.doCalcGains:
521  self.log.info('Computing gains for detector %s' % detNum)
522  gains, nomGains = self.estimateGains(dataRef, visitPairs)
523  dataRef.put(gains, datasetType='brighterFatterGain')
524  self.log.debug('Finished gain estimation for detector %s' % detNum)
525  else:
526  gains = dataRef.get('brighterFatterGain')
527  if not gains:
528  raise RuntimeError('Failed to retrieved gains for detector %s' % detNum)
529  self.log.info('Retrieved stored gain for detector %s' % detNum)
530  self.log.debug('Detector %s has gains %s' % (detNum, gains))
531  else: # we fake the gains as 1 for now, and correct later
532  gains = BrighterFatterGain({}, {})
533  for ampName in ampNames:
534  gains.gains[ampName] = 1.0
535  # We'll use the ptc.py code to calculate the gains, so we set this up
536  ptcConfig = MeasurePhotonTransferCurveTaskConfig()
537  ptcConfig.isrForbiddenSteps = []
538  ptcConfig.doFitBootstrap = True
539  ptcConfig.ptcFitType = 'POLYNOMIAL' # default Astier doesn't work for gain correction
540  ptcConfig.polynomialFitDegree = 3
541  ptcConfig.minMeanSignal = self.config.minMeanSignal
542  ptcConfig.maxMeanSignal = self.config.maxMeanSignal
543  ptcTask = MeasurePhotonTransferCurveTask(config=ptcConfig)
544  ptcDataset = PhotonTransferCurveDataset(ampNames)
545 
546  # Loop over pairs of visits
547  # calculating the cross-correlations at the required level
548  for (v1, v2) in visitPairs:
549  dataRef.dataId['expId'] = v1
550  exp1 = self.isr.runDataRef(dataRef).exposure
551  dataRef.dataId['expId'] = v2
552  exp2 = self.isr.runDataRef(dataRef).exposure
553  del dataRef.dataId['expId']
554  checkExpLengthEqual(exp1, exp2, v1, v2, raiseWithMessage=True)
555 
556  self.log.info('Preparing images for cross-correlation calculation for detector %s' % detNum)
557  # note the shape of these returns depends on level
558  _scaledMaskedIms1, _means1 = self._makeCroppedExposures(exp1, gains, self.config.level)
559  _scaledMaskedIms2, _means2 = self._makeCroppedExposures(exp2, gains, self.config.level)
560 
561  # Compute the cross-correlation and means
562  # at the appropriate config.level:
563  # - "DETECTOR": one key, so compare the two visits to each other
564  # - "AMP": n_amp keys, comparing each amplifier of one visit
565  # to the same amplifier in the visit its paired with
566  for det_object in _scaledMaskedIms1.keys(): # det_object is ampName or detName depending on level
567  self.log.debug("Calculating correlations for %s" % det_object)
568  _xcorr, _mean = self._crossCorrelate(_scaledMaskedIms1[det_object],
569  _scaledMaskedIms2[det_object])
570  xcorrs[det_object].append(_xcorr)
571  means[det_object].append([_means1[det_object], _means2[det_object]])
572  if self.config.level != 'DETECTOR':
573  # Populate the ptcDataset for running fitting in the PTC task
574  expTime = exp1.getInfo().getVisitInfo().getExposureTime()
575  ptcDataset.rawExpTimes[det_object].append(expTime)
576  ptcDataset.rawMeans[det_object].append((_means1[det_object] + _means2[det_object]) / 2.0)
577  ptcDataset.rawVars[det_object].append(_xcorr[0, 0] / 2.0)
578 
579  # TODO: DM-15305 improve debug functionality here.
580  # This is position 1 for the removed code.
581 
582  # Save the raw means and xcorrs so we can look at them before any modifications
583  rawMeans = copy.deepcopy(means)
584  rawXcorrs = copy.deepcopy(xcorrs)
585 
586  # gains are always and only pre-applied for DETECTOR
587  # so for all other levels we now calculate them from the correlations
588  # and apply them
589  if self.config.level != 'DETECTOR':
590  if self.config.doCalcGains: # Run the PTC task for calculating the gains, put results
591  self.log.info('Calculating gains for detector %s using PTC task' % detNum)
592  ptcDataset = ptcTask.fitPtc(ptcDataset, ptcConfig.ptcFitType)
593  dataRef.put(ptcDataset, datasetType='photonTransferCurveDataset')
594  self.log.debug('Finished gain estimation for detector %s' % detNum)
595  else: # load results - confirmed to work much earlier on, so can be relied upon here
596  ptcDataset = dataRef.get('photonTransferCurveDataset')
597 
598  self._applyGains(means, xcorrs, ptcDataset)
599 
600  if self.config.doPlotPtcs:
601  dirname = dataRef.getUri(datasetType='cpPipePlotRoot', write=True)
602  if not os.path.exists(dirname):
603  os.makedirs(dirname)
604  detNum = dataRef.dataId[self.config.ccdKey]
605  filename = f"PTC_det{detNum}.pdf"
606  filenameFull = os.path.join(dirname, filename)
607  with PdfPages(filenameFull) as pdfPages:
608  ptcTask._plotPtc(ptcDataset, ptcConfig.ptcFitType, pdfPages)
609 
610  # having calculated and applied the gains for all code-paths we can now
611  # generate the kernel(s)
612  self.log.info('Generating kernel(s) for %s' % detNum)
613  for det_object in xcorrs.keys(): # looping over either detectors or amps
614  if self.config.level == 'DETECTOR':
615  objId = 'detector %s' % det_object
616  elif self.config.level == 'AMP':
617  objId = 'detector %s AMP %s' % (detNum, det_object)
618 
619  try:
620  meanXcorr, kernel = self.generateKernel(xcorrs[det_object], means[det_object], objId)
621  kernels[det_object] = kernel
622  meanXcorrs[det_object] = meanXcorr
623  except RuntimeError:
624  # bad amps will cause failures here which we want to ignore
625  self.log.warn('RuntimeError during kernel generation for %s' % objId)
626  continue
627 
628  bfKernel = BrighterFatterKernel(self.config.level)
629  bfKernel.means = means
630  bfKernel.rawMeans = rawMeans
631  bfKernel.rawXcorrs = rawXcorrs
632  bfKernel.xCorrs = xcorrs
633  bfKernel.meanXcorrs = meanXcorrs
634  bfKernel.originalLevel = self.config.level
635  try:
636  bfKernel.gain = ptcDataset.gain
637  bfKernel.gainErr = ptcDataset.gainErr
638  bfKernel.noise = ptcDataset.noise
639  bfKernel.noiseErr = ptcDataset.noiseErr
640  except NameError: # we don't have a ptcDataset to store results from
641  pass
642 
643  if self.config.level == 'AMP':
644  bfKernel.ampwiseKernels = kernels
645  ex = self.config.ignoreAmpsForAveraging
646  bfKernel.detectorKernel = bfKernel.makeDetectorKernelFromAmpwiseKernels(detNum, ampsToExclude=ex)
647 
648  elif self.config.level == 'DETECTOR':
649  bfKernel.detectorKernel = kernels
650  else:
651  raise RuntimeError('Invalid level for kernel calculation; this should not be possible.')
652 
653  dataRef.put(bfKernel)
654 
655  self.log.info('Finished generating kernel(s) for %s' % detNum)
656  return pipeBase.Struct(exitStatus=0)
657 
def checkExpLengthEqual(exp1, exp2, v1=None, v2=None, raiseWithMessage=False)
Definition: utils.py:588
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174

◆ successiveOverRelax()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.successiveOverRelax (   self,
  source,
  maxIter = None,
  eLevel = None 
)
An implementation of the successive over relaxation (SOR) method.

A numerical method for solving a system of linear equations
with faster convergence than the Gauss-Seidel method.

Parameters:
-----------
source : `numpy.ndarray`
    The input array.
maxIter : `int`, optional
    Maximum number of iterations to attempt before aborting.
eLevel : `float`, optional
    The target error level at which we deem convergence to have
occurred.

Returns:
--------
output : `numpy.ndarray`
    The solution.

Definition at line 1407 of file makeBrighterFatterKernel.py.

1407  def successiveOverRelax(self, source, maxIter=None, eLevel=None):
1408  """An implementation of the successive over relaxation (SOR) method.
1409 
1410  A numerical method for solving a system of linear equations
1411  with faster convergence than the Gauss-Seidel method.
1412 
1413  Parameters:
1414  -----------
1415  source : `numpy.ndarray`
1416  The input array.
1417  maxIter : `int`, optional
1418  Maximum number of iterations to attempt before aborting.
1419  eLevel : `float`, optional
1420  The target error level at which we deem convergence to have
1421  occurred.
1422 
1423  Returns:
1424  --------
1425  output : `numpy.ndarray`
1426  The solution.
1427  """
1428  if not maxIter:
1429  maxIter = self.config.maxIterSuccessiveOverRelaxation
1430  if not eLevel:
1431  eLevel = self.config.eLevelSuccessiveOverRelaxation
1432 
1433  assert source.shape[0] == source.shape[1], "Input array must be square"
1434  # initialize, and set boundary conditions
1435  func = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
1436  resid = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
1437  rhoSpe = np.cos(np.pi/source.shape[0]) # Here a square grid is assumed
1438 
1439  # Calculate the initial error
1440  for i in range(1, func.shape[0] - 1):
1441  for j in range(1, func.shape[1] - 1):
1442  resid[i, j] = (func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
1443  + func[i + 1, j] - 4*func[i, j] - source[i - 1, j - 1])
1444  inError = np.sum(np.abs(resid))
1445 
1446  # Iterate until convergence
1447  # We perform two sweeps per cycle,
1448  # updating 'odd' and 'even' points separately
1449  nIter = 0
1450  omega = 1.0
1451  dx = 1.0
1452  while nIter < maxIter*2:
1453  outError = 0
1454  if nIter%2 == 0:
1455  for i in range(1, func.shape[0] - 1, 2):
1456  for j in range(1, func.shape[1] - 1, 2):
1457  resid[i, j] = float(func[i, j-1] + func[i, j + 1] + func[i - 1, j]
1458  + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
1459  func[i, j] += omega*resid[i, j]*.25
1460  for i in range(2, func.shape[0] - 1, 2):
1461  for j in range(2, func.shape[1] - 1, 2):
1462  resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
1463  + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
1464  func[i, j] += omega*resid[i, j]*.25
1465  else:
1466  for i in range(1, func.shape[0] - 1, 2):
1467  for j in range(2, func.shape[1] - 1, 2):
1468  resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
1469  + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
1470  func[i, j] += omega*resid[i, j]*.25
1471  for i in range(2, func.shape[0] - 1, 2):
1472  for j in range(1, func.shape[1] - 1, 2):
1473  resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
1474  + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
1475  func[i, j] += omega*resid[i, j]*.25
1476  outError = np.sum(np.abs(resid))
1477  if outError < inError*eLevel:
1478  break
1479  if nIter == 0:
1480  omega = 1.0/(1 - rhoSpe*rhoSpe/2.0)
1481  else:
1482  omega = 1.0/(1 - rhoSpe*rhoSpe*omega/4.0)
1483  nIter += 1
1484 
1485  if nIter >= maxIter*2:
1486  self.log.warn("Failure: SuccessiveOverRelaxation did not converge in %s iterations."
1487  "\noutError: %s, inError: %s," % (nIter//2, outError, inError*eLevel))
1488  else:
1489  self.log.info("Success: SuccessiveOverRelaxation converged in %s iterations."
1490  "\noutError: %s, inError: %s", nIter//2, outError, inError*eLevel)
1491  return func[1: -1, 1: -1]
1492 

◆ validateIsrConfig()

def lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.validateIsrConfig (   self)
Check that appropriate ISR settings are being used
for brighter-fatter kernel calculation.

Definition at line 422 of file makeBrighterFatterKernel.py.

422  def validateIsrConfig(self):
423  """Check that appropriate ISR settings are being used
424  for brighter-fatter kernel calculation."""
425 
426  # How should we handle saturation/bad regions?
427  # 'doSaturationInterpolation': True
428  # 'doNanInterpAfterFlat': False
429  # 'doSaturation': True
430  # 'doSuspect': True
431  # 'doWidenSaturationTrails': True
432  # 'doSetBadRegions': True
433 
434  configDict = self.config.isr.toDict()
435 
436  for configParam in self.config.isrMandatorySteps:
437  if configDict[configParam] is False:
438  raise RuntimeError(f'Must set config.isr.{configParam} to True '
439  'for brighter-fatter kernel calculation')
440 
441  for configParam in self.config.isrForbiddenSteps:
442  if configDict[configParam] is True:
443  raise RuntimeError(f'Must set config.isr.{configParam} to False '
444  'for brighter-fatter kernel calculation')
445 
446  for configParam in self.config.isrDesirableSteps:
447  if configParam not in configDict:
448  self.log.info('Failed to find key %s in the isr config dict. You probably want '
449  'to set the equivalent for your obs_package to True.', configParam)
450  continue
451  if configDict[configParam] is False:
452  self.log.warn('Found config.isr.%s set to False for brighter-fatter kernel calculation. '
453  'It is probably desirable to have this set to True', configParam)
454 
455  # subtask settings
456  if not self.config.isr.assembleCcd.doTrim:
457  raise RuntimeError('Must trim when assembling CCDs. Set config.isr.assembleCcd.doTrim to True')
458 
def validateIsrConfig(isrTask, mandatory=None, forbidden=None, desirable=None, undesirable=None, checkTrim=True, logName=None)
Definition: utils.py:626

Member Data Documentation

◆ ConfigClass

lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.ConfigClass = MakeBrighterFatterKernelTaskConfig
static

Definition at line 377 of file makeBrighterFatterKernel.py.

◆ debug

lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.debug

Definition at line 384 of file makeBrighterFatterKernel.py.

◆ disp1

lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.disp1

Definition at line 395 of file makeBrighterFatterKernel.py.

◆ disp2

lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.disp2

Definition at line 396 of file makeBrighterFatterKernel.py.

◆ RunnerClass

lsst.cp.pipe.makeBrighterFatterKernel.MakeBrighterFatterKernelTask.RunnerClass = PairedVisitListTaskRunner
static

Definition at line 376 of file makeBrighterFatterKernel.py.


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