Coverage for python / lsst / fgcmcal / fgcmFitCycle.py: 16%
685 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-10 01:56 -0800
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 ("SigmaFgcmReservedStarsCrunchedNonPhotom_Plot",
299 "Plot",
300 "Histograms for intrinsic scatter for reserved bright stars (including non-photometric).",
301 bandDims),
302 ("SigmaFgcmPullsAllStars_Plot",
303 "Plot",
304 "Histograms for pulls for all bright stars.",
305 bandDims),
306 ("SigmaFgcmPullsReservedStars_Plot",
307 "Plot",
308 "Histograms for pulls for reserved bright stars.",
309 bandDims),
310 ("SigmaFgcmPullsReservedStarsCrunched_Plot",
311 "Plot",
312 "Histograms for pulls for reserved bright stars (after gray correction).",
313 bandDims),
314 ("SigmaFgcmPullsReservedStarsCrunchedNonPhotom_Plot",
315 "Plot",
316 "Histograms for pulls for reserved bright stars (including non-photometric).",
317 bandDims),
318 ("SigmaCal_Plot",
319 "Plot",
320 "Plot showing scatter as a function of systematic error floor.",
321 bandDims),
322 ("SigmaRef_Plot",
323 "Plot",
324 "Histograms of scatter with respect to reference stars.",
325 bandDims),
326 ("RefResidVsColorAll_Plot",
327 "Plot",
328 "Plot of reference star residuals vs. color (all stars).",
329 bandDims),
330 ("RefResidVsColorCut_Plot",
331 "Plot",
332 "Plot of reference star residuals vs. color (reference star color cuts).",
333 bandDims),
334 ]
335 )
337 plotConnections.extend(
338 [
339 ("I1R1_Plot", "Plot", "Plot of fgcm R1 vs. I1.", filterDims),
340 ("I1_Plot", "Plot", "Focal plane map of fgcm I1.", filterDims),
341 ("R1_Plot", "Plot", "Focal plane map of fgcm R1.", filterDims),
342 ("R1mI1Matchscale_Plot", "Plot", "Focal plane map of fgcm R1 - I1.", filterDims),
343 ("R1mI1_Plot", "Plot", "Focal plane map of fgcm R1 - I1 (rescaled).", filterDims),
344 ("R1mI1_vs_mjd_Plot", "Plot", "R1 - I1 residuals vs. mjd.", filterDims),
345 ("CompareRedblueMirrorchrom_Plot",
346 "Plot",
347 "Comparison of mirror chromaticity residuals for red/blue stars.",
348 filterDims),
349 ("CcdChromaticity_Plot",
350 "Plot",
351 "Focal plane map of fgcm ccd chromaticity.",
352 filterDims),
353 ("EpsilonDetector_Plot",
354 "Plot",
355 "Focal plane map of background over/undersubtraction.",
356 filterDims),
357 ("EpsilonDetectorMatchscale_Plot",
358 "Plot",
359 "Focal plane map of background over/undersubtraction.",
360 filterDims),
361 ]
362 )
363 for epoch in epochs:
364 plotConnections.extend(
365 [
366 (
367 f"Superstar_{epoch}_Plot",
368 "Plot",
369 "Plot of illumination Correction.",
370 filterDims,
371 ),
372 ],
373 )
374 if self.config.superStarPlotCcdResiduals:
375 plotConnections.extend(
376 [
377 (
378 f"SuperstarResidual_{epoch}_Plot",
379 "Plot",
380 "Binned illumination correction residuals.",
381 filterDetectorDims,
382 ),
383 (
384 f"SuperstarResidualStd_{epoch}_Plot",
385 "Plot",
386 "Binned illumination correction residual stdev.",
387 filterDetectorDims,
388 ),
389 ],
390 )
392 if config.doMultipleCycles:
393 # Multiple cycle run.
395 # All but the final cycle get appended here.
396 for cycle in range(config.multipleCyclesFinalCycleNumber):
397 outputConnections = copy.copy(inputAndOutputConnections)
398 outputConnections.extend(multicycleOutputConnections)
399 if config.outputZeropointsBeforeFinalCycle:
400 outputConnections.extend(optionalZpOutputConnections)
401 if config.outputStandardsBeforeFinalCycle:
402 outputConnections.extend(optionalStarOutputConnections)
404 if config.doPlots:
405 # We will make the plots if doPlots is True and either
406 # it is the next-to-last cycle or the
407 # doPlotsBeforeFinalCycles configuration is set which
408 # means we want plots for all cycles.
409 # The final cycle doPlots is configured below.
410 if cycle == (config.multipleCyclesFinalCycleNumber - 1) \
411 or config.doPlotsBeforeFinalCycles:
412 outputConnections.extend(plotConnections)
414 for (name, storageClass, doc, dims) in outputConnections:
415 connectionName = f"fgcm_Cycle{cycle}_{name}"
416 storageName = connectionName
417 outConnection = connectionTypes.Output(
418 name=storageName,
419 storageClass=storageClass,
420 doc=doc,
421 dimensions=dims,
422 multiple=(len(dims) > 1),
423 )
424 setattr(self, connectionName, outConnection)
426 # The final cycle has everything.
427 outputConnections = copy.copy(inputAndOutputConnections)
428 outputConnections.extend(multicycleOutputConnections)
429 outputConnections.extend(optionalZpOutputConnections)
430 outputConnections.extend(optionalStarOutputConnections)
431 if config.doPlots:
432 outputConnections.extend(plotConnections)
433 for (name, storageClass, doc, dims) in outputConnections:
434 connectionName = f"fgcm_Cycle{config.multipleCyclesFinalCycleNumber}_{name}"
435 storageName = connectionName
436 outConnection = connectionTypes.Output(
437 name=storageName,
438 storageClass=storageClass,
439 doc=doc,
440 dimensions=dims,
441 multiple=(len(dims) > 1),
442 )
443 setattr(self, connectionName, outConnection)
444 else:
445 # Single cycle run.
446 if config.cycleNumber > 0:
447 inputConnections = copy.copy(inputAndOutputConnections)
448 else:
449 inputConnections = []
450 outputConnections = copy.copy(inputAndOutputConnections)
451 # The following configurations are also useful for runs
452 # where fit cycles are run one-at-a-time since it is
453 # not typical to look at these outputs for intermediate
454 # steps.
455 if config.isFinalCycle or config.outputZeropointsBeforeFinalCycle:
456 outputConnections.extend(optionalZpOutputConnections)
457 if config.isFinalCycle or config.outputStandardsBeforeFinalCycle:
458 outputConnections.extend(optionalStarOutputConnections)
460 if config.doPlots:
461 outputConnections.extend(plotConnections)
463 for (name, storageClass, doc, dims) in inputConnections:
464 connectionName = f"fgcm{name}Input"
465 storageName = f"fgcm_Cycle{config.cycleNumber - 1}_{name}"
466 inConnection = connectionTypes.PrerequisiteInput(
467 name=storageName,
468 storageClass=storageClass,
469 doc=doc,
470 dimensions=dims,
471 )
472 setattr(self, connectionName, inConnection)
474 for (name, storageClass, doc, dims) in outputConnections:
475 connectionName = f"fgcm{name}"
476 storageName = f"fgcm_Cycle{config.cycleNumber}_{name}"
477 # Plots have unique names as well.
478 if storageClass == "Plot":
479 connectionName = storageName
480 outConnection = connectionTypes.Output(
481 name=storageName,
482 storageClass=storageClass,
483 doc=doc,
484 dimensions=dims,
485 multiple=(len(dims) > 1),
486 )
487 setattr(self, connectionName, outConnection)
489 if not config.doReferenceCalibration:
490 self.inputs.remove("fgcmReferenceStars")
491 self.inputs.remove("fgcmReferenceStarsParquet")
493 if config.useParquetCatalogFormat:
494 self.inputs.remove("fgcmStarObservations")
495 self.inputs.remove("fgcmStarIds")
496 self.inputs.remove("fgcmStarIndices")
497 if config.doReferenceCalibration:
498 self.inputs.remove("fgcmReferenceStars")
499 else:
500 self.inputs.remove("fgcmStarObservationsParquet")
501 self.inputs.remove("fgcmStarIdsParquet")
502 if config.doReferenceCalibration:
503 self.inputs.remove("fgcmReferenceStarsParquet")
506class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
507 pipelineConnections=FgcmFitCycleConnections):
508 """Config for FgcmFitCycle"""
510 doMultipleCycles = pexConfig.Field(
511 doc="Run multiple fit cycles in one task",
512 dtype=bool,
513 default=False,
514 )
515 useParquetCatalogFormat = pexConfig.Field(
516 doc="Use parquet catalog format?",
517 dtype=bool,
518 default=True,
519 )
520 multipleCyclesFinalCycleNumber = pexConfig.RangeField(
521 doc=("Final cycle number in multiple cycle mode. The initial cycle "
522 "is 0, with limited parameters fit. The next cycle is 1 with "
523 "full parameter fit. The final cycle is a clean-up with no "
524 "parameters fit. There will be a total of "
525 "(multipleCycleFinalCycleNumber + 1) cycles run, and the final "
526 "cycle number cannot be less than 2."),
527 dtype=int,
528 default=5,
529 min=2,
530 max=MULTIPLE_CYCLES_MAX,
531 inclusiveMax=True,
532 )
533 bands = pexConfig.ListField(
534 doc="Bands to run calibration",
535 dtype=str,
536 default=[],
537 )
538 fitBands = pexConfig.ListField(
539 doc=("Bands to use in atmospheric fit. The bands not listed here will have "
540 "the atmosphere constrained from the 'fitBands' on the same night. "
541 "Must be a subset of `config.bands`"),
542 dtype=str,
543 default=[],
544 )
545 requiredBands = pexConfig.ListField(
546 doc=("Bands that are required for a star to be considered a calibration star. "
547 "Must be a subset of `config.bands`"),
548 dtype=str,
549 default=[],
550 )
551 physicalFilterMap = pexConfig.DictField(
552 doc="Mapping from 'physicalFilter' to band.",
553 keytype=str,
554 itemtype=str,
555 default={},
556 )
557 doReferenceCalibration = pexConfig.Field(
558 doc="Use reference catalog as additional constraint on calibration",
559 dtype=bool,
560 default=True,
561 )
562 refStarSnMin = pexConfig.Field(
563 doc="Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.",
564 dtype=float,
565 default=50.0,
566 )
567 refStarOutlierNSig = pexConfig.Field(
568 doc=("Number of sigma compared to average mag for reference star to be considered an outlier. "
569 "Computed per-band, and if it is an outlier in any band it is rejected from fits."),
570 dtype=float,
571 default=4.0,
572 )
573 applyRefStarColorCuts = pexConfig.Field(
574 doc=("Apply color cuts defined in ``starColorCuts`` to reference stars? "
575 "These cuts are in addition to any cuts defined in ``refStarColorCuts``"),
576 dtype=bool,
577 default=True,
578 )
579 refStarMaxFracUse = pexConfig.Field(
580 doc=("Maximum fraction of reference stars to use in the fit. Remainder will "
581 "be used only for validation."),
582 dtype=float,
583 default=0.5,
584 )
585 useExposureReferenceOffset = pexConfig.Field(
586 doc=("Use per-exposure (visit) offsets between calibrated stars and reference stars "
587 "for final zeropoints? This may help uniformity for disjoint surveys."),
588 dtype=bool,
589 default=False,
590 )
591 nStarPerRun = pexConfig.Field(
592 doc="Number of stars to run in each chunk. Larger values tend to be faster for large "
593 "datasets, at the expense of some memory overhead.",
594 dtype=int,
595 default=200000,
596 )
597 nStarPerGrayRun = pexConfig.Field(
598 doc="Number of stars to run in each chunk (including gray correction). Increasing this "
599 "value may increase the peak memory use significantly.",
600 dtype=int,
601 default=50000,
602 )
603 nExpPerRun = pexConfig.Field(
604 doc="Number of exposures to run in each chunk",
605 dtype=int,
606 default=1000,
607 )
608 reserveFraction = pexConfig.Field(
609 doc="Fraction of stars to reserve for testing",
610 dtype=float,
611 default=0.1,
612 )
613 freezeStdAtmosphere = pexConfig.Field(
614 doc="Freeze atmosphere parameters to standard (for testing)",
615 dtype=bool,
616 default=False,
617 )
618 precomputeSuperStarInitialCycle = pexConfig.Field(
619 doc="Precompute superstar flat for initial cycle",
620 dtype=bool,
621 default=False,
622 )
623 superStarSubCcdDict = pexConfig.DictField(
624 doc=("Per-band specification on whether to compute superstar flat on sub-ccd scale. "
625 "Must have one entry per band."),
626 keytype=str,
627 itemtype=bool,
628 default={},
629 )
630 superStarSubCcdChebyshevOrder = pexConfig.Field(
631 doc=("Order of the 2D chebyshev polynomials for sub-ccd superstar fit. "
632 "Global default is first-order polynomials, and should be overridden "
633 "on a camera-by-camera basis depending on the ISR."),
634 dtype=int,
635 default=1,
636 )
637 superStarSubCcdTriangular = pexConfig.Field(
638 doc=("Should the sub-ccd superstar chebyshev matrix be triangular to "
639 "suppress high-order cross terms?"),
640 dtype=bool,
641 default=False,
642 )
643 superStarSigmaClip = pexConfig.Field(
644 doc="Number of sigma to clip outliers when selecting for superstar flats",
645 dtype=float,
646 default=5.0,
647 )
648 superStarPlotCcdResiduals = pexConfig.Field(
649 doc="If plotting is enabled, should per-detector residuals be plotted? "
650 "This may produce a lot of output, and should be used only for "
651 "debugging purposes.",
652 dtype=bool,
653 default=False,
654 )
655 superStarForceZeroMean = pexConfig.Field(
656 doc="When computing the super-star flat, force the focal-plane mean to "
657 "zero (per band)? This should only be used when computing stand-alone "
658 "illumination corrections.",
659 dtype=bool,
660 default=False,
661 )
662 focalPlaneSigmaClip = pexConfig.Field(
663 doc="Number of sigma to clip outliers per focal-plane.",
664 dtype=float,
665 default=4.0,
666 )
667 ccdGraySubCcdDict = pexConfig.DictField(
668 doc=("Per-band specification on whether to compute achromatic per-ccd residual "
669 "('ccd gray') on a sub-ccd scale. This will only be used as a fallback "
670 "if ``ccdGrayFocalPlaneDict[band]`` is set to True."),
671 keytype=str,
672 itemtype=bool,
673 default={},
674 )
675 ccdGraySubCcdChebyshevOrder = pexConfig.Field(
676 doc="Order of the 2D chebyshev polynomials for sub-ccd gray fit.",
677 dtype=int,
678 default=1,
679 )
680 ccdGraySubCcdTriangular = pexConfig.Field(
681 doc=("Should the sub-ccd gray chebyshev matrix be triangular to "
682 "suppress high-order cross terms?"),
683 dtype=bool,
684 default=True,
685 )
686 ccdGrayFocalPlaneDict = pexConfig.DictField(
687 doc=("Per-band specification on whether to compute focal-plane residual "
688 "('ccd gray') corrections. The focal-plane residuals will be used "
689 "instead of per-CCD residuals except when there are insufficient "
690 "good CCDs (see ``ccdGrayFocalPlaneFitMinCcd``). Therefore, this "
691 "is not recommended for large focal planes with non-photometric "
692 "exposures."),
693 keytype=str,
694 itemtype=bool,
695 default={},
696 )
697 ccdGrayFocalPlaneFitMinCcd = pexConfig.Field(
698 doc=("Minimum number of 'good' CCDs required to perform focal-plane "
699 "gray corrections. If there are fewer good CCDs then the gray "
700 "correction is computed per-ccd."),
701 dtype=int,
702 default=1,
703 )
704 ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field(
705 doc="Order of the 2D chebyshev polynomials for focal plane fit.",
706 dtype=int,
707 default=3,
708 )
709 ccdGrayFocalPlaneMaxStars = pexConfig.Field(
710 doc="Maximum number of stars to use for focal plane fit. Required to keep "
711 "matrix memory usage from running away. If there are more stars than "
712 "this then they will be down-sampled.",
713 dtype=int,
714 default=50_000,
715 )
716 cycleNumber = pexConfig.Field(
717 doc=("FGCM fit cycle number. This is automatically incremented after each run "
718 "and stage of outlier rejection. See cookbook for details."),
719 dtype=int,
720 default=None,
721 )
722 isFinalCycle = pexConfig.Field(
723 doc=("Is this the final cycle of the fitting? Will automatically compute final "
724 "selection of stars and photometric exposures, and will output zeropoints "
725 "and standard stars for use in fgcmOutputProducts"),
726 dtype=bool,
727 default=False,
728 )
729 maxIterBeforeFinalCycle = pexConfig.Field(
730 doc=("Maximum fit iterations, prior to final cycle. The number of iterations "
731 "will always be 0 in the final cycle for cleanup and final selection."),
732 dtype=int,
733 default=50,
734 )
735 deltaMagBkgOffsetPercentile = pexConfig.Field(
736 doc=("Percentile brightest stars on a visit/ccd to use to compute net "
737 "offset from local background subtraction."),
738 dtype=float,
739 default=0.25,
740 )
741 deltaMagBkgPerCcd = pexConfig.Field(
742 doc=("Compute net offset from local background subtraction per-ccd? "
743 "Otherwise, use computation per visit."),
744 dtype=bool,
745 default=False,
746 )
747 utBoundary = pexConfig.Field(
748 doc="Boundary (in UTC) from day-to-day",
749 dtype=float,
750 default=None,
751 )
752 washMjds = pexConfig.ListField(
753 doc="Mirror wash MJDs",
754 dtype=float,
755 default=(0.0,),
756 )
757 epochMjds = pexConfig.ListField(
758 doc="Epoch boundaries in MJD",
759 dtype=float,
760 default=(0.0,),
761 )
762 minObsPerBand = pexConfig.Field(
763 doc="Minimum good observations per band",
764 dtype=int,
765 default=2,
766 )
767 # TODO: When DM-16511 is done, it will be possible to get the
768 # telescope latitude directly from the camera.
769 latitude = pexConfig.Field(
770 doc="Observatory latitude",
771 dtype=float,
772 default=None,
773 )
774 mirrorArea = pexConfig.Field(
775 doc="Mirror area (square meters) of telescope. If not set, will "
776 "try to estimate from camera.telescopeDiameter.",
777 dtype=float,
778 default=None,
779 optional=True,
780 )
781 cameraGain = pexConfig.Field(
782 doc="Average camera gain. If not set, will use the median of the "
783 "camera model/detector/amplifier gains.",
784 dtype=float,
785 default=None,
786 optional=True,
787 )
788 defaultCameraOrientation = pexConfig.Field(
789 doc="Default camera orientation for QA plots.",
790 dtype=float,
791 default=None,
792 )
793 brightObsGrayMax = pexConfig.Field(
794 doc="Maximum gray extinction to be considered bright observation",
795 dtype=float,
796 default=0.15,
797 )
798 minStarPerCcd = pexConfig.Field(
799 doc=("Minimum number of good stars per CCD to be used in calibration fit. "
800 "CCDs with fewer stars will have their calibration estimated from other "
801 "CCDs in the same visit, with zeropoint error increased accordingly."),
802 dtype=int,
803 default=5,
804 )
805 minCcdPerExp = pexConfig.Field(
806 doc=("Minimum number of good CCDs per exposure/visit to be used in calibration fit. "
807 "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."),
808 dtype=int,
809 default=5,
810 )
811 maxCcdGrayErr = pexConfig.Field(
812 doc="Maximum error on CCD gray offset to be considered photometric",
813 dtype=float,
814 default=0.05,
815 )
816 minStarPerExp = pexConfig.Field(
817 doc=("Minimum number of good stars per exposure/visit to be used in calibration fit. "
818 "Visits with fewer good stars will have CCD zeropoints estimated where possible."),
819 dtype=int,
820 default=600,
821 )
822 minExpPerNight = pexConfig.Field(
823 doc="Minimum number of good exposures/visits to consider a partly photometric night",
824 dtype=int,
825 default=10,
826 )
827 expGrayInitialCut = pexConfig.Field(
828 doc=("Maximum exposure/visit gray value for initial selection of possible photometric "
829 "observations."),
830 dtype=float,
831 default=-0.25,
832 )
833 expFwhmCutDict = pexConfig.DictField(
834 doc=("Per-band specification on maximum exposure FWHM (arcseconds) that will "
835 "be considered for the model fit. Exposures with median FWHM larger "
836 "than this threshold will get zeropoints based on matching to good "
837 "stars."),
838 keytype=str,
839 itemtype=float,
840 default={},
841 )
842 expGrayPhotometricCutDict = pexConfig.DictField(
843 doc=("Per-band specification on maximum (negative) achromatic exposure residual "
844 "('gray term') for a visit to be considered photometric. Must have one "
845 "entry per band. Broad-band filters should be -0.05."),
846 keytype=str,
847 itemtype=float,
848 default={},
849 )
850 expGrayHighCutDict = pexConfig.DictField(
851 doc=("Per-band specification on maximum (positive) achromatic exposure residual "
852 "('gray term') for a visit to be considered photometric. Must have one "
853 "entry per band. Broad-band filters should be 0.2."),
854 keytype=str,
855 itemtype=float,
856 default={},
857 )
858 expGrayRecoverCut = pexConfig.Field(
859 doc=("Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. "
860 "Visits with more gray extinction will only get CCD zeropoints if there are "
861 "sufficient star observations (minStarPerCcd) on that CCD."),
862 dtype=float,
863 default=-1.0,
864 )
865 expVarGrayPhotometricCutDict = pexConfig.DictField(
866 doc=("Per-band specification on maximum exposure variance to be considered possibly "
867 "photometric. Must have one entry per band. Broad-band filters should be "
868 "0.0005."),
869 keytype=str,
870 itemtype=float,
871 default={},
872 )
873 expGrayErrRecoverCut = pexConfig.Field(
874 doc=("Maximum exposure gray error to be able to recover bad ccds via interpolation. "
875 "Visits with more gray variance will only get CCD zeropoints if there are "
876 "sufficient star observations (minStarPerCcd) on that CCD."),
877 dtype=float,
878 default=0.05,
879 )
880 aperCorrUsePsfFwhm = pexConfig.Field(
881 doc="Use PSF FWHM for aperture corrections? Otherwise delta_aper statistic is used.",
882 dtype=bool,
883 default=False,
884 )
885 aperCorrPerCcd = pexConfig.Field(
886 doc="Use aperture corrections per-ccd (detector) instead of per-visit?",
887 dtype=bool,
888 default=False,
889 )
890 aperCorrFitNBins = pexConfig.Field(
891 doc=("Number of aperture bins used in aperture correction fit. When set to 0"
892 "no fit will be performed, and the config.aperCorrInputSlopes will be "
893 "used if available."),
894 dtype=int,
895 default=10,
896 )
897 aperCorrInputSlopeDict = pexConfig.DictField(
898 doc=("Per-band specification of aperture correction input slope parameters. These "
899 "are used on the first fit iteration, and aperture correction parameters will "
900 "be updated from the data if config.aperCorrFitNBins > 0. It is recommended "
901 "to set this when there is insufficient data to fit the parameters (e.g. "
902 "tract mode)."),
903 keytype=str,
904 itemtype=float,
905 default={},
906 )
907 sedboundaryterms = pexConfig.ConfigField(
908 doc="Mapping from bands to SED boundary term names used is sedterms.",
909 dtype=SedboundarytermDict,
910 )
911 sedterms = pexConfig.ConfigField(
912 doc="Mapping from terms to bands for fgcm linear SED approximations.",
913 dtype=SedtermDict,
914 )
915 sigFgcmMaxErr = pexConfig.Field(
916 doc="Maximum mag error for fitting sigma_FGCM",
917 dtype=float,
918 default=0.01,
919 )
920 sigFgcmMaxEGrayDict = pexConfig.DictField(
921 doc=("Per-band specification for maximum (absolute) achromatic residual (gray value) "
922 "for observations in sigma_fgcm (raw repeatability). Broad-band filters "
923 "should be 0.05."),
924 keytype=str,
925 itemtype=float,
926 default={},
927 )
928 ccdGrayMaxStarErr = pexConfig.Field(
929 doc=("Maximum error on a star observation to use in ccd gray (achromatic residual) "
930 "computation"),
931 dtype=float,
932 default=0.10,
933 )
934 approxThroughputDict = pexConfig.DictField(
935 doc=("Per-band specification of the approximate overall throughput at the start of "
936 "calibration observations. Must have one entry per band. Typically should "
937 "be 1.0."),
938 keytype=str,
939 itemtype=float,
940 default={},
941 )
942 sigmaCalRange = pexConfig.ListField(
943 doc="Allowed range for systematic error floor estimation",
944 dtype=float,
945 default=(0.001, 0.003),
946 )
947 sigmaCalFitPercentile = pexConfig.ListField(
948 doc="Magnitude percentile range to fit systematic error floor",
949 dtype=float,
950 default=(0.05, 0.15),
951 )
952 sigmaCalPlotPercentile = pexConfig.ListField(
953 doc="Magnitude percentile range to plot systematic error floor",
954 dtype=float,
955 default=(0.05, 0.95),
956 )
957 sigma0Phot = pexConfig.Field(
958 doc="Systematic error floor for all zeropoints",
959 dtype=float,
960 default=0.003,
961 )
962 mapLongitudeRef = pexConfig.Field(
963 doc="Reference longitude for plotting maps",
964 dtype=float,
965 default=0.0,
966 )
967 mapNSide = pexConfig.Field(
968 doc="Healpix nside for plotting maps",
969 dtype=int,
970 default=256,
971 )
972 outfileBase = pexConfig.Field(
973 doc="Filename start for plot output files",
974 dtype=str,
975 default=None,
976 )
977 starColorCuts = pexConfig.ListField(
978 doc=("Encoded star-color cuts (using calibration star colors). "
979 "This is a list with each entry a string of the format "
980 "``band1,band2,low,high`` such that only stars of color "
981 "low < band1 - band2 < high will be used for calibration."),
982 dtype=str,
983 default=("NO_DATA",),
984 )
985 refStarColorCuts = pexConfig.ListField(
986 doc=("Encoded star color cuts specifically to apply to reference stars. "
987 "This is a list with each entry a string of the format "
988 "``band1,band2,low,high`` such that only stars of color "
989 "low < band1 - band2 < high will be used as reference stars."),
990 dtype=str,
991 default=("NO_DATA",),
992 )
993 colorSplitBands = pexConfig.ListField(
994 doc="Band names to use to split stars by color. Must have 2 entries.",
995 dtype=str,
996 length=2,
997 default=('g', 'i'),
998 )
999 modelMagErrors = pexConfig.Field(
1000 doc="Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)",
1001 dtype=bool,
1002 default=True,
1003 )
1004 useQuadraticPwv = pexConfig.Field(
1005 doc="Model PWV with a quadratic term for variation through the night?",
1006 dtype=bool,
1007 default=False,
1008 )
1009 instrumentParsPerBand = pexConfig.Field(
1010 doc=("Model instrumental parameters per band? "
1011 "Otherwise, instrumental parameters (QE changes with time) are "
1012 "shared among all bands."),
1013 dtype=bool,
1014 default=False,
1015 )
1016 instrumentSlopeMinDeltaT = pexConfig.Field(
1017 doc=("Minimum time change (in days) between observations to use in constraining "
1018 "instrument slope."),
1019 dtype=float,
1020 default=20.0,
1021 )
1022 fitMirrorChromaticity = pexConfig.Field(
1023 doc="Fit (intraband) mirror chromatic term?",
1024 dtype=bool,
1025 default=False,
1026 )
1027 fitCcdChromaticityDict = pexConfig.DictField(
1028 doc="Specification on whether to compute first-order quantum efficiency (QE) "
1029 "adjustments. Key is band, and value will be True or False. Any band "
1030 "not explicitly specified will default to False.",
1031 keytype=str,
1032 itemtype=bool,
1033 default={},
1034 )
1035 coatingMjds = pexConfig.ListField(
1036 doc="Mirror coating dates in MJD",
1037 dtype=float,
1038 default=(0.0,),
1039 )
1040 outputStandardsBeforeFinalCycle = pexConfig.Field(
1041 doc="Output standard stars prior to final cycle? Used in debugging.",
1042 dtype=bool,
1043 default=False,
1044 )
1045 outputZeropointsBeforeFinalCycle = pexConfig.Field(
1046 doc="Output standard stars prior to final cycle? Used in debugging.",
1047 dtype=bool,
1048 default=False,
1049 )
1050 useRepeatabilityForExpGrayCutsDict = pexConfig.DictField(
1051 doc=("Per-band specification on whether to use star repeatability (instead of exposures) "
1052 "for computing photometric cuts. Recommended for tract mode or bands with few visits."),
1053 keytype=str,
1054 itemtype=bool,
1055 default={},
1056 )
1057 autoPhotometricCutNSig = pexConfig.Field(
1058 doc=("Number of sigma for automatic computation of (low) photometric cut. "
1059 "Cut is based on exposure gray width (per band), unless "
1060 "useRepeatabilityForExpGrayCuts is set, in which case the star "
1061 "repeatability is used (also per band)."),
1062 dtype=float,
1063 default=3.0,
1064 )
1065 autoHighCutNSig = pexConfig.Field(
1066 doc=("Number of sigma for automatic computation of (high) outlier cut. "
1067 "Cut is based on exposure gray width (per band), unless "
1068 "useRepeatabilityForExpGrayCuts is set, in which case the star "
1069 "repeatability is used (also per band)."),
1070 dtype=float,
1071 default=4.0,
1072 )
1073 quietMode = pexConfig.Field(
1074 doc="Be less verbose with logging.",
1075 dtype=bool,
1076 default=False,
1077 )
1078 doPlots = pexConfig.Field(
1079 doc="Make fgcm QA plots.",
1080 dtype=bool,
1081 default=True,
1082 )
1083 doPlotsBeforeFinalCycles = pexConfig.Field(
1084 doc="Make fgcm QA plots before the final two fit cycles? This only applies in"
1085 "multi-cycle mode, and if doPlots is True. These are typically the most"
1086 "important QA plots to inspect.",
1087 dtype=bool,
1088 default=False,
1089 )
1090 randomSeed = pexConfig.Field(
1091 doc="Random seed for fgcm for consistency in tests.",
1092 dtype=int,
1093 default=None,
1094 optional=True,
1095 )
1096 deltaAperFitMinNgoodObs = pexConfig.Field(
1097 doc="Minimum number of good observations to use mean delta-aper values in fits.",
1098 dtype=int,
1099 default=2,
1100 )
1101 deltaAperFitPerCcdNx = pexConfig.Field(
1102 doc=("Number of x bins per ccd when computing delta-aper background offsets. "
1103 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
1104 dtype=int,
1105 default=10,
1106 )
1107 deltaAperFitPerCcdNy = pexConfig.Field(
1108 doc=("Number of y bins per ccd when computing delta-aper background offsets. "
1109 "Only used when ``doComputeDeltaAperPerCcd`` is True."),
1110 dtype=int,
1111 default=10,
1112 )
1113 deltaAperFitSpatialNside = pexConfig.Field(
1114 doc="Healpix nside to compute spatial delta-aper background offset maps.",
1115 dtype=int,
1116 default=64,
1117 )
1118 deltaAperInnerRadiusArcsec = pexConfig.Field(
1119 doc=("Inner radius used to compute deltaMagAper (arcseconds). "
1120 "Must be positive and less than ``deltaAperOuterRadiusArcsec`` if "
1121 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
1122 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
1123 dtype=float,
1124 default=0.0,
1125 )
1126 deltaAperOuterRadiusArcsec = pexConfig.Field(
1127 doc=("Outer radius used to compute deltaMagAper (arcseconds). "
1128 "Must be positive and greater than ``deltaAperInnerRadiusArcsec`` if "
1129 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, "
1130 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."),
1131 dtype=float,
1132 default=0.0,
1133 )
1134 doComputeDeltaAperPerVisit = pexConfig.Field(
1135 doc=("Do the computation of delta-aper background offsets per visit? "
1136 "Note: this option can be very slow when there are many visits."),
1137 dtype=bool,
1138 default=False,
1139 )
1140 doComputeDeltaAperPerStar = pexConfig.Field(
1141 doc="Do the computation of delta-aper mean values per star?",
1142 dtype=bool,
1143 default=True,
1144 )
1145 doComputeDeltaAperMap = pexConfig.Field(
1146 doc=("Do the computation of delta-aper spatial maps? "
1147 "This is only used if ``doComputeDeltaAperPerStar`` is True,"),
1148 dtype=bool,
1149 default=False,
1150 )
1151 doComputeDeltaAperPerCcd = pexConfig.Field(
1152 doc="Do the computation of per-ccd delta-aper background offsets?",
1153 dtype=bool,
1154 default=False,
1155 )
1157 def validate(self):
1158 super().validate()
1160 if self.connections.previousCycleNumber != str(self.cycleNumber - 1):
1161 msg = "cycleNumber in template must be connections.previousCycleNumber + 1"
1162 raise RuntimeError(msg)
1163 if self.connections.cycleNumber != str(self.cycleNumber):
1164 msg = "cycleNumber in template must be equal to connections.cycleNumber"
1165 raise RuntimeError(msg)
1167 for band in self.fitBands:
1168 if band not in self.bands:
1169 msg = 'fitBand %s not in bands' % (band)
1170 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg)
1171 for band in self.requiredBands:
1172 if band not in self.bands:
1173 msg = 'requiredBand %s not in bands' % (band)
1174 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg)
1175 for band in self.colorSplitBands:
1176 if band not in self.bands:
1177 msg = 'colorSplitBand %s not in bands' % (band)
1178 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg)
1179 for band in self.bands:
1180 if band not in self.superStarSubCcdDict:
1181 msg = 'band %s not in superStarSubCcdDict' % (band)
1182 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict,
1183 self, msg)
1184 if band not in self.ccdGraySubCcdDict:
1185 msg = 'band %s not in ccdGraySubCcdDict' % (band)
1186 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict,
1187 self, msg)
1188 if band not in self.expGrayPhotometricCutDict:
1189 msg = 'band %s not in expGrayPhotometricCutDict' % (band)
1190 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict,
1191 self, msg)
1192 if band not in self.expGrayHighCutDict:
1193 msg = 'band %s not in expGrayHighCutDict' % (band)
1194 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict,
1195 self, msg)
1196 if band not in self.expVarGrayPhotometricCutDict:
1197 msg = 'band %s not in expVarGrayPhotometricCutDict' % (band)
1198 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict,
1199 self, msg)
1200 if band not in self.sigFgcmMaxEGrayDict:
1201 msg = 'band %s not in sigFgcmMaxEGrayDict' % (band)
1202 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict,
1203 self, msg)
1204 if band not in self.approxThroughputDict:
1205 msg = 'band %s not in approxThroughputDict' % (band)
1206 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict,
1207 self, msg)
1208 if band not in self.useRepeatabilityForExpGrayCutsDict:
1209 msg = 'band %s not in useRepeatabilityForExpGrayCutsDict' % (band)
1210 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
1211 self, msg)
1213 if self.doComputeDeltaAperPerVisit or self.doComputeDeltaAperMap \
1214 or self.doComputeDeltaAperPerCcd:
1215 if self.deltaAperInnerRadiusArcsec <= 0.0:
1216 msg = 'deltaAperInnerRadiusArcsec must be positive if deltaAper computations are turned on.'
1217 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperInnerRadiusArcsec,
1218 self, msg)
1219 if self.deltaAperOuterRadiusArcsec <= 0.0:
1220 msg = 'deltaAperOuterRadiusArcsec must be positive if deltaAper computations are turned on.'
1221 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
1222 self, msg)
1223 if self.deltaAperOuterRadiusArcsec <= self.deltaAperInnerRadiusArcsec:
1224 msg = ('deltaAperOuterRadiusArcsec must be greater than deltaAperInnerRadiusArcsec if '
1225 'deltaAper computations are turned on.')
1226 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec,
1227 self, msg)
1230class FgcmFitCycleTask(pipeBase.PipelineTask):
1231 """
1232 Run Single fit cycle for FGCM global calibration
1233 """
1235 ConfigClass = FgcmFitCycleConfig
1236 _DefaultName = "fgcmFitCycle"
1238 def __init__(self, initInputs=None, **kwargs):
1239 super().__init__(**kwargs)
1241 self.multiCycleLoaded = False
1243 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1244 camera = butlerQC.get(inputRefs.camera)
1246 nCore = butlerQC.resources.num_cores
1248 handleDict = {}
1250 handleDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
1251 handleDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
1253 if self.config.useParquetCatalogFormat:
1254 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservationsParquet)
1255 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIdsParquet)
1256 if self.config.doReferenceCalibration:
1257 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStarsParquet)
1258 else:
1259 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations)
1260 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds)
1261 handleDict['fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices)
1262 if self.config.doReferenceCalibration:
1263 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars)
1264 if self.config.cycleNumber > 0:
1265 handleDict['fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput)
1266 handleDict['fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput)
1268 fgcmDatasetDict = None
1269 fgcmMultiCycleObjects = None
1270 if self.config.doMultipleCycles:
1271 # Run multiple cycles at once.
1272 config = copy.copy(self.config)
1273 config.update(cycleNumber=0)
1274 for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1):
1275 if cycle == self.config.multipleCyclesFinalCycleNumber:
1276 config.update(isFinalCycle=True)
1278 if cycle > 0:
1279 handleDict['fgcmFlaggedStars'] = fgcmDatasetDict['fgcmFlaggedStars']
1280 handleDict['fgcmFitParameters'] = fgcmDatasetDict['fgcmFitParameters']
1282 # Set up plot outputs.
1283 # Note that because the plot outputConnections are skipped in
1284 # the connections class if config.doPlot is False, this will be
1285 # a no-op in that case.
1286 plotHandleDict = {}
1287 for outputRefName in outputRefs.keys():
1288 if outputRefName.endswith("Plot") and f"Cycle{cycle}" in outputRefName:
1289 refs = getattr(outputRefs, outputRefName)
1290 if isinstance(refs, (tuple, list)):
1291 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions:
1292 for ref in refs:
1293 physical_filter = ref.dataId["physical_filter"]
1294 detector = ref.dataId["detector"]
1295 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}"
1296 plotHandleDict[handleDictKey] = ref
1297 elif "physical_filter" in refs[0].dimensions:
1298 for ref in refs:
1299 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}"
1300 plotHandleDict[handleDictKey] = ref
1301 elif "band" in refs[0].dimensions:
1302 for ref in refs:
1303 handleDictKey = f"{outputRefName}_{ref.dataId['band']}"
1304 plotHandleDict[handleDictKey] = ref
1305 else:
1306 plotHandleDict[outputRefName] = refs
1308 fgcmDatasetDict, config, fgcmMultiCycleObjects = self._fgcmFitCycle(
1309 camera,
1310 handleDict,
1311 butlerQC=butlerQC,
1312 plotHandleDict=plotHandleDict,
1313 config=config,
1314 nCore=nCore,
1315 fgcmMultiCycleObjects=fgcmMultiCycleObjects,
1316 )
1317 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'],
1318 getattr(outputRefs, f'fgcm_Cycle{cycle}_FitParameters'))
1319 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'],
1320 getattr(outputRefs, f'fgcm_Cycle{cycle}_FlaggedStars'))
1321 butlerQC.put(config,
1322 getattr(outputRefs, f'fgcm_Cycle{cycle}_OutputConfig'))
1323 if self.outputZeropoints:
1324 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'],
1325 getattr(outputRefs, f'fgcm_Cycle{cycle}_Zeropoints'))
1326 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'],
1327 getattr(outputRefs, f'fgcm_Cycle{cycle}_AtmosphereParameters'))
1328 if self.outputStandards:
1329 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'],
1330 getattr(outputRefs, f'fgcm_Cycle{cycle}_StandardStars'))
1331 else:
1332 # Run a single cycle
1334 # Set up plot outputs.
1335 # Note that nothing will go in the dict if doPlots is False.
1336 plotHandleDict = {}
1337 for outputRefName in outputRefs.keys():
1338 if outputRefName.endswith("Plot") and f"Cycle{self.config.cycleNumber}" in outputRefName:
1339 refs = getattr(outputRefs, outputRefName)
1340 if isinstance(refs, (tuple, list)):
1341 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions:
1342 for ref in refs:
1343 physical_filter = ref.dataId["physical_filter"]
1344 detector = ref.dataId["detector"]
1345 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}"
1346 plotHandleDict[handleDictKey] = ref
1347 elif "physical_filter" in refs[0].dimensions:
1348 for ref in refs:
1349 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}"
1350 plotHandleDict[handleDictKey] = ref
1351 elif "band" in refs[0].dimensions:
1352 for ref in refs:
1353 handleDictKey = f"{outputRefName}_{ref.dataId['band']}"
1354 plotHandleDict[handleDictKey] = ref
1355 else:
1356 plotHandleDict[outputRefName] = refs
1358 fgcmDatasetDict, _ = self._fgcmFitCycle(
1359 camera,
1360 handleDict,
1361 nCore=nCore,
1362 butlerQC=butlerQC,
1363 plotHandleDict=plotHandleDict,
1364 multiCycle=False,
1365 )
1367 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], outputRefs.fgcmFitParameters)
1368 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars)
1369 if self.outputZeropoints:
1370 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], outputRefs.fgcmZeropoints)
1371 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters)
1372 if self.outputStandards:
1373 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], outputRefs.fgcmStandardStars)
1375 def _fgcmFitCycle(
1376 self,
1377 camera,
1378 handleDict,
1379 butlerQC=None,
1380 plotHandleDict=None,
1381 config=None,
1382 nCore=1,
1383 multiCycle=True,
1384 fgcmMultiCycleObjects=None,
1385 ):
1386 """
1387 Run the fit cycle
1389 Parameters
1390 ----------
1391 camera : `lsst.afw.cameraGeom.Camera`
1392 handleDict : `dict`
1393 All handles are `lsst.daf.butler.DeferredDatasetHandle`
1394 handle dictionary with keys:
1396 ``"fgcmLookUpTable"``
1397 handle for the FGCM look-up table.
1398 ``"fgcmVisitCatalog"``
1399 handle for visit summary catalog.
1400 ``"fgcmStarObservations"``
1401 handle for star observation catalog.
1402 ``"fgcmStarIds"``
1403 handle for star id catalog.
1404 ``"fgcmStarIndices"``
1405 handle for star index catalog.
1406 ``"fgcmReferenceStars"``
1407 handle for matched reference star catalog.
1408 ``"fgcmFlaggedStars"``
1409 handle for flagged star catalog.
1410 ``"fgcmFitParameters"``
1411 handle for fit parameter catalog.
1412 butlerQC : `lsst.pipe.base.QuantumContext`, optional
1413 Quantum context used for serializing plots.
1414 plotHandleDict : `dict` [`str`, `lsst.daf.butler.DatasetRef`], optional
1415 Dictionary of plot dataset refs, keyed by plot name.
1416 config : `lsst.pex.config.Config`, optional
1417 Configuration to use to override self.config.
1418 nCore : `int`, optional
1419 Number of cores to use during fitting.
1420 multiCycle : `bool`, optional
1421 Is this part of a multicycle run?
1422 fgcmMultiCycleObjects : `dict`, optional
1423 Dictionary of multi cycle objects from previous cycle.
1424 This should be blank on first call; subsequent calls will use
1425 the dictionary returned from the previous call. The keys include
1426 ``fgcmFitCycle`` (the primary fgcm fit object) and ``fgcmExpInfo``
1427 (information about each exposure/visit).
1429 Returns
1430 -------
1431 fgcmDatasetDict : `dict`
1432 Dictionary of datasets to persist.
1433 config : `lsst.pex.config.Config`
1434 Configuration object for next cycle.
1435 fgcmMultiCycleObjects : `dict`, optional
1436 Dictionary of multi cycle objects; returned if multicycle is True.
1437 To be passed to next call of this method (see Parameters above).
1438 """
1439 if config is not None:
1440 _config = config
1441 else:
1442 _config = self.config
1444 if multiCycle and not self.multiCycleLoaded:
1445 fgcmMultiCycleObjects = {}
1447 # Set defaults on whether to output standards and zeropoints
1448 self.maxIter = _config.maxIterBeforeFinalCycle
1449 self.outputStandards = _config.outputStandardsBeforeFinalCycle
1450 self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle
1451 self.resetFitParameters = True
1453 if _config.isFinalCycle:
1454 # This is the final fit cycle, so we do not want to reset fit
1455 # parameters, we want to run a final "clean-up" with 0 fit iterations,
1456 # and we always want to output standards and zeropoints
1457 self.maxIter = 0
1458 self.outputStandards = True
1459 self.outputZeropoints = True
1460 self.resetFitParameters = False
1462 # Check if we want to do plots.
1463 doPlots = _config.doPlots
1464 if doPlots and multiCycle:
1465 if _config.cycleNumber < (_config.multipleCyclesFinalCycleNumber - 1) \
1466 and not _config.doPlotsBeforeFinalCycles:
1467 doPlots = False
1469 if not self.multiCycleLoaded:
1470 lutCat = handleDict['fgcmLookUpTable'].get()
1471 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
1472 dict(_config.physicalFilterMap))
1473 del lutCat
1475 filterNames = lutIndexVals[0]['FILTERNAMES']
1477 configDict = makeConfigDict(
1478 _config,
1479 self.log,
1480 camera,
1481 self.maxIter,
1482 self.resetFitParameters,
1483 self.outputZeropoints,
1484 filterNames,
1485 nCore=nCore,
1486 doPlots=doPlots,
1487 )
1489 visitCat = handleDict['fgcmVisitCatalog'].get()
1490 fgcmExpInfo = translateVisitCatalog(visitCat)
1491 del visitCat
1493 if len(camera) == fgcmLut.nCCD:
1494 useScienceDetectors = False
1495 else:
1496 # If the LUT has a different number of detectors than
1497 # the camera, then we only want to use science detectors
1498 # in the focal plane projector.
1499 useScienceDetectors = True
1501 self.log.info("Loading focal plane projector")
1502 focalPlaneProjector = FocalPlaneProjector(
1503 camera,
1504 self.config.defaultCameraOrientation,
1505 useScienceDetectors=useScienceDetectors,
1506 )
1508 noFitsDict = {
1509 'lutIndex': lutIndexVals,
1510 'lutStd': lutStd,
1511 'expInfo': fgcmExpInfo,
1512 'focalPlaneProjector': focalPlaneProjector,
1513 }
1515 fgcmFitCycle = fgcm.FgcmFitCycle(
1516 configDict,
1517 useFits=False,
1518 noFitsDict=noFitsDict,
1519 noOutput=True,
1520 butlerQC=butlerQC,
1521 plotHandleDict=plotHandleDict,
1522 )
1524 if fgcmFitCycle.initialCycle:
1525 # cycle = 0, initial cycle
1526 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(
1527 fgcmFitCycle.fgcmConfig,
1528 fgcmLut,
1529 fgcmExpInfo,
1530 butlerQC=butlerQC,
1531 plotHandleDict=plotHandleDict,
1532 )
1533 else:
1534 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog):
1535 parCat = handleDict['fgcmFitParameters']
1536 else:
1537 parCat = handleDict['fgcmFitParameters'].get()
1538 inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1539 del parCat
1541 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(
1542 fgcmFitCycle.fgcmConfig,
1543 fgcmExpInfo,
1544 inParInfo,
1545 inParams,
1546 inSuperStar,
1547 butlerQC=butlerQC,
1548 plotHandleDict=plotHandleDict,
1549 )
1551 fgcmStars = fgcm.FgcmStars(
1552 fgcmFitCycle.fgcmConfig,
1553 butlerQC=butlerQC,
1554 plotHandleDict=plotHandleDict,
1555 )
1557 starObs = handleDict['fgcmStarObservations'].get()
1558 starIds = handleDict['fgcmStarIds'].get()
1559 if not self.config.useParquetCatalogFormat:
1560 starIndices = handleDict['fgcmStarIndices'].get()
1561 else:
1562 starIndices = None
1564 # grab the flagged stars if available
1565 if 'fgcmFlaggedStars' in handleDict:
1566 if isinstance(handleDict['fgcmFlaggedStars'], afwTable.BaseCatalog):
1567 flaggedStars = handleDict['fgcmFlaggedStars']
1568 else:
1569 flaggedStars = handleDict['fgcmFlaggedStars'].get()
1570 flagId = flaggedStars['objId'][:]
1571 flagFlag = flaggedStars['objFlag'][:]
1573 del flaggedStars
1574 elif self.config.useParquetCatalogFormat:
1575 # If we are using the parquet catalog format, then that means that
1576 # reserved stars have already been flagged. We extract the flags here
1577 # to input to fgcm, which will then be persisted (with additional
1578 # quality flags) as the fgcmFlaggedStars datatype in subsequent
1579 # fit cycles.
1580 flagged = (starIds['obj_flag'] > 0)
1581 flagId = starIds['fgcm_id'][flagged]
1582 flagFlag = starIds['obj_flag'][flagged]
1583 else:
1584 flagId = None
1585 flagFlag = None
1587 if _config.doReferenceCalibration:
1588 refStars = handleDict['fgcmReferenceStars'].get()
1590 refMag, refMagErr = extractReferenceMags(
1591 refStars,
1592 _config.bands,
1593 _config.physicalFilterMap,
1594 )
1596 refId = refStars['fgcm_id'][:]
1597 else:
1598 refStars = None
1599 refId = None
1600 refMag = None
1601 refMagErr = None
1603 # match star observations to visits
1604 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
1605 # actually be transferred into fgcm using the indexing below.
1606 if self.config.useParquetCatalogFormat:
1607 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'])
1608 else:
1609 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'][starIndices['obsIndex']])
1611 # The fgcmStars.loadStars method will copy all the star
1612 # information into special shared memory objects that will
1613 # not blow up the memory usage when used with python
1614 # multiprocessing. Once all the numbers are copied, it is
1615 # necessary to release all references to the objects that
1616 # previously stored the data to ensure that the garbage
1617 # collector can clear the memory, and ensure that this
1618 # memory is not copied when multiprocessing kicks in.
1620 if self.config.useParquetCatalogFormat:
1621 # Note that the ra/dec coordinates for the parquet format are in
1622 # degrees, which is what fgcm expects.
1623 self.log.info("Loading stars and such")
1624 fgcmStars.loadStars(
1625 fgcmPars,
1626 starObs['visit'][:],
1627 starObs['detector'][:],
1628 starObs['ra'][:],
1629 starObs['dec'][:],
1630 starObs['inst_mag'][:],
1631 starObs['inst_mag_err'][:],
1632 fgcmExpInfo['FILTERNAME'][visitIndex],
1633 starIds['fgcm_id'][:],
1634 starIds['ra'][:],
1635 starIds['dec'][:],
1636 starIds['obs_arr_index'][:],
1637 starIds['n_obs'][:],
1638 obsX=starObs['x'][:],
1639 obsY=starObs['y'][:],
1640 obsDeltaMagBkg=starObs['delta_mag_bkg'][:],
1641 obsDeltaAper=starObs['delta_mag_aper'][:],
1642 refID=refId,
1643 refMag=refMag,
1644 refMagErr=refMagErr,
1645 flagID=flagId,
1646 flagFlag=flagFlag,
1647 computeNobs=True,
1648 objIDAlternate=starIds['isolated_star_id'],
1649 )
1650 else:
1651 # We determine the conversion from the native units
1652 # (typically radians) to degrees for the first star.
1653 # This allows us to treat coord_ra/coord_dec as
1654 # numpy arrays rather than Angles, which would we
1655 # approximately 600x slower.
1656 conv = starObs[0]['ra'].asDegrees() / float(starObs[0]['ra'])
1658 fgcmStars.loadStars(
1659 fgcmPars,
1660 starObs['visit'][starIndices['obsIndex']],
1661 starObs['ccd'][starIndices['obsIndex']],
1662 starObs['ra'][starIndices['obsIndex']] * conv,
1663 starObs['dec'][starIndices['obsIndex']] * conv,
1664 starObs['instMag'][starIndices['obsIndex']],
1665 starObs['instMagErr'][starIndices['obsIndex']],
1666 fgcmExpInfo['FILTERNAME'][visitIndex],
1667 starIds['fgcm_id'][:],
1668 starIds['ra'][:],
1669 starIds['dec'][:],
1670 starIds['obsArrIndex'][:],
1671 starIds['nObs'][:],
1672 obsX=starObs['x'][starIndices['obsIndex']],
1673 obsY=starObs['y'][starIndices['obsIndex']],
1674 obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']],
1675 obsDeltaAper=starObs['deltaMagAper'][starIndices['obsIndex']],
1676 psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']],
1677 refID=refId,
1678 refMag=refMag,
1679 refMagErr=refMagErr,
1680 flagID=flagId,
1681 flagFlag=flagFlag,
1682 computeNobs=True,
1683 )
1685 # Release all references to temporary objects holding star data
1686 # (see above)
1687 del starObs
1688 del starIds
1689 del starIndices
1690 del flagId
1691 del flagFlag
1692 del refStars
1693 del refId
1694 del refMag
1695 del refMagErr
1697 fgcmFitCycle.setLUT(fgcmLut)
1698 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
1700 if multiCycle:
1701 fgcmMultiCycleObjects["fgcmFitCycle"] = fgcmFitCycle
1702 fgcmMultiCycleObjects["fgcmExpInfo"] = fgcmExpInfo
1704 fgcmFitCycle.setPars(fgcmPars)
1705 fgcmFitCycle.finishSetup()
1706 else:
1707 # We already have an fgcmFitCycle object loaded.
1708 fgcmFitCycle = fgcmMultiCycleObjects["fgcmFitCycle"]
1709 fgcmExpInfo = fgcmMultiCycleObjects["fgcmExpInfo"]
1711 # Update configs.
1712 self.log.info("Updating fgcm configuration for cycle %d", _config.cycleNumber)
1713 fgcmFitCycle.updateConfigNextCycle(
1714 _config.cycleNumber,
1715 maxIter=self.maxIter,
1716 resetParameters=self.resetFitParameters,
1717 outputStandards=self.outputStandards,
1718 outputZeropoints=self.outputZeropoints,
1719 freezeStdAtmosphere=_config.freezeStdAtmosphere,
1720 expGrayPhotometricCutDict=dict(_config.expGrayPhotometricCutDict),
1721 expGrayHighCutDict=dict(_config.expGrayHighCutDict),
1722 doPlots=doPlots,
1723 )
1725 # Reload parameters.
1726 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog):
1727 parCat = handleDict['fgcmFitParameters']
1728 else:
1729 parCat = handleDict['fgcmFitParameters'].get()
1730 inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1731 del parCat
1733 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(
1734 fgcmFitCycle.fgcmConfig,
1735 fgcmExpInfo,
1736 inParInfo,
1737 inParams,
1738 inSuperStar,
1739 butlerQC=butlerQC,
1740 plotHandleDict=plotHandleDict,
1741 )
1743 # Reset star magnitudes.
1744 fgcmFitCycle.fgcmStars.reloadStarMagnitudes()
1745 fgcmFitCycle.fgcmStars.computeAllNobs(fgcmPars)
1747 fgcmFitCycle.setPars(fgcmPars)
1748 fgcmFitCycle.finishReset(butlerQC=butlerQC, plotHandleDict=plotHandleDict)
1750 fgcmFitCycle.run()
1752 ##################
1753 # Persistance
1754 ##################
1756 fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle)
1758 # Output the config for the next cycle
1759 # We need to make a copy since the input one has been frozen
1761 updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i]) for
1762 i, b in enumerate(_config.bands)}
1763 updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i]) for
1764 i, band in enumerate(_config.bands)}
1766 outConfig = copy.copy(_config)
1767 outConfig.update(cycleNumber=(_config.cycleNumber + 1),
1768 precomputeSuperStarInitialCycle=False,
1769 freezeStdAtmosphere=False,
1770 expGrayPhotometricCutDict=updatedPhotometricCutDict,
1771 expGrayHighCutDict=updatedHighCutDict)
1773 outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber),
1774 cycleNumber=str(_config.cycleNumber + 1))
1776 if not multiCycle:
1777 configFileName = '%s_cycle%02d_config.py' % (outConfig.outfileBase,
1778 outConfig.cycleNumber)
1779 outConfig.save(configFileName)
1781 if _config.isFinalCycle == 1:
1782 # We are done, ready to output products
1783 self.log.info("Everything is in place to run fgcmOutputProducts.py")
1784 else:
1785 self.log.info("Saved config for next cycle to %s" % (configFileName))
1786 self.log.info("Be sure to look at:")
1787 self.log.info(" config.expGrayPhotometricCut")
1788 self.log.info(" config.expGrayHighCut")
1789 self.log.info("If you are satisfied with the fit, please set:")
1790 self.log.info(" config.isFinalCycle = True")
1792 if not multiCycle or config.isFinalCycle:
1793 fgcmFitCycle.freeSharedMemory()
1795 if multiCycle:
1796 self.multiCycleLoaded = True
1798 return fgcmDatasetDict, outConfig, fgcmMultiCycleObjects
1799 else:
1800 return fgcmDatasetDict, outConfig
1802 def _loadParameters(self, parCat):
1803 """
1804 Load FGCM parameters from a previous fit cycle
1806 Parameters
1807 ----------
1808 parCat : `lsst.afw.table.BaseCatalog`
1809 Parameter catalog in afw table form.
1811 Returns
1812 -------
1813 inParInfo: `numpy.ndarray`
1814 Numpy array parameter information formatted for input to fgcm
1815 inParameters: `numpy.ndarray`
1816 Numpy array parameter values formatted for input to fgcm
1817 inSuperStar: `numpy.array`
1818 Superstar flat formatted for input to fgcm
1819 """
1820 parLutFilterNames = np.array(parCat[0]['lutFilterNames'].split(','))
1821 parFitBands = np.array(parCat[0]['fitBands'].split(','))
1823 inParInfo = np.zeros(1, dtype=[('NCCD', 'i4'),
1824 ('LUTFILTERNAMES', parLutFilterNames.dtype.str,
1825 (parLutFilterNames.size, )),
1826 ('FITBANDS', parFitBands.dtype.str, (parFitBands.size, )),
1827 ('LNTAUUNIT', 'f8'),
1828 ('LNTAUSLOPEUNIT', 'f8'),
1829 ('ALPHAUNIT', 'f8'),
1830 ('LNPWVUNIT', 'f8'),
1831 ('LNPWVSLOPEUNIT', 'f8'),
1832 ('LNPWVQUADRATICUNIT', 'f8'),
1833 ('LNPWVGLOBALUNIT', 'f8'),
1834 ('O3UNIT', 'f8'),
1835 ('QESYSUNIT', 'f8'),
1836 ('FILTEROFFSETUNIT', 'f8'),
1837 ('HASEXTERNALPWV', 'i2'),
1838 ('HASEXTERNALTAU', 'i2')])
1839 inParInfo['NCCD'] = parCat['nCcd']
1840 inParInfo['LUTFILTERNAMES'][:] = parLutFilterNames
1841 inParInfo['FITBANDS'][:] = parFitBands
1842 inParInfo['HASEXTERNALPWV'] = parCat['hasExternalPwv']
1843 inParInfo['HASEXTERNALTAU'] = parCat['hasExternalTau']
1845 inParams = np.zeros(1, dtype=[('PARALPHA', 'f8', (parCat['parAlpha'].size, )),
1846 ('PARO3', 'f8', (parCat['parO3'].size, )),
1847 ('PARLNTAUINTERCEPT', 'f8',
1848 (parCat['parLnTauIntercept'].size, )),
1849 ('PARLNTAUSLOPE', 'f8',
1850 (parCat['parLnTauSlope'].size, )),
1851 ('PARLNPWVINTERCEPT', 'f8',
1852 (parCat['parLnPwvIntercept'].size, )),
1853 ('PARLNPWVSLOPE', 'f8',
1854 (parCat['parLnPwvSlope'].size, )),
1855 ('PARLNPWVQUADRATIC', 'f8',
1856 (parCat['parLnPwvQuadratic'].size, )),
1857 ('PARQESYSINTERCEPT', 'f8',
1858 (parCat['parQeSysIntercept'].size, )),
1859 ('COMPQESYSSLOPE', 'f8',
1860 (parCat['compQeSysSlope'].size, )),
1861 ('PARFILTEROFFSET', 'f8',
1862 (parCat['parFilterOffset'].size, )),
1863 ('PARFILTEROFFSETFITFLAG', 'i2',
1864 (parCat['parFilterOffsetFitFlag'].size, )),
1865 ('PARRETRIEVEDLNPWVSCALE', 'f8'),
1866 ('PARRETRIEVEDLNPWVOFFSET', 'f8'),
1867 ('PARRETRIEVEDLNPWVNIGHTLYOFFSET', 'f8',
1868 (parCat['parRetrievedLnPwvNightlyOffset'].size, )),
1869 ('COMPABSTHROUGHPUT', 'f8',
1870 (parCat['compAbsThroughput'].size, )),
1871 ('COMPREFOFFSET', 'f8',
1872 (parCat['compRefOffset'].size, )),
1873 ('COMPREFSIGMA', 'f8',
1874 (parCat['compRefSigma'].size, )),
1875 ('COMPMIRRORCHROMATICITY', 'f8',
1876 (parCat['compMirrorChromaticity'].size, )),
1877 ('MIRRORCHROMATICITYPIVOT', 'f8',
1878 (parCat['mirrorChromaticityPivot'].size, )),
1879 ('COMPCCDCHROMATICITY', 'f8',
1880 (parCat['compCcdChromaticity'].size, )),
1881 ('COMPMEDIANSEDSLOPE', 'f8',
1882 (parCat['compMedianSedSlope'].size, )),
1883 ('COMPAPERCORRPIVOT', 'f8',
1884 (parCat['compAperCorrPivot'].size, )),
1885 ('COMPAPERCORRSLOPE', 'f8',
1886 (parCat['compAperCorrSlope'].size, )),
1887 ('COMPAPERCORRSLOPEERR', 'f8',
1888 (parCat['compAperCorrSlopeErr'].size, )),
1889 ('COMPAPERCORRRANGE', 'f8',
1890 (parCat['compAperCorrRange'].size, )),
1891 ('COMPMODELERREXPTIMEPIVOT', 'f8',
1892 (parCat['compModelErrExptimePivot'].size, )),
1893 ('COMPMODELERRFWHMPIVOT', 'f8',
1894 (parCat['compModelErrFwhmPivot'].size, )),
1895 ('COMPMODELERRSKYPIVOT', 'f8',
1896 (parCat['compModelErrSkyPivot'].size, )),
1897 ('COMPMODELERRPARS', 'f8',
1898 (parCat['compModelErrPars'].size, )),
1899 ('COMPEXPGRAY', 'f8',
1900 (parCat['compExpGray'].size, )),
1901 ('COMPVARGRAY', 'f8',
1902 (parCat['compVarGray'].size, )),
1903 ('COMPEXPDELTAMAGBKG', 'f8',
1904 (parCat['compExpDeltaMagBkg'].size, )),
1905 ('COMPNGOODSTARPEREXP', 'i4',
1906 (parCat['compNGoodStarPerExp'].size, )),
1907 ('COMPEXPREFOFFSET', 'f8',
1908 (parCat['compExpRefOffset'].size, )),
1909 ('COMPSIGFGCM', 'f8',
1910 (parCat['compSigFgcm'].size, )),
1911 ('COMPSIGMACAL', 'f8',
1912 (parCat['compSigmaCal'].size, )),
1913 ('COMPRETRIEVEDLNPWV', 'f8',
1914 (parCat['compRetrievedLnPwv'].size, )),
1915 ('COMPRETRIEVEDLNPWVRAW', 'f8',
1916 (parCat['compRetrievedLnPwvRaw'].size, )),
1917 ('COMPRETRIEVEDLNPWVFLAG', 'i2',
1918 (parCat['compRetrievedLnPwvFlag'].size, )),
1919 ('COMPRETRIEVEDTAUNIGHT', 'f8',
1920 (parCat['compRetrievedTauNight'].size, )),
1921 ('COMPEPSILON', 'f8',
1922 (parCat['compEpsilon'].size, )),
1923 ('COMPMEDDELTAAPER', 'f8',
1924 (parCat['compMedDeltaAper'].size, )),
1925 ('COMPGLOBALEPSILON', 'f4',
1926 (parCat['compGlobalEpsilon'].size, )),
1927 ('COMPEPSILONMAP', 'f4',
1928 (parCat['compEpsilonMap'].size, )),
1929 ('COMPEPSILONNSTARMAP', 'i4',
1930 (parCat['compEpsilonNStarMap'].size, )),
1931 ('COMPEPSILONCCDMAP', 'f4',
1932 (parCat['compEpsilonCcdMap'].size, )),
1933 ('COMPEPSILONCCDNSTARMAP', 'i4',
1934 (parCat['compEpsilonCcdNStarMap'].size, ))])
1936 inParams['PARALPHA'][:] = parCat['parAlpha'][0, :]
1937 inParams['PARO3'][:] = parCat['parO3'][0, :]
1938 inParams['PARLNTAUINTERCEPT'][:] = parCat['parLnTauIntercept'][0, :]
1939 inParams['PARLNTAUSLOPE'][:] = parCat['parLnTauSlope'][0, :]
1940 inParams['PARLNPWVINTERCEPT'][:] = parCat['parLnPwvIntercept'][0, :]
1941 inParams['PARLNPWVSLOPE'][:] = parCat['parLnPwvSlope'][0, :]
1942 inParams['PARLNPWVQUADRATIC'][:] = parCat['parLnPwvQuadratic'][0, :]
1943 inParams['PARQESYSINTERCEPT'][:] = parCat['parQeSysIntercept'][0, :]
1944 inParams['COMPQESYSSLOPE'][:] = parCat['compQeSysSlope'][0, :]
1945 inParams['PARFILTEROFFSET'][:] = parCat['parFilterOffset'][0, :]
1946 inParams['PARFILTEROFFSETFITFLAG'][:] = parCat['parFilterOffsetFitFlag'][0, :]
1947 inParams['PARRETRIEVEDLNPWVSCALE'] = parCat['parRetrievedLnPwvScale']
1948 inParams['PARRETRIEVEDLNPWVOFFSET'] = parCat['parRetrievedLnPwvOffset']
1949 inParams['PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat['parRetrievedLnPwvNightlyOffset'][0, :]
1950 inParams['COMPABSTHROUGHPUT'][:] = parCat['compAbsThroughput'][0, :]
1951 inParams['COMPREFOFFSET'][:] = parCat['compRefOffset'][0, :]
1952 inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :]
1953 inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :]
1954 inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :]
1955 inParams['COMPCCDCHROMATICITY'][:] = parCat['compCcdChromaticity'][0, :]
1956 inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :]
1957 inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :]
1958 inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :]
1959 inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :]
1960 inParams['COMPAPERCORRRANGE'][:] = parCat['compAperCorrRange'][0, :]
1961 inParams['COMPMODELERREXPTIMEPIVOT'][:] = parCat['compModelErrExptimePivot'][0, :]
1962 inParams['COMPMODELERRFWHMPIVOT'][:] = parCat['compModelErrFwhmPivot'][0, :]
1963 inParams['COMPMODELERRSKYPIVOT'][:] = parCat['compModelErrSkyPivot'][0, :]
1964 inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :]
1965 inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :]
1966 inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
1967 inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
1968 inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
1969 inParams['COMPEXPREFOFFSET'][:] = parCat['compExpRefOffset'][0, :]
1970 inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
1971 inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
1972 inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :]
1973 inParams['COMPRETRIEVEDLNPWVRAW'][:] = parCat['compRetrievedLnPwvRaw'][0, :]
1974 inParams['COMPRETRIEVEDLNPWVFLAG'][:] = parCat['compRetrievedLnPwvFlag'][0, :]
1975 inParams['COMPRETRIEVEDTAUNIGHT'][:] = parCat['compRetrievedTauNight'][0, :]
1976 inParams['COMPEPSILON'][:] = parCat['compEpsilon'][0, :]
1977 inParams['COMPMEDDELTAAPER'][:] = parCat['compMedDeltaAper'][0, :]
1978 inParams['COMPGLOBALEPSILON'][:] = parCat['compGlobalEpsilon'][0, :]
1979 inParams['COMPEPSILONMAP'][:] = parCat['compEpsilonMap'][0, :]
1980 inParams['COMPEPSILONNSTARMAP'][:] = parCat['compEpsilonNStarMap'][0, :]
1981 inParams['COMPEPSILONCCDMAP'][:] = parCat['compEpsilonCcdMap'][0, :]
1982 inParams['COMPEPSILONCCDNSTARMAP'][:] = parCat['compEpsilonCcdNStarMap'][0, :]
1984 inSuperStar = np.zeros(parCat['superstarSize'][0, :], dtype='f8')
1985 inSuperStar[:, :, :, :] = parCat['superstar'][0, :].reshape(inSuperStar.shape)
1987 return (inParInfo, inParams, inSuperStar)
1989 def _makeFgcmOutputDatasets(self, fgcmFitCycle):
1990 """
1991 Persist FGCM datasets through the butler.
1993 Parameters
1994 ----------
1995 fgcmFitCycle: `lsst.fgcm.FgcmFitCycle`
1996 Fgcm Fit cycle object
1997 """
1998 fgcmDatasetDict = {}
2000 # Save the parameters
2001 parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays()
2003 parSchema = afwTable.Schema()
2005 comma = ','
2006 lutFilterNameString = comma.join([n.decode('utf-8')
2007 for n in parInfo['LUTFILTERNAMES'][0]])
2008 fitBandString = comma.join([n.decode('utf-8')
2009 for n in parInfo['FITBANDS'][0]])
2011 parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat,
2012 lutFilterNameString, fitBandString)
2013 parCat = self._makeParCatalog(parSchema, parInfo, pars,
2014 fgcmFitCycle.fgcmPars.parSuperStarFlat,
2015 lutFilterNameString, fitBandString)
2017 fgcmDatasetDict['fgcmFitParameters'] = parCat
2019 # Save the indices of the flagged stars
2020 # (stars that have been (a) reserved from the fit for testing and
2021 # (b) bad stars that have failed quality checks.)
2022 flagStarSchema = self._makeFlagStarSchema()
2023 flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices()
2024 flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct)
2026 fgcmDatasetDict['fgcmFlaggedStars'] = flagStarCat
2028 # Save the zeropoint information and atmospheres only if desired
2029 if self.outputZeropoints:
2030 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
2031 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
2033 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
2034 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
2036 fgcmDatasetDict['fgcmZeropoints'] = zptCat
2038 # Save atmosphere values
2039 # These are generated by the same code that generates zeropoints
2040 atmSchema = makeAtmSchema()
2041 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
2043 fgcmDatasetDict['fgcmAtmosphereParameters'] = atmCat
2045 # Save the standard stars (if configured)
2046 if self.outputStandards:
2047 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
2048 stdSchema = makeStdSchema(len(goodBands))
2049 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
2051 fgcmDatasetDict['fgcmStandardStars'] = stdCat
2053 return fgcmDatasetDict
2055 def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
2056 lutFilterNameString, fitBandString):
2057 """
2058 Make the parameter persistence schema
2060 Parameters
2061 ----------
2062 parInfo: `numpy.ndarray`
2063 Parameter information returned by fgcm
2064 pars: `numpy.ndarray`
2065 Parameter values returned by fgcm
2066 parSuperStarFlat: `numpy.array`
2067 Superstar flat values returned by fgcm
2068 lutFilterNameString: `str`
2069 Combined string of all the lutFilterNames
2070 fitBandString: `str`
2071 Combined string of all the fitBands
2073 Returns
2074 -------
2075 parSchema: `afwTable.schema`
2076 """
2078 parSchema = afwTable.Schema()
2080 # parameter info section
2081 parSchema.addField('nCcd', type=np.int32, doc='Number of CCDs')
2082 parSchema.addField('lutFilterNames', type=str, doc='LUT Filter names in parameter file',
2083 size=len(lutFilterNameString))
2084 parSchema.addField('fitBands', type=str, doc='Bands that were fit',
2085 size=len(fitBandString))
2086 parSchema.addField('lnTauUnit', type=np.float64, doc='Step units for ln(AOD)')
2087 parSchema.addField('lnTauSlopeUnit', type=np.float64,
2088 doc='Step units for ln(AOD) slope')
2089 parSchema.addField('alphaUnit', type=np.float64, doc='Step units for alpha')
2090 parSchema.addField('lnPwvUnit', type=np.float64, doc='Step units for ln(pwv)')
2091 parSchema.addField('lnPwvSlopeUnit', type=np.float64,
2092 doc='Step units for ln(pwv) slope')
2093 parSchema.addField('lnPwvQuadraticUnit', type=np.float64,
2094 doc='Step units for ln(pwv) quadratic term')
2095 parSchema.addField('lnPwvGlobalUnit', type=np.float64,
2096 doc='Step units for global ln(pwv) parameters')
2097 parSchema.addField('o3Unit', type=np.float64, doc='Step units for O3')
2098 parSchema.addField('qeSysUnit', type=np.float64, doc='Step units for mirror gray')
2099 parSchema.addField('filterOffsetUnit', type=np.float64, doc='Step units for filter offset')
2100 parSchema.addField('hasExternalPwv', type=np.int32, doc='Parameters fit using external pwv')
2101 parSchema.addField('hasExternalTau', type=np.int32, doc='Parameters fit using external tau')
2103 # parameter section
2104 parSchema.addField('parAlpha', type='ArrayD', doc='Alpha parameter vector',
2105 size=pars['PARALPHA'].size)
2106 parSchema.addField('parO3', type='ArrayD', doc='O3 parameter vector',
2107 size=pars['PARO3'].size)
2108 parSchema.addField('parLnTauIntercept', type='ArrayD',
2109 doc='ln(Tau) intercept parameter vector',
2110 size=pars['PARLNTAUINTERCEPT'].size)
2111 parSchema.addField('parLnTauSlope', type='ArrayD',
2112 doc='ln(Tau) slope parameter vector',
2113 size=pars['PARLNTAUSLOPE'].size)
2114 parSchema.addField('parLnPwvIntercept', type='ArrayD', doc='ln(pwv) intercept parameter vector',
2115 size=pars['PARLNPWVINTERCEPT'].size)
2116 parSchema.addField('parLnPwvSlope', type='ArrayD', doc='ln(pwv) slope parameter vector',
2117 size=pars['PARLNPWVSLOPE'].size)
2118 parSchema.addField('parLnPwvQuadratic', type='ArrayD', doc='ln(pwv) quadratic parameter vector',
2119 size=pars['PARLNPWVQUADRATIC'].size)
2120 parSchema.addField('parQeSysIntercept', type='ArrayD', doc='Mirror gray intercept parameter vector',
2121 size=pars['PARQESYSINTERCEPT'].size)
2122 parSchema.addField('compQeSysSlope', type='ArrayD', doc='Mirror gray slope parameter vector',
2123 size=pars[0]['COMPQESYSSLOPE'].size)
2124 parSchema.addField('parFilterOffset', type='ArrayD', doc='Filter offset parameter vector',
2125 size=pars['PARFILTEROFFSET'].size)
2126 parSchema.addField('parFilterOffsetFitFlag', type='ArrayI', doc='Filter offset parameter fit flag',
2127 size=pars['PARFILTEROFFSETFITFLAG'].size)
2128 parSchema.addField('parRetrievedLnPwvScale', type=np.float64,
2129 doc='Global scale for retrieved ln(pwv)')
2130 parSchema.addField('parRetrievedLnPwvOffset', type=np.float64,
2131 doc='Global offset for retrieved ln(pwv)')
2132 parSchema.addField('parRetrievedLnPwvNightlyOffset', type='ArrayD',
2133 doc='Nightly offset for retrieved ln(pwv)',
2134 size=pars['PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size)
2135 parSchema.addField('compAbsThroughput', type='ArrayD',
2136 doc='Absolute throughput (relative to transmission curves)',
2137 size=pars['COMPABSTHROUGHPUT'].size)
2138 parSchema.addField('compRefOffset', type='ArrayD',
2139 doc='Offset between reference stars and calibrated stars',
2140 size=pars['COMPREFOFFSET'].size)
2141 parSchema.addField('compRefSigma', type='ArrayD',
2142 doc='Width of reference star/calibrated star distribution',
2143 size=pars['COMPREFSIGMA'].size)
2144 parSchema.addField('compMirrorChromaticity', type='ArrayD',
2145 doc='Computed mirror chromaticity terms',
2146 size=pars['COMPMIRRORCHROMATICITY'].size)
2147 parSchema.addField('mirrorChromaticityPivot', type='ArrayD',
2148 doc='Mirror chromaticity pivot mjd',
2149 size=pars['MIRRORCHROMATICITYPIVOT'].size)
2150 parSchema.addField('compCcdChromaticity', type='ArrayD',
2151 doc='Computed CCD chromaticity terms',
2152 size=pars['COMPCCDCHROMATICITY'].size)
2153 parSchema.addField('compMedianSedSlope', type='ArrayD',
2154 doc='Computed median SED slope (per band)',
2155 size=pars['COMPMEDIANSEDSLOPE'].size)
2156 parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot',
2157 size=pars['COMPAPERCORRPIVOT'].size)
2158 parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope',
2159 size=pars['COMPAPERCORRSLOPE'].size)
2160 parSchema.addField('compAperCorrSlopeErr', type='ArrayD', doc='Aperture correction slope error',
2161 size=pars['COMPAPERCORRSLOPEERR'].size)
2162 parSchema.addField('compAperCorrRange', type='ArrayD', doc='Aperture correction range',
2163 size=pars['COMPAPERCORRRANGE'].size)
2164 parSchema.addField('compModelErrExptimePivot', type='ArrayD', doc='Model error exptime pivot',
2165 size=pars['COMPMODELERREXPTIMEPIVOT'].size)
2166 parSchema.addField('compModelErrFwhmPivot', type='ArrayD', doc='Model error fwhm pivot',
2167 size=pars['COMPMODELERRFWHMPIVOT'].size)
2168 parSchema.addField('compModelErrSkyPivot', type='ArrayD', doc='Model error sky pivot',
2169 size=pars['COMPMODELERRSKYPIVOT'].size)
2170 parSchema.addField('compModelErrPars', type='ArrayD', doc='Model error parameters',
2171 size=pars['COMPMODELERRPARS'].size)
2172 parSchema.addField('compExpGray', type='ArrayD', doc='Computed exposure gray',
2173 size=pars['COMPEXPGRAY'].size)
2174 parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance',
2175 size=pars['COMPVARGRAY'].size)
2176 parSchema.addField('compExpDeltaMagBkg', type='ArrayD',
2177 doc='Computed exposure offset due to background',
2178 size=pars['COMPEXPDELTAMAGBKG'].size)
2179 parSchema.addField('compNGoodStarPerExp', type='ArrayI',
2180 doc='Computed number of good stars per exposure',
2181 size=pars['COMPNGOODSTARPEREXP'].size)
2182 parSchema.addField('compExpRefOffset', type='ArrayD',
2183 doc='Computed per-visit median offset between standard stars and ref stars.',
2184 size=pars['COMPEXPREFOFFSET'].size)
2185 parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)',
2186 size=pars['COMPSIGFGCM'].size)
2187 parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)',
2188 size=pars['COMPSIGMACAL'].size)
2189 parSchema.addField('compRetrievedLnPwv', type='ArrayD', doc='Retrieved ln(pwv) (smoothed)',
2190 size=pars['COMPRETRIEVEDLNPWV'].size)
2191 parSchema.addField('compRetrievedLnPwvRaw', type='ArrayD', doc='Retrieved ln(pwv) (raw)',
2192 size=pars['COMPRETRIEVEDLNPWVRAW'].size)
2193 parSchema.addField('compRetrievedLnPwvFlag', type='ArrayI', doc='Retrieved ln(pwv) Flag',
2194 size=pars['COMPRETRIEVEDLNPWVFLAG'].size)
2195 parSchema.addField('compRetrievedTauNight', type='ArrayD', doc='Retrieved tau (per night)',
2196 size=pars['COMPRETRIEVEDTAUNIGHT'].size)
2197 parSchema.addField('compEpsilon', type='ArrayD',
2198 doc='Computed epsilon background offset per visit (nJy/arcsec2)',
2199 size=pars['COMPEPSILON'].size)
2200 parSchema.addField('compMedDeltaAper', type='ArrayD',
2201 doc='Median delta mag aper per visit',
2202 size=pars['COMPMEDDELTAAPER'].size)
2203 parSchema.addField('compGlobalEpsilon', type='ArrayD',
2204 doc='Computed epsilon bkg offset (global) (nJy/arcsec2)',
2205 size=pars['COMPGLOBALEPSILON'].size)
2206 parSchema.addField('compEpsilonMap', type='ArrayD',
2207 doc='Computed epsilon maps (nJy/arcsec2)',
2208 size=pars['COMPEPSILONMAP'].size)
2209 parSchema.addField('compEpsilonNStarMap', type='ArrayI',
2210 doc='Number of stars per pixel in computed epsilon maps',
2211 size=pars['COMPEPSILONNSTARMAP'].size)
2212 parSchema.addField('compEpsilonCcdMap', type='ArrayD',
2213 doc='Computed epsilon ccd maps (nJy/arcsec2)',
2214 size=pars['COMPEPSILONCCDMAP'].size)
2215 parSchema.addField('compEpsilonCcdNStarMap', type='ArrayI',
2216 doc='Number of stars per ccd bin in epsilon ccd maps',
2217 size=pars['COMPEPSILONCCDNSTARMAP'].size)
2218 parSchema.addField('epochMjdStart', type='ArrayD',
2219 doc='Epoch MJD start times',
2220 size=pars['EPOCHMJDSTART'].size)
2221 parSchema.addField('epochMjdEnd', type='ArrayD',
2222 doc='EpochMJD end times',
2223 size=pars['EPOCHMJDEND'].size)
2224 # superstarflat section
2225 parSchema.addField('superstarSize', type='ArrayI', doc='Superstar matrix size',
2226 size=4)
2227 parSchema.addField('superstar', type='ArrayD', doc='Superstar matrix (flattened)',
2228 size=parSuperStarFlat.size)
2230 return parSchema
2232 def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
2233 lutFilterNameString, fitBandString):
2234 """
2235 Make the FGCM parameter catalog for persistence
2237 Parameters
2238 ----------
2239 parSchema: `lsst.afw.table.Schema`
2240 Parameter catalog schema
2241 pars: `numpy.ndarray`
2242 FGCM parameters to put into parCat
2243 parSuperStarFlat: `numpy.array`
2244 FGCM superstar flat array to put into parCat
2245 lutFilterNameString: `str`
2246 Combined string of all the lutFilterNames
2247 fitBandString: `str`
2248 Combined string of all the fitBands
2250 Returns
2251 -------
2252 parCat: `afwTable.BasicCatalog`
2253 Atmosphere and instrumental model parameter catalog for persistence
2254 """
2256 parCat = afwTable.BaseCatalog(parSchema)
2257 parCat.reserve(1)
2259 # The parameter catalog just has one row, with many columns for all the
2260 # atmosphere and instrument fit parameters
2261 rec = parCat.addNew()
2263 # info section
2264 rec['nCcd'] = parInfo['NCCD'][0]
2265 rec['lutFilterNames'] = lutFilterNameString
2266 rec['fitBands'] = fitBandString
2267 # note these are not currently supported here.
2268 rec['hasExternalPwv'] = 0
2269 rec['hasExternalTau'] = 0
2271 # parameter section
2273 scalarNames = ['parRetrievedLnPwvScale', 'parRetrievedLnPwvOffset']
2275 arrNames = ['parAlpha', 'parO3', 'parLnTauIntercept', 'parLnTauSlope',
2276 'parLnPwvIntercept', 'parLnPwvSlope', 'parLnPwvQuadratic',
2277 'parQeSysIntercept', 'compQeSysSlope',
2278 'parRetrievedLnPwvNightlyOffset', 'compAperCorrPivot',
2279 'parFilterOffset', 'parFilterOffsetFitFlag',
2280 'compAbsThroughput', 'compRefOffset', 'compRefSigma',
2281 'compMirrorChromaticity', 'mirrorChromaticityPivot', 'compCcdChromaticity',
2282 'compAperCorrSlope', 'compAperCorrSlopeErr', 'compAperCorrRange',
2283 'compModelErrExptimePivot', 'compModelErrFwhmPivot',
2284 'compModelErrSkyPivot', 'compModelErrPars',
2285 'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm',
2286 'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope',
2287 'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
2288 'compRetrievedTauNight', 'compEpsilon', 'compMedDeltaAper',
2289 'compGlobalEpsilon', 'compEpsilonMap', 'compEpsilonNStarMap',
2290 'compEpsilonCcdMap', 'compEpsilonCcdNStarMap', 'compExpRefOffset',
2291 'epochMjdStart', 'epochMjdEnd']
2293 for scalarName in scalarNames:
2294 rec[scalarName] = pars[scalarName.upper()][0]
2296 for arrName in arrNames:
2297 rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:]
2299 # superstar section
2300 rec['superstarSize'][:] = parSuperStarFlat.shape
2301 rec['superstar'][:] = parSuperStarFlat.ravel()
2303 return parCat
2305 def _makeFlagStarSchema(self):
2306 """
2307 Make the flagged-stars schema
2309 Returns
2310 -------
2311 flagStarSchema: `lsst.afw.table.Schema`
2312 """
2314 flagStarSchema = afwTable.Schema()
2316 flagStarSchema.addField('objId', type=np.int32, doc='FGCM object id')
2317 flagStarSchema.addField('objFlag', type=np.int32, doc='FGCM object flag')
2319 return flagStarSchema
2321 def _makeFlagStarCat(self, flagStarSchema, flagStarStruct):
2322 """
2323 Make the flagged star catalog for persistence
2325 Parameters
2326 ----------
2327 flagStarSchema: `lsst.afw.table.Schema`
2328 Flagged star schema
2329 flagStarStruct: `numpy.ndarray`
2330 Flagged star structure from fgcm
2332 Returns
2333 -------
2334 flagStarCat: `lsst.afw.table.BaseCatalog`
2335 Flagged star catalog for persistence
2336 """
2338 flagStarCat = afwTable.BaseCatalog(flagStarSchema)
2339 flagStarCat.resize(flagStarStruct.size)
2341 flagStarCat['objId'][:] = flagStarStruct['OBJID']
2342 flagStarCat['objFlag'][:] = flagStarStruct['OBJFLAG']
2344 return flagStarCat