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