23 __all__ = [
'PairedVisitListTaskRunner',
'SingleVisitListTaskRunner',
24 'NonexistentDatasetTaskDataIdContainer',
'parseCmdlineNumberString',
25 'countMaskedPixels',
'checkExpLengthEqual',
'ddict2dict']
29 from scipy.optimize
import leastsq
30 import numpy.polynomial.polynomial
as poly
31 from scipy.stats
import norm
43 """Correct measured sigma to account for clipping.
45 If we clip our input data and then measure sigma, then the
46 measured sigma is smaller than the true value because real
47 points beyond the clip threshold have been removed. This is a
48 small (1.5% at nSigClip=3) effect when nSigClip >~ 3, but the
49 default parameters for measure crosstalk use nSigClip=2.0.
50 This causes the measured sigma to be about 15% smaller than
51 real. This formula corrects the issue, for the symmetric case
52 (upper clip threshold equal to lower clip threshold).
57 Number of sigma the measurement was clipped by.
62 Scale factor to increase the measured sigma by.
65 varFactor = 1.0 - (2 * nSigClip * norm.pdf(nSigClip)) / (norm.cdf(nSigClip) - norm.cdf(-nSigClip))
66 return 1.0 / np.sqrt(varFactor)
70 """Calculate weighted reduced chi2.
76 List with measured data.
79 List with modeled data.
81 weightsMeasured : `list`
82 List with weights for the measured data.
85 Number of data points.
88 Number of parameters in the model.
93 redWeightedChi2 : `float`
94 Reduced weighted chi2.
97 wRes = (measured - model)*weightsMeasured
98 return ((wRes*wRes).sum())/(nData-nParsModel)
101 def makeMockFlats(expTime, gain=1.0, readNoiseElectrons=5, fluxElectrons=1000,
102 randomSeedFlat1=1984, randomSeedFlat2=666, powerLawBfParams=[],
104 """Create a pair or mock flats with isrMock.
109 Exposure time of the flats.
111 gain : `float`, optional
114 readNoiseElectrons : `float`, optional
115 Read noise rms, in electrons.
117 fluxElectrons : `float`, optional
118 Flux of flats, in electrons per second.
120 randomSeedFlat1 : `int`, optional
121 Random seed for the normal distrubutions for the mean signal and noise (flat1).
123 randomSeedFlat2 : `int`, optional
124 Random seed for the normal distrubutions for the mean signal and noise (flat2).
126 powerLawBfParams : `list`, optional
127 Parameters for `galsim.cdmodel.PowerLawCD` to simulate the brightter-fatter effect.
129 expId1 : `int`, optional
130 Exposure ID for first flat.
132 expId2 : `int`, optional
133 Exposure ID for second flat.
138 flatExp1 : `lsst.afw.image.exposure.exposure.ExposureF`
139 First exposure of flat field pair.
141 flatExp2 : `lsst.afw.image.exposure.exposure.ExposureF`
142 Second exposure of flat field pair.
146 The parameters of `galsim.cdmodel.PowerLawCD` are `n, r0, t0, rx, tx, r, t, alpha`. For more
147 information about their meaning, see the Galsim documentation
148 https://galsim-developers.github.io/GalSim/_build/html/_modules/galsim/cdmodel.html
149 and Gruen+15 (1501.02802).
151 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)
153 flatFlux = fluxElectrons
154 flatMean = flatFlux*expTime
155 readNoise = readNoiseElectrons
157 mockImageConfig = isrMock.IsrMock.ConfigClass()
159 mockImageConfig.flatDrop = 0.99999
160 mockImageConfig.isTrimmed =
True
162 flatExp1 = isrMock.FlatMock(config=mockImageConfig).
run()
163 flatExp2 = flatExp1.clone()
164 (shapeY, shapeX) = flatExp1.getDimensions()
165 flatWidth = np.sqrt(flatMean)
167 rng1 = np.random.RandomState(randomSeedFlat1)
168 flatData1 = rng1.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng1.normal(0.0, readNoise,
170 rng2 = np.random.RandomState(randomSeedFlat2)
171 flatData2 = rng2.normal(flatMean, flatWidth, (shapeX, shapeY)) + rng2.normal(0.0, readNoise,
174 if len(powerLawBfParams):
175 if not len(powerLawBfParams) == 8:
176 raise RuntimeError(
"Wrong number of parameters for `galsim.cdmodel.PowerLawCD`. "
177 f
"Expected 8; passed {len(powerLawBfParams)}.")
178 cd = galsim.cdmodel.PowerLawCD(*powerLawBfParams)
179 tempFlatData1 = galsim.Image(flatData1)
180 temp2FlatData1 = cd.applyForward(tempFlatData1)
182 tempFlatData2 = galsim.Image(flatData2)
183 temp2FlatData2 = cd.applyForward(tempFlatData2)
185 flatExp1.image.array[:] = temp2FlatData1.array/gain
186 flatExp2.image.array[:] = temp2FlatData2.array/gain
188 flatExp1.image.array[:] = flatData1/gain
189 flatExp2.image.array[:] = flatData2/gain
194 flatExp1.getInfo().setVisitInfo(visitInfoExp1)
195 flatExp2.getInfo().setVisitInfo(visitInfoExp2)
197 return flatExp1, flatExp2
201 """Count the number of pixels in a given mask plane."""
202 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
203 nPix = np.where(np.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
208 """Subclass of TaskRunner for handling intrinsically paired visits.
210 This transforms the processed arguments generated by the ArgumentParser
211 into the arguments expected by tasks which take visit pairs for their
214 Such tasks' run() methods tend to take two arguments,
215 one of which is the dataRef (as usual), and the other is the list
216 of visit-pairs, in the form of a list of tuples.
217 This list is supplied on the command line as documented,
218 and this class parses that, and passes the parsed version
221 See pipeBase.TaskRunner for more information.
226 """Parse the visit list and pass through explicitly."""
228 for visitStringPair
in parsedCmd.visitPairs:
229 visitStrings = visitStringPair.split(
",")
230 if len(visitStrings) != 2:
231 raise RuntimeError(
"Found {} visits in {} instead of 2".
format(len(visitStrings),
234 visits = [int(visit)
for visit
in visitStrings]
236 raise RuntimeError(
"Could not parse {} as two integer visit numbers".
format(visitStringPair))
237 visitPairs.append(visits)
239 return pipeBase.TaskRunner.getTargetList(parsedCmd, visitPairs=visitPairs, **kwargs)
243 """Parse command line numerical expression sytax and return as list of int
245 Take an input of the form "'1..5:2^123..126'" as a string, and return
246 a list of ints as [1, 3, 5, 123, 124, 125, 126]
249 for subString
in inputString.split(
"^"):
250 mat = re.search(
r"^(\d+)\.\.(\d+)(?::(\d+))?$", subString)
252 v1 = int(mat.group(1))
253 v2 = int(mat.group(2))
255 v3 = int(v3)
if v3
else 1
256 for v
in range(v1, v2 + 1, v3):
257 outList.append(int(v))
259 outList.append(int(subString))
264 """Subclass of TaskRunner for tasks requiring a list of visits per dataRef.
266 This transforms the processed arguments generated by the ArgumentParser
267 into the arguments expected by tasks which require a list of visits
268 to be supplied for each dataRef, as is common in `lsst.cp.pipe` code.
270 Such tasks' run() methods tend to take two arguments,
271 one of which is the dataRef (as usual), and the other is the list
273 This list is supplied on the command line as documented,
274 and this class parses that, and passes the parsed version
277 See `lsst.pipe.base.TaskRunner` for more information.
282 """Parse the visit list and pass through explicitly."""
285 assert len(parsedCmd.visitList) == 1,
'visitList parsing assumptions violated'
288 return pipeBase.TaskRunner.getTargetList(parsedCmd, visitList=visits, **kwargs)
292 """A DataIdContainer for the tasks for which the output does
296 """Compute refList based on idList.
298 This method must be defined as the dataset does not exist before this
304 Results of parsing the command-line.
308 Not called if ``add_id_argument`` called
309 with ``doMakeDataRefList=False``.
310 Note that this is almost a copy-and-paste of the vanilla
311 implementation, but without checking if the datasets already exist,
312 as this task exists to make them.
314 if self.datasetType
is None:
315 raise RuntimeError(
"Must call setDatasetType first")
316 butler = namespace.butler
317 for dataId
in self.idList:
318 refList =
list(butler.subset(datasetType=self.datasetType, level=self.level, dataId=dataId))
322 namespace.log.warn(
"No data found for dataId=%s", dataId)
324 self.refList += refList
327 def irlsFit(initialParams, dataX, dataY, function, weightsY=None, weightType='Cauchy'):
328 """Iteratively reweighted least squares fit.
330 This uses the `lsst.cp.pipe.utils.fitLeastSq`, but applies weights
331 based on the Cauchy distribution by default. Other weight options
332 are implemented. See e.g. Holland and Welsch, 1977,
333 doi:10.1080/03610927708827533
337 initialParams : `list` [`float`]
339 dataX : `numpy.array` [`float`]
341 dataY : `numpy.array` [`float`]
345 weightsY : `numpy.array` [`float`]
346 Weights to apply to the data.
347 weightType : `str`, optional
348 Type of weighting to use. One of Cauchy, Anderson, bisquare,
349 box, Welsch, Huber, logistic, or Fair.
353 polyFit : `list` [`float`]
354 Final best fit parameters.
355 polyFitErr : `list` [`float`]
356 Final errors on fit parameters.
359 weightsY : `list` [`float`]
360 Final weights used for each point.
365 Raised if an unknown weightType string is passed.
369 weightsY = np.ones_like(dataX)
371 polyFit, polyFitErr, chiSq =
fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
372 for iteration
in range(10):
373 resid = np.abs(dataY - function(polyFit, dataX)) / np.sqrt(dataY)
374 if weightType ==
'Cauchy':
378 weightsY = 1.0 / (1.0 + np.square(Z))
379 elif weightType ==
'Anderson':
382 Z = resid / (1.339 * np.pi)
383 weightsY = np.where(Z < 1.0, np.sinc(Z), 0.0)
384 elif weightType ==
'bisquare':
388 weightsY = np.where(Z < 1.0, 1.0 - np.square(Z), 0.0)
389 elif weightType ==
'box':
392 weightsY = np.where(resid < 2.795, 1.0, 0.0)
393 elif weightType ==
'Welsch':
397 weightsY = np.exp(-1.0 * np.square(Z))
398 elif weightType ==
'Huber':
402 weightsY = np.where(Z < 1.0, 1.0, 1 / Z)
403 elif weightType ==
'logistic':
407 weightsY = np.tanh(Z) / Z
408 elif weightType ==
'Fair':
412 weightsY = (1.0 / (1.0 + (Z)))
414 raise RuntimeError(f
"Unknown weighting type: {weightType}")
415 polyFit, polyFitErr, chiSq =
fitLeastSq(initialParams, dataX, dataY, function, weightsY=weightsY)
417 return polyFit, polyFitErr, chiSq, weightsY
420 def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None):
421 """Do a fit and estimate the parameter errors using using scipy.optimize.leastq.
423 optimize.leastsq returns the fractional covariance matrix. To estimate the
424 standard deviation of the fit parameters, multiply the entries of this matrix
425 by the unweighted reduced chi squared and take the square root of the diagonal elements.
429 initialParams : `list` of `float`
430 initial values for fit parameters. For ptcFitType=POLYNOMIAL, its length
431 determines the degree of the polynomial.
433 dataX : `numpy.array` of `float`
434 Data in the abscissa axis.
436 dataY : `numpy.array` of `float`
437 Data in the ordinate axis.
439 function : callable object (function)
440 Function to fit the data with.
442 weightsY : `numpy.array` of `float`
443 Weights of the data in the ordinate axis.
447 pFitSingleLeastSquares : `list` of `float`
448 List with fitted parameters.
450 pErrSingleLeastSquares : `list` of `float`
451 List with errors for fitted parameters.
453 reducedChiSqSingleLeastSquares : `float`
454 Reduced chi squared, unweighted if weightsY is not provided.
457 weightsY = np.ones(len(dataX))
459 def errFunc(p, x, y, weightsY=None):
461 weightsY = np.ones(len(x))
462 return (function(p, x) - y)*weightsY
464 pFit, pCov, infoDict, errMessage, success = leastsq(errFunc, initialParams,
465 args=(dataX, dataY, weightsY), full_output=1,
468 if (len(dataY) > len(initialParams))
and pCov
is not None:
473 pCov = np.zeros((len(initialParams), len(initialParams)))
475 reducedChiSq = np.nan
478 for i
in range(len(pFit)):
479 errorVec.append(np.fabs(pCov[i][i])**0.5)
481 pFitSingleLeastSquares = pFit
482 pErrSingleLeastSquares = np.array(errorVec)
484 return pFitSingleLeastSquares, pErrSingleLeastSquares, reducedChiSq
487 def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.):
488 """Do a fit using least squares and bootstrap to estimate parameter errors.
490 The bootstrap error bars are calculated by fitting 100 random data sets.
494 initialParams : `list` of `float`
495 initial values for fit parameters. For ptcFitType=POLYNOMIAL, its length
496 determines the degree of the polynomial.
498 dataX : `numpy.array` of `float`
499 Data in the abscissa axis.
501 dataY : `numpy.array` of `float`
502 Data in the ordinate axis.
504 function : callable object (function)
505 Function to fit the data with.
507 weightsY : `numpy.array` of `float`, optional.
508 Weights of the data in the ordinate axis.
510 confidenceSigma : `float`, optional.
511 Number of sigmas that determine confidence interval for the bootstrap errors.
515 pFitBootstrap : `list` of `float`
516 List with fitted parameters.
518 pErrBootstrap : `list` of `float`
519 List with errors for fitted parameters.
521 reducedChiSqBootstrap : `float`
522 Reduced chi squared, unweighted if weightsY is not provided.
525 weightsY = np.ones(len(dataX))
527 def errFunc(p, x, y, weightsY):
529 weightsY = np.ones(len(x))
530 return (function(p, x) - y)*weightsY
533 pFit, _ = leastsq(errFunc, initialParams, args=(dataX, dataY, weightsY), full_output=0)
536 residuals = errFunc(pFit, dataX, dataY, weightsY)
540 randomDelta = np.random.normal(0., np.fabs(residuals), len(dataY))
541 randomDataY = dataY + randomDelta
542 randomFit, _ = leastsq(errFunc, initialParams,
543 args=(dataX, randomDataY, weightsY), full_output=0)
544 pars.append(randomFit)
545 pars = np.array(pars)
546 meanPfit = np.mean(pars, 0)
549 errPfit = confidenceSigma*np.std(pars, 0)
550 pFitBootstrap = meanPfit
551 pErrBootstrap = errPfit
555 return pFitBootstrap, pErrBootstrap, reducedChiSq
559 """Polynomial function definition
563 Polynomial coefficients. Its length determines the polynomial order.
570 Ordinate array after evaluating polynomial of order len(pars)-1 at `x`.
572 return poly.polyval(x, [*pars])
576 """Single brighter-fatter parameter model for PTC; Equation 16 of Astier+19.
581 Parameters of the model: a00 (brightter-fatter), gain (e/ADU), and noise (e^2).
588 C_00 (variance) in ADU^2.
590 a00, gain, noise = pars
591 return 0.5/(a00*gain*gain)*(np.exp(2*a00*x*gain)-1) + noise/(gain*gain)
595 """Arrange exposures by exposure time.
599 exposureList : `list`[`lsst.afw.image.exposure.exposure.ExposureF`]
600 Input list of exposures.
602 exposureIdList : `list`[`int`]
603 List of exposure ids as obtained by dataId[`exposure`].
607 flatsAtExpTime : `dict` [`float`,
608 `list`[(`lsst.afw.image.exposure.exposure.ExposureF`, `int`)]]
609 Dictionary that groups flat-field exposures (and their IDs) that have
610 the same exposure time (seconds).
613 assert len(exposureList) == len(exposureIdList),
"Different lengths for exp. list and exp. ID lists"
614 for exp, expId
in zip(exposureList, exposureIdList):
615 expTime = exp.getInfo().getVisitInfo().getExposureTime()
616 listAtExpTime = flatsAtExpTime.setdefault(expTime, [])
617 listAtExpTime.append((exp, expId))
619 return flatsAtExpTime
623 """Arrange exposures by exposure ID.
625 There is no guarantee that this will properly group exposures, but
626 allows a sequence of flats that have different illumination
627 (despite having the same exposure time) to be processed.
631 exposureList : `list`[`lsst.afw.image.exposure.exposure.ExposureF`]
632 Input list of exposures.
634 exposureIdList : `list`[`int`]
635 List of exposure ids as obtained by dataId[`exposure`].
639 flatsAtExpId : `dict` [`float`,
640 `list`[(`lsst.afw.image.exposure.exposure.ExposureF`, `int`)]]
641 Dictionary that groups flat-field exposures (and their IDs)
642 sequentially by their exposure id.
647 This algorithm sorts the input exposures by their exposure id, and
648 then assigns each pair of exposures (exp_j, exp_{j+1}) to pair k,
649 such that 2*k = j, where j is the python index of one of the
650 exposures (starting from zero). By checking for the IndexError
651 while appending, we can ensure that there will only ever be fully
656 assert len(exposureList) == len(exposureIdList),
"Different lengths for exp. list and exp. ID lists"
658 sortedExposures = sorted(zip(exposureList, exposureIdList), key=
lambda pair: pair[1])
660 for jPair, expTuple
in enumerate(sortedExposures):
663 listAtExpId = flatsAtExpId.setdefault(kPair, [])
665 listAtExpId.append(expTuple)
666 listAtExpId.append(sortedExposures[jPair + 1])
674 """Check the exposure lengths of two exposures are equal.
678 exp1 : `lsst.afw.image.exposure.ExposureF`
679 First exposure to check
680 exp2 : `lsst.afw.image.exposure.ExposureF`
681 Second exposure to check
682 v1 : `int` or `str`, optional
683 First visit of the visit pair
684 v2 : `int` or `str`, optional
685 Second visit of the visit pair
686 raiseWithMessage : `bool`
687 If True, instead of returning a bool, raise a RuntimeError if exposure
688 times are not equal, with a message about which visits mismatch if the
689 information is available.
694 Raised if the exposure lengths of the two exposures are not equal
696 expTime1 = exp1.getInfo().getVisitInfo().getExposureTime()
697 expTime2 = exp2.getInfo().getVisitInfo().getExposureTime()
698 if expTime1 != expTime2:
700 msg =
"Exposure lengths for visit pairs must be equal. " + \
701 "Found %s and %s" % (expTime1, expTime2)
703 msg +=
" for visit pair %s, %s" % (v1, v2)
704 raise RuntimeError(msg)
710 def validateIsrConfig(isrTask, mandatory=None, forbidden=None, desirable=None, undesirable=None,
711 checkTrim=True, logName=None):
712 """Check that appropriate ISR settings have been selected for the task.
714 Note that this checks that the task itself is configured correctly rather
715 than checking a config.
719 isrTask : `lsst.ip.isr.IsrTask`
720 The task whose config is to be validated
722 mandatory : `iterable` of `str`
723 isr steps that must be set to True. Raises if False or missing
725 forbidden : `iterable` of `str`
726 isr steps that must be set to False. Raises if True, warns if missing
728 desirable : `iterable` of `str`
729 isr steps that should probably be set to True. Warns is False, info if
732 undesirable : `iterable` of `str`
733 isr steps that should probably be set to False. Warns is True, info if
737 Check to ensure the isrTask's assembly subtask is trimming the images.
738 This is a separate config as it is very ugly to do this within the
739 normal configuration lists as it is an option of a sub task.
744 Raised if ``mandatory`` config parameters are False,
745 or if ``forbidden`` parameters are True.
748 Raised if parameter ``isrTask`` is an invalid type.
752 Logs warnings using an isrValidation logger for desirable/undesirable
753 options that are of the wrong polarity or if keys are missing.
755 if not isinstance(isrTask, ipIsr.IsrTask):
756 raise TypeError(f
'Must supply an instance of lsst.ip.isr.IsrTask not {type(isrTask)}')
758 configDict = isrTask.config.toDict()
760 if logName
and isinstance(logName, str):
761 log = lsst.log.getLogger(logName)
763 log = lsst.log.getLogger(
"isrValidation")
766 for configParam
in mandatory:
767 if configParam
not in configDict:
768 raise RuntimeError(f
"Mandatory parameter {configParam} not found in the isr configuration.")
769 if configDict[configParam]
is False:
770 raise RuntimeError(f
"Must set config.isr.{configParam} to True for this task.")
773 for configParam
in forbidden:
774 if configParam
not in configDict:
775 log.warn(f
"Failed to find forbidden key {configParam} in the isr config. The keys in the"
776 " forbidden list should each have an associated Field in IsrConfig:"
777 " check that there is not a typo in this case.")
779 if configDict[configParam]
is True:
780 raise RuntimeError(f
"Must set config.isr.{configParam} to False for this task.")
783 for configParam
in desirable:
784 if configParam
not in configDict:
785 log.info(f
"Failed to find key {configParam} in the isr config. You probably want"
786 " to set the equivalent for your obs_package to True.")
788 if configDict[configParam]
is False:
789 log.warn(f
"Found config.isr.{configParam} set to False for this task."
790 " The cp_pipe Config recommends setting this to True.")
792 for configParam
in undesirable:
793 if configParam
not in configDict:
794 log.info(f
"Failed to find key {configParam} in the isr config. You probably want"
795 " to set the equivalent for your obs_package to False.")
797 if configDict[configParam]
is True:
798 log.warn(f
"Found config.isr.{configParam} set to True for this task."
799 " The cp_pipe Config recommends setting this to False.")
802 if not isrTask.assembleCcd.config.doTrim:
803 raise RuntimeError(
"Must trim when assembling CCDs. Set config.isr.assembleCcd.doTrim to True")
807 """Convert nested default dictionaries to regular dictionaries.
809 This is needed to prevent yaml persistence issues.
814 A possibly nested set of `defaultdict`.
819 A possibly nested set of `dict`.
821 for k, v
in d.items():
822 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 sigmaClipCorrection(nSigClip)
def arrangeFlatsByExpTime(exposureList, exposureIdList)
def fitBootstrap(initialParams, dataX, dataY, function, weightsY=None, confidenceSigma=1.)
def fitLeastSq(initialParams, dataX, dataY, function, weightsY=None)
def parseCmdlineNumberString(inputString)
def irlsFit(initialParams, dataX, dataY, function, weightsY=None, weightType='Cauchy')
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 arrangeFlatsByExpId(exposureList, exposureIdList)
def funcPolynomial(pars, x)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)