23 """Perform a single fit cycle of FGCM.
25 This task runs a single "fit cycle" of fgcm. Prior to running this task
26 one must run both fgcmMakeLut (to construct the atmosphere and instrumental
27 look-up-table) and fgcmBuildStars (to extract visits and star observations
30 The fgcmFitCycle is meant to be run multiple times, and is tracked by the
31 'cycleNumber'. After each run of the fit cycle, diagnostic plots should
32 be inspected to set parameters for outlier rejection on the following
33 cycle. Please see the fgcmcal Cookbook for details.
46 from lsst.utils.timer
import timeMethod
48 from .utilities
import makeConfigDict, translateFgcmLut, translateVisitCatalog
49 from .utilities
import extractReferenceMags
50 from .utilities
import makeZptSchema, makeZptCat
51 from .utilities
import makeAtmSchema, makeAtmCat, makeStdSchema, makeStdCat
52 from .sedterms
import SedboundarytermDict, SedtermDict
53 from .utilities
import lookupStaticCalibrations
54 from .focalPlaneProjector
import FocalPlaneProjector
58 __all__ = [
'FgcmFitCycleConfig',
'FgcmFitCycleTask',
'FgcmFitCycleRunner']
60 MULTIPLE_CYCLES_MAX = 10
64 dimensions=(
"instrument",),
65 defaultTemplates={
"previousCycleNumber":
"-1",
67 camera = connectionTypes.PrerequisiteInput(
68 doc=
"Camera instrument",
70 storageClass=
"Camera",
71 dimensions=(
"instrument",),
72 lookupFunction=lookupStaticCalibrations,
76 fgcmLookUpTable = connectionTypes.PrerequisiteInput(
77 doc=(
"Atmosphere + instrument look-up-table for FGCM throughput and "
78 "chromatic corrections."),
79 name=
"fgcmLookUpTable",
80 storageClass=
"Catalog",
81 dimensions=(
"instrument",),
85 fgcmVisitCatalog = connectionTypes.Input(
86 doc=
"Catalog of visit information for fgcm",
87 name=
"fgcmVisitCatalog",
88 storageClass=
"Catalog",
89 dimensions=(
"instrument",),
93 fgcmStarObservations = connectionTypes.Input(
94 doc=
"Catalog of star observations for fgcm",
95 name=
"fgcmStarObservations",
96 storageClass=
"Catalog",
97 dimensions=(
"instrument",),
101 fgcmStarIds = connectionTypes.Input(
102 doc=
"Catalog of fgcm calibration star IDs",
104 storageClass=
"Catalog",
105 dimensions=(
"instrument",),
109 fgcmStarIndices = connectionTypes.Input(
110 doc=
"Catalog of fgcm calibration star indices",
111 name=
"fgcmStarIndices",
112 storageClass=
"Catalog",
113 dimensions=(
"instrument",),
117 fgcmReferenceStars = connectionTypes.Input(
118 doc=
"Catalog of fgcm-matched reference stars",
119 name=
"fgcmReferenceStars",
120 storageClass=
"Catalog",
121 dimensions=(
"instrument",),
125 fgcmFlaggedStarsInput = connectionTypes.PrerequisiteInput(
126 doc=
"Catalog of flagged stars for fgcm calibration from previous fit cycle",
127 name=
"fgcmFlaggedStars{previousCycleNumber}",
128 storageClass=
"Catalog",
129 dimensions=(
"instrument",),
133 fgcmFitParametersInput = connectionTypes.PrerequisiteInput(
134 doc=
"Catalog of fgcm fit parameters from previous fit cycle",
135 name=
"fgcmFitParameters{previousCycleNumber}",
136 storageClass=
"Catalog",
137 dimensions=(
"instrument",),
141 fgcmFitParameters = connectionTypes.Output(
142 doc=
"Catalog of fgcm fit parameters from current fit cycle",
143 name=
"fgcmFitParameters{cycleNumber}",
144 storageClass=
"Catalog",
145 dimensions=(
"instrument",),
148 fgcmFlaggedStars = connectionTypes.Output(
149 doc=
"Catalog of flagged stars for fgcm calibration from current fit cycle",
150 name=
"fgcmFlaggedStars{cycleNumber}",
151 storageClass=
"Catalog",
152 dimensions=(
"instrument",),
155 fgcmZeropoints = connectionTypes.Output(
156 doc=
"Catalog of fgcm zeropoint data from current fit cycle",
157 name=
"fgcmZeropoints{cycleNumber}",
158 storageClass=
"Catalog",
159 dimensions=(
"instrument",),
162 fgcmAtmosphereParameters = connectionTypes.Output(
163 doc=
"Catalog of atmospheric fit parameters from current fit cycle",
164 name=
"fgcmAtmosphereParameters{cycleNumber}",
165 storageClass=
"Catalog",
166 dimensions=(
"instrument",),
169 fgcmStandardStars = connectionTypes.Output(
170 doc=
"Catalog of standard star magnitudes from current fit cycle",
171 name=
"fgcmStandardStars{cycleNumber}",
172 storageClass=
"SimpleCatalog",
173 dimensions=(
"instrument",),
179 for cycle
in range(MULTIPLE_CYCLES_MAX):
180 vars()[f
"fgcmFitParameters{cycle}"] = connectionTypes.Output(
181 doc=f
"Catalog of fgcm fit parameters from fit cycle {cycle}",
182 name=f
"fgcmFitParameters{cycle}",
183 storageClass=
"Catalog",
184 dimensions=(
"instrument",),
186 vars()[f
"fgcmFlaggedStars{cycle}"] = connectionTypes.Output(
187 doc=f
"Catalog of flagged stars for fgcm calibration from fit cycle {cycle}",
188 name=f
"fgcmFlaggedStars{cycle}",
189 storageClass=
"Catalog",
190 dimensions=(
"instrument",),
192 vars()[f
"fgcmZeropoints{cycle}"] = connectionTypes.Output(
193 doc=f
"Catalog of fgcm zeropoint data from fit cycle {cycle}",
194 name=f
"fgcmZeropoints{cycle}",
195 storageClass=
"Catalog",
196 dimensions=(
"instrument",),
198 vars()[f
"fgcmAtmosphereParameters{cycle}"] = connectionTypes.Output(
199 doc=f
"Catalog of atmospheric fit parameters from fit cycle {cycle}",
200 name=f
"fgcmAtmosphereParameters{cycle}",
201 storageClass=
"Catalog",
202 dimensions=(
"instrument",),
204 vars()[f
"fgcmStandardStars{cycle}"] = connectionTypes.Output(
205 doc=f
"Catalog of standard star magnitudes from fit cycle {cycle}",
206 name=f
"fgcmStandardStars{cycle}",
207 storageClass=
"SimpleCatalog",
208 dimensions=(
"instrument",),
211 def __init__(self, *, config=None):
212 super().__init__(config=config)
214 if not config.doReferenceCalibration:
215 self.inputs.remove(
"fgcmReferenceStars")
217 if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber:
218 raise ValueError(
"cycleNumber must be of integer format")
219 if str(int(config.connections.previousCycleNumber)) != config.connections.previousCycleNumber:
220 raise ValueError(
"previousCycleNumber must be of integer format")
221 if int(config.connections.previousCycleNumber) != (int(config.connections.cycleNumber) - 1):
222 raise ValueError(
"previousCycleNumber must be 1 less than cycleNumber")
224 if int(config.connections.cycleNumber) == 0:
225 self.prerequisiteInputs.remove(
"fgcmFlaggedStarsInput")
226 self.prerequisiteInputs.remove(
"fgcmFitParametersInput")
228 if not self.config.doMultipleCycles:
230 if not self.config.isFinalCycle
and not self.config.outputStandardsBeforeFinalCycle:
231 self.outputs.remove(
"fgcmStandardStars")
233 if not self.config.isFinalCycle
and not self.config.outputZeropointsBeforeFinalCycle:
234 self.outputs.remove(
"fgcmZeropoints")
235 self.outputs.remove(
"fgcmAtmosphereParameters")
238 for cycle
in range(0, MULTIPLE_CYCLES_MAX):
239 self.outputs.remove(f
"fgcmFitParameters{cycle}")
240 self.outputs.remove(f
"fgcmFlaggedStars{cycle}")
241 self.outputs.remove(f
"fgcmZeropoints{cycle}")
242 self.outputs.remove(f
"fgcmAtmosphereParameters{cycle}")
243 self.outputs.remove(f
"fgcmStandardStars{cycle}")
248 self.outputs.remove(
"fgcmFitParameters")
249 self.outputs.remove(
"fgcmFlaggedStars")
250 self.outputs.remove(
"fgcmZeropoints")
251 self.outputs.remove(
"fgcmAtmosphereParameters")
252 self.outputs.remove(
"fgcmStandardStars")
255 for cycle
in range(self.config.multipleCyclesFinalCycleNumber + 1,
256 MULTIPLE_CYCLES_MAX):
257 self.outputs.remove(f
"fgcmFitParameters{cycle}")
258 self.outputs.remove(f
"fgcmFlaggedStars{cycle}")
259 self.outputs.remove(f
"fgcmZeropoints{cycle}")
260 self.outputs.remove(f
"fgcmAtmosphereParameters{cycle}")
261 self.outputs.remove(f
"fgcmStandardStars{cycle}")
264 for cycle
in range(self.config.multipleCyclesFinalCycleNumber):
265 if not self.config.outputZeropointsBeforeFinalCycle:
266 self.outputs.remove(f
"fgcmZeropoints{cycle}")
267 self.outputs.remove(f
"fgcmAtmosphereParameters{cycle}")
268 if not self.config.outputStandardsBeforeFinalCycle:
269 self.outputs.remove(f
"fgcmStandardStars{cycle}")
272 class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
273 pipelineConnections=FgcmFitCycleConnections):
274 """Config for FgcmFitCycle"""
276 doMultipleCycles = pexConfig.Field(
277 doc=
"Run multiple fit cycles in one task",
281 multipleCyclesFinalCycleNumber = pexConfig.RangeField(
282 doc=(
"Final cycle number in multiple cycle mode. The initial cycle "
283 "is 0, with limited parameters fit. The next cycle is 1 with "
284 "full parameter fit. The final cycle is a clean-up with no "
285 "parameters fit. There will be a total of "
286 "(multipleCycleFinalCycleNumber + 1) cycles run, and the final "
287 "cycle number cannot be less than 2."),
291 max=MULTIPLE_CYCLES_MAX,
294 bands = pexConfig.ListField(
295 doc=
"Bands to run calibration",
299 fitBands = pexConfig.ListField(
300 doc=(
"Bands to use in atmospheric fit. The bands not listed here will have "
301 "the atmosphere constrained from the 'fitBands' on the same night. "
302 "Must be a subset of `config.bands`"),
306 requiredBands = pexConfig.ListField(
307 doc=(
"Bands that are required for a star to be considered a calibration star. "
308 "Must be a subset of `config.bands`"),
315 physicalFilterMap = pexConfig.DictField(
316 doc=
"Mapping from 'physicalFilter' to band.",
321 doReferenceCalibration = pexConfig.Field(
322 doc=
"Use reference catalog as additional constraint on calibration",
326 refStarSnMin = pexConfig.Field(
327 doc=
"Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.",
331 refStarOutlierNSig = pexConfig.Field(
332 doc=(
"Number of sigma compared to average mag for reference star to be considered an outlier. "
333 "Computed per-band, and if it is an outlier in any band it is rejected from fits."),
337 applyRefStarColorCuts = pexConfig.Field(
338 doc=
"Apply color cuts to reference stars?",
342 nCore = pexConfig.Field(
343 doc=
"Number of cores to use",
347 nStarPerRun = pexConfig.Field(
348 doc=
"Number of stars to run in each chunk",
352 nExpPerRun = pexConfig.Field(
353 doc=
"Number of exposures to run in each chunk",
357 reserveFraction = pexConfig.Field(
358 doc=
"Fraction of stars to reserve for testing",
362 freezeStdAtmosphere = pexConfig.Field(
363 doc=
"Freeze atmosphere parameters to standard (for testing)",
367 precomputeSuperStarInitialCycle = pexConfig.Field(
368 doc=
"Precompute superstar flat for initial cycle",
372 superStarSubCcdDict = pexConfig.DictField(
373 doc=(
"Per-band specification on whether to compute superstar flat on sub-ccd scale. "
374 "Must have one entry per band."),
379 superStarSubCcdChebyshevOrder = pexConfig.Field(
380 doc=(
"Order of the 2D chebyshev polynomials for sub-ccd superstar fit. "
381 "Global default is first-order polynomials, and should be overridden "
382 "on a camera-by-camera basis depending on the ISR."),
386 superStarSubCcdTriangular = pexConfig.Field(
387 doc=(
"Should the sub-ccd superstar chebyshev matrix be triangular to "
388 "suppress high-order cross terms?"),
392 superStarSigmaClip = pexConfig.Field(
393 doc=
"Number of sigma to clip outliers when selecting for superstar flats",
397 focalPlaneSigmaClip = pexConfig.Field(
398 doc=
"Number of sigma to clip outliers per focal-plane.",
402 ccdGraySubCcdDict = pexConfig.DictField(
403 doc=(
"Per-band specification on whether to compute achromatic per-ccd residual "
404 "('ccd gray') on a sub-ccd scale."),
409 ccdGraySubCcdChebyshevOrder = pexConfig.Field(
410 doc=
"Order of the 2D chebyshev polynomials for sub-ccd gray fit.",
414 ccdGraySubCcdTriangular = pexConfig.Field(
415 doc=(
"Should the sub-ccd gray chebyshev matrix be triangular to "
416 "suppress high-order cross terms?"),
420 ccdGrayFocalPlaneDict = pexConfig.DictField(
421 doc=(
"Per-band specification on whether to compute focal-plane residual "
422 "('ccd gray') corrections."),
427 ccdGrayFocalPlaneFitMinCcd = pexConfig.Field(
428 doc=(
"Minimum number of 'good' CCDs required to perform focal-plane "
429 "gray corrections. If there are fewer good CCDs then the gray "
430 "correction is computed per-ccd."),
434 ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field(
435 doc=
"Order of the 2D chebyshev polynomials for focal plane fit.",
439 cycleNumber = pexConfig.Field(
440 doc=(
"FGCM fit cycle number. This is automatically incremented after each run "
441 "and stage of outlier rejection. See cookbook for details."),
445 isFinalCycle = pexConfig.Field(
446 doc=(
"Is this the final cycle of the fitting? Will automatically compute final "
447 "selection of stars and photometric exposures, and will output zeropoints "
448 "and standard stars for use in fgcmOutputProducts"),
452 maxIterBeforeFinalCycle = pexConfig.Field(
453 doc=(
"Maximum fit iterations, prior to final cycle. The number of iterations "
454 "will always be 0 in the final cycle for cleanup and final selection."),
458 deltaMagBkgOffsetPercentile = pexConfig.Field(
459 doc=(
"Percentile brightest stars on a visit/ccd to use to compute net "
460 "offset from local background subtraction."),
464 deltaMagBkgPerCcd = pexConfig.Field(
465 doc=(
"Compute net offset from local background subtraction per-ccd? "
466 "Otherwise, use computation per visit."),
470 utBoundary = pexConfig.Field(
471 doc=
"Boundary (in UTC) from day-to-day",
475 washMjds = pexConfig.ListField(
476 doc=
"Mirror wash MJDs",
480 epochMjds = pexConfig.ListField(
481 doc=
"Epoch boundaries in MJD",
485 minObsPerBand = pexConfig.Field(
486 doc=
"Minimum good observations per band",
492 latitude = pexConfig.Field(
493 doc=
"Observatory latitude",
497 defaultCameraOrientation = pexConfig.Field(
498 doc=
"Default camera orientation for QA plots.",
502 brightObsGrayMax = pexConfig.Field(
503 doc=
"Maximum gray extinction to be considered bright observation",
507 minStarPerCcd = pexConfig.Field(
508 doc=(
"Minimum number of good stars per CCD to be used in calibration fit. "
509 "CCDs with fewer stars will have their calibration estimated from other "
510 "CCDs in the same visit, with zeropoint error increased accordingly."),
514 minCcdPerExp = pexConfig.Field(
515 doc=(
"Minimum number of good CCDs per exposure/visit to be used in calibration fit. "
516 "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."),
520 maxCcdGrayErr = pexConfig.Field(
521 doc=
"Maximum error on CCD gray offset to be considered photometric",
525 minStarPerExp = pexConfig.Field(
526 doc=(
"Minimum number of good stars per exposure/visit to be used in calibration fit. "
527 "Visits with fewer good stars will have CCD zeropoints estimated where possible."),
531 minExpPerNight = pexConfig.Field(
532 doc=
"Minimum number of good exposures/visits to consider a partly photometric night",
536 expGrayInitialCut = pexConfig.Field(
537 doc=(
"Maximum exposure/visit gray value for initial selection of possible photometric "
542 expGrayPhotometricCutDict = pexConfig.DictField(
543 doc=(
"Per-band specification on maximum (negative) achromatic exposure residual "
544 "('gray term') for a visit to be considered photometric. Must have one "
545 "entry per band. Broad-band filters should be -0.05."),
550 expGrayHighCutDict = pexConfig.DictField(
551 doc=(
"Per-band specification on maximum (positive) achromatic exposure residual "
552 "('gray term') for a visit to be considered photometric. Must have one "
553 "entry per band. Broad-band filters should be 0.2."),
558 expGrayRecoverCut = pexConfig.Field(
559 doc=(
"Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. "
560 "Visits with more gray extinction will only get CCD zeropoints if there are "
561 "sufficient star observations (minStarPerCcd) on that CCD."),
565 expVarGrayPhotometricCutDict = pexConfig.DictField(
566 doc=(
"Per-band specification on maximum exposure variance to be considered possibly "
567 "photometric. Must have one entry per band. Broad-band filters should be "
573 expGrayErrRecoverCut = pexConfig.Field(
574 doc=(
"Maximum exposure gray error to be able to recover bad ccds via interpolation. "
575 "Visits with more gray variance will only get CCD zeropoints if there are "
576 "sufficient star observations (minStarPerCcd) on that CCD."),
580 aperCorrFitNBins = pexConfig.Field(
581 doc=(
"Number of aperture bins used in aperture correction fit. When set to 0"
582 "no fit will be performed, and the config.aperCorrInputSlopes will be "
583 "used if available."),
587 aperCorrInputSlopeDict = pexConfig.DictField(
588 doc=(
"Per-band specification of aperture correction input slope parameters. These "
589 "are used on the first fit iteration, and aperture correction parameters will "
590 "be updated from the data if config.aperCorrFitNBins > 0. It is recommended "
591 "to set this when there is insufficient data to fit the parameters (e.g. "
597 sedboundaryterms = pexConfig.ConfigField(
598 doc=
"Mapping from bands to SED boundary term names used is sedterms.",
599 dtype=SedboundarytermDict,
601 sedterms = pexConfig.ConfigField(
602 doc=
"Mapping from terms to bands for fgcm linear SED approximations.",
605 sigFgcmMaxErr = pexConfig.Field(
606 doc=
"Maximum mag error for fitting sigma_FGCM",
610 sigFgcmMaxEGrayDict = pexConfig.DictField(
611 doc=(
"Per-band specification for maximum (absolute) achromatic residual (gray value) "
612 "for observations in sigma_fgcm (raw repeatability). Broad-band filters "
618 ccdGrayMaxStarErr = pexConfig.Field(
619 doc=(
"Maximum error on a star observation to use in ccd gray (achromatic residual) "
624 approxThroughputDict = pexConfig.DictField(
625 doc=(
"Per-band specification of the approximate overall throughput at the start of "
626 "calibration observations. Must have one entry per band. Typically should "
632 sigmaCalRange = pexConfig.ListField(
633 doc=
"Allowed range for systematic error floor estimation",
635 default=(0.001, 0.003),
637 sigmaCalFitPercentile = pexConfig.ListField(
638 doc=
"Magnitude percentile range to fit systematic error floor",
640 default=(0.05, 0.15),
642 sigmaCalPlotPercentile = pexConfig.ListField(
643 doc=
"Magnitude percentile range to plot systematic error floor",
645 default=(0.05, 0.95),
647 sigma0Phot = pexConfig.Field(
648 doc=
"Systematic error floor for all zeropoints",
652 mapLongitudeRef = pexConfig.Field(
653 doc=
"Reference longitude for plotting maps",
657 mapNSide = pexConfig.Field(
658 doc=
"Healpix nside for plotting maps",
662 outfileBase = pexConfig.Field(
663 doc=
"Filename start for plot output files",
667 starColorCuts = pexConfig.ListField(
668 doc=
"Encoded star-color cuts (to be cleaned up)",
670 default=(
"NO_DATA",),
672 colorSplitBands = pexConfig.ListField(
673 doc=
"Band names to use to split stars by color. Must have 2 entries.",
678 modelMagErrors = pexConfig.Field(
679 doc=
"Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)",
683 useQuadraticPwv = pexConfig.Field(
684 doc=
"Model PWV with a quadratic term for variation through the night?",
688 instrumentParsPerBand = pexConfig.Field(
689 doc=(
"Model instrumental parameters per band? "
690 "Otherwise, instrumental parameters (QE changes with time) are "
691 "shared among all bands."),
695 instrumentSlopeMinDeltaT = pexConfig.Field(
696 doc=(
"Minimum time change (in days) between observations to use in constraining "
697 "instrument slope."),
701 fitMirrorChromaticity = pexConfig.Field(
702 doc=
"Fit (intraband) mirror chromatic term?",
706 coatingMjds = pexConfig.ListField(
707 doc=
"Mirror coating dates in MJD",
711 outputStandardsBeforeFinalCycle = pexConfig.Field(
712 doc=
"Output standard stars prior to final cycle? Used in debugging.",
716 outputZeropointsBeforeFinalCycle = pexConfig.Field(
717 doc=
"Output standard stars prior to final cycle? Used in debugging.",
721 useRepeatabilityForExpGrayCutsDict = pexConfig.DictField(
722 doc=(
"Per-band specification on whether to use star repeatability (instead of exposures) "
723 "for computing photometric cuts. Recommended for tract mode or bands with few visits."),
728 autoPhotometricCutNSig = pexConfig.Field(
729 doc=(
"Number of sigma for automatic computation of (low) photometric cut. "
730 "Cut is based on exposure gray width (per band), unless "
731 "useRepeatabilityForExpGrayCuts is set, in which case the star "
732 "repeatability is used (also per band)."),
736 autoHighCutNSig = pexConfig.Field(
737 doc=(
"Number of sigma for automatic computation of (high) outlier cut. "
738 "Cut is based on exposure gray width (per band), unless "
739 "useRepeatabilityForExpGrayCuts is set, in which case the star "
740 "repeatability is used (also per band)."),
744 quietMode = pexConfig.Field(
745 doc=
"Be less verbose with logging.",
749 doPlots = pexConfig.Field(
750 doc=
"Make fgcm QA plots.",
754 randomSeed = pexConfig.Field(
755 doc=
"Random seed for fgcm for consistency in tests.",
764 if self.connections.previousCycleNumber != str(self.cycleNumber - 1):
765 msg =
"cycleNumber in template must be connections.previousCycleNumber + 1"
766 raise RuntimeError(msg)
767 if self.connections.cycleNumber != str(self.cycleNumber):
768 msg =
"cycleNumber in template must be equal to connections.cycleNumber"
769 raise RuntimeError(msg)
771 for band
in self.fitBands:
772 if band
not in self.bands:
773 msg =
'fitBand %s not in bands' % (band)
774 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg)
775 for band
in self.requiredBands:
776 if band
not in self.bands:
777 msg =
'requiredBand %s not in bands' % (band)
778 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg)
779 for band
in self.colorSplitBands:
780 if band
not in self.bands:
781 msg =
'colorSplitBand %s not in bands' % (band)
782 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg)
783 for band
in self.bands:
784 if band
not in self.superStarSubCcdDict:
785 msg =
'band %s not in superStarSubCcdDict' % (band)
786 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict,
788 if band
not in self.ccdGraySubCcdDict:
789 msg =
'band %s not in ccdGraySubCcdDict' % (band)
790 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict,
792 if band
not in self.expGrayPhotometricCutDict:
793 msg =
'band %s not in expGrayPhotometricCutDict' % (band)
794 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict,
796 if band
not in self.expGrayHighCutDict:
797 msg =
'band %s not in expGrayHighCutDict' % (band)
798 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict,
800 if band
not in self.expVarGrayPhotometricCutDict:
801 msg =
'band %s not in expVarGrayPhotometricCutDict' % (band)
802 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict,
804 if band
not in self.sigFgcmMaxEGrayDict:
805 msg =
'band %s not in sigFgcmMaxEGrayDict' % (band)
806 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict,
808 if band
not in self.approxThroughputDict:
809 msg =
'band %s not in approxThroughputDict' % (band)
810 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict,
812 if band
not in self.useRepeatabilityForExpGrayCutsDict:
813 msg =
'band %s not in useRepeatabilityForExpGrayCutsDict' % (band)
814 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
818 class FgcmFitCycleRunner(pipeBase.ButlerInitializedTaskRunner):
819 """Subclass of TaskRunner for fgcmFitCycleTask
821 fgcmFitCycleTask.run() takes one argument, the butler, and uses
822 stars and visits previously extracted from dataRefs by
824 This Runner does not perform any dataRef parallelization, but the FGCM
825 code called by the Task uses python multiprocessing (see the "ncores"
830 def getTargetList(parsedCmd):
832 Return a list with one element, the butler.
834 return [parsedCmd.butler]
836 def __call__(self, butler):
840 butler: `lsst.daf.persistence.Butler`
844 exitStatus: `list` with `pipeBase.Struct`
845 exitStatus (0: success; 1: failure)
848 task = self.TaskClass(config=self.config, log=self.log)
852 task.runDataRef(butler)
855 task.runDataRef(butler)
856 except Exception
as e:
858 task.log.fatal(
"Failed: %s" % e)
859 if not isinstance(e, pipeBase.TaskError):
860 traceback.print_exc(file=sys.stderr)
862 task.writeMetadata(butler)
865 return [pipeBase.Struct(exitStatus=exitStatus)]
867 def run(self, parsedCmd):
869 Run the task, with no multiprocessing
873 parsedCmd: ArgumentParser parsed command line
878 if self.precall(parsedCmd):
879 targetList = self.getTargetList(parsedCmd)
881 resultList = self(targetList[0])
886 class FgcmFitCycleTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
888 Run Single fit cycle for FGCM global calibration
891 ConfigClass = FgcmFitCycleConfig
892 RunnerClass = FgcmFitCycleRunner
893 _DefaultName =
"fgcmFitCycle"
895 def __init__(self, butler=None, initInputs=None, **kwargs):
896 super().__init__(**kwargs)
899 def _getMetadataName(self):
902 def runQuantum(self, butlerQC, inputRefs, outputRefs):
903 camera = butlerQC.get(inputRefs.camera)
907 dataRefDict[
'fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
908 dataRefDict[
'fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
909 dataRefDict[
'fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations)
910 dataRefDict[
'fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds)
911 dataRefDict[
'fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices)
912 if self.config.doReferenceCalibration:
913 dataRefDict[
'fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars)
914 if self.config.cycleNumber > 0:
915 dataRefDict[
'fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput)
916 dataRefDict[
'fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput)
918 fgcmDatasetDict =
None
919 if self.config.doMultipleCycles:
921 config = copy.copy(self.config)
922 config.update(cycleNumber=0)
923 for cycle
in range(self.config.multipleCyclesFinalCycleNumber + 1):
924 if cycle == self.config.multipleCyclesFinalCycleNumber:
925 config.update(isFinalCycle=
True)
928 dataRefDict[
'fgcmFlaggedStars'] = fgcmDatasetDict[
'fgcmFlaggedStars']
929 dataRefDict[
'fgcmFitParameters'] = fgcmDatasetDict[
'fgcmFitParameters']
931 fgcmDatasetDict, config = self._fgcmFitCycle(camera, dataRefDict, config=config)
932 butlerQC.put(fgcmDatasetDict[
'fgcmFitParameters'],
933 getattr(outputRefs, f
'fgcmFitParameters{cycle}'))
934 butlerQC.put(fgcmDatasetDict[
'fgcmFlaggedStars'],
935 getattr(outputRefs, f
'fgcmFlaggedStars{cycle}'))
936 if self.outputZeropoints:
937 butlerQC.put(fgcmDatasetDict[
'fgcmZeropoints'],
938 getattr(outputRefs, f
'fgcmZeropoints{cycle}'))
939 butlerQC.put(fgcmDatasetDict[
'fgcmAtmosphereParameters'],
940 getattr(outputRefs, f
'fgcmAtmosphereParameters{cycle}'))
941 if self.outputStandards:
942 butlerQC.put(fgcmDatasetDict[
'fgcmStandardStars'],
943 getattr(outputRefs, f
'fgcmStandardStars{cycle}'))
946 fgcmDatasetDict, _ = self._fgcmFitCycle(camera, dataRefDict)
948 butlerQC.put(fgcmDatasetDict[
'fgcmFitParameters'], outputRefs.fgcmFitParameters)
949 butlerQC.put(fgcmDatasetDict[
'fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars)
950 if self.outputZeropoints:
951 butlerQC.put(fgcmDatasetDict[
'fgcmZeropoints'], outputRefs.fgcmZeropoints)
952 butlerQC.put(fgcmDatasetDict[
'fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters)
953 if self.outputStandards:
954 butlerQC.put(fgcmDatasetDict[
'fgcmStandardStars'], outputRefs.fgcmStandardStars)
957 def runDataRef(self, butler):
959 Run a single fit cycle for FGCM
963 butler: `lsst.daf.persistence.Butler`
965 self._checkDatasetsExist(butler)
968 dataRefDict[
'fgcmLookUpTable'] = butler.dataRef(
'fgcmLookUpTable')
969 dataRefDict[
'fgcmVisitCatalog'] = butler.dataRef(
'fgcmVisitCatalog')
970 dataRefDict[
'fgcmStarObservations'] = butler.dataRef(
'fgcmStarObservations')
971 dataRefDict[
'fgcmStarIds'] = butler.dataRef(
'fgcmStarIds')
972 dataRefDict[
'fgcmStarIndices'] = butler.dataRef(
'fgcmStarIndices')
973 if self.config.doReferenceCalibration:
974 dataRefDict[
'fgcmReferenceStars'] = butler.dataRef(
'fgcmReferenceStars')
975 if self.config.cycleNumber > 0:
976 lastCycle = self.config.cycleNumber - 1
977 dataRefDict[
'fgcmFlaggedStars'] = butler.dataRef(
'fgcmFlaggedStars',
979 dataRefDict[
'fgcmFitParameters'] = butler.dataRef(
'fgcmFitParameters',
982 camera = butler.get(
'camera')
983 fgcmDatasetDict, _ = self._fgcmFitCycle(camera, dataRefDict)
985 butler.put(fgcmDatasetDict[
'fgcmFitParameters'],
'fgcmFitParameters',
986 fgcmcycle=self.config.cycleNumber)
987 butler.put(fgcmDatasetDict[
'fgcmFlaggedStars'],
'fgcmFlaggedStars',
988 fgcmcycle=self.config.cycleNumber)
989 if self.outputZeropoints:
990 butler.put(fgcmDatasetDict[
'fgcmZeropoints'],
'fgcmZeropoints',
991 fgcmcycle=self.config.cycleNumber)
992 butler.put(fgcmDatasetDict[
'fgcmAtmosphereParameters'],
'fgcmAtmosphereParameters',
993 fgcmcycle=self.config.cycleNumber)
994 if self.outputStandards:
995 butler.put(fgcmDatasetDict[
'fgcmStandardStars'],
'fgcmStandardStars',
996 fgcmcycle=self.config.cycleNumber)
998 def writeConfig(self, butler, clobber=False, doBackup=True):
999 """Write the configuration used for processing the data, or check that an existing
1000 one is equal to the new one if present. This is an override of the regular
1001 version from pipe_base that knows about fgcmcycle.
1005 butler : `lsst.daf.persistence.Butler`
1006 Data butler used to write the config. The config is written to dataset type
1007 `CmdLineTask._getConfigName`.
1008 clobber : `bool`, optional
1009 A boolean flag that controls what happens if a config already has been saved:
1010 - `True`: overwrite or rename the existing config, depending on ``doBackup``.
1011 - `False`: raise `TaskError` if this config does not match the existing config.
1012 doBackup : `bool`, optional
1013 Set to `True` to backup the config files if clobbering.
1015 configName = self._getConfigName()
1016 if configName
is None:
1019 butler.put(self.config, configName, doBackup=doBackup, fgcmcycle=self.config.cycleNumber)
1020 elif butler.datasetExists(configName, write=
True, fgcmcycle=self.config.cycleNumber):
1023 oldConfig = butler.get(configName, immediate=
True, fgcmcycle=self.config.cycleNumber)
1024 except Exception
as exc:
1025 raise type(exc)(
"Unable to read stored config file %s (%s); consider using --clobber-config" %
1028 def logConfigMismatch(msg):
1029 self.log.
fatal(
"Comparing configuration: %s", msg)
1031 if not self.config.compare(oldConfig, shortcut=
False, output=logConfigMismatch):
1032 raise pipeBase.TaskError(
1033 f
"Config does not match existing task config {configName!r} on disk; tasks configurations"
1034 " must be consistent within the same output repo (override with --clobber-config)")
1036 butler.put(self.config, configName, fgcmcycle=self.config.cycleNumber)
1038 def _fgcmFitCycle(self, camera, dataRefDict, config=None):
1044 camera : `lsst.afw.cameraGeom.Camera`
1045 dataRefDict : `dict`
1046 All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
1047 `lsst.daf.butler.DeferredDatasetHandle` (gen3)
1048 dataRef dictionary with keys:
1050 ``"fgcmLookUpTable"``
1051 dataRef for the FGCM look-up table.
1052 ``"fgcmVisitCatalog"``
1053 dataRef for visit summary catalog.
1054 ``"fgcmStarObservations"``
1055 dataRef for star observation catalog.
1057 dataRef for star id catalog.
1058 ``"fgcmStarIndices"``
1059 dataRef for star index catalog.
1060 ``"fgcmReferenceStars"``
1061 dataRef for matched reference star catalog.
1062 ``"fgcmFlaggedStars"``
1063 dataRef for flagged star catalog.
1064 ``"fgcmFitParameters"``
1065 dataRef for fit parameter catalog.
1066 config : `lsst.pex.config.Config`, optional
1067 Configuration to use to override self.config.
1071 fgcmDatasetDict : `dict`
1072 Dictionary of datasets to persist.
1074 if config
is not None:
1077 _config = self.config
1080 self.maxIter = _config.maxIterBeforeFinalCycle
1081 self.outputStandards = _config.outputStandardsBeforeFinalCycle
1082 self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle
1083 self.resetFitParameters =
True
1085 if _config.isFinalCycle:
1090 self.outputStandards =
True
1091 self.outputZeropoints =
True
1092 self.resetFitParameters =
False
1094 lutCat = dataRefDict[
'fgcmLookUpTable'].get()
1096 dict(_config.physicalFilterMap))
1100 self.maxIter, self.resetFitParameters,
1101 self.outputZeropoints,
1102 lutIndexVals[0][
'FILTERNAMES'])
1105 visitCat = dataRefDict[
'fgcmVisitCatalog'].get()
1110 self.config.defaultCameraOrientation)
1112 noFitsDict = {
'lutIndex': lutIndexVals,
1114 'expInfo': fgcmExpInfo,
1115 'focalPlaneProjector': focalPlaneProjector}
1118 fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=
False,
1119 noFitsDict=noFitsDict, noOutput=
True)
1122 if (fgcmFitCycle.initialCycle):
1124 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
1129 parCat = dataRefDict[
'fgcmFitParameters']
1131 parCat = dataRefDict[
'fgcmFitParameters'].get()
1132 inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1134 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
1141 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
1143 starObs = dataRefDict[
'fgcmStarObservations'].get()
1144 starIds = dataRefDict[
'fgcmStarIds'].get()
1145 starIndices = dataRefDict[
'fgcmStarIndices'].get()
1148 if 'fgcmFlaggedStars' in dataRefDict:
1150 flaggedStars = dataRefDict[
'fgcmFlaggedStars']
1152 flaggedStars = dataRefDict[
'fgcmFlaggedStars'].get()
1153 flagId = flaggedStars[
'objId'][:]
1154 flagFlag = flaggedStars[
'objFlag'][:]
1160 if _config.doReferenceCalibration:
1161 refStars = dataRefDict[
'fgcmReferenceStars'].get()
1165 _config.physicalFilterMap)
1166 refId = refStars[
'fgcm_id'][:]
1176 visitIndex = np.searchsorted(fgcmExpInfo[
'VISIT'], starObs[
'visit'][starIndices[
'obsIndex']])
1188 conv = starObs[0][
'ra'].asDegrees() / float(starObs[0][
'ra'])
1190 fgcmStars.loadStars(fgcmPars,
1191 starObs[
'visit'][starIndices[
'obsIndex']],
1192 starObs[
'ccd'][starIndices[
'obsIndex']],
1193 starObs[
'ra'][starIndices[
'obsIndex']] * conv,
1194 starObs[
'dec'][starIndices[
'obsIndex']] * conv,
1195 starObs[
'instMag'][starIndices[
'obsIndex']],
1196 starObs[
'instMagErr'][starIndices[
'obsIndex']],
1197 fgcmExpInfo[
'FILTERNAME'][visitIndex],
1198 starIds[
'fgcm_id'][:],
1201 starIds[
'obsArrIndex'][:],
1203 obsX=starObs[
'x'][starIndices[
'obsIndex']],
1204 obsY=starObs[
'y'][starIndices[
'obsIndex']],
1205 obsDeltaMagBkg=starObs[
'deltaMagBkg'][starIndices[
'obsIndex']],
1206 psfCandidate=starObs[
'psf_candidate'][starIndices[
'obsIndex']],
1209 refMagErr=refMagErr,
1227 fgcmFitCycle.setLUT(fgcmLut)
1228 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
1229 fgcmFitCycle.setPars(fgcmPars)
1232 fgcmFitCycle.finishSetup()
1241 fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle)
1246 updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i])
for
1247 i, b
in enumerate(_config.bands)}
1248 updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i])
for
1249 i, band
in enumerate(_config.bands)}
1251 outConfig = copy.copy(_config)
1252 outConfig.update(cycleNumber=(_config.cycleNumber + 1),
1253 precomputeSuperStarInitialCycle=
False,
1254 freezeStdAtmosphere=
False,
1255 expGrayPhotometricCutDict=updatedPhotometricCutDict,
1256 expGrayHighCutDict=updatedHighCutDict)
1258 outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber),
1259 cycleNumber=str(_config.cycleNumber + 1))
1261 configFileName =
'%s_cycle%02d_config.py' % (outConfig.outfileBase,
1262 outConfig.cycleNumber)
1263 outConfig.save(configFileName)
1265 if _config.isFinalCycle == 1:
1267 self.log.
info(
"Everything is in place to run fgcmOutputProducts.py")
1269 self.log.
info(
"Saved config for next cycle to %s" % (configFileName))
1270 self.log.
info(
"Be sure to look at:")
1271 self.log.
info(
" config.expGrayPhotometricCut")
1272 self.log.
info(
" config.expGrayHighCut")
1273 self.log.
info(
"If you are satisfied with the fit, please set:")
1274 self.log.
info(
" config.isFinalCycle = True")
1276 fgcmFitCycle.freeSharedMemory()
1278 return fgcmDatasetDict, outConfig
1280 def _checkDatasetsExist(self, butler):
1282 Check if necessary datasets exist to run fgcmFitCycle
1286 butler: `lsst.daf.persistence.Butler`
1291 If any of fgcmVisitCatalog, fgcmStarObservations, fgcmStarIds,
1292 fgcmStarIndices, fgcmLookUpTable datasets do not exist.
1293 If cycleNumber > 0, then also checks for fgcmFitParameters,
1297 if not butler.datasetExists(
'fgcmVisitCatalog'):
1298 raise RuntimeError(
"Could not find fgcmVisitCatalog in repo!")
1299 if not butler.datasetExists(
'fgcmStarObservations'):
1300 raise RuntimeError(
"Could not find fgcmStarObservations in repo!")
1301 if not butler.datasetExists(
'fgcmStarIds'):
1302 raise RuntimeError(
"Could not find fgcmStarIds in repo!")
1303 if not butler.datasetExists(
'fgcmStarIndices'):
1304 raise RuntimeError(
"Could not find fgcmStarIndices in repo!")
1305 if not butler.datasetExists(
'fgcmLookUpTable'):
1306 raise RuntimeError(
"Could not find fgcmLookUpTable in repo!")
1309 if (self.config.cycleNumber > 0):
1310 if not butler.datasetExists(
'fgcmFitParameters',
1311 fgcmcycle=self.config.cycleNumber-1):
1312 raise RuntimeError(
"Could not find fgcmFitParameters for previous cycle (%d) in repo!" %
1313 (self.config.cycleNumber-1))
1314 if not butler.datasetExists(
'fgcmFlaggedStars',
1315 fgcmcycle=self.config.cycleNumber-1):
1316 raise RuntimeError(
"Could not find fgcmFlaggedStars for previous cycle (%d) in repo!" %
1317 (self.config.cycleNumber-1))
1320 if self.config.doReferenceCalibration:
1321 if not butler.datasetExists(
'fgcmReferenceStars'):
1322 raise RuntimeError(
"Could not find fgcmReferenceStars in repo, and "
1323 "doReferenceCalibration is True.")
1325 def _loadParameters(self, parCat):
1327 Load FGCM parameters from a previous fit cycle
1331 parCat : `lsst.afw.table.BaseCatalog`
1332 Parameter catalog in afw table form.
1336 inParInfo: `numpy.ndarray`
1337 Numpy array parameter information formatted for input to fgcm
1338 inParameters: `numpy.ndarray`
1339 Numpy array parameter values formatted for input to fgcm
1340 inSuperStar: `numpy.array`
1341 Superstar flat formatted for input to fgcm
1343 parLutFilterNames = np.array(parCat[0][
'lutFilterNames'].split(
','))
1344 parFitBands = np.array(parCat[0][
'fitBands'].split(
','))
1346 inParInfo = np.zeros(1, dtype=[(
'NCCD',
'i4'),
1347 (
'LUTFILTERNAMES', parLutFilterNames.dtype.str,
1348 (parLutFilterNames.size, )),
1349 (
'FITBANDS', parFitBands.dtype.str, (parFitBands.size, )),
1350 (
'LNTAUUNIT',
'f8'),
1351 (
'LNTAUSLOPEUNIT',
'f8'),
1352 (
'ALPHAUNIT',
'f8'),
1353 (
'LNPWVUNIT',
'f8'),
1354 (
'LNPWVSLOPEUNIT',
'f8'),
1355 (
'LNPWVQUADRATICUNIT',
'f8'),
1356 (
'LNPWVGLOBALUNIT',
'f8'),
1358 (
'QESYSUNIT',
'f8'),
1359 (
'FILTEROFFSETUNIT',
'f8'),
1360 (
'HASEXTERNALPWV',
'i2'),
1361 (
'HASEXTERNALTAU',
'i2')])
1362 inParInfo[
'NCCD'] = parCat[
'nCcd']
1363 inParInfo[
'LUTFILTERNAMES'][:] = parLutFilterNames
1364 inParInfo[
'FITBANDS'][:] = parFitBands
1365 inParInfo[
'HASEXTERNALPWV'] = parCat[
'hasExternalPwv']
1366 inParInfo[
'HASEXTERNALTAU'] = parCat[
'hasExternalTau']
1368 inParams = np.zeros(1, dtype=[(
'PARALPHA',
'f8', (parCat[
'parAlpha'].size, )),
1369 (
'PARO3',
'f8', (parCat[
'parO3'].size, )),
1370 (
'PARLNTAUINTERCEPT',
'f8',
1371 (parCat[
'parLnTauIntercept'].size, )),
1372 (
'PARLNTAUSLOPE',
'f8',
1373 (parCat[
'parLnTauSlope'].size, )),
1374 (
'PARLNPWVINTERCEPT',
'f8',
1375 (parCat[
'parLnPwvIntercept'].size, )),
1376 (
'PARLNPWVSLOPE',
'f8',
1377 (parCat[
'parLnPwvSlope'].size, )),
1378 (
'PARLNPWVQUADRATIC',
'f8',
1379 (parCat[
'parLnPwvQuadratic'].size, )),
1380 (
'PARQESYSINTERCEPT',
'f8',
1381 (parCat[
'parQeSysIntercept'].size, )),
1382 (
'COMPQESYSSLOPE',
'f8',
1383 (parCat[
'compQeSysSlope'].size, )),
1384 (
'PARFILTEROFFSET',
'f8',
1385 (parCat[
'parFilterOffset'].size, )),
1386 (
'PARFILTEROFFSETFITFLAG',
'i2',
1387 (parCat[
'parFilterOffsetFitFlag'].size, )),
1388 (
'PARRETRIEVEDLNPWVSCALE',
'f8'),
1389 (
'PARRETRIEVEDLNPWVOFFSET',
'f8'),
1390 (
'PARRETRIEVEDLNPWVNIGHTLYOFFSET',
'f8',
1391 (parCat[
'parRetrievedLnPwvNightlyOffset'].size, )),
1392 (
'COMPABSTHROUGHPUT',
'f8',
1393 (parCat[
'compAbsThroughput'].size, )),
1394 (
'COMPREFOFFSET',
'f8',
1395 (parCat[
'compRefOffset'].size, )),
1396 (
'COMPREFSIGMA',
'f8',
1397 (parCat[
'compRefSigma'].size, )),
1398 (
'COMPMIRRORCHROMATICITY',
'f8',
1399 (parCat[
'compMirrorChromaticity'].size, )),
1400 (
'MIRRORCHROMATICITYPIVOT',
'f8',
1401 (parCat[
'mirrorChromaticityPivot'].size, )),
1402 (
'COMPMEDIANSEDSLOPE',
'f8',
1403 (parCat[
'compMedianSedSlope'].size, )),
1404 (
'COMPAPERCORRPIVOT',
'f8',
1405 (parCat[
'compAperCorrPivot'].size, )),
1406 (
'COMPAPERCORRSLOPE',
'f8',
1407 (parCat[
'compAperCorrSlope'].size, )),
1408 (
'COMPAPERCORRSLOPEERR',
'f8',
1409 (parCat[
'compAperCorrSlopeErr'].size, )),
1410 (
'COMPAPERCORRRANGE',
'f8',
1411 (parCat[
'compAperCorrRange'].size, )),
1412 (
'COMPMODELERREXPTIMEPIVOT',
'f8',
1413 (parCat[
'compModelErrExptimePivot'].size, )),
1414 (
'COMPMODELERRFWHMPIVOT',
'f8',
1415 (parCat[
'compModelErrFwhmPivot'].size, )),
1416 (
'COMPMODELERRSKYPIVOT',
'f8',
1417 (parCat[
'compModelErrSkyPivot'].size, )),
1418 (
'COMPMODELERRPARS',
'f8',
1419 (parCat[
'compModelErrPars'].size, )),
1420 (
'COMPEXPGRAY',
'f8',
1421 (parCat[
'compExpGray'].size, )),
1422 (
'COMPVARGRAY',
'f8',
1423 (parCat[
'compVarGray'].size, )),
1424 (
'COMPEXPDELTAMAGBKG',
'f8',
1425 (parCat[
'compExpDeltaMagBkg'].size, )),
1426 (
'COMPNGOODSTARPEREXP',
'i4',
1427 (parCat[
'compNGoodStarPerExp'].size, )),
1428 (
'COMPSIGFGCM',
'f8',
1429 (parCat[
'compSigFgcm'].size, )),
1430 (
'COMPSIGMACAL',
'f8',
1431 (parCat[
'compSigmaCal'].size, )),
1432 (
'COMPRETRIEVEDLNPWV',
'f8',
1433 (parCat[
'compRetrievedLnPwv'].size, )),
1434 (
'COMPRETRIEVEDLNPWVRAW',
'f8',
1435 (parCat[
'compRetrievedLnPwvRaw'].size, )),
1436 (
'COMPRETRIEVEDLNPWVFLAG',
'i2',
1437 (parCat[
'compRetrievedLnPwvFlag'].size, )),
1438 (
'COMPRETRIEVEDTAUNIGHT',
'f8',
1439 (parCat[
'compRetrievedTauNight'].size, ))])
1441 inParams[
'PARALPHA'][:] = parCat[
'parAlpha'][0, :]
1442 inParams[
'PARO3'][:] = parCat[
'parO3'][0, :]
1443 inParams[
'PARLNTAUINTERCEPT'][:] = parCat[
'parLnTauIntercept'][0, :]
1444 inParams[
'PARLNTAUSLOPE'][:] = parCat[
'parLnTauSlope'][0, :]
1445 inParams[
'PARLNPWVINTERCEPT'][:] = parCat[
'parLnPwvIntercept'][0, :]
1446 inParams[
'PARLNPWVSLOPE'][:] = parCat[
'parLnPwvSlope'][0, :]
1447 inParams[
'PARLNPWVQUADRATIC'][:] = parCat[
'parLnPwvQuadratic'][0, :]
1448 inParams[
'PARQESYSINTERCEPT'][:] = parCat[
'parQeSysIntercept'][0, :]
1449 inParams[
'COMPQESYSSLOPE'][:] = parCat[
'compQeSysSlope'][0, :]
1450 inParams[
'PARFILTEROFFSET'][:] = parCat[
'parFilterOffset'][0, :]
1451 inParams[
'PARFILTEROFFSETFITFLAG'][:] = parCat[
'parFilterOffsetFitFlag'][0, :]
1452 inParams[
'PARRETRIEVEDLNPWVSCALE'] = parCat[
'parRetrievedLnPwvScale']
1453 inParams[
'PARRETRIEVEDLNPWVOFFSET'] = parCat[
'parRetrievedLnPwvOffset']
1454 inParams[
'PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat[
'parRetrievedLnPwvNightlyOffset'][0, :]
1455 inParams[
'COMPABSTHROUGHPUT'][:] = parCat[
'compAbsThroughput'][0, :]
1456 inParams[
'COMPREFOFFSET'][:] = parCat[
'compRefOffset'][0, :]
1457 inParams[
'COMPREFSIGMA'][:] = parCat[
'compRefSigma'][0, :]
1458 inParams[
'COMPMIRRORCHROMATICITY'][:] = parCat[
'compMirrorChromaticity'][0, :]
1459 inParams[
'MIRRORCHROMATICITYPIVOT'][:] = parCat[
'mirrorChromaticityPivot'][0, :]
1460 inParams[
'COMPMEDIANSEDSLOPE'][:] = parCat[
'compMedianSedSlope'][0, :]
1461 inParams[
'COMPAPERCORRPIVOT'][:] = parCat[
'compAperCorrPivot'][0, :]
1462 inParams[
'COMPAPERCORRSLOPE'][:] = parCat[
'compAperCorrSlope'][0, :]
1463 inParams[
'COMPAPERCORRSLOPEERR'][:] = parCat[
'compAperCorrSlopeErr'][0, :]
1464 inParams[
'COMPAPERCORRRANGE'][:] = parCat[
'compAperCorrRange'][0, :]
1465 inParams[
'COMPMODELERREXPTIMEPIVOT'][:] = parCat[
'compModelErrExptimePivot'][0, :]
1466 inParams[
'COMPMODELERRFWHMPIVOT'][:] = parCat[
'compModelErrFwhmPivot'][0, :]
1467 inParams[
'COMPMODELERRSKYPIVOT'][:] = parCat[
'compModelErrSkyPivot'][0, :]
1468 inParams[
'COMPMODELERRPARS'][:] = parCat[
'compModelErrPars'][0, :]
1469 inParams[
'COMPEXPGRAY'][:] = parCat[
'compExpGray'][0, :]
1470 inParams[
'COMPVARGRAY'][:] = parCat[
'compVarGray'][0, :]
1471 inParams[
'COMPEXPDELTAMAGBKG'][:] = parCat[
'compExpDeltaMagBkg'][0, :]
1472 inParams[
'COMPNGOODSTARPEREXP'][:] = parCat[
'compNGoodStarPerExp'][0, :]
1473 inParams[
'COMPSIGFGCM'][:] = parCat[
'compSigFgcm'][0, :]
1474 inParams[
'COMPSIGMACAL'][:] = parCat[
'compSigmaCal'][0, :]
1475 inParams[
'COMPRETRIEVEDLNPWV'][:] = parCat[
'compRetrievedLnPwv'][0, :]
1476 inParams[
'COMPRETRIEVEDLNPWVRAW'][:] = parCat[
'compRetrievedLnPwvRaw'][0, :]
1477 inParams[
'COMPRETRIEVEDLNPWVFLAG'][:] = parCat[
'compRetrievedLnPwvFlag'][0, :]
1478 inParams[
'COMPRETRIEVEDTAUNIGHT'][:] = parCat[
'compRetrievedTauNight'][0, :]
1480 inSuperStar = np.zeros(parCat[
'superstarSize'][0, :], dtype=
'f8')
1481 inSuperStar[:, :, :, :] = parCat[
'superstar'][0, :].reshape(inSuperStar.shape)
1483 return (inParInfo, inParams, inSuperStar)
1485 def _makeFgcmOutputDatasets(self, fgcmFitCycle):
1487 Persist FGCM datasets through the butler.
1491 fgcmFitCycle: `lsst.fgcm.FgcmFitCycle`
1492 Fgcm Fit cycle object
1494 fgcmDatasetDict = {}
1497 parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays()
1502 lutFilterNameString = comma.join([n.decode(
'utf-8')
1503 for n
in parInfo[
'LUTFILTERNAMES'][0]])
1504 fitBandString = comma.join([n.decode(
'utf-8')
1505 for n
in parInfo[
'FITBANDS'][0]])
1507 parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat,
1508 lutFilterNameString, fitBandString)
1509 parCat = self._makeParCatalog(parSchema, parInfo, pars,
1510 fgcmFitCycle.fgcmPars.parSuperStarFlat,
1511 lutFilterNameString, fitBandString)
1513 fgcmDatasetDict[
'fgcmFitParameters'] = parCat
1518 flagStarSchema = self._makeFlagStarSchema()
1519 flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices()
1520 flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct)
1522 fgcmDatasetDict[
'fgcmFlaggedStars'] = flagStarCat
1525 if self.outputZeropoints:
1526 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct[
'FGCM_FZPT_SSTAR_CHEB'].shape[1]
1527 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct[
'FGCM_FZPT_CHEB'].shape[1]
1530 zptCat =
makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
1532 fgcmDatasetDict[
'fgcmZeropoints'] = zptCat
1537 atmCat =
makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
1539 fgcmDatasetDict[
'fgcmAtmosphereParameters'] = atmCat
1542 if self.outputStandards:
1543 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
1545 stdCat =
makeStdCat(stdSchema, stdStruct, goodBands)
1547 fgcmDatasetDict[
'fgcmStandardStars'] = stdCat
1549 return fgcmDatasetDict
1551 def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
1552 lutFilterNameString, fitBandString):
1554 Make the parameter persistence schema
1558 parInfo: `numpy.ndarray`
1559 Parameter information returned by fgcm
1560 pars: `numpy.ndarray`
1561 Parameter values returned by fgcm
1562 parSuperStarFlat: `numpy.array`
1563 Superstar flat values returned by fgcm
1564 lutFilterNameString: `str`
1565 Combined string of all the lutFilterNames
1566 fitBandString: `str`
1567 Combined string of all the fitBands
1571 parSchema: `afwTable.schema`
1577 parSchema.addField(
'nCcd', type=np.int32, doc=
'Number of CCDs')
1578 parSchema.addField(
'lutFilterNames', type=str, doc=
'LUT Filter names in parameter file',
1579 size=len(lutFilterNameString))
1580 parSchema.addField(
'fitBands', type=str, doc=
'Bands that were fit',
1581 size=len(fitBandString))
1582 parSchema.addField(
'lnTauUnit', type=np.float64, doc=
'Step units for ln(AOD)')
1583 parSchema.addField(
'lnTauSlopeUnit', type=np.float64,
1584 doc=
'Step units for ln(AOD) slope')
1585 parSchema.addField(
'alphaUnit', type=np.float64, doc=
'Step units for alpha')
1586 parSchema.addField(
'lnPwvUnit', type=np.float64, doc=
'Step units for ln(pwv)')
1587 parSchema.addField(
'lnPwvSlopeUnit', type=np.float64,
1588 doc=
'Step units for ln(pwv) slope')
1589 parSchema.addField(
'lnPwvQuadraticUnit', type=np.float64,
1590 doc=
'Step units for ln(pwv) quadratic term')
1591 parSchema.addField(
'lnPwvGlobalUnit', type=np.float64,
1592 doc=
'Step units for global ln(pwv) parameters')
1593 parSchema.addField(
'o3Unit', type=np.float64, doc=
'Step units for O3')
1594 parSchema.addField(
'qeSysUnit', type=np.float64, doc=
'Step units for mirror gray')
1595 parSchema.addField(
'filterOffsetUnit', type=np.float64, doc=
'Step units for filter offset')
1596 parSchema.addField(
'hasExternalPwv', type=np.int32, doc=
'Parameters fit using external pwv')
1597 parSchema.addField(
'hasExternalTau', type=np.int32, doc=
'Parameters fit using external tau')
1600 parSchema.addField(
'parAlpha', type=
'ArrayD', doc=
'Alpha parameter vector',
1601 size=pars[
'PARALPHA'].size)
1602 parSchema.addField(
'parO3', type=
'ArrayD', doc=
'O3 parameter vector',
1603 size=pars[
'PARO3'].size)
1604 parSchema.addField(
'parLnTauIntercept', type=
'ArrayD',
1605 doc=
'ln(Tau) intercept parameter vector',
1606 size=pars[
'PARLNTAUINTERCEPT'].size)
1607 parSchema.addField(
'parLnTauSlope', type=
'ArrayD',
1608 doc=
'ln(Tau) slope parameter vector',
1609 size=pars[
'PARLNTAUSLOPE'].size)
1610 parSchema.addField(
'parLnPwvIntercept', type=
'ArrayD', doc=
'ln(pwv) intercept parameter vector',
1611 size=pars[
'PARLNPWVINTERCEPT'].size)
1612 parSchema.addField(
'parLnPwvSlope', type=
'ArrayD', doc=
'ln(pwv) slope parameter vector',
1613 size=pars[
'PARLNPWVSLOPE'].size)
1614 parSchema.addField(
'parLnPwvQuadratic', type=
'ArrayD', doc=
'ln(pwv) quadratic parameter vector',
1615 size=pars[
'PARLNPWVQUADRATIC'].size)
1616 parSchema.addField(
'parQeSysIntercept', type=
'ArrayD', doc=
'Mirror gray intercept parameter vector',
1617 size=pars[
'PARQESYSINTERCEPT'].size)
1618 parSchema.addField(
'compQeSysSlope', type=
'ArrayD', doc=
'Mirror gray slope parameter vector',
1619 size=pars[0][
'COMPQESYSSLOPE'].size)
1620 parSchema.addField(
'parFilterOffset', type=
'ArrayD', doc=
'Filter offset parameter vector',
1621 size=pars[
'PARFILTEROFFSET'].size)
1622 parSchema.addField(
'parFilterOffsetFitFlag', type=
'ArrayI', doc=
'Filter offset parameter fit flag',
1623 size=pars[
'PARFILTEROFFSETFITFLAG'].size)
1624 parSchema.addField(
'parRetrievedLnPwvScale', type=np.float64,
1625 doc=
'Global scale for retrieved ln(pwv)')
1626 parSchema.addField(
'parRetrievedLnPwvOffset', type=np.float64,
1627 doc=
'Global offset for retrieved ln(pwv)')
1628 parSchema.addField(
'parRetrievedLnPwvNightlyOffset', type=
'ArrayD',
1629 doc=
'Nightly offset for retrieved ln(pwv)',
1630 size=pars[
'PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size)
1631 parSchema.addField(
'compAbsThroughput', type=
'ArrayD',
1632 doc=
'Absolute throughput (relative to transmission curves)',
1633 size=pars[
'COMPABSTHROUGHPUT'].size)
1634 parSchema.addField(
'compRefOffset', type=
'ArrayD',
1635 doc=
'Offset between reference stars and calibrated stars',
1636 size=pars[
'COMPREFOFFSET'].size)
1637 parSchema.addField(
'compRefSigma', type=
'ArrayD',
1638 doc=
'Width of reference star/calibrated star distribution',
1639 size=pars[
'COMPREFSIGMA'].size)
1640 parSchema.addField(
'compMirrorChromaticity', type=
'ArrayD',
1641 doc=
'Computed mirror chromaticity terms',
1642 size=pars[
'COMPMIRRORCHROMATICITY'].size)
1643 parSchema.addField(
'mirrorChromaticityPivot', type=
'ArrayD',
1644 doc=
'Mirror chromaticity pivot mjd',
1645 size=pars[
'MIRRORCHROMATICITYPIVOT'].size)
1646 parSchema.addField(
'compMedianSedSlope', type=
'ArrayD',
1647 doc=
'Computed median SED slope (per band)',
1648 size=pars[
'COMPMEDIANSEDSLOPE'].size)
1649 parSchema.addField(
'compAperCorrPivot', type=
'ArrayD', doc=
'Aperture correction pivot',
1650 size=pars[
'COMPAPERCORRPIVOT'].size)
1651 parSchema.addField(
'compAperCorrSlope', type=
'ArrayD', doc=
'Aperture correction slope',
1652 size=pars[
'COMPAPERCORRSLOPE'].size)
1653 parSchema.addField(
'compAperCorrSlopeErr', type=
'ArrayD', doc=
'Aperture correction slope error',
1654 size=pars[
'COMPAPERCORRSLOPEERR'].size)
1655 parSchema.addField(
'compAperCorrRange', type=
'ArrayD', doc=
'Aperture correction range',
1656 size=pars[
'COMPAPERCORRRANGE'].size)
1657 parSchema.addField(
'compModelErrExptimePivot', type=
'ArrayD', doc=
'Model error exptime pivot',
1658 size=pars[
'COMPMODELERREXPTIMEPIVOT'].size)
1659 parSchema.addField(
'compModelErrFwhmPivot', type=
'ArrayD', doc=
'Model error fwhm pivot',
1660 size=pars[
'COMPMODELERRFWHMPIVOT'].size)
1661 parSchema.addField(
'compModelErrSkyPivot', type=
'ArrayD', doc=
'Model error sky pivot',
1662 size=pars[
'COMPMODELERRSKYPIVOT'].size)
1663 parSchema.addField(
'compModelErrPars', type=
'ArrayD', doc=
'Model error parameters',
1664 size=pars[
'COMPMODELERRPARS'].size)
1665 parSchema.addField(
'compExpGray', type=
'ArrayD', doc=
'Computed exposure gray',
1666 size=pars[
'COMPEXPGRAY'].size)
1667 parSchema.addField(
'compVarGray', type=
'ArrayD', doc=
'Computed exposure variance',
1668 size=pars[
'COMPVARGRAY'].size)
1669 parSchema.addField(
'compExpDeltaMagBkg', type=
'ArrayD',
1670 doc=
'Computed exposure offset due to background',
1671 size=pars[
'COMPEXPDELTAMAGBKG'].size)
1672 parSchema.addField(
'compNGoodStarPerExp', type=
'ArrayI',
1673 doc=
'Computed number of good stars per exposure',
1674 size=pars[
'COMPNGOODSTARPEREXP'].size)
1675 parSchema.addField(
'compSigFgcm', type=
'ArrayD', doc=
'Computed sigma_fgcm (intrinsic repeatability)',
1676 size=pars[
'COMPSIGFGCM'].size)
1677 parSchema.addField(
'compSigmaCal', type=
'ArrayD', doc=
'Computed sigma_cal (systematic error floor)',
1678 size=pars[
'COMPSIGMACAL'].size)
1679 parSchema.addField(
'compRetrievedLnPwv', type=
'ArrayD', doc=
'Retrieved ln(pwv) (smoothed)',
1680 size=pars[
'COMPRETRIEVEDLNPWV'].size)
1681 parSchema.addField(
'compRetrievedLnPwvRaw', type=
'ArrayD', doc=
'Retrieved ln(pwv) (raw)',
1682 size=pars[
'COMPRETRIEVEDLNPWVRAW'].size)
1683 parSchema.addField(
'compRetrievedLnPwvFlag', type=
'ArrayI', doc=
'Retrieved ln(pwv) Flag',
1684 size=pars[
'COMPRETRIEVEDLNPWVFLAG'].size)
1685 parSchema.addField(
'compRetrievedTauNight', type=
'ArrayD', doc=
'Retrieved tau (per night)',
1686 size=pars[
'COMPRETRIEVEDTAUNIGHT'].size)
1688 parSchema.addField(
'superstarSize', type=
'ArrayI', doc=
'Superstar matrix size',
1690 parSchema.addField(
'superstar', type=
'ArrayD', doc=
'Superstar matrix (flattened)',
1691 size=parSuperStarFlat.size)
1695 def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
1696 lutFilterNameString, fitBandString):
1698 Make the FGCM parameter catalog for persistence
1702 parSchema: `lsst.afw.table.Schema`
1703 Parameter catalog schema
1704 pars: `numpy.ndarray`
1705 FGCM parameters to put into parCat
1706 parSuperStarFlat: `numpy.array`
1707 FGCM superstar flat array to put into parCat
1708 lutFilterNameString: `str`
1709 Combined string of all the lutFilterNames
1710 fitBandString: `str`
1711 Combined string of all the fitBands
1715 parCat: `afwTable.BasicCatalog`
1716 Atmosphere and instrumental model parameter catalog for persistence
1724 rec = parCat.addNew()
1727 rec[
'nCcd'] = parInfo[
'NCCD']
1728 rec[
'lutFilterNames'] = lutFilterNameString
1729 rec[
'fitBands'] = fitBandString
1731 rec[
'hasExternalPwv'] = 0
1732 rec[
'hasExternalTau'] = 0
1736 scalarNames = [
'parRetrievedLnPwvScale',
'parRetrievedLnPwvOffset']
1738 arrNames = [
'parAlpha',
'parO3',
'parLnTauIntercept',
'parLnTauSlope',
1739 'parLnPwvIntercept',
'parLnPwvSlope',
'parLnPwvQuadratic',
1740 'parQeSysIntercept',
'compQeSysSlope',
1741 'parRetrievedLnPwvNightlyOffset',
'compAperCorrPivot',
1742 'parFilterOffset',
'parFilterOffsetFitFlag',
1743 'compAbsThroughput',
'compRefOffset',
'compRefSigma',
1744 'compMirrorChromaticity',
'mirrorChromaticityPivot',
1745 'compAperCorrSlope',
'compAperCorrSlopeErr',
'compAperCorrRange',
1746 'compModelErrExptimePivot',
'compModelErrFwhmPivot',
1747 'compModelErrSkyPivot',
'compModelErrPars',
1748 'compExpGray',
'compVarGray',
'compNGoodStarPerExp',
'compSigFgcm',
1749 'compSigmaCal',
'compExpDeltaMagBkg',
'compMedianSedSlope',
1750 'compRetrievedLnPwv',
'compRetrievedLnPwvRaw',
'compRetrievedLnPwvFlag',
1751 'compRetrievedTauNight']
1753 for scalarName
in scalarNames:
1754 rec[scalarName] = pars[scalarName.upper()]
1756 for arrName
in arrNames:
1757 rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:]
1760 rec[
'superstarSize'][:] = parSuperStarFlat.shape
1761 rec[
'superstar'][:] = parSuperStarFlat.flatten()
1765 def _makeFlagStarSchema(self):
1767 Make the flagged-stars schema
1771 flagStarSchema: `lsst.afw.table.Schema`
1776 flagStarSchema.addField(
'objId', type=np.int32, doc=
'FGCM object id')
1777 flagStarSchema.addField(
'objFlag', type=np.int32, doc=
'FGCM object flag')
1779 return flagStarSchema
1781 def _makeFlagStarCat(self, flagStarSchema, flagStarStruct):
1783 Make the flagged star catalog for persistence
1787 flagStarSchema: `lsst.afw.table.Schema`
1789 flagStarStruct: `numpy.ndarray`
1790 Flagged star structure from fgcm
1794 flagStarCat: `lsst.afw.table.BaseCatalog`
1795 Flagged star catalog for persistence
1799 flagStarCat.resize(flagStarStruct.size)
1801 flagStarCat[
'objId'][:] = flagStarStruct[
'OBJID']
1802 flagStarCat[
'objFlag'][:] = flagStarStruct[
'OBJFLAG']
Defines the fields and offsets for a table.
def extractReferenceMags(refStars, bands, filterMap)
def makeStdSchema(nBands)
def makeAtmCat(atmSchema, atmStruct)
def makeConfigDict(config, log, camera, maxIter, resetFitParameters, outputZeropoints, lutFilterNames, tract=None)
def translateFgcmLut(lutCat, physicalFilterMap)
def makeZptCat(zptSchema, zpStruct)
def makeStdCat(stdSchema, stdStruct, goodBands)
def makeZptSchema(superStarChebyshevSize, zptChebyshevSize)
def translateVisitCatalog(visitCat)
def run(self, coaddExposures, bbox, wcs)