23 __all__ = [
'PairedVisitListTaskRunner',
'SingleVisitListTaskRunner',
24 'NonexistentDatasetTaskDataIdContainer',
'parseCmdlineNumberString',
25 'countMaskedPixels',
'checkExpLengthEqual',
'ddict2dict']
29 from scipy.optimize
import leastsq
30 import numpy.polynomial.polynomial
as poly
42 """Calculate weighted reduced chi2.
48 List with measured data.
51 List with modeled data.
53 weightsMeasured : `list`
54 List with weights for the measured data.
57 Number of data points.
60 Number of parameters in the model.
65 redWeightedChi2 : `float`
66 Reduced weighted chi2.
69 wRes = (measured - model)*weightsMeasured
70 return ((wRes*wRes).sum())/(nData-nParsModel)
73 def makeMockFlats(expTime, gain=1.0, readNoiseElectrons=5, fluxElectrons=1000,
74 randomSeedFlat1=1984, randomSeedFlat2=666, powerLawBfParams=[],
76 """Create a pair or mock flats with isrMock.
81 Exposure time of the flats.
83 gain : `float`, optional
86 readNoiseElectrons : `float`, optional
87 Read noise rms, in electrons.
89 fluxElectrons : `float`, optional
90 Flux of flats, in electrons per second.
92 randomSeedFlat1 : `int`, optional
93 Random seed for the normal distrubutions for the mean signal and noise (flat1).
95 randomSeedFlat2 : `int`, optional
96 Random seed for the normal distrubutions for the mean signal and noise (flat2).
98 powerLawBfParams : `list`, optional
99 Parameters for `galsim.cdmodel.PowerLawCD` to simulate the brightter-fatter effect.
101 expId1 : `int`, optional
102 Exposure ID for first flat.
104 expId2 : `int`, optional
105 Exposure ID for second flat.
110 flatExp1 : `lsst.afw.image.exposure.exposure.ExposureF`
111 First exposure of flat field pair.
113 flatExp2 : `lsst.afw.image.exposure.exposure.ExposureF`
114 Second exposure of flat field pair.
118 The parameters of `galsim.cdmodel.PowerLawCD` are `n, r0, t0, rx, tx, r, t, alpha`. For more
119 information about their meaning, see the Galsim documentation
120 https://galsim-developers.github.io/GalSim/_build/html/_modules/galsim/cdmodel.html
121 and Gruen+15 (1501.02802).
123 Example: galsim.cdmodel.PowerLawCD(8, 1.1e-7, 1.1e-7, 1.0e-8, 1.0e-8, 1.0e-9, 1.0e-9, 2.0)
125 flatFlux = fluxElectrons
126 flatMean = flatFlux*expTime
127 readNoise = readNoiseElectrons
129 mockImageConfig = isrMock.IsrMock.ConfigClass()
131 mockImageConfig.flatDrop = 0.99999
132 mockImageConfig.isTrimmed =
True
134 flatExp1 = isrMock.FlatMock(config=mockImageConfig).
run()
135 flatExp2 = flatExp1.clone()
136 (shapeY, shapeX) = flatExp1.getDimensions()
137 flatWidth = np.sqrt(flatMean)
139 rng1 = np.random.RandomState(randomSeedFlat1)
140 flatData1 = rng1.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng1.normal(0.0, readNoise,
142 rng2 = np.random.RandomState(randomSeedFlat2)
143 flatData2 = rng2.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng2.normal(0.0, readNoise,
146 if len(powerLawBfParams):
147 if not len(powerLawBfParams) == 8:
148 raise RuntimeError(
"Wrong number of parameters for `galsim.cdmodel.PowerLawCD`. "
149 f
"Expected 8; passed {len(powerLawBfParams)}.")
150 cd = galsim.cdmodel.PowerLawCD(*powerLawBfParams)
151 tempFlatData1 = galsim.Image(flatData1)
152 temp2FlatData1 = cd.applyForward(tempFlatData1)
154 tempFlatData2 = galsim.Image(flatData2)
155 temp2FlatData2 = cd.applyForward(tempFlatData2)
157 flatExp1.image.array[:] = temp2FlatData1.array/gain
158 flatExp2.image.array[:] = temp2FlatData2.array/gain
160 flatExp1.image.array[:] = flatData1/gain
161 flatExp2.image.array[:] = flatData2/gain
166 flatExp1.getInfo().setVisitInfo(visitInfoExp1)
167 flatExp2.getInfo().setVisitInfo(visitInfoExp2)
169 return flatExp1, flatExp2
173 """Count the number of pixels in a given mask plane."""
174 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
175 nPix = np.where(np.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
180 """Subclass of TaskRunner for handling intrinsically paired visits.
182 This transforms the processed arguments generated by the ArgumentParser
183 into the arguments expected by tasks which take visit pairs for their
186 Such tasks' run() methods tend to take two arguments,
187 one of which is the dataRef (as usual), and the other is the list
188 of visit-pairs, in the form of a list of tuples.
189 This list is supplied on the command line as documented,
190 and this class parses that, and passes the parsed version
193 See pipeBase.TaskRunner for more information.
198 """Parse the visit list and pass through explicitly."""
200 for visitStringPair
in parsedCmd.visitPairs:
201 visitStrings = visitStringPair.split(
",")
202 if len(visitStrings) != 2:
203 raise RuntimeError(
"Found {} visits in {} instead of 2".
format(len(visitStrings),
206 visits = [int(visit)
for visit
in visitStrings]
208 raise RuntimeError(
"Could not parse {} as two integer visit numbers".
format(visitStringPair))
209 visitPairs.append(visits)
211 return pipeBase.TaskRunner.getTargetList(parsedCmd, visitPairs=visitPairs, **kwargs)
215 """Parse command line numerical expression sytax and return as list of int
217 Take an input of the form "'1..5:2^123..126'" as a string, and return
218 a list of ints as [1, 3, 5, 123, 124, 125, 126]
221 for subString
in inputString.split(
"^"):
222 mat = re.search(
r"^(\d+)\.\.(\d+)(?::(\d+))?$", subString)
224 v1 = int(mat.group(1))
225 v2 = int(mat.group(2))
227 v3 = int(v3)
if v3
else 1
228 for v
in range(v1, v2 + 1, v3):
229 outList.append(int(v))
231 outList.append(int(subString))
236 """Subclass of TaskRunner for tasks requiring a list of visits per dataRef.
238 This transforms the processed arguments generated by the ArgumentParser
239 into the arguments expected by tasks which require a list of visits
240 to be supplied for each dataRef, as is common in `lsst.cp.pipe` code.
242 Such tasks' run() methods tend to take two arguments,
243 one of which is the dataRef (as usual), and the other is the list
245 This list is supplied on the command line as documented,
246 and this class parses that, and passes the parsed version
249 See `lsst.pipe.base.TaskRunner` for more information.
254 """Parse the visit list and pass through explicitly."""
257 assert len(parsedCmd.visitList) == 1,
'visitList parsing assumptions violated'
260 return pipeBase.TaskRunner.getTargetList(parsedCmd, visitList=visits, **kwargs)
264 """A DataIdContainer for the tasks for which the output does
268 """Compute refList based on idList.
270 This method must be defined as the dataset does not exist before this
276 Results of parsing the command-line.
280 Not called if ``add_id_argument`` called
281 with ``doMakeDataRefList=False``.
282 Note that this is almost a copy-and-paste of the vanilla
283 implementation, but without checking if the datasets already exist,
284 as this task exists to make them.
286 if self.datasetType
is None:
287 raise RuntimeError(
"Must call setDatasetType first")
288 butler = namespace.butler
289 for dataId
in self.idList:
290 refList =
list(butler.subset(datasetType=self.datasetType, level=self.level, dataId=dataId))
294 namespace.log.warn(
"No data found for dataId=%s", dataId)
296 self.refList += refList
299 def irlsFit(initialParams, dataX, dataY, function, weightsY=None):
300 """Iteratively reweighted least squares fit.
302 This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies
303 weights based on the Cauchy distribution to the fitter. See
304 e.g. Holland and Welsch, 1977, doi:10.1080/03610927708827533
308 initialParams : `list` [`float`]
310 dataX : `numpy.array` [`float`]
312 dataY : `numpy.array` [`float`]
316 weightsY : `numpy.array` [`float`]
317 Weights to apply to the data.
321 polyFit : `list` [`float`]
322 Final best fit parameters.
323 polyFitErr : `list` [`float`]
324 Final errors on fit parameters.
327 weightsY : `list` [`float`]
328 Final weights used for each point.
332 weightsY = np.ones_like(dataX)
334 polyFit, polyFitErr, chiSq =
fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
335 for iteration
in range(10):
337 resid = np.abs(dataY - function(polyFit, dataX)) / np.sqrt(dataY)
338 weightsY = 1.0 / (1.0 + np.sqrt(resid / 2.385))
339 polyFit, polyFitErr, chiSq =
fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
341 return polyFit, polyFitErr, chiSq, weightsY
344 def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):
345 """Do a fit and estimate the parameter errors using using scipy.optimize.leastq.
347 optimize.leastsq returns the fractional covariance matrix. To estimate the
348 standard deviation of the fit parameters, multiply the entries of this matrix
349 by the unweighted reduced chi squared and take the square root of the diagonal elements.
353 initialParams : `list` of `float`
354 initial values for fit parameters. For ptcFitType=POLYNOMIAL, its length
355 determines the degree of the polynomial.
357 dataX : `numpy.array` of `float`
358 Data in the abscissa axis.
360 dataY : `numpy.array` of `float`
361 Data in the ordinate axis.
363 function : callable object (function)
364 Function to fit the data with.
366 weightsY : `numpy.array` of `float`
367 Weights of the data in the ordinate axis.
371 pFitSingleLeastSquares : `list` of `float`
372 List with fitted parameters.
374 pErrSingleLeastSquares : `list` of `float`
375 List with errors for fitted parameters.
377 reducedChiSqSingleLeastSquares : `float`
378 Reduced chi squared, unweighted if weightsY is not provided.
381 weightsY = np.ones(len(dataX))
383 def errFunc(p, x, y, weightsY=None):
385 weightsY = np.ones(len(x))
386 return (function(p, x) - y)*weightsY
388 pFit, pCov, infoDict, errMessage, success = leastsq(errFunc, initialParams,
389 args=(dataX, dataY, weightsY), full_output=1,
392 if (len(dataY) > len(initialParams))
and pCov
is not None:
397 pCov = np.zeros((len(initialParams), len(initialParams)))
399 reducedChiSq = np.nan
402 for i
in range(len(pFit)):
403 errorVec.append(np.fabs(pCov[i][i])**0.5)
405 pFitSingleLeastSquares = pFit
406 pErrSingleLeastSquares = np.array(errorVec)
408 return pFitSingleLeastSquares, pErrSingleLeastSquares, reducedChiSq
411 def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.):
412 """Do a fit using least squares and bootstrap to estimate parameter errors.
414 The bootstrap error bars are calculated by fitting 100 random data sets.
418 initialParams : `list` of `float`
419 initial values for fit parameters. For ptcFitType=POLYNOMIAL, its length
420 determines the degree of the polynomial.
422 dataX : `numpy.array` of `float`
423 Data in the abscissa axis.
425 dataY : `numpy.array` of `float`
426 Data in the ordinate axis.
428 function : callable object (function)
429 Function to fit the data with.
431 weightsY : `numpy.array` of `float`, optional.
432 Weights of the data in the ordinate axis.
434 confidenceSigma : `float`, optional.
435 Number of sigmas that determine confidence interval for the bootstrap errors.
439 pFitBootstrap : `list` of `float`
440 List with fitted parameters.
442 pErrBootstrap : `list` of `float`
443 List with errors for fitted parameters.
445 reducedChiSqBootstrap : `float`
446 Reduced chi squared, unweighted if weightsY is not provided.
449 weightsY = np.ones(len(dataX))
451 def errFunc(p, x, y, weightsY):
453 weightsY = np.ones(len(x))
454 return (function(p, x) - y)*weightsY
457 pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY, weightsY), full_output=0)
460 residuals = errFunc(pFit, dataX, dataY, weightsY)
464 randomDelta = np.random.normal(0., np.fabs(residuals), len(dataY))
465 randomDataY = dataY + randomDelta
466 randomFit, _ = leastsq(errFunc, initialParams,
467 args=(dataX, randomDataY, weightsY), full_output=0)
468 pars.append(randomFit)
469 pars = np.array(pars)
470 meanPfit = np.mean(pars, 0)
473 errPfit = confidenceSigma*np.std(pars, 0)
474 pFitBootstrap = meanPfit
475 pErrBootstrap = errPfit
479 return pFitBootstrap, pErrBootstrap, reducedChiSq
483 """Polynomial function definition
487 Polynomial coefficients. Its length determines the polynomial order.
494 Ordinate array after evaluating polynomial of order len(pars)-1 at `x`.
496 return poly.polyval(x, [*pars])
500 """Single brighter-fatter parameter model for PTC; Equation 16 of Astier+19.
505 Parameters of the model: a00 (brightter-fatter), gain (e/ADU), and noise (e^2).
512 C_00 (variance) in ADU^2.
514 a00, gain, noise = pars
515 return 0.5/(a00*gain*gain)*(np.exp(2*a00*x*gain)-1) + noise/(gain*gain)
519 """Arrange exposures by exposure time.
523 exposureList : `list`[`lsst.afw.image.exposure.exposure.ExposureF`]
524 Input list of exposures.
528 flatsAtExpTime : `dict` [`float`,
529 `list`[`lsst.afw.image.exposure.exposure.ExposureF`]]
530 Dictionary that groups flat-field exposures that have the same
531 exposure time (seconds).
534 for exp
in exposureList:
536 expTime = tempFlat.getInfo().getVisitInfo().getExposureTime()
537 listAtExpTime = flatsAtExpTime.setdefault(expTime, [])
538 listAtExpTime.append(tempFlat)
540 return flatsAtExpTime
544 """Arrange exposures by exposure ID.
546 There is no guarantee that this will properly group exposures, but
547 allows a sequence of flats that have different illumination
548 (despite having the same exposure time) to be processed.
552 exposureList : `list`[`lsst.afw.image.exposure.exposure.ExposureF`]
553 Input list of exposures.
557 flatsAtExpId : `dict` [`float`,
558 `list`[`lsst.afw.image.exposure.exposure.ExposureF`]]
559 Dictionary that groups flat-field exposures sequentially by
565 This algorithm sorts the input exposures by their exposure id, and
566 then assigns each pair of exposures (exp_j, exp_{j+1}) to pair k,
567 such that 2*k = j, where j is the python index of one of the
568 exposures (starting from zero). By checking for the IndexError
569 while appending, we can ensure that there will only ever be fully
573 sortedExposures = sorted(exposureList, key=
lambda exp: exp.getInfo().getVisitInfo().
getExposureId())
575 for jPair, exp
in enumerate(sortedExposures):
578 listAtExpId = flatsAtExpId.setdefault(kPair, [])
580 listAtExpId.append(exp)
581 listAtExpId.append(sortedExposures[jPair + 1])
589 """Check the exposure lengths of two exposures are equal.
593 exp1 : `lsst.afw.image.exposure.ExposureF`
594 First exposure to check
595 exp2 : `lsst.afw.image.exposure.ExposureF`
596 Second exposure to check
597 v1 : `int` or `str`, optional
598 First visit of the visit pair
599 v2 : `int` or `str`, optional
600 Second visit of the visit pair
601 raiseWithMessage : `bool`
602 If True, instead of returning a bool, raise a RuntimeError if exposure
603 times are not equal, with a message about which visits mismatch if the
604 information is available.
609 Raised if the exposure lengths of the two exposures are not equal
611 expTime1 = exp1.getInfo().getVisitInfo().getExposureTime()
612 expTime2 = exp2.getInfo().getVisitInfo().getExposureTime()
613 if expTime1 != expTime2:
615 msg =
"Exposure lengths for visit pairs must be equal. " + \
616 "Found %s and %s" % (expTime1, expTime2)
618 msg +=
" for visit pair %s, %s" % (v1, v2)
619 raise RuntimeError(msg)
625 def validateIsrConfig(isrTask, mandatory=None, forbidden=None, desirable=None, undesirable=None,
626 checkTrim=True, logName=None):
627 """Check that appropriate ISR settings have been selected for the task.
629 Note that this checks that the task itself is configured correctly rather
630 than checking a config.
634 isrTask : `lsst.ip.isr.IsrTask`
635 The task whose config is to be validated
637 mandatory : `iterable` of `str`
638 isr steps that must be set to True. Raises if False or missing
640 forbidden : `iterable` of `str`
641 isr steps that must be set to False. Raises if True, warns if missing
643 desirable : `iterable` of `str`
644 isr steps that should probably be set to True. Warns is False, info if
647 undesirable : `iterable` of `str`
648 isr steps that should probably be set to False. Warns is True, info if
652 Check to ensure the isrTask's assembly subtask is trimming the images.
653 This is a separate config as it is very ugly to do this within the
654 normal configuration lists as it is an option of a sub task.
659 Raised if ``mandatory`` config parameters are False,
660 or if ``forbidden`` parameters are True.
663 Raised if parameter ``isrTask`` is an invalid type.
667 Logs warnings using an isrValidation logger for desirable/undesirable
668 options that are of the wrong polarity or if keys are missing.
670 if not isinstance(isrTask, ipIsr.IsrTask):
671 raise TypeError(f
'Must supply an instance of lsst.ip.isr.IsrTask not {type(isrTask)}')
673 configDict = isrTask.config.toDict()
675 if logName
and isinstance(logName, str):
676 log = lsst.log.getLogger(logName)
678 log = lsst.log.getLogger(
"isrValidation")
681 for configParam
in mandatory:
682 if configParam
not in configDict:
683 raise RuntimeError(f
"Mandatory parameter {configParam} not found in the isr configuration.")
684 if configDict[configParam]
is False:
685 raise RuntimeError(f
"Must set config.isr.{configParam} to True for this task.")
688 for configParam
in forbidden:
689 if configParam
not in configDict:
690 log.warn(f
"Failed to find forbidden key {configParam} in the isr config. The keys in the"
691 " forbidden list should each have an associated Field in IsrConfig:"
692 " check that there is not a typo in this case.")
694 if configDict[configParam]
is True:
695 raise RuntimeError(f
"Must set config.isr.{configParam} to False for this task.")
698 for configParam
in desirable:
699 if configParam
not in configDict:
700 log.info(f
"Failed to find key {configParam} in the isr config. You probably want"
701 " to set the equivalent for your obs_package to True.")
703 if configDict[configParam]
is False:
704 log.warn(f
"Found config.isr.{configParam} set to False for this task."
705 " The cp_pipe Config recommends setting this to True.")
707 for configParam
in undesirable:
708 if configParam
not in configDict:
709 log.info(f
"Failed to find key {configParam} in the isr config. You probably want"
710 " to set the equivalent for your obs_package to False.")
712 if configDict[configParam]
is True:
713 log.warn(f
"Found config.isr.{configParam} set to True for this task."
714 " The cp_pipe Config recommends setting this to False.")
717 if not isrTask.assembleCcd.config.doTrim:
718 raise RuntimeError(
"Must trim when assembling CCDs. Set config.isr.assembleCcd.doTrim to True")
722 """Convert nested default dictionaries to regular dictionaries.
724 This is needed to prevent yaml persistence issues.
729 A possibly nested set of `defaultdict`.
734 A possibly nested set of `dict`.
736 for k, v
in d.items():
737 if isinstance(v, dict):
Information about a single exposure of an imaging camera.
def makeDataRefList(self, namespace)
def getTargetList(parsedCmd, **kwargs)
def getTargetList(parsedCmd, **kwargs)
daf::base::PropertyList * list
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def countMaskedPixels(maskedIm, maskPlane)
def validateIsrConfig(isrTask, mandatory=None, forbidden=None, desirable=None, undesirable=None, checkTrim=True, logName=None)
def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.)
def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None)
def arrangeFlatsByExpId(exposureList)
def parseCmdlineNumberString(inputString)
def arrangeFlatsByExpTime(exposureList)
def checkExpLengthEqual(exp1, exp2, v1=None, v2=None, raiseWithMessage=False)
def makeMockFlats(expTime, gain=1.0, readNoiseElectrons=5, fluxElectrons=1000, randomSeedFlat1=1984, randomSeedFlat2=666, powerLawBfParams=[], expId1=0, expId2=1)
def calculateWeightedReducedChi2(measured, model, weightsMeasured, nData, nParsModel)
def irlsFit(initialParams, dataX, dataY, function, weightsY=None)
def funcPolynomial(pars, x)
def getExposureId(self, dataRef)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)