diff --git a/doc/versionHistory.rst b/doc/versionHistory.rst index 6174d181..2661371b 100644 --- a/doc/versionHistory.rst +++ b/doc/versionHistory.rst @@ -6,6 +6,33 @@ Version History ################## +.. _lsst.ts.wep-13.3.0: + +------------- +13.3.0 +------------- + +* Add donut quality tables to outputs even when there are no donuts that pass so that we can match it up to the donut stamps and understand rejections. +* Change default pipeline setting to false for rubinTV upload. + +.. _lsst.ts.wep-13.2.0: + +------------- +13.2.0 +------------- + +* Implemented joint-fitting of donut pairs with Danish. + +.. _lsst.ts.wep-13.1.0: + +------------- +13.1.0 +------------- + +* Set saveHistory=True and loosen convergence criteria in the Danish production pipeline +* Upgrades to the forward modeling util, including specifying flux ratios for blends, miscentering donuts, and simulating "flat" donuts without intensity patterns +* Fixed bug in forward modeling util when adding noise to large flux values + .. _lsst.ts.wep-13.0.4: ------------- @@ -23,7 +50,7 @@ Version History * Task plotPsfFromZern added in comCamRapidAnalysisPipeline and comCamRapidAnalysisDanishPipeline. -.. _lsst.ts.wep-13.0.3: +.. _lsst.ts.wep-13.0.2: ------------- 13.0.2 diff --git a/pipelines/_ingredients/donutVizGroupPipeline.yaml b/pipelines/_ingredients/donutVizGroupPipeline.yaml index 96d5c03d..36ea171b 100644 --- a/pipelines/_ingredients/donutVizGroupPipeline.yaml +++ b/pipelines/_ingredients/donutVizGroupPipeline.yaml @@ -14,17 +14,17 @@ tasks: plotAOSTask: class: lsst.donut.viz.PlotAOSTask config: - doRubinTVUpload: true + doRubinTVUpload: false aggregateDonutStampsTask: class: lsst.donut.viz.AggregateDonutStampsTask plotDonutTask: class: lsst.donut.viz.PlotDonutTask config: - doRubinTVUpload: true + doRubinTVUpload: false plotPsfZernTask: class: lsst.donut.viz.PlotPsfZernTask config: - doRubinTVUpload: true + doRubinTVUpload: false subsets: donutVizGroups: diff --git a/pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml b/pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml index c02dc4d8..b470bbf7 100644 --- a/pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml +++ b/pipelines/_ingredients/wepDirectDetectScienceGroupPipeline.yaml @@ -26,19 +26,6 @@ tasks: estimateZernikes.saveHistory: False estimateZernikes.maskKwargs: { "doMaskBlends": False } donutStampSelector.maxSelect: 5 - isr: - class: lsst.ip.isr.IsrTaskLSST - config: - # Although we don't have to apply the amp offset corrections, we do want - # to compute them for analyzeAmpOffsetMetadata to report on as metrics. - doAmpOffset: true - ampOffset.doApplyAmpOffset: false - # Turn off slow steps in ISR - doBrighterFatter: false - # Mask saturated pixels, - # but turn off quadratic crosstalk because it's currently broken - doSaturation: True - crosstalk.doQuadraticCrosstalkCorrection: False subsets: wepDirectDetect: diff --git a/pipelines/production/comCamDailyProcessing.yaml b/pipelines/production/comCamDailyProcessing.yaml index ee208b76..9edf3eae 100644 --- a/pipelines/production/comCamDailyProcessing.yaml +++ b/pipelines/production/comCamDailyProcessing.yaml @@ -9,22 +9,11 @@ tasks: class: lsst.ts.wep.task.generateDonutDirectDetectTask.GenerateDonutDirectDetectTask config: donutSelector.useCustomMagLimit: True - plotAOSTask: - class: lsst.donut.viz.PlotAOSTask - config: - doRubinTVUpload: false - aggregateDonutStampsTask: - class: lsst.donut.viz.AggregateDonutStampsTask - plotDonutTask: - class: lsst.donut.viz.PlotDonutTask - config: - doRubinTVUpload: false # Define pipeline steps subsets: step1: subset: - - isr - generateDonutDirectDetectTask - cutOutDonutsScienceSensorGroupTask - calcZernikesTask diff --git a/pipelines/production/comCamRapidAnalysisDanishPipeline.yaml b/pipelines/production/comCamRapidAnalysisDanishPipeline.yaml index 9763b4bf..676f8117 100644 --- a/pipelines/production/comCamRapidAnalysisDanishPipeline.yaml +++ b/pipelines/production/comCamRapidAnalysisDanishPipeline.yaml @@ -5,6 +5,19 @@ imports: - $TS_WEP_DIR/pipelines/_ingredients/donutVizGroupPipeline.yaml tasks: + isr: + class: lsst.ip.isr.IsrTaskLSST + config: + # Although we don't have to apply the amp offset corrections, we do want + # to compute them for analyzeAmpOffsetMetadata to report on as metrics. + doAmpOffset: true + ampOffset.doApplyAmpOffset: false + # Turn off slow steps in ISR + doBrighterFatter: false + # Mask saturated pixels, + # but turn off quadratic crosstalk because it's currently broken + doSaturation: True + crosstalk.doQuadraticCrosstalkCorrection: False calcZernikesTask: class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask config: @@ -16,6 +29,23 @@ tasks: estimateZernikes.binning: 2 estimateZernikes.nollIndices: [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 20, 21, 22, 27, 28] + estimateZernikes.saveHistory: true + estimateZernikes.lstsqKwargs: + ftol: 1.0e-3 + xtol: 1.0e-3 + gtol: 1.0e-3 + plotAOSTask: + class: lsst.donut.viz.PlotAOSTask + config: + doRubinTVUpload: true + plotDonutTask: + class: lsst.donut.viz.PlotDonutTask + config: + doRubinTVUpload: true + plotPsfZernTask: + class: lsst.donut.viz.PlotPsfZernTask + config: + doRubinTVUpload: true # Define pipeline steps subsets: diff --git a/pipelines/production/comCamRapidAnalysisPipeline.yaml b/pipelines/production/comCamRapidAnalysisPipeline.yaml index 5f86060a..a51a6169 100644 --- a/pipelines/production/comCamRapidAnalysisPipeline.yaml +++ b/pipelines/production/comCamRapidAnalysisPipeline.yaml @@ -5,6 +5,19 @@ imports: - $TS_WEP_DIR/pipelines/_ingredients/donutVizGroupPipeline.yaml tasks: + isr: + class: lsst.ip.isr.IsrTaskLSST + config: + # Although we don't have to apply the amp offset corrections, we do want + # to compute them for analyzeAmpOffsetMetadata to report on as metrics. + doAmpOffset: true + ampOffset.doApplyAmpOffset: false + # Turn off slow steps in ISR + doBrighterFatter: false + # Mask saturated pixels, + # but turn off quadratic crosstalk because it's currently broken + doSaturation: True + crosstalk.doQuadraticCrosstalkCorrection: False calcZernikesTask: class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask config: @@ -18,6 +31,18 @@ tasks: estimateZernikes.maskKwargs: { "doMaskBlends": False } donutStampSelector.maxSelect: 5 donutStampSelector.maxFracBadPixels: 2.0e-4 + plotAOSTask: + class: lsst.donut.viz.PlotAOSTask + config: + doRubinTVUpload: true + plotDonutTask: + class: lsst.donut.viz.PlotDonutTask + config: + doRubinTVUpload: true + plotPsfZernTask: + class: lsst.donut.viz.PlotPsfZernTask + config: + doRubinTVUpload: true # Define pipeline steps subsets: diff --git a/python/lsst/ts/wep/estimation/danish.py b/python/lsst/ts/wep/estimation/danish.py index 988a8514..ff651da7 100644 --- a/python/lsst/ts/wep/estimation/danish.py +++ b/python/lsst/ts/wep/estimation/danish.py @@ -50,11 +50,21 @@ class DanishAlgorithm(WfAlgorithm): binning : int, optional Binning factor to apply to the donut stamps before estimating Zernike coefficients. The default value of 1 means no binning. + jointFitPair : bool, optional + Whether to jointly fit intra/extra pairs, when a pair is provided. + If False, Zernikes are estimated for each individually, then + averaged. (the default is True) """ - def __init__(self, lstsqKwargs: Optional[dict] = None, binning: int = 1) -> None: + def __init__( + self, + lstsqKwargs: Optional[dict] = None, + binning: int = 1, + jointFitPair: bool = True, + ) -> None: self.binning = binning self.lstsqKwargs = lstsqKwargs + self.jointFitPair = jointFitPair @property def requiresPairs(self) -> bool: @@ -112,6 +122,26 @@ def lstsqKwargs(self, value: Union[dict, None]) -> None: self._lstsqKwargs = value + @property + def jointFitPair(self) -> bool: + """Whether to jointly fit intra/extra pairs.""" + return self._jointFitPair + + @jointFitPair.setter + def jointFitPair(self, value: bool) -> None: + """Whether to jointly fit intra/extra pairs, when a pair is provided. + + Parameters + ---------- + value : bool + Whether to jointly fit intra/extra pairs, when a pair is provided. + If False, Zernikes are estimated for each individually, then + averaged. + """ + if not isinstance(value, bool): + raise TypeError("jointFitPair must be a bool.") + self._jointFitPair = value + @property def history(self) -> dict: """The algorithm history. @@ -137,39 +167,13 @@ def history(self) -> dict: """ return super().history - def _estimateSingleZk( + def _prepDanish( self, image: Image, zkStart: np.ndarray, nollIndices: np.ndarray, instrument: Instrument, - factory: danish.DonutFactory, - saveHistory: bool, - ) -> Tuple[np.ndarray, dict]: - """Estimate Zernikes (in meters) for a single donut stamp. - - Parameters - ---------- - image : Image - The ts_wep image of the donut stamp - zkStart : np.ndarray - The starting point for the Zernikes - nollIndices : np.ndarray - Noll indices for which to estimate Zernikes - instrument : Instrument - The ts_wep Instrument - factory : danish.DonutFactory - The Danish donut factory - saveHistory : bool - Whether to create a history to be saved - - Returns - ------- - np.ndarray - The Zernike coefficients (in meters) for Noll indices >= 4 - dict - The single-stamp history. This is empty if saveHistory is False. - """ + ): # Warn about using Danish for blended donuts if np.isfinite(image.blendOffsets).sum() > 0: warnings.warn("Danish is currently only setup for non-blended donuts.") @@ -207,13 +211,58 @@ def _estimateSingleZk( if img.shape[0] % 2 == 0: img = img[:-1, :-1] + # Convert field angle from degrees to radians + angle = np.deg2rad(image.fieldAngle) + + return img, angle, zkRef, backgroundStd + + def _estimateSingleZk( + self, + image: Image, + zkStart: np.ndarray, + nollIndices: np.ndarray, + instrument: Instrument, + factory: danish.DonutFactory, + saveHistory: bool, + ) -> Tuple[np.ndarray, dict]: + """Estimate Zernikes (in meters) for a single donut stamp. + + Parameters + ---------- + image : Image + The ts_wep image of the donut stamp + zkStart : np.ndarray + The starting point for the Zernikes + nollIndices : np.ndarray + Noll indices for which to estimate Zernikes + instrument : Instrument + The ts_wep Instrument + factory : danish.DonutFactory + The Danish donut factory + saveHistory : bool + Whether to create a history to be saved + + Returns + ------- + np.ndarray + The Zernike coefficients (in meters) for Noll indices >= 4 + dict + The single-stamp history. This is empty if saveHistory is False. + """ + img, angle, zkRef, backgroundStd = self._prepDanish( + image=image, + zkStart=zkStart, + nollIndices=nollIndices, + instrument=instrument, + ) + # Create the Danish donut model model = danish.SingleDonutModel( factory, z_ref=zkRef, z_terms=nollIndices, - thx=np.deg2rad(image.fieldAngle[0]), - thy=np.deg2rad(image.fieldAngle[1]), + thx=angle[0], + thy=angle[1], npix=img.shape[0], ) @@ -248,7 +297,7 @@ def _estimateSingleZk( dy, fwhm, zkFit, - sky_level=backgroundStd, + sky_level=backgroundStd**2, flux=img.sum(), ) @@ -283,6 +332,172 @@ def _estimateSingleZk( return zkSum, hist + def _estimatePairZk( + self, + I1: Image, + I2: Optional[Image], + zkStartI1: np.ndarray, + zkStartI2: Optional[np.ndarray], + nollIndices: np.ndarray, + instrument: Instrument, + factory: danish.DonutFactory, + saveHistory: bool, + ) -> Tuple[np.ndarray, dict]: + """Estimate Zernikes (in meters) for pairs of donut stamps. + + Parameters + ---------- + I1 : Image + An Image object containing an intra- or extra-focal donut image. + I2 : Image or None + A second image, on the opposite side of focus from I1. Can be None. + zkStartI1 : np.ndarray + The starting Zernikes for I1 (in meters, for Noll indices >= 4) + zkStartI2 : np.ndarray or None + The starting Zernikes for I2 (in meters, for Noll indices >= 4) + nollIndices : np.ndarray + Noll indices for which you wish to estimate Zernike coefficients. + instrument : Instrument + The Instrument object associated with the DonutStamps. + factory : danish.DonutFactory + The Danish donut factory + saveHistory : bool + Whether to save the algorithm history in the self.history + attribute. If True, then self.history contains information + about the most recent time the algorithm was run. + + Returns + ------- + Returns + ------- + np.ndarray + The Zernike coefficients (in meters) for Noll indices >= 4 + dict + The single-stamp history. This is empty if saveHistory is False. + """ + # Prep quantities for both images + img1, angle1, zkRef1, backgroundStd1 = self._prepDanish( + image=I1, + zkStart=zkStartI1, + nollIndices=nollIndices, + instrument=instrument, + ) + img2, angle2, zkRef2, backgroundStd2 = self._prepDanish( + image=I2, + zkStart=zkStartI2, + nollIndices=nollIndices, + instrument=instrument, + ) + + # Package these into lists for Danish + imgs = [img1, img2] + thxs = [angle1[0], angle2[0]] + thys = [angle1[1], angle2[1]] + zkRefs = [zkRef1, zkRef2] + skyLevels = [backgroundStd1**2, backgroundStd2**2] + + # Create Double Zernike tuples + dzTerms = [(1, j) for j in nollIndices] + + # Set field radius to max value from mask params + fieldRadius = np.deg2rad( + np.max( + [ + edge["thetaMax"] + for item in instrument.maskParams.values() + for edge in item.values() + ] + ) + ) + + # Create model + model = danish.MultiDonutModel( + factory, + z_refs=zkRefs, + dz_terms=dzTerms, + field_radius=fieldRadius, + thxs=thxs, + thys=thys, + npix=imgs[0].shape[0], + ) + + # Initial guess + x0 = [0.0] * 2 + [0.0] * 2 + [0.7] + [0.0] * len(dzTerms) + + # Use scipy to optimize the parameters + try: + result = least_squares( + model.chi, + jac=model.jac, + x0=x0, + args=(imgs, skyLevels), + **self.lstsqKwargs, + ) + result = dict(result) + + # Unpack the parameters + dxs, dys, fwhm, zkFit = model.unpack_params(result["x"]) + + # Add the starting zernikes back into the result + zkSum = zkFit + np.nanmean([zkStartI1, zkStartI2], axis=0) + + # Flag that we didn't hit GalSimFFTSizeError + galSimFFTSizeError = False + + # If we're saving the history, compute the model image + if saveHistory: + modelImages = model.model( + dxs, + dys, + fwhm, + zkFit, + sky_levels=skyLevels, + fluxes=np.sum(imgs, axis=(1, 2)), + ) + + # Sometimes this happens with Danish :( + except GalSimFFTSizeError: + # Fill dummy objects + result = None + zkFit = np.full_like(zkStartI1, np.nan) + zkSum = np.full_like(zkStartI1, np.nan) + if saveHistory: + modelImages = np.full_like(imgs, np.nan) + + # Flag the error + galSimFFTSizeError = True + + if saveHistory: + # Save the history + hist = {} + hist[I1.defocalType.value] = { + "image": imgs[0].copy(), + "variance": skyLevels[0], + "nollIndices": nollIndices.copy(), + "zkStart": zkStartI1.copy(), + "lstsqResult": result, + "zkFit": zkFit.copy(), + "zkSum": zkSum.copy(), + "model": modelImages[0], + "GalSimFFTSizeError": galSimFFTSizeError, + } + hist[I2.defocalType.value] = { + "image": imgs[1].copy(), + "variance": skyLevels[1], + "nollIndices": nollIndices.copy(), + "zkStart": zkStartI2.copy(), + "lstsqResult": result, + "zkFit": zkFit.copy(), + "zkSum": zkSum.copy(), + "model": modelImages[1], + "GalSimFFTSizeError": galSimFFTSizeError, + } + + else: + hist = {} + + return zkSum, hist + def _estimateZk( self, I1: Image, @@ -317,7 +532,7 @@ def _estimateZk( Returns ------- np.ndarray - Zernike coefficients (for Noll indices >= 4) estimated from + Zernike coefficients for the provided Noll indices, estimated from the images, in meters. """ # Create the Danish donut factory @@ -329,23 +544,43 @@ def _estimateZk( pixel_scale=instrument.pixelSize * self.binning, ) - # Create an empty history - hist = {} + if I2 is None or not self.jointFitPair: + # Create an empty history + hist = {} - # Estimate for I1 - zk1, hist[I1.defocalType.value] = self._estimateSingleZk( - I1, - zkStartI1, - nollIndices, - instrument, - factory, - saveHistory, - ) + # Estimate for I1 + zk1, hist[I1.defocalType.value] = self._estimateSingleZk( + I1, + zkStartI1, + nollIndices, + instrument, + factory, + saveHistory, + ) + + if I2 is not None: + # If I2 provided, estimate for that donut as well + zk2, hist[I2.defocalType.value] = self._estimateSingleZk( + I2, + zkStartI2, + nollIndices, + instrument, + factory, + saveHistory, + ) + + # Average the Zernikes + zk = np.mean([zk1, zk2], axis=0) + else: + zk = zk1 - if I2 is not None: - # If I2 provided, estimate for that donut as well - zk2, hist[I2.defocalType.value] = self._estimateSingleZk( + hist["zk"] = zk + + else: + zk, hist = self._estimatePairZk( + I1, I2, + zkStartI1, zkStartI2, nollIndices, instrument, @@ -353,13 +588,7 @@ def _estimateZk( saveHistory, ) - # Average the Zernikes - zk = np.mean([zk1, zk2], axis=0) - else: - zk = zk1 - if saveHistory: - hist["zk"] = zk self._history = hist return zk diff --git a/python/lsst/ts/wep/task/calcZernikesTask.py b/python/lsst/ts/wep/task/calcZernikesTask.py index b0990fc4..6c3fc0b3 100644 --- a/python/lsst/ts/wep/task/calcZernikesTask.py +++ b/python/lsst/ts/wep/task/calcZernikesTask.py @@ -346,8 +346,23 @@ def createZkTable( return zkTable - def empty(self) -> pipeBase.Struct: - """Return empty results if no donuts are available.""" + def empty(self, qualityTable=None) -> pipeBase.Struct: + """Return empty results if no donuts are available. If + it is a result of no quality donuts we still include the + quality table results instead of an empty quality table. + + Parameters + ---------- + qualityTable : astropy.table.QTable + Quality table created with donut stamp input. + + Returns + ------- + lsst.pipe.base.Struct + Empty output tables for zernikes. Empty quality table + if no donuts. Otherwise contains quality table + with donuts that all failed to pass quality check. + """ qualityTableCols = [ "SN", "ENTROPY", @@ -356,11 +371,15 @@ def empty(self) -> pipeBase.Struct: "FINAL_SELECT", "DEFOCAL_TYPE", ] + if qualityTable is None: + donutQualityTable = QTable({name: [] for name in qualityTableCols}) + else: + donutQualityTable = qualityTable return pipeBase.Struct( outputZernikesRaw=np.atleast_2d(np.full(len(self.nollIndices), np.nan)), outputZernikesAvg=np.atleast_2d(np.full(len(self.nollIndices), np.nan)), zernikes=self.initZkTable(), - donutQualityTable=QTable({name: [] for name in qualityTableCols}), + donutQualityTable=donutQualityTable, ) @timeMethod @@ -396,7 +415,7 @@ def run( or len(selectionIntra.donutStampsSelect) == 0 ): self.log.info("No donut stamps were selected.") - return self.empty() + return self.empty(qualityTable=donutQualityTable) else: donutQualityTable = QTable([]) diff --git a/python/lsst/ts/wep/task/estimateZernikesDanishTask.py b/python/lsst/ts/wep/task/estimateZernikesDanishTask.py index e68ff073..de2df13c 100644 --- a/python/lsst/ts/wep/task/estimateZernikesDanishTask.py +++ b/python/lsst/ts/wep/task/estimateZernikesDanishTask.py @@ -44,6 +44,12 @@ class EstimateZernikesDanishConfig(EstimateZernikesBaseConfig): doc="Binning factor to apply to the donut stamps before estimating " + "Zernike coefficients. A value of 1 means no binning.", ) + jointFitPair = pexConfig.Field( + dtype=bool, + default=True, + doc="Whether to jointly fit intra/extra pairs, when a pair is provided. " + + "If False, Zernikes are estimated for each individually, then averaged. ", + ) class EstimateZernikesDanishTask(EstimateZernikesBaseTask): diff --git a/python/lsst/ts/wep/utils/modelUtils.py b/python/lsst/ts/wep/utils/modelUtils.py index 19a6f154..01fed252 100644 --- a/python/lsst/ts/wep/utils/modelUtils.py +++ b/python/lsst/ts/wep/utils/modelUtils.py @@ -45,10 +45,16 @@ def forwardModelPair( bandLabel: str = "r", fieldAngleIntra: Union[tuple, None] = None, fieldAngleExtra: Union[tuple, None] = None, + miscenterIntra: Union[tuple, None] = None, + miscenterExtra: Union[tuple, None] = None, blendOffsetsIntra: Union[np.ndarray, tuple, list, None] = None, blendOffsetsExtra: Union[np.ndarray, tuple, list, None] = None, + blendRatiosIntra: Union[np.ndarray, tuple, list, None] = None, + blendRatiosExtra: Union[np.ndarray, tuple, list, None] = None, instConfig: Union[str, dict, Instrument] = "policy:instruments/LsstCam.yaml", nPix: int = 180, + opticalModel: str = "offAxis", + flat: bool = False, ) -> Tuple[np.ndarray, Image, Image]: """Forward model a pair of donuts. @@ -97,17 +103,41 @@ def forwardModelPair( only specified for the intra or extrafocal donut, both donuts use that angle. If neither is specified, the same random angle is used for both. (the default is None) + miscenterIntra : tuple, optional + The amount by which the intrafocal donut is miscentered. A tuple of + (dx, dy) in pixels. If None, a random value between +/- 2 is used + for each. Note the random non-zero offsets are default to avoid + pixel aliasing effects when examining ensembles of wavefront + estimation errors. (the default is None) + miscenterExtra : tuple, optional + The amount by which the extrafocal donut is miscentered. A tuple of + (dx, dy) in pixels. If None, a random value between +/- 2 is used + for each. Note the random non-zero offsets are default to avoid + pixel aliasing effects when examining ensembles of wavefront + estimation errors. (the default is None) blendOffsetsIntra : Iterable or None, optional The blend offsets of the intrafocal donut. (the default is None) blendOffsetsExtra : Iterable or None, optional The blend offsets of the extrafocal donut. (the default is None) + blendRatiosIntra : Iterable or None, optional + Flux ratios of the blends to the central star. If None, 1 is assumed + for all. (the default is None) + blendRatiosExtra : Iterable or None, optional + Flux ratios of the blends to the central star. If None, 1 is assumed + for all. (the default is None) instConfig : str, dict, or Instrument, optional The instrument config for the image mapper. (the default is "policy:instruments/LsstCam.yaml") nPix : int, optional The size of the images. (the default is 180) + opticalModel : str, optional + Which optical model to use for the ImageMapper. Can be "offAxis", + "onAxis", or "paraxial". (the default is "offAxis") + flat : bool, optional + Whether to remove surface brightness fluctuations from the donut. + (the default is False) Returns ------- @@ -119,7 +149,7 @@ def forwardModelPair( The extrafocal image """ # Create the ImageMapper that will forward model the images - mapper = ImageMapper(instConfig=instConfig) + mapper = ImageMapper(instConfig=instConfig, opticalModel=opticalModel) # And the random number generators rng = np.random.default_rng(seed) @@ -181,22 +211,34 @@ def forwardModelPair( ) # Add blends if blendOffsetsIntra is None: - nBlends = 0 + blendRatiosIntra = [0] else: offsets = np.atleast_2d(blendOffsetsIntra) - nBlends = len(offsets) + blendRatiosIntra = ( + np.ones(offsets.shape[0]) if blendRatiosIntra is None else blendRatiosIntra + ) centralDonut = intraStamp.image.copy() - for blendShift in offsets: - intraStamp.image += shift(centralDonut, blendShift[::-1]) + for blendShift, ratio in zip(offsets, blendRatiosIntra): + intraStamp.image += ratio * shift(centralDonut, blendShift[::-1]) + # Flatten surface brightness? + if flat: + mapper.createImageMasks(intraStamp, zkCoeff, isBinary=False) + intraStamp.image = np.clip(intraStamp.mask + intraStamp.maskBlends, 0, 1) + # Normalize the flux + intraStamp.image *= fluxIntra * (1 + sum(blendRatiosIntra)) / intraStamp.image.sum() + # Miscenter image + _miscenterIntra = rng.uniform(-2, 2, size=2) + miscenterIntra = _miscenterIntra if miscenterIntra is None else miscenterIntra + intraStamp.image = shift(intraStamp.image, miscenterIntra[::-1]) # Convolve with Kolmogorov atmosphere if seeing > 0: intraStamp.image = convolve(intraStamp.image, atmKernel, mode="same") - # Normalize the flux - intraStamp.image *= fluxIntra * (1 + nBlends) / intraStamp.image.sum() - # Poisson noise + # Pseudo-Poissonian noise intraStamp.image += background intraStamp.image = np.clip(intraStamp.image, 0, None) - intraStamp.image = rng.poisson(intraStamp.image).astype(float) + intraStamp.image += rng.normal(size=intraStamp.image.shape) * np.sqrt( + intraStamp.image + ) intraStamp.image -= background # Now the extrafocal image @@ -212,22 +254,34 @@ def forwardModelPair( ) # Add blends if blendOffsetsExtra is None: - nBlends = 0 + blendRatiosExtra = [0] else: offsets = np.atleast_2d(blendOffsetsExtra) - nBlends = len(offsets) + blendRatiosExtra = ( + np.ones(offsets.shape[0]) if blendRatiosExtra is None else blendRatiosExtra + ) centralDonut = extraStamp.image.copy() - for blendShift in offsets: - extraStamp.image += shift(centralDonut, blendShift[::-1]) + for blendShift, ratio in zip(offsets, blendRatiosExtra): + extraStamp.image += ratio * shift(centralDonut, blendShift[::-1]) + # Flatten surface brightness? + if flat: + mapper.createImageMasks(extraStamp, zkCoeff, isBinary=False) + extraStamp.image = np.clip(extraStamp.mask + extraStamp.maskBlends, 0, 1) + # Normalize the flux + extraStamp.image *= fluxExtra * (1 + sum(blendRatiosExtra)) / extraStamp.image.sum() + # Miscenter image + _miscenterExtra = rng.uniform(-2, 2, size=2) + miscenterExtra = _miscenterExtra if miscenterExtra is None else miscenterExtra + extraStamp.image = shift(extraStamp.image, miscenterExtra[::-1]) # Convolve with Kolmogorov atmosphere if seeing > 0: extraStamp.image = convolve(extraStamp.image, atmKernel, mode="same") - # Normalize the flux - extraStamp.image *= fluxExtra / extraStamp.image.sum() - # Poisson noise + # Pseudo-Poissonian noise extraStamp.image += background extraStamp.image = np.clip(extraStamp.image, 0, None) - extraStamp.image = rng.poisson(extraStamp.image).astype(float) + extraStamp.image += rng.normal(size=extraStamp.image.shape) * np.sqrt( + extraStamp.image + ) extraStamp.image -= background return zkCoeff, intraStamp, extraStamp diff --git a/tests/estimation/test_danish.py b/tests/estimation/test_danish.py index 6d264393..344474af 100644 --- a/tests/estimation/test_danish.py +++ b/tests/estimation/test_danish.py @@ -50,45 +50,46 @@ def testGoodLstsqKwargs(self): self.assertEqual(hist["lstsqResult"]["nfev"], 1) def testAccuracy(self): - # Create estimator - dan = DanishAlgorithm() - danBin = DanishAlgorithm(binning=2) - - # Try several different random seeds - for seed in [12345, 23451, 34512, 45123, 51234]: - # Get the test data - zkTrue, intra, extra = forwardModelPair(seed=seed) - - # Compute shape of binned images - shapex, shapey = intra.image.shape - binned_shapex = shapex // 2 - binned_shapey = shapey // 2 - - # Ensure odd - if binned_shapex % 2 == 0: - binned_shapex -= 1 - if binned_shapey % 2 == 0: - binned_shapey -= 1 - binned_shape = (binned_shapex, binned_shapey) - - # Test estimation with pairs and single donuts: - for images in [[intra, extra], [intra], [extra]]: - # Estimate Zernikes (in meters) - zkEst = dan.estimateZk(*images) - - # Check that results are fairly accurate - self.assertLess(np.sqrt(np.sum((zkEst - zkTrue) ** 2)), 0.35e-6) - - # Test with binning - zkEst = danBin.estimateZk(*images, saveHistory=True) - self.assertLess(np.sqrt(np.sum((zkEst - zkTrue) ** 2)), 0.35e-6) - - # Test that we binned the images. - if "intra" in danBin.history: - self.assertEqual( - danBin.history["intra"]["image"].shape, binned_shape - ) - if "extra" in danBin.history: - self.assertEqual( - danBin.history["extra"]["image"].shape, binned_shape - ) + for jointFitPair in [True, False]: + # Create estimator + dan = DanishAlgorithm(jointFitPair=jointFitPair) + danBin = DanishAlgorithm(binning=2, jointFitPair=jointFitPair) + + # Try several different random seeds + for seed in [12345, 23451, 34512, 45123]: + # Get the test data + zkTrue, intra, extra = forwardModelPair(seed=seed) + + # Compute shape of binned images + shapex, shapey = intra.image.shape + binned_shapex = shapex // 2 + binned_shapey = shapey // 2 + + # Ensure odd + if binned_shapex % 2 == 0: + binned_shapex -= 1 + if binned_shapey % 2 == 0: + binned_shapey -= 1 + binned_shape = (binned_shapex, binned_shapey) + + # Test estimation with pairs and single donuts: + for images in [[intra, extra], [intra], [extra]]: + # Estimate Zernikes (in meters) + zkEst = dan.estimateZk(*images) + + # Check that results are fairly accurate + self.assertLess(np.sqrt(np.sum((zkEst - zkTrue) ** 2)), 0.35e-6) + + # Test with binning + zkEst = danBin.estimateZk(*images, saveHistory=True) + self.assertLess(np.sqrt(np.sum((zkEst - zkTrue) ** 2)), 0.35e-6) + + # Test that we binned the images. + if "intra" in danBin.history: + self.assertEqual( + danBin.history["intra"]["image"].shape, binned_shape + ) + if "extra" in danBin.history: + self.assertEqual( + danBin.history["extra"]["image"].shape, binned_shape + ) diff --git a/tests/task/test_calcZernikesTieTaskScienceSensor.py b/tests/task/test_calcZernikesTieTaskScienceSensor.py index fa6575dc..3c1d4c1c 100644 --- a/tests/task/test_calcZernikesTieTaskScienceSensor.py +++ b/tests/task/test_calcZernikesTieTaskScienceSensor.py @@ -225,6 +225,13 @@ def testCalcZernikes(self): self.assertIsInstance(struct.outputZernikesRaw, np.ndarray) self.assertIsInstance(struct.outputZernikesAvg, np.ndarray) self.assertIsInstance(struct.zernikes, QTable) + self.assertEqual(len(structNormal.donutQualityTable), 6) + self.assertEqual(len(structNull.donutQualityTable), 0) + + self.config.donutStampSelector.maxSelect = 0 + self.task = CalcZernikesTask(config=self.config) + structAllDonutsFail = self.task.run(donutStampsIntra, donutStampsExtra) + self.assertEqual(len(structAllDonutsFail.donutQualityTable), 6) def testGetCombinedZernikes(self): testArr = np.zeros((2, 19)) diff --git a/tests/utils/test_modelUtils.py b/tests/utils/test_modelUtils.py index 66b4754d..3487099c 100644 --- a/tests/utils/test_modelUtils.py +++ b/tests/utils/test_modelUtils.py @@ -82,8 +82,9 @@ def testZernikes(self): self.assertTrue(np.allclose(0, zk2)) # Divergent Zernikes - with self.assertRaises(ValueError): - forwardModelPair(zkNorm=1e9, zkMax=1e9) + _, intra, extra = forwardModelPair(zkNorm=1e9, zkMax=1e9) + self.assertFalse(np.isfinite(intra.image).any()) + self.assertFalse(np.isfinite(extra.image).any()) def testBlends(self): # Just going to test that blends get added by looking at flux @@ -96,6 +97,17 @@ def testBlends(self): self.assertTrue(intra1.image.sum() < intra2.image.sum()) self.assertTrue(extra1.image.sum() < extra2.image.sum()) + # Ensure that flux ratios work + _, intra3, extra3 = forwardModelPair( + blendOffsetsIntra=((40, 30),), + blendOffsetsExtra=((40, 30),), + blendRatiosIntra=[0.1], + blendRatiosExtra=[2.0], + ) + self.assertTrue(intra1.image.sum() < intra3.image.sum()) + self.assertTrue(intra3.image.sum() < intra2.image.sum()) + self.assertTrue(extra2.image.sum() < extra3.image.sum()) + # Make sure that blendOffsets is robust to 1D specification _, intra1, extra1 = forwardModelPair(blendOffsetsIntra=(40, 30)) _, intra2, extra2 = forwardModelPair(blendOffsetsIntra=((40, 30),)) @@ -107,3 +119,9 @@ def testSeeing(self): _, intra2, extra2 = forwardModelPair(seeing=1, skyLevel=0) self.assertTrue(np.sum(intra1.image > 0) < np.sum(intra2.image > 0)) self.assertTrue(np.sum(extra1.image > 0) < np.sum(extra2.image > 0)) + + def testMiscenter(self): + # Move extrafocal donut all the way out of frame + _, intra, extra = forwardModelPair(skyLevel=0, miscenterExtra=(500, 0)) + self.assertTrue(intra.image.sum() > 0) + self.assertTrue(np.isclose(extra.image.sum(), 0))