Coverage for python/lsst/fgcmcal/fgcmFitCycle.py: 17%
655 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-10 02:27 -0700
1# See COPYRIGHT file at the top of the source tree.
2#
3# This file is part of fgcmcal.
4#
5# Developed for the LSST Data Management System.
6# This product includes software developed by the LSST Project
7# (https://www.lsst.org).
8# See the COPYRIGHT file at the top-level directory of this distribution
9# for details of code ownership.
10#
11# This program is free software: you can redistribute it and/or modify
12# it under the terms of the GNU General Public License as published by
13# the Free Software Foundation, either version 3 of the License, or
14# (at your option) any later version.
15#
16# This program is distributed in the hope that it will be useful,
17# but WITHOUT ANY WARRANTY; without even the implied warranty of
18# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19# GNU General Public License for more details.
20#
21# You should have received a copy of the GNU General Public License
22# along with this program. If not, see <https://www.gnu.org/licenses/>.
23"""Perform a single fit cycle of FGCM.
25This task runs a single "fit cycle" of fgcm. Prior to running this task
26one must run both fgcmMakeLut (to construct the atmosphere and instrumental
27look-up-table) and fgcmBuildStars (to extract visits and star observations
28for the global fit).
30The 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
32be inspected to set parameters for outlier rejection on the following
33cycle. Please see the fgcmcal Cookbook for details.
34"""
36import copy
38import numpy as np
40import lsst.pex.config as pexConfig
41import lsst.pipe.base as pipeBase
42from lsst.pipe.base import connectionTypes
43import lsst.afw.table as afwTable
45from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
46from .utilities import extractReferenceMags
47from .utilities import makeZptSchema, makeZptCat
48from .utilities import makeAtmSchema, makeAtmCat, makeStdSchema, makeStdCat
49from .sedterms import SedboundarytermDict, SedtermDict
50from .focalPlaneProjector import FocalPlaneProjector
52import fgcm
54__all__ = ['FgcmFitCycleConfig', 'FgcmFitCycleTask']
56MULTIPLE_CYCLES_MAX = 10
59class FgcmFitCycleConnections(pipeBase.PipelineTaskConnections,
60 dimensions=("instrument",),
61 defaultTemplates={"previousCycleNumber": "-1",
62 "cycleNumber": "0"}):
63 camera = connectionTypes.PrerequisiteInput(
64 doc="Camera instrument",
65 name="camera",
66 storageClass="Camera",
67 dimensions=("instrument",),
68 isCalibration=True,
69 )
71 fgcmLookUpTable = connectionTypes.PrerequisiteInput(
72 doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
73 "chromatic corrections."),
74 name="fgcmLookUpTable",
75 storageClass="Catalog",
76 dimensions=("instrument",),
77 deferLoad=True,
78 )
80 fgcmVisitCatalog = connectionTypes.Input(
81 doc="Catalog of visit information for fgcm",
82 name="fgcmVisitCatalog",
83 storageClass="Catalog",
84 dimensions=("instrument",),
85 deferLoad=True,
86 )
88 fgcmStarObservationsParquet = connectionTypes.Input(
89 doc=("Catalog of star observations for fgcm, in parquet format. "
90 "Used if useParquetCatalogFormat is True."),
91 name="fgcm_star_observations",
92 storageClass="ArrowAstropy",
93 dimensions=("instrument",),
94 deferLoad=True,
95 )
97 fgcmStarIdsParquet = connectionTypes.Input(
98 doc=("Catalog of fgcm calibration star IDs, in parquet format. "
99 "Used if useParquetCatalogFormat is True."),
100 name="fgcm_star_ids",
101 storageClass="ArrowAstropy",
102 dimensions=("instrument",),
103 deferLoad=True,
104 )
106 fgcmReferenceStarsParquet = connectionTypes.Input(
107 doc=("Catalog of fgcm-matched reference stars, in parquet format. "
108 "Used if useParquetCatalogFormat is True."),
109 name="fgcm_reference_stars",
110 storageClass="ArrowAstropy",
111 dimensions=("instrument",),
112 deferLoad=True,
113 )
115 fgcmStarObservations = connectionTypes.Input(
116 doc=("Catalog of star observations for fgcm; old format. "
117 "Used if useParquetCatalogFormat is False."),
118 name="fgcmStarObservations",
119 storageClass="Catalog",
120 dimensions=("instrument",),
121 deferLoad=True,
122 )
124 fgcmStarIds = connectionTypes.Input(
125 doc=("Catalog of fgcm calibration star IDs. "
126 "Used if useParquetCatalogFormat is False."),
127 name="fgcmStarIds",
128 storageClass="Catalog",
129 dimensions=("instrument",),
130 deferLoad=True,
131 )
133 fgcmStarIndices = connectionTypes.Input(
134 doc=("Catalog of fgcm calibration star indices; old format."
135 "Used if useParquetCatalogFormat is False."),
136 name="fgcmStarIndices",
137 storageClass="Catalog",
138 dimensions=("instrument",),
139 deferLoad=True,
140 )
142 fgcmReferenceStars = connectionTypes.Input(
143 doc=("Catalog of fgcm-matched reference stars; old format."
144 "Used if useParquetCatalogFormat is False."),
145 name="fgcmReferenceStars",
146 storageClass="Catalog",
147 dimensions=("instrument",),
148 deferLoad=True,
149 )
151 def __init__(self, *, config=None):
152 super().__init__(config=config)
154 # The connection parameters ``cycleNumber`` and ``previousCycleNumber``
155 # always comes in as strings, and doing a round-trip via
156 # integer-to-string ensures that they actually are integers.
157 if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber:
158 raise ValueError("cycleNumber must be of integer format")
159 if str(int(config.connections.previousCycleNumber)) != config.connections.previousCycleNumber:
160 raise ValueError("previousCycleNumber must be of integer format")
161 # We additionally confirm that there is a delta of 1 between them.
162 # Note that this math cannot be done in the parameters because they are
163 # strings.
164 if int(config.connections.previousCycleNumber) != (int(config.connections.cycleNumber) - 1):
165 raise ValueError("previousCycleNumber must be 1 less than cycleNumber")
167 instDims = ("instrument",)
168 bandDims = ("instrument", "band")
169 filterDims = ("instrument", "physical_filter")
170 filterDetectorDims = ("instrument", "physical_filter", "detector")
172 inputAndOutputConnections = [
173 ("FitParameters", "Catalog", "Catalog of fgcm fit parameters.", instDims),
174 ("FlaggedStars", "Catalog", "Catalog of flagged stars for fgcm calibration.", instDims),
175 ]
176 multicycleOutputConnections = [
177 ("OutputConfig", "Config", "Configuration for next fgcm fit cycle.", instDims),
178 ]
179 optionalZpOutputConnections = [
180 ("Zeropoints", "Catalog", "Catalog of fgcm zeropoint data.", instDims),
181 ("AtmosphereParameters", "Catalog", "Catalog of atmospheric fit parameters.", instDims),
182 ]
183 optionalStarOutputConnections = [
184 ("StandardStars", "SimpleCatalog", "Catalog of standard star magnitudes.", instDims),
185 ]
187 epochs = [f"epoch{i}" for i in range(len(config.epochMjds))]
189 # All names will be preceeded with ``fgcm_CycleN_``.
190 plotConnections = [
191 ("Zeropoints_Plot", "Plot", "Plot of fgcm zeropoints.", instDims),
192 ("ExpgrayDeep_Plot",
193 "Plot",
194 "Plot of gray term per exposure per time for deep fields.",
195 instDims),
196 ("NightlyAlpha_Plot", "Plot", "Plot of nightly AOD alpha term.", instDims),
197 ("NightlyTau_Plot", "Plot", "Plot of nightly aerosol optical depth (tau).", instDims),
198 ("NightlyPwv_Plot", "Plot", "Plot of nightly water vapor.", instDims),
199 ("NightlyO3_Plot", "Plot", "Plot of nightly ozone.", instDims),
200 ("FilterOffsets_Plot", "Plot", "Plot of in-band filter offsets.", instDims),
201 ("AbsThroughputs_Plot", "Plot", "Plot of absolute throughput fractions.", instDims),
202 ("QESysWashesInitial_Plot", "Plot", "Plot of initial system QE with mirror washes.", instDims),
203 ("QESysWashesFinal_Plot", "Plot", "Plot of final system QE with mirror washes.", instDims),
204 ("RpwvVsRpwvInput_Plot",
205 "Plot",
206 "Plot of change in per-visit ``retrieved`` PWV from previous fit cycle.",
207 instDims),
208 ("RpwvVsRpwvSmooth_Plot",
209 "Plot",
210 "Plot of per-visit ``retrieved`` PWV vs. smoothed PWV.",
211 instDims),
212 ("ModelPwvVsRpwv_Plot",
213 "Plot",
214 "Plot of model PWV vs. per-visit ``retrieved`` PWV.",
215 instDims),
216 ("ChisqFit_Plot",
217 "Plot",
218 "Plot of chisq as a function of iteration.",
219 instDims),
220 ]
222 plotConnections.extend(
223 [
224 ("Apercorr_Plot", "Plot", "Plot of fgcm aperture corrections.", bandDims),
225 ("EpsilonGlobal_Plot",
226 "Plot",
227 "Plot of global background over/undersubtraction.",
228 bandDims),
229 ("EpsilonMap_Plot",
230 "Plot",
231 "Map of spatially varying background over/undersubtraction.",
232 bandDims),
233 ("ExpgrayInitial_Plot",
234 "Plot",
235 "Histogram of initial gray term per exposure.",
236 bandDims),
237 ("CompareRedblueExpgray_Plot",
238 "Plot",
239 "Plot of red/blue split gray term per exposure",
240 bandDims),
241 ("Expgray_Plot", "Plot", "Histogram of gray term per exposure.", bandDims),
242 ("ExpgrayAirmass_Plot",
243 "Plot",
244 "Plot of exposure gray term as a function of airmass.",
245 bandDims),
246 ("ExpgrayCompareMjdRedblue_Plot",
247 "Plot",
248 "Plot of red/blue split gray term per exposure as a function of time.",
249 bandDims),
250 ("ExpgrayUT_Plot",
251 "Plot",
252 "Plot of grey term per exposure as a function of time of night.",
253 bandDims),
254 ("ExpgrayCompareBands_Plot",
255 "Plot",
256 "Plot of gray term per exposure between bands nearby in time.",
257 bandDims),
258 ("ExpgrayReference_Plot",
259 "Plot",
260 "Histogram of gray term per exposure compared to reference mags.",
261 bandDims),
262 ("QESysRefstarsStdInitial_Plot",
263 "Plot",
264 "Plot of reference mag - calibrated (standard) mag vs. time (before fit).",
265 bandDims),
266 ("QESysRefstarsStdFinal_Plot",
267 "Plot",
268 "Plot of reference mag - calibrated (standard) mag vs. time (after fit).",
269 bandDims),
270 ("QESysRefstarsObsInitial_Plot",
271 "Plot",
272 "Plot of reference mag - observed (instrumental) mag vs. time (before fit).",
273 bandDims),
274 ("QESysRefstarsObsFinal_Plot",
275 "Plot",
276 "Plot of reference mag - observed (instrumental) mag vs. time (after fit).",
277 bandDims),
278 ("ModelMagerrInitial_Plot",
279 "Plot",
280 "Plots for magnitude error model, initial estimate.",
281 bandDims),
282 ("ModelMagerrPostfit_Plot",
283 "Plot",
284 "Plots for magnitude error model, after fitting.",
285 bandDims),
286 ("SigmaFgcmAllStars_Plot",
287 "Plot",
288 "Histograms for intrinsic scatter for all bright stars.",
289 bandDims),
290 ("SigmaFgcmReservedStars_Plot",
291 "Plot",
292 "Histograms for intrinsic scatter for reserved bright stars.",
293 bandDims),
294 ("SigmaFgcmReservedStarsCrunched_Plot",
295 "Plot",
296 "Histograms for intrinsic scatter for reserved bright stars (after gray correction).",
297 bandDims),
298 ("SigmaFgcmPullsAllStars_Plot",
299 "Plot",
300 "Histograms for pulls for all bright stars.",
301 bandDims),
302 ("SigmaFgcmPullsReservedStars_Plot",
303 "Plot",
304 "Histograms for pulls for reserved bright stars.",
305 bandDims),
306 ("SigmaFgcmPullsReservedStarsCrunched_Plot",
307 "Plot",
308 "Histograms for pulls for reserved bright stars (after gray correction).",
309 bandDims),
310 ("SigmaCal_Plot",
311 "Plot",
312 "Plot showing scatter as a function of systematic error floor.",
313 bandDims),
314 ("SigmaRef_Plot",
315 "Plot",
316 "Histograms of scatter with respect to reference stars.",
317 bandDims),
318 ("RefResidVsColorAll_Plot",
319 "Plot",
320 "Plot of reference star residuals vs. color (all stars).",
321 bandDims),
322 ("RefResidVsColorCut_Plot",
323 "Plot",
324 "Plot of reference star residuals vs. color (reference star color cuts).",
325 bandDims),
326 ]
327 )
329 plotConnections.extend(
330 [
331 ("I1R1_Plot", "Plot", "Plot of fgcm R1 vs. I1.", filterDims),
332 ("I1_Plot", "Plot", "Focal plane map of fgcm I1.", filterDims),
333 ("R1_Plot", "Plot", "Focal plane map of fgcm R1.", filterDims),
334 ("R1mI1Matchscale_Plot", "Plot", "Focal plane map of fgcm R1 - I1.", filterDims),
335 ("R1mI1_Plot", "Plot", "Focal plane map of fgcm R1 - I1 (rescaled).", filterDims),
336 ("R1mI1_vs_mjd_Plot", "Plot", "R1 - I1 residuals vs. mjd.", filterDims),
337 ("CompareRedblueMirrorchrom_Plot",
338 "Plot",
339 "Comparison of mirror chromaticity residuals for red/blue stars.",
340 filterDims),
341 ("CcdChromaticity_Plot",
342 "Plot",
343 "Focal plane map of fgcm ccd chromaticity.",
344 filterDims),
345 ("EpsilonDetector_Plot",
346 "Plot",
347 "Focal plane map of background over/undersubtraction.",
348 filterDims),
349 ("EpsilonDetectorMatchscale_Plot",
350 "Plot",
351 "Focal plane map of background over/undersubtraction.",
352 filterDims),
353 ]
354 )
355 for epoch in epochs:
356 plotConnections.extend(
357 [
358 (
359 f"Superstar_{epoch}_Plot",
360 "Plot",
361 "Plot of illumination Correction.",
362 filterDims,
363 ),
364 ],
365 )
366 if self.config.superStarPlotCcdResiduals:
367 plotConnections.extend(
368 [
369 (
370 f"SuperstarResidual_{epoch}_Plot",
371 "Plot",
372 "Binned illumination correction residuals.",
373 filterDetectorDims,
374 ),
375 (
376 f"SuperstarResidualStd_{epoch}_Plot",
377 "Plot",
378 "Binned illumination correction residual stdev.",
379 filterDetectorDims,
380 ),
381 ],
382 )
384 if config.doMultipleCycles:
385 # Multiple cycle run.
387 # All but the final cycle get appended here.
388 for cycle in range(config.multipleCyclesFinalCycleNumber):
389 outputConnections = copy.copy(inputAndOutputConnections)
390 outputConnections.extend(multicycleOutputConnections)
391 if config.outputZeropointsBeforeFinalCycle:
392 outputConnections.extend(optionalZpOutputConnections)
393 if config.outputStandardsBeforeFinalCycle:
394 outputConnections.extend(optionalStarOutputConnections)
396 if config.doPlots:
397 # We will make the plots if doPlots is True and either
398 # it is the next-to-last cycle or the
399 # doPlotsBeforeFinalCycles configuration is set which
400 # means we want plots for all cycles.
401 # The final cycle doPlots is configured below.
402 if cycle == (config.multipleCyclesFinalCycleNumber - 1) \
403 or config.doPlotsBeforeFinalCycles:
404 outputConnections.extend(plotConnections)
406 for (name, storageClass, doc, dims) in outputConnections:
407 connectionName = f"fgcm_Cycle{cycle}_{name}"
408 storageName = connectionName
409 outConnection = connectionTypes.Output(
410 name=storageName,
411 storageClass=storageClass,
412 doc=doc,
413 dimensions=dims,
414 multiple=(len(dims) > 1),
415 )
416 setattr(self, connectionName, outConnection)
418 # The final cycle has everything.
419 outputConnections = copy.copy(inputAndOutputConnections)
420 outputConnections.extend(multicycleOutputConnections)
421 outputConnections.extend(optionalZpOutputConnections)
422 outputConnections.extend(optionalStarOutputConnections)
423 if config.doPlots:
424 outputConnections.extend(plotConnections)
425 for (name, storageClass, doc, dims) in outputConnections:
426 connectionName = f"fgcm_Cycle{config.multipleCyclesFinalCycleNumber}_{name}"
427 storageName = connectionName
428 outConnection = connectionTypes.Output(
429 name=storageName,
430 storageClass=storageClass,
431 doc=doc,
432 dimensions=dims,
433 multiple=(len(dims) > 1),
434 )
435 setattr(self, connectionName, outConnection)
436 else:
437 # Single cycle run.
438 if config.cycleNumber > 0:
439 inputConnections = copy.copy(inputAndOutputConnections)
440 else:
441 inputConnections = []
442 outputConnections = copy.copy(inputAndOutputConnections)
443 # The following configurations are also useful for runs
444 # where fit cycles are run one-at-a-time since it is
445 # not typical to look at these outputs for intermediate
446 # steps.
447 if config.isFinalCycle or config.outputZeropointsBeforeFinalCycle:
448 outputConnections.extend(optionalZpOutputConnections)
449 if config.isFinalCycle or config.outputStandardsBeforeFinalCycle:
450 outputConnections.extend(optionalStarOutputConnections)
452 if config.doPlots:
453 outputConnections.extend(plotConnections)
455 for (name, storageClass, doc, dims) in inputConnections:
456 connectionName = f"fgcm{name}Input"
457 storageName = f"fgcm_Cycle{config.cycleNumber - 1}_{name}"
458 inConnection = connectionTypes.PrerequisiteInput(
459 name=storageName,
460 storageClass=storageClass,
461 doc=doc,
462 dimensions=dims,
463 )
464 setattr(self, connectionName, inConnection)
466 for (name, storageClass, doc, dims) in outputConnections:
467 connectionName = f"fgcm{name}"
468 storageName = f"fgcm_Cycle{config.cycleNumber}_{name}"
469 # Plots have unique names as well.
470 if storageClass == "Plot":
471 connectionName = storageName
472 outConnection = connectionTypes.Output(
473 name=storageName,
474 storageClass=storageClass,
475 doc=doc,
476 dimensions=dims,
477 multiple=(len(dims) > 1),
478 )
479 setattr(self, connectionName, outConnection)
481 if not config.doReferenceCalibration:
482 self.inputs.remove("fgcmReferenceStars")
483 self.inputs.remove("fgcmReferenceStarsParquet")
485 if config.useParquetCatalogFormat:
486 self.inputs.remove("fgcmStarObservations")
487 self.inputs.remove("fgcmStarIds")
488 self.inputs.remove("fgcmStarIndices")
489 if config.doReferenceCalibration:
490 self.inputs.remove("fgcmReferenceStars")
491 else:
492 self.inputs.remove("fgcmStarObservationsParquet")
493 self.inputs.remove("fgcmStarIdsParquet")
494 if config.doReferenceCalibration:
495 self.inputs.remove("fgcmReferenceStarsParquet")
498class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
499 pipelineConnections=FgcmFitCycleConnections):
500 """Config for FgcmFitCycle"""
502 doMultipleCycles = pexConfig.Field(
503 doc="Run multiple fit cycles in one task",
504 dtype=bool,
505 default=False,
506 )
507 useParquetCatalogFormat = pexConfig.Field(
508 doc="Use parquet catalog format?",
509 dtype=bool,
510 default=True,
511 )
512 multipleCyclesFinalCycleNumber = pexConfig.RangeField(
513 doc=("Final cycle number in multiple cycle mode. The initial cycle "
514 "is 0, with limited parameters fit. The next cycle is 1 with "
515 "full parameter fit. The final cycle is a clean-up with no "
516 "parameters fit. There will be a total of "
517 "(multipleCycleFinalCycleNumber + 1) cycles run, and the final "
518 "cycle number cannot be less than 2."),
519 dtype=int,
520 default=5,
521 min=2,
522 max=MULTIPLE_CYCLES_MAX,
523 inclusiveMax=True,
524 )
525 bands = pexConfig.ListField(
526 doc="Bands to run calibration",
527 dtype=str,
528 default=[],
529 )
530 fitBands = pexConfig.ListField(
531 doc=("Bands to use in atmospheric fit. The bands not listed here will have "
532 "the atmosphere constrained from the 'fitBands' on the same night. "
533 "Must be a subset of `config.bands`"),
534 dtype=str,
535 default=[],
536 )
537 requiredBands = pexConfig.ListField(
538 doc=("Bands that are required for a star to be considered a calibration star. "
539 "Must be a subset of `config.bands`"),
540 dtype=str,
541 default=[],
542 )
543 physicalFilterMap = pexConfig.DictField(
544 doc="Mapping from 'physicalFilter' to band.",
545 keytype=str,
546 itemtype=str,
547 default={},
548 )
549 doReferenceCalibration = pexConfig.Field(
550 doc="Use reference catalog as additional constraint on calibration",
551 dtype=bool,
552 default=True,
553 )
554 refStarSnMin = pexConfig.Field(
555 doc="Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.",
556 dtype=float,
557 default=50.0,
558 )
559 refStarOutlierNSig = pexConfig.Field(
560 doc=("Number of sigma compared to average mag for reference star to be considered an outlier. "
561 "Computed per-band, and if it is an outlier in any band it is rejected from fits."),
562 dtype=float,
563 default=4.0,
564 )
565 applyRefStarColorCuts = pexConfig.Field(
566 doc=("Apply color cuts defined in ``starColorCuts`` to reference stars? "
567 "These cuts are in addition to any cuts defined in ``refStarColorCuts``"),
568 dtype=bool,
569 default=True,
570 )
571 refStarMaxFracUse = pexConfig.Field(
572 doc=("Maximum fraction of reference stars to use in the fit. Remainder will "
573 "be used only for validation."),
574 dtype=float,
575 default=0.5,
576 )
577 useExposureReferenceOffset = pexConfig.Field(
578 doc=("Use per-exposure (visit) offsets between calibrated stars and reference stars "
579 "for final zeropoints? This may help uniformity for disjoint surveys."),
580 dtype=bool,
581 default=False,
582 )
583 nCore = pexConfig.Field(
584 doc="Number of cores to use",
585 dtype=int,
586 default=4,
587 deprecated="Number of cores is deprecated as a config, and will be removed after v27. "
588 "Please use ``pipetask run --cores-per-quantum`` instead.",
589 )
590 nStarPerRun = pexConfig.Field(
591 doc="Number of stars to run in each chunk",
592 dtype=int,
593 default=200000,
594 )
595 nExpPerRun = pexConfig.Field(
596 doc="Number of exposures to run in each chunk",
597 dtype=int,
598 default=1000,
599 )
600 reserveFraction = pexConfig.Field(
601 doc="Fraction of stars to reserve for testing",
602 dtype=float,
603 default=0.1,
604 )
605 freezeStdAtmosphere = pexConfig.Field(
606 doc="Freeze atmosphere parameters to standard (for testing)",
607 dtype=bool,
608 default=False,
609 )
610 precomputeSuperStarInitialCycle = pexConfig.Field(
611 doc="Precompute superstar flat for initial cycle",
612 dtype=bool,
613 default=False,
614 )
615 superStarSubCcdDict = pexConfig.DictField(
616 doc=("Per-band specification on whether to compute superstar flat on sub-ccd scale. "
617 "Must have one entry per band."),
618 keytype=str,
619 itemtype=bool,
620 default={},
621 )
622 superStarSubCcdChebyshevOrder = pexConfig.Field(
623 doc=("Order of the 2D chebyshev polynomials for sub-ccd superstar fit. "
624 "Global default is first-order polynomials, and should be overridden "
625 "on a camera-by-camera basis depending on the ISR."),
626 dtype=int,
627 default=1,
628 )
629 superStarSubCcdTriangular = pexConfig.Field(
630 doc=("Should the sub-ccd superstar chebyshev matrix be triangular to "
631 "suppress high-order cross terms?"),
632 dtype=bool,
633 default=False,
634 )
635 superStarSigmaClip = pexConfig.Field(
636 doc="Number of sigma to clip outliers when selecting for superstar flats",
637 dtype=float,
638 default=5.0,
639 )
640 superStarPlotCcdResiduals = pexConfig.Field(
641 doc="If plotting is enabled, should per-detector residuals be plotted? "
642 "This may produce a lot of output, and should be used only for "
643 "debugging purposes.",
644 dtype=bool,
645 default=False,
646 )
647 superStarForceZeroMean = pexConfig.Field(
648 doc="When computing the super-star flat, force the focal-plane mean to "
649 "zero (per band)? This should only be used when computing stand-alone "
650 "illumination corrections.",
651 dtype=bool,
652 default=False,
653 )
654 focalPlaneSigmaClip = pexConfig.Field(
655 doc="Number of sigma to clip outliers per focal-plane.",
656 dtype=float,
657 default=4.0,
658 )
659 ccdGraySubCcdDict = pexConfig.DictField(
660 doc=("Per-band specification on whether to compute achromatic per-ccd residual "
661 "('ccd gray') on a sub-ccd scale."),
662 keytype=str,
663 itemtype=bool,
664 default={},
665 )
666 ccdGraySubCcdChebyshevOrder = pexConfig.Field(
667 doc="Order of the 2D chebyshev polynomials for sub-ccd gray fit.",
668 dtype=int,
669 default=1,
670 )
671 ccdGraySubCcdTriangular = pexConfig.Field(
672 doc=("Should the sub-ccd gray chebyshev matrix be triangular to "
673 "suppress high-order cross terms?"),
674 dtype=bool,
675 default=True,
676 )
677 ccdGrayFocalPlaneDict = pexConfig.DictField(
678 doc=("Per-band specification on whether to compute focal-plane residual "
679 "('ccd gray') corrections."),
680 keytype=str,
681 itemtype=bool,
682 default={},
683 )
684 ccdGrayFocalPlaneFitMinCcd = pexConfig.Field(
685 doc=("Minimum number of 'good' CCDs required to perform focal-plane "
686 "gray corrections. If there are fewer good CCDs then the gray "
687 "correction is computed per-ccd."),
688 dtype=int,
689 default=1,
690 )
691 ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field(
692 doc="Order of the 2D chebyshev polynomials for focal plane fit.",
693 dtype=int,
694 default=3,
695 )
696 ccdGrayFocalPlaneMaxStars = pexConfig.Field(
697 doc="Maximum number of stars to use for focal plane fit. Required to keep "
698 "matrix memory usage from running away. If there are more stars than "
699 "this then they will be down-sampled.",
700 dtype=int,
701 default=50_000,
702 )
703 cycleNumber = pexConfig.Field(
704 doc=("FGCM fit cycle number. This is automatically incremented after each run "
705 "and stage of outlier rejection. See cookbook for details."),
706 dtype=int,
707 default=None,
708 )
709 isFinalCycle = pexConfig.Field(
710 doc=("Is this the final cycle of the fitting? Will automatically compute final "
711 "selection of stars and photometric exposures, and will output zeropoints "
712 "and standard stars for use in fgcmOutputProducts"),
713 dtype=bool,
714 default=False,
715 )
716 maxIterBeforeFinalCycle = pexConfig.Field(
717 doc=("Maximum fit iterations, prior to final cycle. The number of iterations "
718 "will always be 0 in the final cycle for cleanup and final selection."),
719 dtype=int,
720 default=50,
721 )
722 deltaMagBkgOffsetPercentile = pexConfig.Field(
723 doc=("Percentile brightest stars on a visit/ccd to use to compute net "
724 "offset from local background subtraction."),
725 dtype=float,
726 default=0.25,
727 )
728 deltaMagBkgPerCcd = pexConfig.Field(
729 doc=("Compute net offset from local background subtraction per-ccd? "
730 "Otherwise, use computation per visit."),
731 dtype=bool,
732 default=False,
733 )
734 utBoundary = pexConfig.Field(
735 doc="Boundary (in UTC) from day-to-day",
736 dtype=float,
737 default=None,
738 )
739 washMjds = pexConfig.ListField(
740 doc="Mirror wash MJDs",
741 dtype=float,
742 default=(0.0,),
743 )
744 epochMjds = pexConfig.ListField(
745 doc="Epoch boundaries in MJD",
746 dtype=float,
747 default=(0.0,),
748 )
749 minObsPerBand = pexConfig.Field(
750 doc="Minimum good observations per band",
751 dtype=int,
752 default=2,
753 )
754 # TODO: When DM-16511 is done, it will be possible to get the
755 # telescope latitude directly from the camera.
756 latitude = pexConfig.Field(
757 doc="Observatory latitude",
758 dtype=float,
759 default=None,
760 )
761 mirrorArea = pexConfig.Field(
762 doc="Mirror area (square meters) of telescope. If not set, will "
763 "try to estimate from camera.telescopeDiameter.",
764 dtype=float,
765 default=None,
766 optional=True,
767 )
768 cameraGain = pexConfig.Field(
769 doc="Average camera gain. If not set, will use the median of the "
770 "camera model/detector/amplifier gains.",
771 dtype=float,
772 default=None,
773 optional=True,
774 )
775 defaultCameraOrientation = pexConfig.Field(
776 doc="Default camera orientation for QA plots.",
777 dtype=float,
778 default=None,
779 )
780 brightObsGrayMax = pexConfig.Field(
781 doc="Maximum gray extinction to be considered bright observation",
782 dtype=float,
783 default=0.15,
784 )
785 minStarPerCcd = pexConfig.Field(
786 doc=("Minimum number of good stars per CCD to be used in calibration fit. "
787 "CCDs with fewer stars will have their calibration estimated from other "
788 "CCDs in the same visit, with zeropoint error increased accordingly."),
789 dtype=int,
790 default=5,
791 )
792 minCcdPerExp = pexConfig.Field(
793 doc=("Minimum number of good CCDs per exposure/visit to be used in calibration fit. "
794 "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."),
795 dtype=int,
796 default=5,
797 )
798 maxCcdGrayErr = pexConfig.Field(
799 doc="Maximum error on CCD gray offset to be considered photometric",
800 dtype=float,
801 default=0.05,
802 )
803 minStarPerExp = pexConfig.Field(
804 doc=("Minimum number of good stars per exposure/visit to be used in calibration fit. "
805 "Visits with fewer good stars will have CCD zeropoints estimated where possible."),
806 dtype=int,
807 default=600,
808 )
809 minExpPerNight = pexConfig.Field(
810 doc="Minimum number of good exposures/visits to consider a partly photometric night",
811 dtype=int,
812 default=10,
813 )
814 expGrayInitialCut = pexConfig.Field(
815 doc=("Maximum exposure/visit gray value for initial selection of possible photometric "
816 "observations."),
817 dtype=float,
818 default=-0.25,
819 )
820 expFwhmCutDict = pexConfig.DictField(
821 doc=("Per-band specification on maximum exposure FWHM (arcseconds) that will "
822 "be considered for the model fit. Exposures with median FWHM larger "
823 "than this threshold will get zeropoints based on matching to good "
824 "stars."),
825 keytype=str,
826 itemtype=float,
827 default={},
828 )
829 expGrayPhotometricCutDict = pexConfig.DictField(
830 doc=("Per-band specification on maximum (negative) achromatic exposure residual "
831 "('gray term') for a visit to be considered photometric. Must have one "
832 "entry per band. Broad-band filters should be -0.05."),
833 keytype=str,
834 itemtype=float,
835 default={},
836 )
837 expGrayHighCutDict = pexConfig.DictField(
838 doc=("Per-band specification on maximum (positive) achromatic exposure residual "
839 "('gray term') for a visit to be considered photometric. Must have one "
840 "entry per band. Broad-band filters should be 0.2."),
841 keytype=str,
842 itemtype=float,
843 default={},
844 )
845 expGrayRecoverCut = pexConfig.Field(
846 doc=("Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. "
847 "Visits with more gray extinction will only get CCD zeropoints if there are "
848 "sufficient star observations (minStarPerCcd) on that CCD."),
849 dtype=float,
850 default=-1.0,
851 )
852 expVarGrayPhotometricCutDict = pexConfig.DictField(
853 doc=("Per-band specification on maximum exposure variance to be considered possibly "
854 "photometric. Must have one entry per band. Broad-band filters should be "
855 "0.0005."),
856 keytype=str,
857 itemtype=float,
858 default={},
859 )
860 expGrayErrRecoverCut = pexConfig.Field(
861 doc=("Maximum exposure gray error to be able to recover bad ccds via interpolation. "
862 "Visits with more gray variance will only get CCD zeropoints if there are "
863 "sufficient star observations (minStarPerCcd) on that CCD."),
864 dtype=float,
865 default=0.05,
866 )
867 aperCorrPerCcd = pexConfig.Field(
868 doc="Use aperture corrections per-ccd (detector) instead of per-visit?",
869 dtype=bool,
870 default=False,
871 )
872 aperCorrFitNBins = pexConfig.Field(
873 doc=("Number of aperture bins used in aperture correction fit. When set to 0"
874 "no fit will be performed, and the config.aperCorrInputSlopes will be "
875 "used if available."),
876 dtype=int,
877 default=10,
878 )
879 aperCorrInputSlopeDict = pexConfig.DictField(
880 doc=("Per-band specification of aperture correction input slope parameters. These "
881 "are used on the first fit iteration, and aperture correction parameters will "
882 "be updated from the data if config.aperCorrFitNBins > 0. It is recommended "
883 "to set this when there is insufficient data to fit the parameters (e.g. "
884 "tract mode)."),
885 keytype=str,
886 itemtype=float,
887 default={},
888 )
889 sedboundaryterms = pexConfig.ConfigField(
890 doc="Mapping from bands to SED boundary term names used is sedterms.",
891 dtype=SedboundarytermDict,
892 )
893 sedterms = pexConfig.ConfigField(
894 doc="Mapping from terms to bands for fgcm linear SED approximations.",
895 dtype=SedtermDict,
896 )
897 sigFgcmMaxErr = pexConfig.Field(
898 doc="Maximum mag error for fitting sigma_FGCM",
899 dtype=float,
900 default=0.01,
901 )
902 sigFgcmMaxEGrayDict = pexConfig.DictField(
903 doc=("Per-band specification for maximum (absolute) achromatic residual (gray value) "
904 "for observations in sigma_fgcm (raw repeatability). Broad-band filters "
905 "should be 0.05."),
906 keytype=str,
907 itemtype=float,
908 default={},
909 )
910 ccdGrayMaxStarErr = pexConfig.Field(
911 doc=("Maximum error on a star observation to use in ccd gray (achromatic residual) "
912 "computation"),
913 dtype=float,
914 default=0.10,
915 )
916 approxThroughputDict = pexConfig.DictField(
917 doc=("Per-band specification of the approximate overall throughput at the start of "
918 "calibration observations. Must have one entry per band. Typically should "
919 "be 1.0."),
920 keytype=str,
921 itemtype=float,
922 default={},
923 )
924 sigmaCalRange = pexConfig.ListField(
925 doc="Allowed range for systematic error floor estimation",
926 dtype=float,
927 default=(0.001, 0.003),
928 )
929 sigmaCalFitPercentile = pexConfig.ListField(
930 doc="Magnitude percentile range to fit systematic error floor",
931 dtype=float,
932 default=(0.05, 0.15),
933 )
934 sigmaCalPlotPercentile = pexConfig.ListField(
935 doc="Magnitude percentile range to plot systematic error floor",
936 dtype=float,
937 default=(0.05, 0.95),
938 )
939 sigma0Phot = pexConfig.Field(
940 doc="Systematic error floor for all zeropoints",
941 dtype=float,
942 default=0.003,
943 )
944 mapLongitudeRef = pexConfig.Field(
945 doc="Reference longitude for plotting maps",
946 dtype=float,
947 default=0.0,
948 )
949 mapNSide = pexConfig.Field(
950 doc="Healpix nside for plotting maps",
951 dtype=int,
952 default=256,
953 )
954 outfileBase = pexConfig.Field(
955 doc="Filename start for plot output files",
956 dtype=str,
957 default=None,
958 )
959 starColorCuts = pexConfig.ListField(
960 doc=("Encoded star-color cuts (using calibration star colors). "
961 "This is a list with each entry a string of the format "
962 "``band1,band2,low,high`` such that only stars of color "
963 "low < band1 - band2 < high will be used for calibration."),
964 dtype=str,
965 default=("NO_DATA",),
966 )
967 refStarColorCuts = pexConfig.ListField(
968 doc=("Encoded star color cuts specifically to apply to reference stars. "
969 "This is a list with each entry a string of the format "
970 "``band1,band2,low,high`` such that only stars of color "
971 "low < band1 - band2 < high will be used as reference stars."),
972 dtype=str,
973 default=("NO_DATA",),
974 )
975 colorSplitBands = pexConfig.ListField(
976 doc="Band names to use to split stars by color. Must have 2 entries.",
977 dtype=str,
978 length=2,
979 default=('g', 'i'),
980 )
981 modelMagErrors = pexConfig.Field(
982 doc="Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)",
983 dtype=bool,
984 default=True,
985 )
986 useQuadraticPwv = pexConfig.Field(
987 doc="Model PWV with a quadratic term for variation through the night?",
988 dtype=bool,
989 default=False,
990 )
991 instrumentParsPerBand = pexConfig.Field(
992 doc=("Model instrumental parameters per band? "
993 "Otherwise, instrumental parameters (QE changes with time) are "
994 "shared among all bands."),
995 dtype=bool,
996 default=False,
997 )
998 instrumentSlopeMinDeltaT = pexConfig.Field(
999 doc=("Minimum time change (in days) between observations to use in constraining "
1000 "instrument slope."),
1001 dtype=float,
1002 default=20.0,
1003 )
1004 fitMirrorChromaticity = pexConfig.Field(
1005 doc="Fit (intraband) mirror chromatic term?",
1006 dtype=bool,
1007 default=False,
1008 )
1009 fitCcdChromaticityDict = pexConfig.DictField(
1010 doc="Specification on whether to compute first-order quantum efficiency (QE) "
1011 "adjustments. Key is band, and value will be True or False. Any band "
1012 "not explicitly specified will default to False.",
1013 keytype=str,
1014 itemtype=bool,
1015 default={},
1016 )
1017 coatingMjds = pexConfig.ListField(
1018 doc="Mirror coating dates in MJD",
1019 dtype=float,
1020 default=(0.0,),
1021 )
1022 outputStandardsBeforeFinalCycle = pexConfig.Field(
1023 doc="Output standard stars prior to final cycle? Used in debugging.",
1024 dtype=bool,
1025 default=False,
1026 )
1027 outputZeropointsBeforeFinalCycle = pexConfig.Field(
1028 doc="Output standard stars prior to final cycle? Used in debugging.",
1029 dtype=bool,
1030 default=False,
1031 )
1032 useRepeatabilityForExpGrayCutsDict = pexConfig.DictField(
1033 doc=("Per-band specification on whether to use star repeatability (instead of exposures) "
1034 "for computing photometric cuts. Recommended for tract mode or bands with few visits."),
1035 keytype=str,
1036 itemtype=bool,
1037 default={},
1038 )
1039 autoPhotometricCutNSig = pexConfig.Field(
1040 doc=("Number of sigma for automatic computation of (low) photometric cut. "
1041 "Cut is based on exposure gray width (per band), unless "
1042 "useRepeatabilityForExpGrayCuts is set, in which case the star "
1043 "repeatability is used (also per band)."),
1044 dtype=float,
1045 default=3.0,
1046 )
1047 autoHighCutNSig = pexConfig.Field(
1048 doc=("Number of sigma for automatic computation of (high) outlier cut. "
1049 "Cut is based on exposure gray width (per band), unless "
1050 "useRepeatabilityForExpGrayCuts is set, in which case the star "
1051 "repeatability is used (also per band)."),
1052 dtype=float,
1053 default=4.0,
1054 )
1055 quietMode = pexConfig.Field(
1056 doc="Be less verbose with logging.",
1057 dtype=bool,
1058 default=False,
1059 )
1060 doPlots = pexConfig.Field(
1061 doc="Make fgcm QA plots.",
1062 dtype=bool,
1063 default=True,
1064 )
1065 doPlotsBeforeFinalCycles = pexConfig.Field(
1066 doc="Make fgcm QA plots before the final two fit cycles? This only applies in"
1067 "multi-cycle mode, and if doPlots is True. These are typically the most"
1068 "important QA plots to inspect.",
1069 dtype=bool,
1070 default=False,
1071 )
1072 randomSeed = pexConfig.Field(
1073 doc="Random seed for fgcm for consistency in tests.",
1074 dtype=int,
1075 default=None,
1076 optional=True,
1077 )
1078 deltaAperFitMinNgoodObs = pexConfig.Field(
1079 doc="Minimum number of good observations to use mean delta-aper values in fits.",
1080 dtype=int,
1081 default=2,
1082 )
1083 deltaAperFitPerCcdNx = pexConfig.Field(
1084 doc=("Number of x bins per ccd when computing delta-aper background offsets. "
1085 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
1086 dtype=int,
1087 default=10,
1088 )
1089 deltaAperFitPerCcdNy = pexConfig.Field(
1090 doc=("Number of y bins per ccd when computing delta-aper background offsets. "
1091 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
1092 dtype=int,
1093 default=10,
1094 )
1095 deltaAperFitSpatialNside = pexConfig.Field(
1096 doc="Healpix nside to compute spatial delta-aper background offset maps.",
1097 dtype=int,
1098 default=64,
1099 )
1100 deltaAperInnerRadiusArcsec = pexConfig.Field(
1101 doc=("Inner radius used to compute deltaMagAper (arcseconds). "
1102 "Must be positive and less than ``deltaAperOuterRadiusArcsec`` if "
1103 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
1104 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
1105 dtype=float,
1106 default=0.0,
1107 )
1108 deltaAperOuterRadiusArcsec = pexConfig.Field(
1109 doc=("Outer radius used to compute deltaMagAper (arcseconds). "
1110 "Must be positive and greater than ``deltaAperInnerRadiusArcsec`` if "
1111 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
1112 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
1113 dtype=float,
1114 default=0.0,
1115 )
1116 doComputeDeltaAperPerVisit = pexConfig.Field(
1117 doc=("Do the computation of delta-aper background offsets per visit? "
1118 "Note: this option can be very slow when there are many visits."),
1119 dtype=bool,
1120 default=False,
1121 )
1122 doComputeDeltaAperPerStar = pexConfig.Field(
1123 doc="Do the computation of delta-aper mean values per star?",
1124 dtype=bool,
1125 default=True,
1126 )
1127 doComputeDeltaAperMap = pexConfig.Field(
1128 doc=("Do the computation of delta-aper spatial maps? "
1129 "This is only used if ``doComputeDeltaAperPerStar`` is True,"),
1130 dtype=bool,
1131 default=False,
1132 )
1133 doComputeDeltaAperPerCcd = pexConfig.Field(
1134 doc="Do the computation of per-ccd delta-aper background offsets?",
1135 dtype=bool,
1136 default=False,
1137 )
1139 def validate(self):
1140 super().validate()
1142 if self.connections.previousCycleNumber != str(self.cycleNumber - 1):
1143 msg = "cycleNumber in template must be connections.previousCycleNumber + 1"
1144 raise RuntimeError(msg)
1145 if self.connections.cycleNumber != str(self.cycleNumber):
1146 msg = "cycleNumber in template must be equal to connections.cycleNumber"
1147 raise RuntimeError(msg)
1149 for band in self.fitBands:
1150 if band not in self.bands:
1151 msg = 'fitBand %s not in bands' % (band)
1152 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg)
1153 for band in self.requiredBands:
1154 if band not in self.bands:
1155 msg = 'requiredBand %s not in bands' % (band)
1156 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg)
1157 for band in self.colorSplitBands:
1158 if band not in self.bands:
1159 msg = 'colorSplitBand %s not in bands' % (band)
1160 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg)
1161 for band in self.bands:
1162 if band not in self.superStarSubCcdDict:
1163 msg = 'band %s not in superStarSubCcdDict' % (band)
1164 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict,
1165 self, msg)
1166 if band not in self.ccdGraySubCcdDict:
1167 msg = 'band %s not in ccdGraySubCcdDict' % (band)
1168 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict,
1169 self, msg)
1170 if band not in self.expGrayPhotometricCutDict:
1171 msg = 'band %s not in expGrayPhotometricCutDict' % (band)
1172 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict,
1173 self, msg)
1174 if band not in self.expGrayHighCutDict:
1175 msg = 'band %s not in expGrayHighCutDict' % (band)
1176 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict,
1177 self, msg)
1178 if band not in self.expVarGrayPhotometricCutDict:
1179 msg = 'band %s not in expVarGrayPhotometricCutDict' % (band)
1180 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict,
1181 self, msg)
1182 if band not in self.sigFgcmMaxEGrayDict:
1183 msg = 'band %s not in sigFgcmMaxEGrayDict' % (band)
1184 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict,
1185 self, msg)
1186 if band not in self.approxThroughputDict:
1187 msg = 'band %s not in approxThroughputDict' % (band)
1188 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict,
1189 self, msg)
1190 if band not in self.useRepeatabilityForExpGrayCutsDict:
1191 msg = 'band %s not in useRepeatabilityForExpGrayCutsDict' % (band)
1192 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
1193 self, msg)
1195 if self.doComputeDeltaAperPerVisit or self.doComputeDeltaAperMap \
1196 or self.doComputeDeltaAperPerCcd:
1197 if self.deltaAperInnerRadiusArcsec <= 0.0:
1198 msg = 'deltaAperInnerRadiusArcsec must be positive if deltaAper computations are turned on.'
1199 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperInnerRadiusArcsec,
1200 self, msg)
1201 if self.deltaAperOuterRadiusArcsec <= 0.0:
1202 msg = 'deltaAperOuterRadiusArcsec must be positive if deltaAper computations are turned on.'
1203 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
1204 self, msg)
1205 if self.deltaAperOuterRadiusArcsec <= self.deltaAperInnerRadiusArcsec:
1206 msg = ('deltaAperOuterRadiusArcsec must be greater than deltaAperInnerRadiusArcsec if '
1207 'deltaAper computations are turned on.')
1208 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
1209 self, msg)
1212class FgcmFitCycleTask(pipeBase.PipelineTask):
1213 """
1214 Run Single fit cycle for FGCM global calibration
1215 """
1217 ConfigClass = FgcmFitCycleConfig
1218 _DefaultName = "fgcmFitCycle"
1220 def __init__(self, initInputs=None, **kwargs):
1221 super().__init__(**kwargs)
1223 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1224 camera = butlerQC.get(inputRefs.camera)
1226 nCore = butlerQC.resources.num_cores
1228 handleDict = {}
1230 handleDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
1231 handleDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
1233 if self.config.useParquetCatalogFormat:
1234 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservationsParquet)
1235 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIdsParquet)
1236 if self.config.doReferenceCalibration:
1237 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStarsParquet)
1238 else:
1239 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations)
1240 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds)
1241 handleDict['fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices)
1242 if self.config.doReferenceCalibration:
1243 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars)
1244 if self.config.cycleNumber > 0:
1245 handleDict['fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput)
1246 handleDict['fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput)
1248 fgcmDatasetDict = None
1249 if self.config.doMultipleCycles:
1250 # Run multiple cycles at once.
1251 config = copy.copy(self.config)
1252 config.update(cycleNumber=0)
1253 for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1):
1254 if cycle == self.config.multipleCyclesFinalCycleNumber:
1255 config.update(isFinalCycle=True)
1257 if cycle > 0:
1258 handleDict['fgcmFlaggedStars'] = fgcmDatasetDict['fgcmFlaggedStars']
1259 handleDict['fgcmFitParameters'] = fgcmDatasetDict['fgcmFitParameters']
1261 # Set up plot outputs.
1262 # Note that because the plot outputConnections are skipped in
1263 # the connections class if config.doPlot is False, this will be
1264 # a no-op in that case.
1265 plotHandleDict = {}
1266 for outputRefName in outputRefs.keys():
1267 if outputRefName.endswith("Plot") and f"Cycle{cycle}" in outputRefName:
1268 refs = getattr(outputRefs, outputRefName)
1269 if isinstance(refs, (tuple, list)):
1270 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions:
1271 for ref in refs:
1272 physical_filter = ref.dataId["physical_filter"]
1273 detector = ref.dataId["detector"]
1274 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}"
1275 plotHandleDict[handleDictKey] = ref
1276 elif "physical_filter" in refs[0].dimensions:
1277 for ref in refs:
1278 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}"
1279 plotHandleDict[handleDictKey] = ref
1280 elif "band" in refs[0].dimensions:
1281 for ref in refs:
1282 handleDictKey = f"{outputRefName}_{ref.dataId['band']}"
1283 plotHandleDict[handleDictKey] = ref
1284 else:
1285 plotHandleDict[outputRefName] = refs
1287 fgcmDatasetDict, config = self._fgcmFitCycle(
1288 camera,
1289 handleDict,
1290 butlerQC=butlerQC,
1291 plotHandleDict=plotHandleDict,
1292 config=config,
1293 nCore=nCore,
1294 )
1295 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'],
1296 getattr(outputRefs, f'fgcm_Cycle{cycle}_FitParameters'))
1297 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'],
1298 getattr(outputRefs, f'fgcm_Cycle{cycle}_FlaggedStars'))
1299 butlerQC.put(config,
1300 getattr(outputRefs, f'fgcm_Cycle{cycle}_OutputConfig'))
1301 if self.outputZeropoints:
1302 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'],
1303 getattr(outputRefs, f'fgcm_Cycle{cycle}_Zeropoints'))
1304 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'],
1305 getattr(outputRefs, f'fgcm_Cycle{cycle}_AtmosphereParameters'))
1306 if self.outputStandards:
1307 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'],
1308 getattr(outputRefs, f'fgcm_Cycle{cycle}_StandardStars'))
1309 else:
1310 # Run a single cycle
1312 # Set up plot outputs.
1313 # Note that nothing will go in the dict if doPlots is False.
1314 plotHandleDict = {}
1315 for outputRefName in outputRefs.keys():
1316 if outputRefName.endswith("Plot") and f"Cycle{self.config.cycleNumber}" in outputRefName:
1317 refs = getattr(outputRefs, outputRefName)
1318 if isinstance(refs, (tuple, list)):
1319 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions:
1320 for ref in refs:
1321 physical_filter = ref.dataId["physical_filter"]
1322 detector = ref.dataId["detector"]
1323 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}"
1324 plotHandleDict[handleDictKey] = ref
1325 elif "physical_filter" in refs[0].dimensions:
1326 for ref in refs:
1327 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}"
1328 plotHandleDict[handleDictKey] = ref
1329 elif "band" in refs[0].dimensions:
1330 for ref in refs:
1331 handleDictKey = f"{outputRefName}_{ref.dataId['band']}"
1332 plotHandleDict[handleDictKey] = ref
1333 else:
1334 plotHandleDict[outputRefName] = refs
1336 fgcmDatasetDict, _ = self._fgcmFitCycle(
1337 camera,
1338 handleDict,
1339 nCore=nCore,
1340 butlerQC=butlerQC,
1341 plotHandleDict=plotHandleDict,
1342 multiCycle=False,
1343 )
1345 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], outputRefs.fgcmFitParameters)
1346 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars)
1347 if self.outputZeropoints:
1348 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], outputRefs.fgcmZeropoints)
1349 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters)
1350 if self.outputStandards:
1351 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], outputRefs.fgcmStandardStars)
1353 def _fgcmFitCycle(
1354 self,
1355 camera,
1356 handleDict,
1357 butlerQC=None,
1358 plotHandleDict=None,
1359 config=None,
1360 nCore=1,
1361 multiCycle=True,
1362 ):
1363 """
1364 Run the fit cycle
1366 Parameters
1367 ----------
1368 camera : `lsst.afw.cameraGeom.Camera`
1369 handleDict : `dict`
1370 All handles are `lsst.daf.butler.DeferredDatasetHandle`
1371 handle dictionary with keys:
1373 ``"fgcmLookUpTable"``
1374 handle for the FGCM look-up table.
1375 ``"fgcmVisitCatalog"``
1376 handle for visit summary catalog.
1377 ``"fgcmStarObservations"``
1378 handle for star observation catalog.
1379 ``"fgcmStarIds"``
1380 handle for star id catalog.
1381 ``"fgcmStarIndices"``
1382 handle for star index catalog.
1383 ``"fgcmReferenceStars"``
1384 handle for matched reference star catalog.
1385 ``"fgcmFlaggedStars"``
1386 handle for flagged star catalog.
1387 ``"fgcmFitParameters"``
1388 handle for fit parameter catalog.
1389 butlerQC : `lsst.pipe.base.QuantumContext`, optional
1390 Quantum context used for serializing plots.
1391 plotHandleDict : `dict` [`str`, `lsst.daf.butler.DatasetRef`], optional
1392 Dictionary of plot dataset refs, keyed by plot name.
1393 config : `lsst.pex.config.Config`, optional
1394 Configuration to use to override self.config.
1395 nCore : `int`, optional
1396 Number of cores to use during fitting.
1397 multiCycle : `bool`, optional
1398 Is this part of a multicycle run?
1400 Returns
1401 -------
1402 fgcmDatasetDict : `dict`
1403 Dictionary of datasets to persist.
1404 """
1405 if config is not None:
1406 _config = config
1407 else:
1408 _config = self.config
1410 # Set defaults on whether to output standards and zeropoints
1411 self.maxIter = _config.maxIterBeforeFinalCycle
1412 self.outputStandards = _config.outputStandardsBeforeFinalCycle
1413 self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle
1414 self.resetFitParameters = True
1416 if _config.isFinalCycle:
1417 # This is the final fit cycle, so we do not want to reset fit
1418 # parameters, we want to run a final "clean-up" with 0 fit iterations,
1419 # and we always want to output standards and zeropoints
1420 self.maxIter = 0
1421 self.outputStandards = True
1422 self.outputZeropoints = True
1423 self.resetFitParameters = False
1425 lutCat = handleDict['fgcmLookUpTable'].get()
1426 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
1427 dict(_config.physicalFilterMap))
1428 del lutCat
1430 # Check if we want to do plots.
1431 doPlots = _config.doPlots
1432 if doPlots and multiCycle:
1433 if _config.cycleNumber < (_config.multipleCyclesFinalCycleNumber - 1) \
1434 and not _config.doPlotsBeforeFinalCycles:
1435 doPlots = False
1437 configDict = makeConfigDict(_config, self.log, camera,
1438 self.maxIter, self.resetFitParameters,
1439 self.outputZeropoints,
1440 lutIndexVals[0]['FILTERNAMES'],
1441 nCore=nCore,
1442 doPlots=doPlots)
1444 # next we need the exposure/visit information
1445 visitCat = handleDict['fgcmVisitCatalog'].get()
1446 fgcmExpInfo = translateVisitCatalog(visitCat)
1447 del visitCat
1449 if len(camera) == fgcmLut.nCCD:
1450 useScienceDetectors = False
1451 else:
1452 # If the LUT has a different number of detectors than
1453 # the camera, then we only want to use science detectors
1454 # in the focal plane projector.
1455 useScienceDetectors = True
1457 focalPlaneProjector = FocalPlaneProjector(
1458 camera,
1459 self.config.defaultCameraOrientation,
1460 useScienceDetectors=useScienceDetectors,
1461 )
1463 noFitsDict = {'lutIndex': lutIndexVals,
1464 'lutStd': lutStd,
1465 'expInfo': fgcmExpInfo,
1466 'focalPlaneProjector': focalPlaneProjector}
1468 # set up the fitter object
1469 fgcmFitCycle = fgcm.FgcmFitCycle(
1470 configDict,
1471 useFits=False,
1472 noFitsDict=noFitsDict,
1473 noOutput=True,
1474 butlerQC=butlerQC,
1475 plotHandleDict=plotHandleDict,
1476 )
1478 # create the parameter object
1479 if (fgcmFitCycle.initialCycle):
1480 # cycle = 0, initial cycle
1481 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
1482 fgcmLut,
1483 fgcmExpInfo,
1484 butlerQC=butlerQC,
1485 plotHandleDict=plotHandleDict)
1486 else:
1487 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog):
1488 parCat = handleDict['fgcmFitParameters']
1489 else:
1490 parCat = handleDict['fgcmFitParameters'].get()
1491 inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1492 del parCat
1493 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
1494 fgcmExpInfo,
1495 inParInfo,
1496 inParams,
1497 inSuperStar,
1498 butlerQC=butlerQC,
1499 plotHandleDict=plotHandleDict)
1501 # set up the stars...
1502 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig, butlerQC=butlerQC, plotHandleDict=plotHandleDict)
1504 starObs = handleDict['fgcmStarObservations'].get()
1505 starIds = handleDict['fgcmStarIds'].get()
1506 if not self.config.useParquetCatalogFormat:
1507 starIndices = handleDict['fgcmStarIndices'].get()
1508 else:
1509 starIndices = None
1511 # grab the flagged stars if available
1512 if 'fgcmFlaggedStars' in handleDict:
1513 if isinstance(handleDict['fgcmFlaggedStars'], afwTable.BaseCatalog):
1514 flaggedStars = handleDict['fgcmFlaggedStars']
1515 else:
1516 flaggedStars = handleDict['fgcmFlaggedStars'].get()
1517 flagId = flaggedStars['objId'][:]
1518 flagFlag = flaggedStars['objFlag'][:]
1520 del flaggedStars
1521 elif self.config.useParquetCatalogFormat:
1522 # If we are using the parquet catalog format, then that means that
1523 # reserved stars have already been flagged. We extract the flags here
1524 # to input to fgcm, which will then be persisted (with additional
1525 # quality flags) as the fgcmFlaggedStars datatype in subsequent
1526 # fit cycles.
1527 (flagged,) = (starIds['obj_flag'] > 0).nonzero()
1528 flagId = starIds['fgcm_id'][flagged]
1529 flagFlag = starIds['obj_flag'][flagged]
1530 else:
1531 flagId = None
1532 flagFlag = None
1534 if _config.doReferenceCalibration:
1535 refStars = handleDict['fgcmReferenceStars'].get()
1537 refMag, refMagErr = extractReferenceMags(refStars,
1538 _config.bands,
1539 _config.physicalFilterMap)
1541 refId = refStars['fgcm_id'][:]
1542 else:
1543 refStars = None
1544 refId = None
1545 refMag = None
1546 refMagErr = None
1548 # match star observations to visits
1549 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
1550 # actually be transferred into fgcm using the indexing below.
1551 if self.config.useParquetCatalogFormat:
1552 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'])
1553 else:
1554 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'][starIndices['obsIndex']])
1556 # The fgcmStars.loadStars method will copy all the star information into
1557 # special shared memory objects that will not blow up the memory usage when
1558 # used with python multiprocessing. Once all the numbers are copied,
1559 # it is necessary to release all references to the objects that previously
1560 # stored the data to ensure that the garbage collector can clear the memory,
1561 # and ensure that this memory is not copied when multiprocessing kicks in.
1563 if self.config.useParquetCatalogFormat:
1564 # Note that the ra/dec coordinates for the parquet format are in
1565 # degrees, which is what fgcm expects.
1566 fgcmStars.loadStars(fgcmPars,
1567 starObs['visit'][:],
1568 starObs['detector'][:],
1569 starObs['ra'][:],
1570 starObs['dec'][:],
1571 starObs['inst_mag'][:],
1572 starObs['inst_mag_err'][:],
1573 fgcmExpInfo['FILTERNAME'][visitIndex],
1574 starIds['fgcm_id'][:],
1575 starIds['ra'][:],
1576 starIds['dec'][:],
1577 starIds['obs_arr_index'][:],
1578 starIds['n_obs'][:],
1579 obsX=starObs['x'][:],
1580 obsY=starObs['y'][:],
1581 obsDeltaMagBkg=starObs['delta_mag_bkg'][:],
1582 obsDeltaAper=starObs['delta_mag_aper'][:],
1583 refID=refId,
1584 refMag=refMag,
1585 refMagErr=refMagErr,
1586 flagID=flagId,
1587 flagFlag=flagFlag,
1588 computeNobs=True,
1589 objIDAlternate=starIds['isolated_star_id'])
1590 else:
1591 # We determine the conversion from the native units (typically radians) to
1592 # degrees for the first star. This allows us to treat coord_ra/coord_dec as
1593 # numpy arrays rather than Angles, which would we approximately 600x slower.
1594 conv = starObs[0]['ra'].asDegrees() / float(starObs[0]['ra'])
1596 fgcmStars.loadStars(fgcmPars,
1597 starObs['visit'][starIndices['obsIndex']],
1598 starObs['ccd'][starIndices['obsIndex']],
1599 starObs['ra'][starIndices['obsIndex']] * conv,
1600 starObs['dec'][starIndices['obsIndex']] * conv,
1601 starObs['instMag'][starIndices['obsIndex']],
1602 starObs['instMagErr'][starIndices['obsIndex']],
1603 fgcmExpInfo['FILTERNAME'][visitIndex],
1604 starIds['fgcm_id'][:],
1605 starIds['ra'][:],
1606 starIds['dec'][:],
1607 starIds['obsArrIndex'][:],
1608 starIds['nObs'][:],
1609 obsX=starObs['x'][starIndices['obsIndex']],
1610 obsY=starObs['y'][starIndices['obsIndex']],
1611 obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']],
1612 obsDeltaAper=starObs['deltaMagAper'][starIndices['obsIndex']],
1613 psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']],
1614 refID=refId,
1615 refMag=refMag,
1616 refMagErr=refMagErr,
1617 flagID=flagId,
1618 flagFlag=flagFlag,
1619 computeNobs=True)
1621 # Release all references to temporary objects holding star data (see above)
1622 del starObs
1623 del starIds
1624 del starIndices
1625 del flagId
1626 del flagFlag
1627 del refStars
1628 del refId
1629 del refMag
1630 del refMagErr
1632 # and set the bits in the cycle object
1633 fgcmFitCycle.setLUT(fgcmLut)
1634 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
1635 fgcmFitCycle.setPars(fgcmPars)
1637 # finish the setup
1638 fgcmFitCycle.finishSetup()
1640 # and run
1641 fgcmFitCycle.run()
1643 ##################
1644 # Persistance
1645 ##################
1647 fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle)
1649 # Output the config for the next cycle
1650 # We need to make a copy since the input one has been frozen
1652 updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i]) for
1653 i, b in enumerate(_config.bands)}
1654 updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i]) for
1655 i, band in enumerate(_config.bands)}
1657 outConfig = copy.copy(_config)
1658 outConfig.update(cycleNumber=(_config.cycleNumber + 1),
1659 precomputeSuperStarInitialCycle=False,
1660 freezeStdAtmosphere=False,
1661 expGrayPhotometricCutDict=updatedPhotometricCutDict,
1662 expGrayHighCutDict=updatedHighCutDict)
1664 outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber),
1665 cycleNumber=str(_config.cycleNumber + 1))
1667 if not multiCycle:
1668 configFileName = '%s_cycle%02d_config.py' % (outConfig.outfileBase,
1669 outConfig.cycleNumber)
1670 outConfig.save(configFileName)
1672 if _config.isFinalCycle == 1:
1673 # We are done, ready to output products
1674 self.log.info("Everything is in place to run fgcmOutputProducts.py")
1675 else:
1676 self.log.info("Saved config for next cycle to %s" % (configFileName))
1677 self.log.info("Be sure to look at:")
1678 self.log.info(" config.expGrayPhotometricCut")
1679 self.log.info(" config.expGrayHighCut")
1680 self.log.info("If you are satisfied with the fit, please set:")
1681 self.log.info(" config.isFinalCycle = True")
1683 fgcmFitCycle.freeSharedMemory()
1685 return fgcmDatasetDict, outConfig
1687 def _loadParameters(self, parCat):
1688 """
1689 Load FGCM parameters from a previous fit cycle
1691 Parameters
1692 ----------
1693 parCat : `lsst.afw.table.BaseCatalog`
1694 Parameter catalog in afw table form.
1696 Returns
1697 -------
1698 inParInfo: `numpy.ndarray`
1699 Numpy array parameter information formatted for input to fgcm
1700 inParameters: `numpy.ndarray`
1701 Numpy array parameter values formatted for input to fgcm
1702 inSuperStar: `numpy.array`
1703 Superstar flat formatted for input to fgcm
1704 """
1705 parLutFilterNames = np.array(parCat[0]['lutFilterNames'].split(','))
1706 parFitBands = np.array(parCat[0]['fitBands'].split(','))
1708 inParInfo = np.zeros(1, dtype=[('NCCD', 'i4'),
1709 ('LUTFILTERNAMES', parLutFilterNames.dtype.str,
1710 (parLutFilterNames.size, )),
1711 ('FITBANDS', parFitBands.dtype.str, (parFitBands.size, )),
1712 ('LNTAUUNIT', 'f8'),
1713 ('LNTAUSLOPEUNIT', 'f8'),
1714 ('ALPHAUNIT', 'f8'),
1715 ('LNPWVUNIT', 'f8'),
1716 ('LNPWVSLOPEUNIT', 'f8'),
1717 ('LNPWVQUADRATICUNIT', 'f8'),
1718 ('LNPWVGLOBALUNIT', 'f8'),
1719 ('O3UNIT', 'f8'),
1720 ('QESYSUNIT', 'f8'),
1721 ('FILTEROFFSETUNIT', 'f8'),
1722 ('HASEXTERNALPWV', 'i2'),
1723 ('HASEXTERNALTAU', 'i2')])
1724 inParInfo['NCCD'] = parCat['nCcd']
1725 inParInfo['LUTFILTERNAMES'][:] = parLutFilterNames
1726 inParInfo['FITBANDS'][:] = parFitBands
1727 inParInfo['HASEXTERNALPWV'] = parCat['hasExternalPwv']
1728 inParInfo['HASEXTERNALTAU'] = parCat['hasExternalTau']
1730 inParams = np.zeros(1, dtype=[('PARALPHA', 'f8', (parCat['parAlpha'].size, )),
1731 ('PARO3', 'f8', (parCat['parO3'].size, )),
1732 ('PARLNTAUINTERCEPT', 'f8',
1733 (parCat['parLnTauIntercept'].size, )),
1734 ('PARLNTAUSLOPE', 'f8',
1735 (parCat['parLnTauSlope'].size, )),
1736 ('PARLNPWVINTERCEPT', 'f8',
1737 (parCat['parLnPwvIntercept'].size, )),
1738 ('PARLNPWVSLOPE', 'f8',
1739 (parCat['parLnPwvSlope'].size, )),
1740 ('PARLNPWVQUADRATIC', 'f8',
1741 (parCat['parLnPwvQuadratic'].size, )),
1742 ('PARQESYSINTERCEPT', 'f8',
1743 (parCat['parQeSysIntercept'].size, )),
1744 ('COMPQESYSSLOPE', 'f8',
1745 (parCat['compQeSysSlope'].size, )),
1746 ('PARFILTEROFFSET', 'f8',
1747 (parCat['parFilterOffset'].size, )),
1748 ('PARFILTEROFFSETFITFLAG', 'i2',
1749 (parCat['parFilterOffsetFitFlag'].size, )),
1750 ('PARRETRIEVEDLNPWVSCALE', 'f8'),
1751 ('PARRETRIEVEDLNPWVOFFSET', 'f8'),
1752 ('PARRETRIEVEDLNPWVNIGHTLYOFFSET', 'f8',
1753 (parCat['parRetrievedLnPwvNightlyOffset'].size, )),
1754 ('COMPABSTHROUGHPUT', 'f8',
1755 (parCat['compAbsThroughput'].size, )),
1756 ('COMPREFOFFSET', 'f8',
1757 (parCat['compRefOffset'].size, )),
1758 ('COMPREFSIGMA', 'f8',
1759 (parCat['compRefSigma'].size, )),
1760 ('COMPMIRRORCHROMATICITY', 'f8',
1761 (parCat['compMirrorChromaticity'].size, )),
1762 ('MIRRORCHROMATICITYPIVOT', 'f8',
1763 (parCat['mirrorChromaticityPivot'].size, )),
1764 ('COMPCCDCHROMATICITY', 'f8',
1765 (parCat['compCcdChromaticity'].size, )),
1766 ('COMPMEDIANSEDSLOPE', 'f8',
1767 (parCat['compMedianSedSlope'].size, )),
1768 ('COMPAPERCORRPIVOT', 'f8',
1769 (parCat['compAperCorrPivot'].size, )),
1770 ('COMPAPERCORRSLOPE', 'f8',
1771 (parCat['compAperCorrSlope'].size, )),
1772 ('COMPAPERCORRSLOPEERR', 'f8',
1773 (parCat['compAperCorrSlopeErr'].size, )),
1774 ('COMPAPERCORRRANGE', 'f8',
1775 (parCat['compAperCorrRange'].size, )),
1776 ('COMPMODELERREXPTIMEPIVOT', 'f8',
1777 (parCat['compModelErrExptimePivot'].size, )),
1778 ('COMPMODELERRFWHMPIVOT', 'f8',
1779 (parCat['compModelErrFwhmPivot'].size, )),
1780 ('COMPMODELERRSKYPIVOT', 'f8',
1781 (parCat['compModelErrSkyPivot'].size, )),
1782 ('COMPMODELERRPARS', 'f8',
1783 (parCat['compModelErrPars'].size, )),
1784 ('COMPEXPGRAY', 'f8',
1785 (parCat['compExpGray'].size, )),
1786 ('COMPVARGRAY', 'f8',
1787 (parCat['compVarGray'].size, )),
1788 ('COMPEXPDELTAMAGBKG', 'f8',
1789 (parCat['compExpDeltaMagBkg'].size, )),
1790 ('COMPNGOODSTARPEREXP', 'i4',
1791 (parCat['compNGoodStarPerExp'].size, )),
1792 ('COMPEXPREFOFFSET', 'f8',
1793 (parCat['compExpRefOffset'].size, )),
1794 ('COMPSIGFGCM', 'f8',
1795 (parCat['compSigFgcm'].size, )),
1796 ('COMPSIGMACAL', 'f8',
1797 (parCat['compSigmaCal'].size, )),
1798 ('COMPRETRIEVEDLNPWV', 'f8',
1799 (parCat['compRetrievedLnPwv'].size, )),
1800 ('COMPRETRIEVEDLNPWVRAW', 'f8',
1801 (parCat['compRetrievedLnPwvRaw'].size, )),
1802 ('COMPRETRIEVEDLNPWVFLAG', 'i2',
1803 (parCat['compRetrievedLnPwvFlag'].size, )),
1804 ('COMPRETRIEVEDTAUNIGHT', 'f8',
1805 (parCat['compRetrievedTauNight'].size, )),
1806 ('COMPEPSILON', 'f8',
1807 (parCat['compEpsilon'].size, )),
1808 ('COMPMEDDELTAAPER', 'f8',
1809 (parCat['compMedDeltaAper'].size, )),
1810 ('COMPGLOBALEPSILON', 'f4',
1811 (parCat['compGlobalEpsilon'].size, )),
1812 ('COMPEPSILONMAP', 'f4',
1813 (parCat['compEpsilonMap'].size, )),
1814 ('COMPEPSILONNSTARMAP', 'i4',
1815 (parCat['compEpsilonNStarMap'].size, )),
1816 ('COMPEPSILONCCDMAP', 'f4',
1817 (parCat['compEpsilonCcdMap'].size, )),
1818 ('COMPEPSILONCCDNSTARMAP', 'i4',
1819 (parCat['compEpsilonCcdNStarMap'].size, ))])
1821 inParams['PARALPHA'][:] = parCat['parAlpha'][0, :]
1822 inParams['PARO3'][:] = parCat['parO3'][0, :]
1823 inParams['PARLNTAUINTERCEPT'][:] = parCat['parLnTauIntercept'][0, :]
1824 inParams['PARLNTAUSLOPE'][:] = parCat['parLnTauSlope'][0, :]
1825 inParams['PARLNPWVINTERCEPT'][:] = parCat['parLnPwvIntercept'][0, :]
1826 inParams['PARLNPWVSLOPE'][:] = parCat['parLnPwvSlope'][0, :]
1827 inParams['PARLNPWVQUADRATIC'][:] = parCat['parLnPwvQuadratic'][0, :]
1828 inParams['PARQESYSINTERCEPT'][:] = parCat['parQeSysIntercept'][0, :]
1829 inParams['COMPQESYSSLOPE'][:] = parCat['compQeSysSlope'][0, :]
1830 inParams['PARFILTEROFFSET'][:] = parCat['parFilterOffset'][0, :]
1831 inParams['PARFILTEROFFSETFITFLAG'][:] = parCat['parFilterOffsetFitFlag'][0, :]
1832 inParams['PARRETRIEVEDLNPWVSCALE'] = parCat['parRetrievedLnPwvScale']
1833 inParams['PARRETRIEVEDLNPWVOFFSET'] = parCat['parRetrievedLnPwvOffset']
1834 inParams['PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat['parRetrievedLnPwvNightlyOffset'][0, :]
1835 inParams['COMPABSTHROUGHPUT'][:] = parCat['compAbsThroughput'][0, :]
1836 inParams['COMPREFOFFSET'][:] = parCat['compRefOffset'][0, :]
1837 inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :]
1838 inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :]
1839 inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :]
1840 inParams['COMPCCDCHROMATICITY'][:] = parCat['compCcdChromaticity'][0, :]
1841 inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :]
1842 inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :]
1843 inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :]
1844 inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :]
1845 inParams['COMPAPERCORRRANGE'][:] = parCat['compAperCorrRange'][0, :]
1846 inParams['COMPMODELERREXPTIMEPIVOT'][:] = parCat['compModelErrExptimePivot'][0, :]
1847 inParams['COMPMODELERRFWHMPIVOT'][:] = parCat['compModelErrFwhmPivot'][0, :]
1848 inParams['COMPMODELERRSKYPIVOT'][:] = parCat['compModelErrSkyPivot'][0, :]
1849 inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :]
1850 inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :]
1851 inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
1852 inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
1853 inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
1854 inParams['COMPEXPREFOFFSET'][:] = parCat['compExpRefOffset'][0, :]
1855 inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
1856 inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
1857 inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :]
1858 inParams['COMPRETRIEVEDLNPWVRAW'][:] = parCat['compRetrievedLnPwvRaw'][0, :]
1859 inParams['COMPRETRIEVEDLNPWVFLAG'][:] = parCat['compRetrievedLnPwvFlag'][0, :]
1860 inParams['COMPRETRIEVEDTAUNIGHT'][:] = parCat['compRetrievedTauNight'][0, :]
1861 inParams['COMPEPSILON'][:] = parCat['compEpsilon'][0, :]
1862 inParams['COMPMEDDELTAAPER'][:] = parCat['compMedDeltaAper'][0, :]
1863 inParams['COMPGLOBALEPSILON'][:] = parCat['compGlobalEpsilon'][0, :]
1864 inParams['COMPEPSILONMAP'][:] = parCat['compEpsilonMap'][0, :]
1865 inParams['COMPEPSILONNSTARMAP'][:] = parCat['compEpsilonNStarMap'][0, :]
1866 inParams['COMPEPSILONCCDMAP'][:] = parCat['compEpsilonCcdMap'][0, :]
1867 inParams['COMPEPSILONCCDNSTARMAP'][:] = parCat['compEpsilonCcdNStarMap'][0, :]
1869 inSuperStar = np.zeros(parCat['superstarSize'][0, :], dtype='f8')
1870 inSuperStar[:, :, :, :] = parCat['superstar'][0, :].reshape(inSuperStar.shape)
1872 return (inParInfo, inParams, inSuperStar)
1874 def _makeFgcmOutputDatasets(self, fgcmFitCycle):
1875 """
1876 Persist FGCM datasets through the butler.
1878 Parameters
1879 ----------
1880 fgcmFitCycle: `lsst.fgcm.FgcmFitCycle`
1881 Fgcm Fit cycle object
1882 """
1883 fgcmDatasetDict = {}
1885 # Save the parameters
1886 parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays()
1888 parSchema = afwTable.Schema()
1890 comma = ','
1891 lutFilterNameString = comma.join([n.decode('utf-8')
1892 for n in parInfo['LUTFILTERNAMES'][0]])
1893 fitBandString = comma.join([n.decode('utf-8')
1894 for n in parInfo['FITBANDS'][0]])
1896 parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat,
1897 lutFilterNameString, fitBandString)
1898 parCat = self._makeParCatalog(parSchema, parInfo, pars,
1899 fgcmFitCycle.fgcmPars.parSuperStarFlat,
1900 lutFilterNameString, fitBandString)
1902 fgcmDatasetDict['fgcmFitParameters'] = parCat
1904 # Save the indices of the flagged stars
1905 # (stars that have been (a) reserved from the fit for testing and
1906 # (b) bad stars that have failed quality checks.)
1907 flagStarSchema = self._makeFlagStarSchema()
1908 flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices()
1909 flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct)
1911 fgcmDatasetDict['fgcmFlaggedStars'] = flagStarCat
1913 # Save the zeropoint information and atmospheres only if desired
1914 if self.outputZeropoints:
1915 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
1916 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
1918 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
1919 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
1921 fgcmDatasetDict['fgcmZeropoints'] = zptCat
1923 # Save atmosphere values
1924 # These are generated by the same code that generates zeropoints
1925 atmSchema = makeAtmSchema()
1926 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
1928 fgcmDatasetDict['fgcmAtmosphereParameters'] = atmCat
1930 # Save the standard stars (if configured)
1931 if self.outputStandards:
1932 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
1933 stdSchema = makeStdSchema(len(goodBands))
1934 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
1936 fgcmDatasetDict['fgcmStandardStars'] = stdCat
1938 return fgcmDatasetDict
1940 def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
1941 lutFilterNameString, fitBandString):
1942 """
1943 Make the parameter persistence schema
1945 Parameters
1946 ----------
1947 parInfo: `numpy.ndarray`
1948 Parameter information returned by fgcm
1949 pars: `numpy.ndarray`
1950 Parameter values returned by fgcm
1951 parSuperStarFlat: `numpy.array`
1952 Superstar flat values returned by fgcm
1953 lutFilterNameString: `str`
1954 Combined string of all the lutFilterNames
1955 fitBandString: `str`
1956 Combined string of all the fitBands
1958 Returns
1959 -------
1960 parSchema: `afwTable.schema`
1961 """
1963 parSchema = afwTable.Schema()
1965 # parameter info section
1966 parSchema.addField('nCcd', type=np.int32, doc='Number of CCDs')
1967 parSchema.addField('lutFilterNames', type=str, doc='LUT Filter names in parameter file',
1968 size=len(lutFilterNameString))
1969 parSchema.addField('fitBands', type=str, doc='Bands that were fit',
1970 size=len(fitBandString))
1971 parSchema.addField('lnTauUnit', type=np.float64, doc='Step units for ln(AOD)')
1972 parSchema.addField('lnTauSlopeUnit', type=np.float64,
1973 doc='Step units for ln(AOD) slope')
1974 parSchema.addField('alphaUnit', type=np.float64, doc='Step units for alpha')
1975 parSchema.addField('lnPwvUnit', type=np.float64, doc='Step units for ln(pwv)')
1976 parSchema.addField('lnPwvSlopeUnit', type=np.float64,
1977 doc='Step units for ln(pwv) slope')
1978 parSchema.addField('lnPwvQuadraticUnit', type=np.float64,
1979 doc='Step units for ln(pwv) quadratic term')
1980 parSchema.addField('lnPwvGlobalUnit', type=np.float64,
1981 doc='Step units for global ln(pwv) parameters')
1982 parSchema.addField('o3Unit', type=np.float64, doc='Step units for O3')
1983 parSchema.addField('qeSysUnit', type=np.float64, doc='Step units for mirror gray')
1984 parSchema.addField('filterOffsetUnit', type=np.float64, doc='Step units for filter offset')
1985 parSchema.addField('hasExternalPwv', type=np.int32, doc='Parameters fit using external pwv')
1986 parSchema.addField('hasExternalTau', type=np.int32, doc='Parameters fit using external tau')
1988 # parameter section
1989 parSchema.addField('parAlpha', type='ArrayD', doc='Alpha parameter vector',
1990 size=pars['PARALPHA'].size)
1991 parSchema.addField('parO3', type='ArrayD', doc='O3 parameter vector',
1992 size=pars['PARO3'].size)
1993 parSchema.addField('parLnTauIntercept', type='ArrayD',
1994 doc='ln(Tau) intercept parameter vector',
1995 size=pars['PARLNTAUINTERCEPT'].size)
1996 parSchema.addField('parLnTauSlope', type='ArrayD',
1997 doc='ln(Tau) slope parameter vector',
1998 size=pars['PARLNTAUSLOPE'].size)
1999 parSchema.addField('parLnPwvIntercept', type='ArrayD', doc='ln(pwv) intercept parameter vector',
2000 size=pars['PARLNPWVINTERCEPT'].size)
2001 parSchema.addField('parLnPwvSlope', type='ArrayD', doc='ln(pwv) slope parameter vector',
2002 size=pars['PARLNPWVSLOPE'].size)
2003 parSchema.addField('parLnPwvQuadratic', type='ArrayD', doc='ln(pwv) quadratic parameter vector',
2004 size=pars['PARLNPWVQUADRATIC'].size)
2005 parSchema.addField('parQeSysIntercept', type='ArrayD', doc='Mirror gray intercept parameter vector',
2006 size=pars['PARQESYSINTERCEPT'].size)
2007 parSchema.addField('compQeSysSlope', type='ArrayD', doc='Mirror gray slope parameter vector',
2008 size=pars[0]['COMPQESYSSLOPE'].size)
2009 parSchema.addField('parFilterOffset', type='ArrayD', doc='Filter offset parameter vector',
2010 size=pars['PARFILTEROFFSET'].size)
2011 parSchema.addField('parFilterOffsetFitFlag', type='ArrayI', doc='Filter offset parameter fit flag',
2012 size=pars['PARFILTEROFFSETFITFLAG'].size)
2013 parSchema.addField('parRetrievedLnPwvScale', type=np.float64,
2014 doc='Global scale for retrieved ln(pwv)')
2015 parSchema.addField('parRetrievedLnPwvOffset', type=np.float64,
2016 doc='Global offset for retrieved ln(pwv)')
2017 parSchema.addField('parRetrievedLnPwvNightlyOffset', type='ArrayD',
2018 doc='Nightly offset for retrieved ln(pwv)',
2019 size=pars['PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size)
2020 parSchema.addField('compAbsThroughput', type='ArrayD',
2021 doc='Absolute throughput (relative to transmission curves)',
2022 size=pars['COMPABSTHROUGHPUT'].size)
2023 parSchema.addField('compRefOffset', type='ArrayD',
2024 doc='Offset between reference stars and calibrated stars',
2025 size=pars['COMPREFOFFSET'].size)
2026 parSchema.addField('compRefSigma', type='ArrayD',
2027 doc='Width of reference star/calibrated star distribution',
2028 size=pars['COMPREFSIGMA'].size)
2029 parSchema.addField('compMirrorChromaticity', type='ArrayD',
2030 doc='Computed mirror chromaticity terms',
2031 size=pars['COMPMIRRORCHROMATICITY'].size)
2032 parSchema.addField('mirrorChromaticityPivot', type='ArrayD',
2033 doc='Mirror chromaticity pivot mjd',
2034 size=pars['MIRRORCHROMATICITYPIVOT'].size)
2035 parSchema.addField('compCcdChromaticity', type='ArrayD',
2036 doc='Computed CCD chromaticity terms',
2037 size=pars['COMPCCDCHROMATICITY'].size)
2038 parSchema.addField('compMedianSedSlope', type='ArrayD',
2039 doc='Computed median SED slope (per band)',
2040 size=pars['COMPMEDIANSEDSLOPE'].size)
2041 parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot',
2042 size=pars['COMPAPERCORRPIVOT'].size)
2043 parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope',
2044 size=pars['COMPAPERCORRSLOPE'].size)
2045 parSchema.addField('compAperCorrSlopeErr', type='ArrayD', doc='Aperture correction slope error',
2046 size=pars['COMPAPERCORRSLOPEERR'].size)
2047 parSchema.addField('compAperCorrRange', type='ArrayD', doc='Aperture correction range',
2048 size=pars['COMPAPERCORRRANGE'].size)
2049 parSchema.addField('compModelErrExptimePivot', type='ArrayD', doc='Model error exptime pivot',
2050 size=pars['COMPMODELERREXPTIMEPIVOT'].size)
2051 parSchema.addField('compModelErrFwhmPivot', type='ArrayD', doc='Model error fwhm pivot',
2052 size=pars['COMPMODELERRFWHMPIVOT'].size)
2053 parSchema.addField('compModelErrSkyPivot', type='ArrayD', doc='Model error sky pivot',
2054 size=pars['COMPMODELERRSKYPIVOT'].size)
2055 parSchema.addField('compModelErrPars', type='ArrayD', doc='Model error parameters',
2056 size=pars['COMPMODELERRPARS'].size)
2057 parSchema.addField('compExpGray', type='ArrayD', doc='Computed exposure gray',
2058 size=pars['COMPEXPGRAY'].size)
2059 parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance',
2060 size=pars['COMPVARGRAY'].size)
2061 parSchema.addField('compExpDeltaMagBkg', type='ArrayD',
2062 doc='Computed exposure offset due to background',
2063 size=pars['COMPEXPDELTAMAGBKG'].size)
2064 parSchema.addField('compNGoodStarPerExp', type='ArrayI',
2065 doc='Computed number of good stars per exposure',
2066 size=pars['COMPNGOODSTARPEREXP'].size)
2067 parSchema.addField('compExpRefOffset', type='ArrayD',
2068 doc='Computed per-visit median offset between standard stars and ref stars.',
2069 size=pars['COMPEXPREFOFFSET'].size)
2070 parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)',
2071 size=pars['COMPSIGFGCM'].size)
2072 parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)',
2073 size=pars['COMPSIGMACAL'].size)
2074 parSchema.addField('compRetrievedLnPwv', type='ArrayD', doc='Retrieved ln(pwv) (smoothed)',
2075 size=pars['COMPRETRIEVEDLNPWV'].size)
2076 parSchema.addField('compRetrievedLnPwvRaw', type='ArrayD', doc='Retrieved ln(pwv) (raw)',
2077 size=pars['COMPRETRIEVEDLNPWVRAW'].size)
2078 parSchema.addField('compRetrievedLnPwvFlag', type='ArrayI', doc='Retrieved ln(pwv) Flag',
2079 size=pars['COMPRETRIEVEDLNPWVFLAG'].size)
2080 parSchema.addField('compRetrievedTauNight', type='ArrayD', doc='Retrieved tau (per night)',
2081 size=pars['COMPRETRIEVEDTAUNIGHT'].size)
2082 parSchema.addField('compEpsilon', type='ArrayD',
2083 doc='Computed epsilon background offset per visit (nJy/arcsec2)',
2084 size=pars['COMPEPSILON'].size)
2085 parSchema.addField('compMedDeltaAper', type='ArrayD',
2086 doc='Median delta mag aper per visit',
2087 size=pars['COMPMEDDELTAAPER'].size)
2088 parSchema.addField('compGlobalEpsilon', type='ArrayD',
2089 doc='Computed epsilon bkg offset (global) (nJy/arcsec2)',
2090 size=pars['COMPGLOBALEPSILON'].size)
2091 parSchema.addField('compEpsilonMap', type='ArrayD',
2092 doc='Computed epsilon maps (nJy/arcsec2)',
2093 size=pars['COMPEPSILONMAP'].size)
2094 parSchema.addField('compEpsilonNStarMap', type='ArrayI',
2095 doc='Number of stars per pixel in computed epsilon maps',
2096 size=pars['COMPEPSILONNSTARMAP'].size)
2097 parSchema.addField('compEpsilonCcdMap', type='ArrayD',
2098 doc='Computed epsilon ccd maps (nJy/arcsec2)',
2099 size=pars['COMPEPSILONCCDMAP'].size)
2100 parSchema.addField('compEpsilonCcdNStarMap', type='ArrayI',
2101 doc='Number of stars per ccd bin in epsilon ccd maps',
2102 size=pars['COMPEPSILONCCDNSTARMAP'].size)
2103 parSchema.addField('epochMjdStart', type='ArrayD',
2104 doc='Epoch MJD start times',
2105 size=pars['EPOCHMJDSTART'].size)
2106 parSchema.addField('epochMjdEnd', type='ArrayD',
2107 doc='EpochMJD end times',
2108 size=pars['EPOCHMJDEND'].size)
2109 # superstarflat section
2110 parSchema.addField('superstarSize', type='ArrayI', doc='Superstar matrix size',
2111 size=4)
2112 parSchema.addField('superstar', type='ArrayD', doc='Superstar matrix (flattened)',
2113 size=parSuperStarFlat.size)
2115 return parSchema
2117 def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
2118 lutFilterNameString, fitBandString):
2119 """
2120 Make the FGCM parameter catalog for persistence
2122 Parameters
2123 ----------
2124 parSchema: `lsst.afw.table.Schema`
2125 Parameter catalog schema
2126 pars: `numpy.ndarray`
2127 FGCM parameters to put into parCat
2128 parSuperStarFlat: `numpy.array`
2129 FGCM superstar flat array to put into parCat
2130 lutFilterNameString: `str`
2131 Combined string of all the lutFilterNames
2132 fitBandString: `str`
2133 Combined string of all the fitBands
2135 Returns
2136 -------
2137 parCat: `afwTable.BasicCatalog`
2138 Atmosphere and instrumental model parameter catalog for persistence
2139 """
2141 parCat = afwTable.BaseCatalog(parSchema)
2142 parCat.reserve(1)
2144 # The parameter catalog just has one row, with many columns for all the
2145 # atmosphere and instrument fit parameters
2146 rec = parCat.addNew()
2148 # info section
2149 rec['nCcd'] = parInfo['NCCD'][0]
2150 rec['lutFilterNames'] = lutFilterNameString
2151 rec['fitBands'] = fitBandString
2152 # note these are not currently supported here.
2153 rec['hasExternalPwv'] = 0
2154 rec['hasExternalTau'] = 0
2156 # parameter section
2158 scalarNames = ['parRetrievedLnPwvScale', 'parRetrievedLnPwvOffset']
2160 arrNames = ['parAlpha', 'parO3', 'parLnTauIntercept', 'parLnTauSlope',
2161 'parLnPwvIntercept', 'parLnPwvSlope', 'parLnPwvQuadratic',
2162 'parQeSysIntercept', 'compQeSysSlope',
2163 'parRetrievedLnPwvNightlyOffset', 'compAperCorrPivot',
2164 'parFilterOffset', 'parFilterOffsetFitFlag',
2165 'compAbsThroughput', 'compRefOffset', 'compRefSigma',
2166 'compMirrorChromaticity', 'mirrorChromaticityPivot', 'compCcdChromaticity',
2167 'compAperCorrSlope', 'compAperCorrSlopeErr', 'compAperCorrRange',
2168 'compModelErrExptimePivot', 'compModelErrFwhmPivot',
2169 'compModelErrSkyPivot', 'compModelErrPars',
2170 'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm',
2171 'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope',
2172 'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
2173 'compRetrievedTauNight', 'compEpsilon', 'compMedDeltaAper',
2174 'compGlobalEpsilon', 'compEpsilonMap', 'compEpsilonNStarMap',
2175 'compEpsilonCcdMap', 'compEpsilonCcdNStarMap', 'compExpRefOffset',
2176 'epochMjdStart', 'epochMjdEnd']
2178 for scalarName in scalarNames:
2179 rec[scalarName] = pars[scalarName.upper()][0]
2181 for arrName in arrNames:
2182 rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:]
2184 # superstar section
2185 rec['superstarSize'][:] = parSuperStarFlat.shape
2186 rec['superstar'][:] = parSuperStarFlat.ravel()
2188 return parCat
2190 def _makeFlagStarSchema(self):
2191 """
2192 Make the flagged-stars schema
2194 Returns
2195 -------
2196 flagStarSchema: `lsst.afw.table.Schema`
2197 """
2199 flagStarSchema = afwTable.Schema()
2201 flagStarSchema.addField('objId', type=np.int32, doc='FGCM object id')
2202 flagStarSchema.addField('objFlag', type=np.int32, doc='FGCM object flag')
2204 return flagStarSchema
2206 def _makeFlagStarCat(self, flagStarSchema, flagStarStruct):
2207 """
2208 Make the flagged star catalog for persistence
2210 Parameters
2211 ----------
2212 flagStarSchema: `lsst.afw.table.Schema`
2213 Flagged star schema
2214 flagStarStruct: `numpy.ndarray`
2215 Flagged star structure from fgcm
2217 Returns
2218 -------
2219 flagStarCat: `lsst.afw.table.BaseCatalog`
2220 Flagged star catalog for persistence
2221 """
2223 flagStarCat = afwTable.BaseCatalog(flagStarSchema)
2224 flagStarCat.resize(flagStarStruct.size)
2226 flagStarCat['objId'][:] = flagStarStruct['OBJID']
2227 flagStarCat['objFlag'][:] = flagStarStruct['OBJFLAG']
2229 return flagStarCat