LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
fgcmFitCycle.py
Go to the documentation of this file.
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.
24 
25 This task runs a single "fit cycle" of fgcm. Prior to running this task
26 one must run both fgcmMakeLut (to construct the atmosphere and instrumental
27 look-up-table) and fgcmBuildStars (to extract visits and star observations
28 for the global fit).
29 
30 The fgcmFitCycle is meant to be run multiple times, and is tracked by the
31 'cycleNumber'. After each run of the fit cycle, diagnostic plots should
32 be inspected to set parameters for outlier rejection on the following
33 cycle. Please see the fgcmcal Cookbook for details.
34 """
35 
36 import sys
37 import traceback
38 import copy
39 
40 import numpy as np
41 
42 import lsst.pex.config as pexConfig
43 import lsst.pipe.base as pipeBase
44 from lsst.pipe.base import connectionTypes
45 import lsst.afw.table as afwTable
46 from lsst.utils.timer import timeMethod
47 
48 from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
49 from .utilities import extractReferenceMags
50 from .utilities import makeZptSchema, makeZptCat
51 from .utilities import makeAtmSchema, makeAtmCat, makeStdSchema, makeStdCat
52 from .sedterms import SedboundarytermDict, SedtermDict
53 from .utilities import lookupStaticCalibrations
54 from .focalPlaneProjector import FocalPlaneProjector
55 
56 import fgcm
57 
58 __all__ = ['FgcmFitCycleConfig', 'FgcmFitCycleTask', 'FgcmFitCycleRunner']
59 
60 MULTIPLE_CYCLES_MAX = 10
61 
62 
63 class FgcmFitCycleConnections(pipeBase.PipelineTaskConnections,
64  dimensions=("instrument",),
65  defaultTemplates={"previousCycleNumber": "-1",
66  "cycleNumber": "0"}):
67  camera = connectionTypes.PrerequisiteInput(
68  doc="Camera instrument",
69  name="camera",
70  storageClass="Camera",
71  dimensions=("instrument",),
72  lookupFunction=lookupStaticCalibrations,
73  isCalibration=True,
74  )
75 
76  fgcmLookUpTable = connectionTypes.PrerequisiteInput(
77  doc=("Atmosphere + instrument look-up-table for FGCM throughput and "
78  "chromatic corrections."),
79  name="fgcmLookUpTable",
80  storageClass="Catalog",
81  dimensions=("instrument",),
82  deferLoad=True,
83  )
84 
85  fgcmVisitCatalog = connectionTypes.Input(
86  doc="Catalog of visit information for fgcm",
87  name="fgcmVisitCatalog",
88  storageClass="Catalog",
89  dimensions=("instrument",),
90  deferLoad=True,
91  )
92 
93  fgcmStarObservations = connectionTypes.Input(
94  doc="Catalog of star observations for fgcm",
95  name="fgcmStarObservations",
96  storageClass="Catalog",
97  dimensions=("instrument",),
98  deferLoad=True,
99  )
100 
101  fgcmStarIds = connectionTypes.Input(
102  doc="Catalog of fgcm calibration star IDs",
103  name="fgcmStarIds",
104  storageClass="Catalog",
105  dimensions=("instrument",),
106  deferLoad=True,
107  )
108 
109  fgcmStarIndices = connectionTypes.Input(
110  doc="Catalog of fgcm calibration star indices",
111  name="fgcmStarIndices",
112  storageClass="Catalog",
113  dimensions=("instrument",),
114  deferLoad=True,
115  )
116 
117  fgcmReferenceStars = connectionTypes.Input(
118  doc="Catalog of fgcm-matched reference stars",
119  name="fgcmReferenceStars",
120  storageClass="Catalog",
121  dimensions=("instrument",),
122  deferLoad=True,
123  )
124 
125  fgcmFlaggedStarsInput = connectionTypes.PrerequisiteInput(
126  doc="Catalog of flagged stars for fgcm calibration from previous fit cycle",
127  name="fgcmFlaggedStars{previousCycleNumber}",
128  storageClass="Catalog",
129  dimensions=("instrument",),
130  deferLoad=True,
131  )
132 
133  fgcmFitParametersInput = connectionTypes.PrerequisiteInput(
134  doc="Catalog of fgcm fit parameters from previous fit cycle",
135  name="fgcmFitParameters{previousCycleNumber}",
136  storageClass="Catalog",
137  dimensions=("instrument",),
138  deferLoad=True,
139  )
140 
141  fgcmFitParameters = connectionTypes.Output(
142  doc="Catalog of fgcm fit parameters from current fit cycle",
143  name="fgcmFitParameters{cycleNumber}",
144  storageClass="Catalog",
145  dimensions=("instrument",),
146  )
147 
148  fgcmFlaggedStars = connectionTypes.Output(
149  doc="Catalog of flagged stars for fgcm calibration from current fit cycle",
150  name="fgcmFlaggedStars{cycleNumber}",
151  storageClass="Catalog",
152  dimensions=("instrument",),
153  )
154 
155  fgcmZeropoints = connectionTypes.Output(
156  doc="Catalog of fgcm zeropoint data from current fit cycle",
157  name="fgcmZeropoints{cycleNumber}",
158  storageClass="Catalog",
159  dimensions=("instrument",),
160  )
161 
162  fgcmAtmosphereParameters = connectionTypes.Output(
163  doc="Catalog of atmospheric fit parameters from current fit cycle",
164  name="fgcmAtmosphereParameters{cycleNumber}",
165  storageClass="Catalog",
166  dimensions=("instrument",),
167  )
168 
169  fgcmStandardStars = connectionTypes.Output(
170  doc="Catalog of standard star magnitudes from current fit cycle",
171  name="fgcmStandardStars{cycleNumber}",
172  storageClass="SimpleCatalog",
173  dimensions=("instrument",),
174  )
175 
176  # Add connections for running multiple cycles
177  # This uses vars() to alter the class __dict__ to programmatically
178  # write many similar outputs.
179  for cycle in range(MULTIPLE_CYCLES_MAX):
180  vars()[f"fgcmFitParameters{cycle}"] = connectionTypes.Output(
181  doc=f"Catalog of fgcm fit parameters from fit cycle {cycle}",
182  name=f"fgcmFitParameters{cycle}",
183  storageClass="Catalog",
184  dimensions=("instrument",),
185  )
186  vars()[f"fgcmFlaggedStars{cycle}"] = connectionTypes.Output(
187  doc=f"Catalog of flagged stars for fgcm calibration from fit cycle {cycle}",
188  name=f"fgcmFlaggedStars{cycle}",
189  storageClass="Catalog",
190  dimensions=("instrument",),
191  )
192  vars()[f"fgcmZeropoints{cycle}"] = connectionTypes.Output(
193  doc=f"Catalog of fgcm zeropoint data from fit cycle {cycle}",
194  name=f"fgcmZeropoints{cycle}",
195  storageClass="Catalog",
196  dimensions=("instrument",),
197  )
198  vars()[f"fgcmAtmosphereParameters{cycle}"] = connectionTypes.Output(
199  doc=f"Catalog of atmospheric fit parameters from fit cycle {cycle}",
200  name=f"fgcmAtmosphereParameters{cycle}",
201  storageClass="Catalog",
202  dimensions=("instrument",),
203  )
204  vars()[f"fgcmStandardStars{cycle}"] = connectionTypes.Output(
205  doc=f"Catalog of standard star magnitudes from fit cycle {cycle}",
206  name=f"fgcmStandardStars{cycle}",
207  storageClass="SimpleCatalog",
208  dimensions=("instrument",),
209  )
210 
211  def __init__(self, *, config=None):
212  super().__init__(config=config)
213 
214  if not config.doReferenceCalibration:
215  self.inputs.remove("fgcmReferenceStars")
216 
217  if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber:
218  raise ValueError("cycleNumber must be of integer format")
219  if str(int(config.connections.previousCycleNumber)) != config.connections.previousCycleNumber:
220  raise ValueError("previousCycleNumber must be of integer format")
221  if int(config.connections.previousCycleNumber) != (int(config.connections.cycleNumber) - 1):
222  raise ValueError("previousCycleNumber must be 1 less than cycleNumber")
223 
224  if int(config.connections.cycleNumber) == 0:
225  self.prerequisiteInputs.remove("fgcmFlaggedStarsInput")
226  self.prerequisiteInputs.remove("fgcmFitParametersInput")
227 
228  if not self.config.doMultipleCycles:
229  # Single-cycle run
230  if not self.config.isFinalCycle and not self.config.outputStandardsBeforeFinalCycle:
231  self.outputs.remove("fgcmStandardStars")
232 
233  if not self.config.isFinalCycle and not self.config.outputZeropointsBeforeFinalCycle:
234  self.outputs.remove("fgcmZeropoints")
235  self.outputs.remove("fgcmAtmosphereParameters")
236 
237  # Remove all multiple cycle outputs
238  for cycle in range(0, MULTIPLE_CYCLES_MAX):
239  self.outputs.remove(f"fgcmFitParameters{cycle}")
240  self.outputs.remove(f"fgcmFlaggedStars{cycle}")
241  self.outputs.remove(f"fgcmZeropoints{cycle}")
242  self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
243  self.outputs.remove(f"fgcmStandardStars{cycle}")
244 
245  else:
246  # Multiple-cycle run
247  # Remove single-cycle outputs
248  self.outputs.remove("fgcmFitParameters")
249  self.outputs.remove("fgcmFlaggedStars")
250  self.outputs.remove("fgcmZeropoints")
251  self.outputs.remove("fgcmAtmosphereParameters")
252  self.outputs.remove("fgcmStandardStars")
253 
254  # Remove outputs from cycles that are not used
255  for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1,
256  MULTIPLE_CYCLES_MAX):
257  self.outputs.remove(f"fgcmFitParameters{cycle}")
258  self.outputs.remove(f"fgcmFlaggedStars{cycle}")
259  self.outputs.remove(f"fgcmZeropoints{cycle}")
260  self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
261  self.outputs.remove(f"fgcmStandardStars{cycle}")
262 
263  # Remove non-final-cycle outputs if necessary
264  for cycle in range(self.config.multipleCyclesFinalCycleNumber):
265  if not self.config.outputZeropointsBeforeFinalCycle:
266  self.outputs.remove(f"fgcmZeropoints{cycle}")
267  self.outputs.remove(f"fgcmAtmosphereParameters{cycle}")
268  if not self.config.outputStandardsBeforeFinalCycle:
269  self.outputs.remove(f"fgcmStandardStars{cycle}")
270 
271 
272 class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
273  pipelineConnections=FgcmFitCycleConnections):
274  """Config for FgcmFitCycle"""
275 
276  doMultipleCycles = pexConfig.Field(
277  doc="Run multiple fit cycles in one task",
278  dtype=bool,
279  default=False,
280  )
281  multipleCyclesFinalCycleNumber = pexConfig.RangeField(
282  doc=("Final cycle number in multiple cycle mode. The initial cycle "
283  "is 0, with limited parameters fit. The next cycle is 1 with "
284  "full parameter fit. The final cycle is a clean-up with no "
285  "parameters fit. There will be a total of "
286  "(multipleCycleFinalCycleNumber + 1) cycles run, and the final "
287  "cycle number cannot be less than 2."),
288  dtype=int,
289  default=5,
290  min=2,
291  max=MULTIPLE_CYCLES_MAX,
292  inclusiveMax=True,
293  )
294  bands = pexConfig.ListField(
295  doc="Bands to run calibration",
296  dtype=str,
297  default=[],
298  )
299  fitBands = pexConfig.ListField(
300  doc=("Bands to use in atmospheric fit. The bands not listed here will have "
301  "the atmosphere constrained from the 'fitBands' on the same night. "
302  "Must be a subset of `config.bands`"),
303  dtype=str,
304  default=[],
305  )
306  requiredBands = pexConfig.ListField(
307  doc=("Bands that are required for a star to be considered a calibration star. "
308  "Must be a subset of `config.bands`"),
309  dtype=str,
310  default=[],
311  )
312  # The following config will not be necessary after Gen2 retirement.
313  # In the meantime, it is set to 'filterDefinitions.filter_to_band' which
314  # is easiest to access in the config file.
315  physicalFilterMap = pexConfig.DictField(
316  doc="Mapping from 'physicalFilter' to band.",
317  keytype=str,
318  itemtype=str,
319  default={},
320  )
321  doReferenceCalibration = pexConfig.Field(
322  doc="Use reference catalog as additional constraint on calibration",
323  dtype=bool,
324  default=True,
325  )
326  refStarSnMin = pexConfig.Field(
327  doc="Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.",
328  dtype=float,
329  default=50.0,
330  )
331  refStarOutlierNSig = pexConfig.Field(
332  doc=("Number of sigma compared to average mag for reference star to be considered an outlier. "
333  "Computed per-band, and if it is an outlier in any band it is rejected from fits."),
334  dtype=float,
335  default=4.0,
336  )
337  applyRefStarColorCuts = pexConfig.Field(
338  doc="Apply color cuts to reference stars?",
339  dtype=bool,
340  default=True,
341  )
342  nCore = pexConfig.Field(
343  doc="Number of cores to use",
344  dtype=int,
345  default=4,
346  )
347  nStarPerRun = pexConfig.Field(
348  doc="Number of stars to run in each chunk",
349  dtype=int,
350  default=200000,
351  )
352  nExpPerRun = pexConfig.Field(
353  doc="Number of exposures to run in each chunk",
354  dtype=int,
355  default=1000,
356  )
357  reserveFraction = pexConfig.Field(
358  doc="Fraction of stars to reserve for testing",
359  dtype=float,
360  default=0.1,
361  )
362  freezeStdAtmosphere = pexConfig.Field(
363  doc="Freeze atmosphere parameters to standard (for testing)",
364  dtype=bool,
365  default=False,
366  )
367  precomputeSuperStarInitialCycle = pexConfig.Field(
368  doc="Precompute superstar flat for initial cycle",
369  dtype=bool,
370  default=False,
371  )
372  superStarSubCcdDict = pexConfig.DictField(
373  doc=("Per-band specification on whether to compute superstar flat on sub-ccd scale. "
374  "Must have one entry per band."),
375  keytype=str,
376  itemtype=bool,
377  default={},
378  )
379  superStarSubCcdChebyshevOrder = pexConfig.Field(
380  doc=("Order of the 2D chebyshev polynomials for sub-ccd superstar fit. "
381  "Global default is first-order polynomials, and should be overridden "
382  "on a camera-by-camera basis depending on the ISR."),
383  dtype=int,
384  default=1,
385  )
386  superStarSubCcdTriangular = pexConfig.Field(
387  doc=("Should the sub-ccd superstar chebyshev matrix be triangular to "
388  "suppress high-order cross terms?"),
389  dtype=bool,
390  default=False,
391  )
392  superStarSigmaClip = pexConfig.Field(
393  doc="Number of sigma to clip outliers when selecting for superstar flats",
394  dtype=float,
395  default=5.0,
396  )
397  focalPlaneSigmaClip = pexConfig.Field(
398  doc="Number of sigma to clip outliers per focal-plane.",
399  dtype=float,
400  default=4.0,
401  )
402  ccdGraySubCcdDict = pexConfig.DictField(
403  doc=("Per-band specification on whether to compute achromatic per-ccd residual "
404  "('ccd gray') on a sub-ccd scale."),
405  keytype=str,
406  itemtype=bool,
407  default={},
408  )
409  ccdGraySubCcdChebyshevOrder = pexConfig.Field(
410  doc="Order of the 2D chebyshev polynomials for sub-ccd gray fit.",
411  dtype=int,
412  default=1,
413  )
414  ccdGraySubCcdTriangular = pexConfig.Field(
415  doc=("Should the sub-ccd gray chebyshev matrix be triangular to "
416  "suppress high-order cross terms?"),
417  dtype=bool,
418  default=True,
419  )
420  ccdGrayFocalPlaneDict = pexConfig.DictField(
421  doc=("Per-band specification on whether to compute focal-plane residual "
422  "('ccd gray') corrections."),
423  keytype=str,
424  itemtype=bool,
425  default={},
426  )
427  ccdGrayFocalPlaneFitMinCcd = pexConfig.Field(
428  doc=("Minimum number of 'good' CCDs required to perform focal-plane "
429  "gray corrections. If there are fewer good CCDs then the gray "
430  "correction is computed per-ccd."),
431  dtype=int,
432  default=1,
433  )
434  ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field(
435  doc="Order of the 2D chebyshev polynomials for focal plane fit.",
436  dtype=int,
437  default=3,
438  )
439  cycleNumber = pexConfig.Field(
440  doc=("FGCM fit cycle number. This is automatically incremented after each run "
441  "and stage of outlier rejection. See cookbook for details."),
442  dtype=int,
443  default=None,
444  )
445  isFinalCycle = pexConfig.Field(
446  doc=("Is this the final cycle of the fitting? Will automatically compute final "
447  "selection of stars and photometric exposures, and will output zeropoints "
448  "and standard stars for use in fgcmOutputProducts"),
449  dtype=bool,
450  default=False,
451  )
452  maxIterBeforeFinalCycle = pexConfig.Field(
453  doc=("Maximum fit iterations, prior to final cycle. The number of iterations "
454  "will always be 0 in the final cycle for cleanup and final selection."),
455  dtype=int,
456  default=50,
457  )
458  deltaMagBkgOffsetPercentile = pexConfig.Field(
459  doc=("Percentile brightest stars on a visit/ccd to use to compute net "
460  "offset from local background subtraction."),
461  dtype=float,
462  default=0.25,
463  )
464  deltaMagBkgPerCcd = pexConfig.Field(
465  doc=("Compute net offset from local background subtraction per-ccd? "
466  "Otherwise, use computation per visit."),
467  dtype=bool,
468  default=False,
469  )
470  utBoundary = pexConfig.Field(
471  doc="Boundary (in UTC) from day-to-day",
472  dtype=float,
473  default=None,
474  )
475  washMjds = pexConfig.ListField(
476  doc="Mirror wash MJDs",
477  dtype=float,
478  default=(0.0,),
479  )
480  epochMjds = pexConfig.ListField(
481  doc="Epoch boundaries in MJD",
482  dtype=float,
483  default=(0.0,),
484  )
485  minObsPerBand = pexConfig.Field(
486  doc="Minimum good observations per band",
487  dtype=int,
488  default=2,
489  )
490  # TODO: When DM-16511 is done, it will be possible to get the
491  # telescope latitude directly from the camera.
492  latitude = pexConfig.Field(
493  doc="Observatory latitude",
494  dtype=float,
495  default=None,
496  )
497  defaultCameraOrientation = pexConfig.Field(
498  doc="Default camera orientation for QA plots.",
499  dtype=float,
500  default=None,
501  )
502  brightObsGrayMax = pexConfig.Field(
503  doc="Maximum gray extinction to be considered bright observation",
504  dtype=float,
505  default=0.15,
506  )
507  minStarPerCcd = pexConfig.Field(
508  doc=("Minimum number of good stars per CCD to be used in calibration fit. "
509  "CCDs with fewer stars will have their calibration estimated from other "
510  "CCDs in the same visit, with zeropoint error increased accordingly."),
511  dtype=int,
512  default=5,
513  )
514  minCcdPerExp = pexConfig.Field(
515  doc=("Minimum number of good CCDs per exposure/visit to be used in calibration fit. "
516  "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."),
517  dtype=int,
518  default=5,
519  )
520  maxCcdGrayErr = pexConfig.Field(
521  doc="Maximum error on CCD gray offset to be considered photometric",
522  dtype=float,
523  default=0.05,
524  )
525  minStarPerExp = pexConfig.Field(
526  doc=("Minimum number of good stars per exposure/visit to be used in calibration fit. "
527  "Visits with fewer good stars will have CCD zeropoints estimated where possible."),
528  dtype=int,
529  default=600,
530  )
531  minExpPerNight = pexConfig.Field(
532  doc="Minimum number of good exposures/visits to consider a partly photometric night",
533  dtype=int,
534  default=10,
535  )
536  expGrayInitialCut = pexConfig.Field(
537  doc=("Maximum exposure/visit gray value for initial selection of possible photometric "
538  "observations."),
539  dtype=float,
540  default=-0.25,
541  )
542  expGrayPhotometricCutDict = pexConfig.DictField(
543  doc=("Per-band specification on maximum (negative) achromatic exposure residual "
544  "('gray term') for a visit to be considered photometric. Must have one "
545  "entry per band. Broad-band filters should be -0.05."),
546  keytype=str,
547  itemtype=float,
548  default={},
549  )
550  expGrayHighCutDict = pexConfig.DictField(
551  doc=("Per-band specification on maximum (positive) achromatic exposure residual "
552  "('gray term') for a visit to be considered photometric. Must have one "
553  "entry per band. Broad-band filters should be 0.2."),
554  keytype=str,
555  itemtype=float,
556  default={},
557  )
558  expGrayRecoverCut = pexConfig.Field(
559  doc=("Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. "
560  "Visits with more gray extinction will only get CCD zeropoints if there are "
561  "sufficient star observations (minStarPerCcd) on that CCD."),
562  dtype=float,
563  default=-1.0,
564  )
565  expVarGrayPhotometricCutDict = pexConfig.DictField(
566  doc=("Per-band specification on maximum exposure variance to be considered possibly "
567  "photometric. Must have one entry per band. Broad-band filters should be "
568  "0.0005."),
569  keytype=str,
570  itemtype=float,
571  default={},
572  )
573  expGrayErrRecoverCut = pexConfig.Field(
574  doc=("Maximum exposure gray error to be able to recover bad ccds via interpolation. "
575  "Visits with more gray variance will only get CCD zeropoints if there are "
576  "sufficient star observations (minStarPerCcd) on that CCD."),
577  dtype=float,
578  default=0.05,
579  )
580  aperCorrFitNBins = pexConfig.Field(
581  doc=("Number of aperture bins used in aperture correction fit. When set to 0"
582  "no fit will be performed, and the config.aperCorrInputSlopes will be "
583  "used if available."),
584  dtype=int,
585  default=10,
586  )
587  aperCorrInputSlopeDict = pexConfig.DictField(
588  doc=("Per-band specification of aperture correction input slope parameters. These "
589  "are used on the first fit iteration, and aperture correction parameters will "
590  "be updated from the data if config.aperCorrFitNBins > 0. It is recommended "
591  "to set this when there is insufficient data to fit the parameters (e.g. "
592  "tract mode)."),
593  keytype=str,
594  itemtype=float,
595  default={},
596  )
597  sedboundaryterms = pexConfig.ConfigField(
598  doc="Mapping from bands to SED boundary term names used is sedterms.",
599  dtype=SedboundarytermDict,
600  )
601  sedterms = pexConfig.ConfigField(
602  doc="Mapping from terms to bands for fgcm linear SED approximations.",
603  dtype=SedtermDict,
604  )
605  sigFgcmMaxErr = pexConfig.Field(
606  doc="Maximum mag error for fitting sigma_FGCM",
607  dtype=float,
608  default=0.01,
609  )
610  sigFgcmMaxEGrayDict = pexConfig.DictField(
611  doc=("Per-band specification for maximum (absolute) achromatic residual (gray value) "
612  "for observations in sigma_fgcm (raw repeatability). Broad-band filters "
613  "should be 0.05."),
614  keytype=str,
615  itemtype=float,
616  default={},
617  )
618  ccdGrayMaxStarErr = pexConfig.Field(
619  doc=("Maximum error on a star observation to use in ccd gray (achromatic residual) "
620  "computation"),
621  dtype=float,
622  default=0.10,
623  )
624  approxThroughputDict = pexConfig.DictField(
625  doc=("Per-band specification of the approximate overall throughput at the start of "
626  "calibration observations. Must have one entry per band. Typically should "
627  "be 1.0."),
628  keytype=str,
629  itemtype=float,
630  default={},
631  )
632  sigmaCalRange = pexConfig.ListField(
633  doc="Allowed range for systematic error floor estimation",
634  dtype=float,
635  default=(0.001, 0.003),
636  )
637  sigmaCalFitPercentile = pexConfig.ListField(
638  doc="Magnitude percentile range to fit systematic error floor",
639  dtype=float,
640  default=(0.05, 0.15),
641  )
642  sigmaCalPlotPercentile = pexConfig.ListField(
643  doc="Magnitude percentile range to plot systematic error floor",
644  dtype=float,
645  default=(0.05, 0.95),
646  )
647  sigma0Phot = pexConfig.Field(
648  doc="Systematic error floor for all zeropoints",
649  dtype=float,
650  default=0.003,
651  )
652  mapLongitudeRef = pexConfig.Field(
653  doc="Reference longitude for plotting maps",
654  dtype=float,
655  default=0.0,
656  )
657  mapNSide = pexConfig.Field(
658  doc="Healpix nside for plotting maps",
659  dtype=int,
660  default=256,
661  )
662  outfileBase = pexConfig.Field(
663  doc="Filename start for plot output files",
664  dtype=str,
665  default=None,
666  )
667  starColorCuts = pexConfig.ListField(
668  doc="Encoded star-color cuts (to be cleaned up)",
669  dtype=str,
670  default=("NO_DATA",),
671  )
672  colorSplitBands = pexConfig.ListField(
673  doc="Band names to use to split stars by color. Must have 2 entries.",
674  dtype=str,
675  length=2,
676  default=('g', 'i'),
677  )
678  modelMagErrors = pexConfig.Field(
679  doc="Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)",
680  dtype=bool,
681  default=True,
682  )
683  useQuadraticPwv = pexConfig.Field(
684  doc="Model PWV with a quadratic term for variation through the night?",
685  dtype=bool,
686  default=False,
687  )
688  instrumentParsPerBand = pexConfig.Field(
689  doc=("Model instrumental parameters per band? "
690  "Otherwise, instrumental parameters (QE changes with time) are "
691  "shared among all bands."),
692  dtype=bool,
693  default=False,
694  )
695  instrumentSlopeMinDeltaT = pexConfig.Field(
696  doc=("Minimum time change (in days) between observations to use in constraining "
697  "instrument slope."),
698  dtype=float,
699  default=20.0,
700  )
701  fitMirrorChromaticity = pexConfig.Field(
702  doc="Fit (intraband) mirror chromatic term?",
703  dtype=bool,
704  default=False,
705  )
706  coatingMjds = pexConfig.ListField(
707  doc="Mirror coating dates in MJD",
708  dtype=float,
709  default=(0.0,),
710  )
711  outputStandardsBeforeFinalCycle = pexConfig.Field(
712  doc="Output standard stars prior to final cycle? Used in debugging.",
713  dtype=bool,
714  default=False,
715  )
716  outputZeropointsBeforeFinalCycle = pexConfig.Field(
717  doc="Output standard stars prior to final cycle? Used in debugging.",
718  dtype=bool,
719  default=False,
720  )
721  useRepeatabilityForExpGrayCutsDict = pexConfig.DictField(
722  doc=("Per-band specification on whether to use star repeatability (instead of exposures) "
723  "for computing photometric cuts. Recommended for tract mode or bands with few visits."),
724  keytype=str,
725  itemtype=bool,
726  default={},
727  )
728  autoPhotometricCutNSig = pexConfig.Field(
729  doc=("Number of sigma for automatic computation of (low) photometric cut. "
730  "Cut is based on exposure gray width (per band), unless "
731  "useRepeatabilityForExpGrayCuts is set, in which case the star "
732  "repeatability is used (also per band)."),
733  dtype=float,
734  default=3.0,
735  )
736  autoHighCutNSig = pexConfig.Field(
737  doc=("Number of sigma for automatic computation of (high) outlier cut. "
738  "Cut is based on exposure gray width (per band), unless "
739  "useRepeatabilityForExpGrayCuts is set, in which case the star "
740  "repeatability is used (also per band)."),
741  dtype=float,
742  default=4.0,
743  )
744  quietMode = pexConfig.Field(
745  doc="Be less verbose with logging.",
746  dtype=bool,
747  default=False,
748  )
749  doPlots = pexConfig.Field(
750  doc="Make fgcm QA plots.",
751  dtype=bool,
752  default=True,
753  )
754  randomSeed = pexConfig.Field(
755  doc="Random seed for fgcm for consistency in tests.",
756  dtype=int,
757  default=None,
758  optional=True,
759  )
760 
761  def validate(self):
762  super().validate()
763 
764  if self.connections.previousCycleNumber != str(self.cycleNumber - 1):
765  msg = "cycleNumber in template must be connections.previousCycleNumber + 1"
766  raise RuntimeError(msg)
767  if self.connections.cycleNumber != str(self.cycleNumber):
768  msg = "cycleNumber in template must be equal to connections.cycleNumber"
769  raise RuntimeError(msg)
770 
771  for band in self.fitBands:
772  if band not in self.bands:
773  msg = 'fitBand %s not in bands' % (band)
774  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg)
775  for band in self.requiredBands:
776  if band not in self.bands:
777  msg = 'requiredBand %s not in bands' % (band)
778  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg)
779  for band in self.colorSplitBands:
780  if band not in self.bands:
781  msg = 'colorSplitBand %s not in bands' % (band)
782  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg)
783  for band in self.bands:
784  if band not in self.superStarSubCcdDict:
785  msg = 'band %s not in superStarSubCcdDict' % (band)
786  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict,
787  self, msg)
788  if band not in self.ccdGraySubCcdDict:
789  msg = 'band %s not in ccdGraySubCcdDict' % (band)
790  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict,
791  self, msg)
792  if band not in self.expGrayPhotometricCutDict:
793  msg = 'band %s not in expGrayPhotometricCutDict' % (band)
794  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict,
795  self, msg)
796  if band not in self.expGrayHighCutDict:
797  msg = 'band %s not in expGrayHighCutDict' % (band)
798  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict,
799  self, msg)
800  if band not in self.expVarGrayPhotometricCutDict:
801  msg = 'band %s not in expVarGrayPhotometricCutDict' % (band)
802  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict,
803  self, msg)
804  if band not in self.sigFgcmMaxEGrayDict:
805  msg = 'band %s not in sigFgcmMaxEGrayDict' % (band)
806  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict,
807  self, msg)
808  if band not in self.approxThroughputDict:
809  msg = 'band %s not in approxThroughputDict' % (band)
810  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict,
811  self, msg)
812  if band not in self.useRepeatabilityForExpGrayCutsDict:
813  msg = 'band %s not in useRepeatabilityForExpGrayCutsDict' % (band)
814  raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
815  self, msg)
816 
817 
818 class FgcmFitCycleRunner(pipeBase.ButlerInitializedTaskRunner):
819  """Subclass of TaskRunner for fgcmFitCycleTask
820 
821  fgcmFitCycleTask.run() takes one argument, the butler, and uses
822  stars and visits previously extracted from dataRefs by
823  fgcmBuildStars.
824  This Runner does not perform any dataRef parallelization, but the FGCM
825  code called by the Task uses python multiprocessing (see the "ncores"
826  config option).
827  """
828 
829  @staticmethod
830  def getTargetList(parsedCmd):
831  """
832  Return a list with one element, the butler.
833  """
834  return [parsedCmd.butler]
835 
836  def __call__(self, butler):
837  """
838  Parameters
839  ----------
840  butler: `lsst.daf.persistence.Butler`
841 
842  Returns
843  -------
844  exitStatus: `list` with `pipeBase.Struct`
845  exitStatus (0: success; 1: failure)
846  """
847 
848  task = self.TaskClass(config=self.config, log=self.log)
849 
850  exitStatus = 0
851  if self.doRaise:
852  task.runDataRef(butler)
853  else:
854  try:
855  task.runDataRef(butler)
856  except Exception as e:
857  exitStatus = 1
858  task.log.fatal("Failed: %s" % e)
859  if not isinstance(e, pipeBase.TaskError):
860  traceback.print_exc(file=sys.stderr)
861 
862  task.writeMetadata(butler)
863 
864  # The task does not return any results:
865  return [pipeBase.Struct(exitStatus=exitStatus)]
866 
867  def run(self, parsedCmd):
868  """
869  Run the task, with no multiprocessing
870 
871  Parameters
872  ----------
873  parsedCmd: ArgumentParser parsed command line
874  """
875 
876  resultList = []
877 
878  if self.precall(parsedCmd):
879  targetList = self.getTargetList(parsedCmd)
880  # make sure that we only get 1
881  resultList = self(targetList[0])
882 
883  return resultList
884 
885 
886 class FgcmFitCycleTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
887  """
888  Run Single fit cycle for FGCM global calibration
889  """
890 
891  ConfigClass = FgcmFitCycleConfig
892  RunnerClass = FgcmFitCycleRunner
893  _DefaultName = "fgcmFitCycle"
894 
895  def __init__(self, butler=None, initInputs=None, **kwargs):
896  super().__init__(**kwargs)
897 
898  # no saving of metadata for now
899  def _getMetadataName(self):
900  return None
901 
902  def runQuantum(self, butlerQC, inputRefs, outputRefs):
903  camera = butlerQC.get(inputRefs.camera)
904 
905  dataRefDict = {}
906 
907  dataRefDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable)
908  dataRefDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog)
909  dataRefDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations)
910  dataRefDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds)
911  dataRefDict['fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices)
912  if self.config.doReferenceCalibration:
913  dataRefDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars)
914  if self.config.cycleNumber > 0:
915  dataRefDict['fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput)
916  dataRefDict['fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput)
917 
918  fgcmDatasetDict = None
919  if self.config.doMultipleCycles:
920  # Run multiple cycles at once.
921  config = copy.copy(self.config)
922  config.update(cycleNumber=0)
923  for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1):
924  if cycle == self.config.multipleCyclesFinalCycleNumber:
925  config.update(isFinalCycle=True)
926 
927  if cycle > 0:
928  dataRefDict['fgcmFlaggedStars'] = fgcmDatasetDict['fgcmFlaggedStars']
929  dataRefDict['fgcmFitParameters'] = fgcmDatasetDict['fgcmFitParameters']
930 
931  fgcmDatasetDict, config = self._fgcmFitCycle(camera, dataRefDict, config=config)
932  butlerQC.put(fgcmDatasetDict['fgcmFitParameters'],
933  getattr(outputRefs, f'fgcmFitParameters{cycle}'))
934  butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'],
935  getattr(outputRefs, f'fgcmFlaggedStars{cycle}'))
936  if self.outputZeropoints:
937  butlerQC.put(fgcmDatasetDict['fgcmZeropoints'],
938  getattr(outputRefs, f'fgcmZeropoints{cycle}'))
939  butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'],
940  getattr(outputRefs, f'fgcmAtmosphereParameters{cycle}'))
941  if self.outputStandards:
942  butlerQC.put(fgcmDatasetDict['fgcmStandardStars'],
943  getattr(outputRefs, f'fgcmStandardStars{cycle}'))
944  else:
945  # Run a single cycle
946  fgcmDatasetDict, _ = self._fgcmFitCycle(camera, dataRefDict)
947 
948  butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], outputRefs.fgcmFitParameters)
949  butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars)
950  if self.outputZeropoints:
951  butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], outputRefs.fgcmZeropoints)
952  butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters)
953  if self.outputStandards:
954  butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], outputRefs.fgcmStandardStars)
955 
956  @timeMethod
957  def runDataRef(self, butler):
958  """
959  Run a single fit cycle for FGCM
960 
961  Parameters
962  ----------
963  butler: `lsst.daf.persistence.Butler`
964  """
965  self._checkDatasetsExist(butler)
966 
967  dataRefDict = {}
968  dataRefDict['fgcmLookUpTable'] = butler.dataRef('fgcmLookUpTable')
969  dataRefDict['fgcmVisitCatalog'] = butler.dataRef('fgcmVisitCatalog')
970  dataRefDict['fgcmStarObservations'] = butler.dataRef('fgcmStarObservations')
971  dataRefDict['fgcmStarIds'] = butler.dataRef('fgcmStarIds')
972  dataRefDict['fgcmStarIndices'] = butler.dataRef('fgcmStarIndices')
973  if self.config.doReferenceCalibration:
974  dataRefDict['fgcmReferenceStars'] = butler.dataRef('fgcmReferenceStars')
975  if self.config.cycleNumber > 0:
976  lastCycle = self.config.cycleNumber - 1
977  dataRefDict['fgcmFlaggedStars'] = butler.dataRef('fgcmFlaggedStars',
978  fgcmcycle=lastCycle)
979  dataRefDict['fgcmFitParameters'] = butler.dataRef('fgcmFitParameters',
980  fgcmcycle=lastCycle)
981 
982  camera = butler.get('camera')
983  fgcmDatasetDict, _ = self._fgcmFitCycle(camera, dataRefDict)
984 
985  butler.put(fgcmDatasetDict['fgcmFitParameters'], 'fgcmFitParameters',
986  fgcmcycle=self.config.cycleNumber)
987  butler.put(fgcmDatasetDict['fgcmFlaggedStars'], 'fgcmFlaggedStars',
988  fgcmcycle=self.config.cycleNumber)
989  if self.outputZeropoints:
990  butler.put(fgcmDatasetDict['fgcmZeropoints'], 'fgcmZeropoints',
991  fgcmcycle=self.config.cycleNumber)
992  butler.put(fgcmDatasetDict['fgcmAtmosphereParameters'], 'fgcmAtmosphereParameters',
993  fgcmcycle=self.config.cycleNumber)
994  if self.outputStandards:
995  butler.put(fgcmDatasetDict['fgcmStandardStars'], 'fgcmStandardStars',
996  fgcmcycle=self.config.cycleNumber)
997 
998  def writeConfig(self, butler, clobber=False, doBackup=True):
999  """Write the configuration used for processing the data, or check that an existing
1000  one is equal to the new one if present. This is an override of the regular
1001  version from pipe_base that knows about fgcmcycle.
1002 
1003  Parameters
1004  ----------
1005  butler : `lsst.daf.persistence.Butler`
1006  Data butler used to write the config. The config is written to dataset type
1007  `CmdLineTask._getConfigName`.
1008  clobber : `bool`, optional
1009  A boolean flag that controls what happens if a config already has been saved:
1010  - `True`: overwrite or rename the existing config, depending on ``doBackup``.
1011  - `False`: raise `TaskError` if this config does not match the existing config.
1012  doBackup : `bool`, optional
1013  Set to `True` to backup the config files if clobbering.
1014  """
1015  configName = self._getConfigName()
1016  if configName is None:
1017  return
1018  if clobber:
1019  butler.put(self.config, configName, doBackup=doBackup, fgcmcycle=self.config.cycleNumber)
1020  elif butler.datasetExists(configName, write=True, fgcmcycle=self.config.cycleNumber):
1021  # this may be subject to a race condition; see #2789
1022  try:
1023  oldConfig = butler.get(configName, immediate=True, fgcmcycle=self.config.cycleNumber)
1024  except Exception as exc:
1025  raise type(exc)("Unable to read stored config file %s (%s); consider using --clobber-config" %
1026  (configName, exc))
1027 
1028  def logConfigMismatch(msg):
1029  self.log.fatal("Comparing configuration: %s", msg)
1030 
1031  if not self.config.compare(oldConfig, shortcut=False, output=logConfigMismatch):
1032  raise pipeBase.TaskError(
1033  f"Config does not match existing task config {configName!r} on disk; tasks configurations"
1034  " must be consistent within the same output repo (override with --clobber-config)")
1035  else:
1036  butler.put(self.config, configName, fgcmcycle=self.config.cycleNumber)
1037 
1038  def _fgcmFitCycle(self, camera, dataRefDict, config=None):
1039  """
1040  Run the fit cycle
1041 
1042  Parameters
1043  ----------
1044  camera : `lsst.afw.cameraGeom.Camera`
1045  dataRefDict : `dict`
1046  All dataRefs are `lsst.daf.persistence.ButlerDataRef` (gen2) or
1047  `lsst.daf.butler.DeferredDatasetHandle` (gen3)
1048  dataRef dictionary with keys:
1049 
1050  ``"fgcmLookUpTable"``
1051  dataRef for the FGCM look-up table.
1052  ``"fgcmVisitCatalog"``
1053  dataRef for visit summary catalog.
1054  ``"fgcmStarObservations"``
1055  dataRef for star observation catalog.
1056  ``"fgcmStarIds"``
1057  dataRef for star id catalog.
1058  ``"fgcmStarIndices"``
1059  dataRef for star index catalog.
1060  ``"fgcmReferenceStars"``
1061  dataRef for matched reference star catalog.
1062  ``"fgcmFlaggedStars"``
1063  dataRef for flagged star catalog.
1064  ``"fgcmFitParameters"``
1065  dataRef for fit parameter catalog.
1066  config : `lsst.pex.config.Config`, optional
1067  Configuration to use to override self.config.
1068 
1069  Returns
1070  -------
1071  fgcmDatasetDict : `dict`
1072  Dictionary of datasets to persist.
1073  """
1074  if config is not None:
1075  _config = config
1076  else:
1077  _config = self.config
1078 
1079  # Set defaults on whether to output standards and zeropoints
1080  self.maxIter = _config.maxIterBeforeFinalCycle
1081  self.outputStandards = _config.outputStandardsBeforeFinalCycle
1082  self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle
1083  self.resetFitParameters = True
1084 
1085  if _config.isFinalCycle:
1086  # This is the final fit cycle, so we do not want to reset fit
1087  # parameters, we want to run a final "clean-up" with 0 fit iterations,
1088  # and we always want to output standards and zeropoints
1089  self.maxIter = 0
1090  self.outputStandards = True
1091  self.outputZeropoints = True
1092  self.resetFitParameters = False
1093 
1094  lutCat = dataRefDict['fgcmLookUpTable'].get()
1095  fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
1096  dict(_config.physicalFilterMap))
1097  del lutCat
1098 
1099  configDict = makeConfigDict(_config, self.log, camera,
1100  self.maxIter, self.resetFitParameters,
1101  self.outputZeropoints,
1102  lutIndexVals[0]['FILTERNAMES'])
1103 
1104  # next we need the exposure/visit information
1105  visitCat = dataRefDict['fgcmVisitCatalog'].get()
1106  fgcmExpInfo = translateVisitCatalog(visitCat)
1107  del visitCat
1108 
1109  focalPlaneProjector = FocalPlaneProjector(camera,
1110  self.config.defaultCameraOrientation)
1111 
1112  noFitsDict = {'lutIndex': lutIndexVals,
1113  'lutStd': lutStd,
1114  'expInfo': fgcmExpInfo,
1115  'focalPlaneProjector': focalPlaneProjector}
1116 
1117  # set up the fitter object
1118  fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
1119  noFitsDict=noFitsDict, noOutput=True)
1120 
1121  # create the parameter object
1122  if (fgcmFitCycle.initialCycle):
1123  # cycle = 0, initial cycle
1124  fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
1125  fgcmLut,
1126  fgcmExpInfo)
1127  else:
1128  if isinstance(dataRefDict['fgcmFitParameters'], afwTable.BaseCatalog):
1129  parCat = dataRefDict['fgcmFitParameters']
1130  else:
1131  parCat = dataRefDict['fgcmFitParameters'].get()
1132  inParInfo, inParams, inSuperStar = self._loadParameters(parCat)
1133  del parCat
1134  fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
1135  fgcmExpInfo,
1136  inParInfo,
1137  inParams,
1138  inSuperStar)
1139 
1140  # set up the stars...
1141  fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
1142 
1143  starObs = dataRefDict['fgcmStarObservations'].get()
1144  starIds = dataRefDict['fgcmStarIds'].get()
1145  starIndices = dataRefDict['fgcmStarIndices'].get()
1146 
1147  # grab the flagged stars if available
1148  if 'fgcmFlaggedStars' in dataRefDict:
1149  if isinstance(dataRefDict['fgcmFlaggedStars'], afwTable.BaseCatalog):
1150  flaggedStars = dataRefDict['fgcmFlaggedStars']
1151  else:
1152  flaggedStars = dataRefDict['fgcmFlaggedStars'].get()
1153  flagId = flaggedStars['objId'][:]
1154  flagFlag = flaggedStars['objFlag'][:]
1155  else:
1156  flaggedStars = None
1157  flagId = None
1158  flagFlag = None
1159 
1160  if _config.doReferenceCalibration:
1161  refStars = dataRefDict['fgcmReferenceStars'].get()
1162 
1163  refMag, refMagErr = extractReferenceMags(refStars,
1164  _config.bands,
1165  _config.physicalFilterMap)
1166  refId = refStars['fgcm_id'][:]
1167  else:
1168  refStars = None
1169  refId = None
1170  refMag = None
1171  refMagErr = None
1172 
1173  # match star observations to visits
1174  # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
1175  # actually be transferred into fgcm using the indexing below.
1176  visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'][starIndices['obsIndex']])
1177 
1178  # The fgcmStars.loadStars method will copy all the star information into
1179  # special shared memory objects that will not blow up the memory usage when
1180  # used with python multiprocessing. Once all the numbers are copied,
1181  # it is necessary to release all references to the objects that previously
1182  # stored the data to ensure that the garbage collector can clear the memory,
1183  # and ensure that this memory is not copied when multiprocessing kicks in.
1184 
1185  # We determine the conversion from the native units (typically radians) to
1186  # degrees for the first star. This allows us to treat coord_ra/coord_dec as
1187  # numpy arrays rather than Angles, which would we approximately 600x slower.
1188  conv = starObs[0]['ra'].asDegrees() / float(starObs[0]['ra'])
1189 
1190  fgcmStars.loadStars(fgcmPars,
1191  starObs['visit'][starIndices['obsIndex']],
1192  starObs['ccd'][starIndices['obsIndex']],
1193  starObs['ra'][starIndices['obsIndex']] * conv,
1194  starObs['dec'][starIndices['obsIndex']] * conv,
1195  starObs['instMag'][starIndices['obsIndex']],
1196  starObs['instMagErr'][starIndices['obsIndex']],
1197  fgcmExpInfo['FILTERNAME'][visitIndex],
1198  starIds['fgcm_id'][:],
1199  starIds['ra'][:],
1200  starIds['dec'][:],
1201  starIds['obsArrIndex'][:],
1202  starIds['nObs'][:],
1203  obsX=starObs['x'][starIndices['obsIndex']],
1204  obsY=starObs['y'][starIndices['obsIndex']],
1205  obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']],
1206  psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']],
1207  refID=refId,
1208  refMag=refMag,
1209  refMagErr=refMagErr,
1210  flagID=flagId,
1211  flagFlag=flagFlag,
1212  computeNobs=True)
1213 
1214  # Release all references to temporary objects holding star data (see above)
1215  del starObs
1216  del starIds
1217  del starIndices
1218  del flagId
1219  del flagFlag
1220  del flaggedStars
1221  del refStars
1222  del refId
1223  del refMag
1224  del refMagErr
1225 
1226  # and set the bits in the cycle object
1227  fgcmFitCycle.setLUT(fgcmLut)
1228  fgcmFitCycle.setStars(fgcmStars, fgcmPars)
1229  fgcmFitCycle.setPars(fgcmPars)
1230 
1231  # finish the setup
1232  fgcmFitCycle.finishSetup()
1233 
1234  # and run
1235  fgcmFitCycle.run()
1236 
1237 
1240 
1241  fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle)
1242 
1243  # Output the config for the next cycle
1244  # We need to make a copy since the input one has been frozen
1245 
1246  updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i]) for
1247  i, b in enumerate(_config.bands)}
1248  updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i]) for
1249  i, band in enumerate(_config.bands)}
1250 
1251  outConfig = copy.copy(_config)
1252  outConfig.update(cycleNumber=(_config.cycleNumber + 1),
1253  precomputeSuperStarInitialCycle=False,
1254  freezeStdAtmosphere=False,
1255  expGrayPhotometricCutDict=updatedPhotometricCutDict,
1256  expGrayHighCutDict=updatedHighCutDict)
1257 
1258  outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber),
1259  cycleNumber=str(_config.cycleNumber + 1))
1260 
1261  configFileName = '%s_cycle%02d_config.py' % (outConfig.outfileBase,
1262  outConfig.cycleNumber)
1263  outConfig.save(configFileName)
1264 
1265  if _config.isFinalCycle == 1:
1266  # We are done, ready to output products
1267  self.log.info("Everything is in place to run fgcmOutputProducts.py")
1268  else:
1269  self.log.info("Saved config for next cycle to %s" % (configFileName))
1270  self.log.info("Be sure to look at:")
1271  self.log.info(" config.expGrayPhotometricCut")
1272  self.log.info(" config.expGrayHighCut")
1273  self.log.info("If you are satisfied with the fit, please set:")
1274  self.log.info(" config.isFinalCycle = True")
1275 
1276  fgcmFitCycle.freeSharedMemory()
1277 
1278  return fgcmDatasetDict, outConfig
1279 
1280  def _checkDatasetsExist(self, butler):
1281  """
1282  Check if necessary datasets exist to run fgcmFitCycle
1283 
1284  Parameters
1285  ----------
1286  butler: `lsst.daf.persistence.Butler`
1287 
1288  Raises
1289  ------
1290  RuntimeError
1291  If any of fgcmVisitCatalog, fgcmStarObservations, fgcmStarIds,
1292  fgcmStarIndices, fgcmLookUpTable datasets do not exist.
1293  If cycleNumber > 0, then also checks for fgcmFitParameters,
1294  fgcmFlaggedStars.
1295  """
1296 
1297  if not butler.datasetExists('fgcmVisitCatalog'):
1298  raise RuntimeError("Could not find fgcmVisitCatalog in repo!")
1299  if not butler.datasetExists('fgcmStarObservations'):
1300  raise RuntimeError("Could not find fgcmStarObservations in repo!")
1301  if not butler.datasetExists('fgcmStarIds'):
1302  raise RuntimeError("Could not find fgcmStarIds in repo!")
1303  if not butler.datasetExists('fgcmStarIndices'):
1304  raise RuntimeError("Could not find fgcmStarIndices in repo!")
1305  if not butler.datasetExists('fgcmLookUpTable'):
1306  raise RuntimeError("Could not find fgcmLookUpTable in repo!")
1307 
1308  # Need additional datasets if we are not the initial cycle
1309  if (self.config.cycleNumber > 0):
1310  if not butler.datasetExists('fgcmFitParameters',
1311  fgcmcycle=self.config.cycleNumber-1):
1312  raise RuntimeError("Could not find fgcmFitParameters for previous cycle (%d) in repo!" %
1313  (self.config.cycleNumber-1))
1314  if not butler.datasetExists('fgcmFlaggedStars',
1315  fgcmcycle=self.config.cycleNumber-1):
1316  raise RuntimeError("Could not find fgcmFlaggedStars for previous cycle (%d) in repo!" %
1317  (self.config.cycleNumber-1))
1318 
1319  # And additional dataset if we want reference calibration
1320  if self.config.doReferenceCalibration:
1321  if not butler.datasetExists('fgcmReferenceStars'):
1322  raise RuntimeError("Could not find fgcmReferenceStars in repo, and "
1323  "doReferenceCalibration is True.")
1324 
1325  def _loadParameters(self, parCat):
1326  """
1327  Load FGCM parameters from a previous fit cycle
1328 
1329  Parameters
1330  ----------
1331  parCat : `lsst.afw.table.BaseCatalog`
1332  Parameter catalog in afw table form.
1333 
1334  Returns
1335  -------
1336  inParInfo: `numpy.ndarray`
1337  Numpy array parameter information formatted for input to fgcm
1338  inParameters: `numpy.ndarray`
1339  Numpy array parameter values formatted for input to fgcm
1340  inSuperStar: `numpy.array`
1341  Superstar flat formatted for input to fgcm
1342  """
1343  parLutFilterNames = np.array(parCat[0]['lutFilterNames'].split(','))
1344  parFitBands = np.array(parCat[0]['fitBands'].split(','))
1345 
1346  inParInfo = np.zeros(1, dtype=[('NCCD', 'i4'),
1347  ('LUTFILTERNAMES', parLutFilterNames.dtype.str,
1348  (parLutFilterNames.size, )),
1349  ('FITBANDS', parFitBands.dtype.str, (parFitBands.size, )),
1350  ('LNTAUUNIT', 'f8'),
1351  ('LNTAUSLOPEUNIT', 'f8'),
1352  ('ALPHAUNIT', 'f8'),
1353  ('LNPWVUNIT', 'f8'),
1354  ('LNPWVSLOPEUNIT', 'f8'),
1355  ('LNPWVQUADRATICUNIT', 'f8'),
1356  ('LNPWVGLOBALUNIT', 'f8'),
1357  ('O3UNIT', 'f8'),
1358  ('QESYSUNIT', 'f8'),
1359  ('FILTEROFFSETUNIT', 'f8'),
1360  ('HASEXTERNALPWV', 'i2'),
1361  ('HASEXTERNALTAU', 'i2')])
1362  inParInfo['NCCD'] = parCat['nCcd']
1363  inParInfo['LUTFILTERNAMES'][:] = parLutFilterNames
1364  inParInfo['FITBANDS'][:] = parFitBands
1365  inParInfo['HASEXTERNALPWV'] = parCat['hasExternalPwv']
1366  inParInfo['HASEXTERNALTAU'] = parCat['hasExternalTau']
1367 
1368  inParams = np.zeros(1, dtype=[('PARALPHA', 'f8', (parCat['parAlpha'].size, )),
1369  ('PARO3', 'f8', (parCat['parO3'].size, )),
1370  ('PARLNTAUINTERCEPT', 'f8',
1371  (parCat['parLnTauIntercept'].size, )),
1372  ('PARLNTAUSLOPE', 'f8',
1373  (parCat['parLnTauSlope'].size, )),
1374  ('PARLNPWVINTERCEPT', 'f8',
1375  (parCat['parLnPwvIntercept'].size, )),
1376  ('PARLNPWVSLOPE', 'f8',
1377  (parCat['parLnPwvSlope'].size, )),
1378  ('PARLNPWVQUADRATIC', 'f8',
1379  (parCat['parLnPwvQuadratic'].size, )),
1380  ('PARQESYSINTERCEPT', 'f8',
1381  (parCat['parQeSysIntercept'].size, )),
1382  ('COMPQESYSSLOPE', 'f8',
1383  (parCat['compQeSysSlope'].size, )),
1384  ('PARFILTEROFFSET', 'f8',
1385  (parCat['parFilterOffset'].size, )),
1386  ('PARFILTEROFFSETFITFLAG', 'i2',
1387  (parCat['parFilterOffsetFitFlag'].size, )),
1388  ('PARRETRIEVEDLNPWVSCALE', 'f8'),
1389  ('PARRETRIEVEDLNPWVOFFSET', 'f8'),
1390  ('PARRETRIEVEDLNPWVNIGHTLYOFFSET', 'f8',
1391  (parCat['parRetrievedLnPwvNightlyOffset'].size, )),
1392  ('COMPABSTHROUGHPUT', 'f8',
1393  (parCat['compAbsThroughput'].size, )),
1394  ('COMPREFOFFSET', 'f8',
1395  (parCat['compRefOffset'].size, )),
1396  ('COMPREFSIGMA', 'f8',
1397  (parCat['compRefSigma'].size, )),
1398  ('COMPMIRRORCHROMATICITY', 'f8',
1399  (parCat['compMirrorChromaticity'].size, )),
1400  ('MIRRORCHROMATICITYPIVOT', 'f8',
1401  (parCat['mirrorChromaticityPivot'].size, )),
1402  ('COMPMEDIANSEDSLOPE', 'f8',
1403  (parCat['compMedianSedSlope'].size, )),
1404  ('COMPAPERCORRPIVOT', 'f8',
1405  (parCat['compAperCorrPivot'].size, )),
1406  ('COMPAPERCORRSLOPE', 'f8',
1407  (parCat['compAperCorrSlope'].size, )),
1408  ('COMPAPERCORRSLOPEERR', 'f8',
1409  (parCat['compAperCorrSlopeErr'].size, )),
1410  ('COMPAPERCORRRANGE', 'f8',
1411  (parCat['compAperCorrRange'].size, )),
1412  ('COMPMODELERREXPTIMEPIVOT', 'f8',
1413  (parCat['compModelErrExptimePivot'].size, )),
1414  ('COMPMODELERRFWHMPIVOT', 'f8',
1415  (parCat['compModelErrFwhmPivot'].size, )),
1416  ('COMPMODELERRSKYPIVOT', 'f8',
1417  (parCat['compModelErrSkyPivot'].size, )),
1418  ('COMPMODELERRPARS', 'f8',
1419  (parCat['compModelErrPars'].size, )),
1420  ('COMPEXPGRAY', 'f8',
1421  (parCat['compExpGray'].size, )),
1422  ('COMPVARGRAY', 'f8',
1423  (parCat['compVarGray'].size, )),
1424  ('COMPEXPDELTAMAGBKG', 'f8',
1425  (parCat['compExpDeltaMagBkg'].size, )),
1426  ('COMPNGOODSTARPEREXP', 'i4',
1427  (parCat['compNGoodStarPerExp'].size, )),
1428  ('COMPSIGFGCM', 'f8',
1429  (parCat['compSigFgcm'].size, )),
1430  ('COMPSIGMACAL', 'f8',
1431  (parCat['compSigmaCal'].size, )),
1432  ('COMPRETRIEVEDLNPWV', 'f8',
1433  (parCat['compRetrievedLnPwv'].size, )),
1434  ('COMPRETRIEVEDLNPWVRAW', 'f8',
1435  (parCat['compRetrievedLnPwvRaw'].size, )),
1436  ('COMPRETRIEVEDLNPWVFLAG', 'i2',
1437  (parCat['compRetrievedLnPwvFlag'].size, )),
1438  ('COMPRETRIEVEDTAUNIGHT', 'f8',
1439  (parCat['compRetrievedTauNight'].size, ))])
1440 
1441  inParams['PARALPHA'][:] = parCat['parAlpha'][0, :]
1442  inParams['PARO3'][:] = parCat['parO3'][0, :]
1443  inParams['PARLNTAUINTERCEPT'][:] = parCat['parLnTauIntercept'][0, :]
1444  inParams['PARLNTAUSLOPE'][:] = parCat['parLnTauSlope'][0, :]
1445  inParams['PARLNPWVINTERCEPT'][:] = parCat['parLnPwvIntercept'][0, :]
1446  inParams['PARLNPWVSLOPE'][:] = parCat['parLnPwvSlope'][0, :]
1447  inParams['PARLNPWVQUADRATIC'][:] = parCat['parLnPwvQuadratic'][0, :]
1448  inParams['PARQESYSINTERCEPT'][:] = parCat['parQeSysIntercept'][0, :]
1449  inParams['COMPQESYSSLOPE'][:] = parCat['compQeSysSlope'][0, :]
1450  inParams['PARFILTEROFFSET'][:] = parCat['parFilterOffset'][0, :]
1451  inParams['PARFILTEROFFSETFITFLAG'][:] = parCat['parFilterOffsetFitFlag'][0, :]
1452  inParams['PARRETRIEVEDLNPWVSCALE'] = parCat['parRetrievedLnPwvScale']
1453  inParams['PARRETRIEVEDLNPWVOFFSET'] = parCat['parRetrievedLnPwvOffset']
1454  inParams['PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat['parRetrievedLnPwvNightlyOffset'][0, :]
1455  inParams['COMPABSTHROUGHPUT'][:] = parCat['compAbsThroughput'][0, :]
1456  inParams['COMPREFOFFSET'][:] = parCat['compRefOffset'][0, :]
1457  inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :]
1458  inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :]
1459  inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :]
1460  inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :]
1461  inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :]
1462  inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :]
1463  inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :]
1464  inParams['COMPAPERCORRRANGE'][:] = parCat['compAperCorrRange'][0, :]
1465  inParams['COMPMODELERREXPTIMEPIVOT'][:] = parCat['compModelErrExptimePivot'][0, :]
1466  inParams['COMPMODELERRFWHMPIVOT'][:] = parCat['compModelErrFwhmPivot'][0, :]
1467  inParams['COMPMODELERRSKYPIVOT'][:] = parCat['compModelErrSkyPivot'][0, :]
1468  inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :]
1469  inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :]
1470  inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
1471  inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
1472  inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
1473  inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
1474  inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
1475  inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :]
1476  inParams['COMPRETRIEVEDLNPWVRAW'][:] = parCat['compRetrievedLnPwvRaw'][0, :]
1477  inParams['COMPRETRIEVEDLNPWVFLAG'][:] = parCat['compRetrievedLnPwvFlag'][0, :]
1478  inParams['COMPRETRIEVEDTAUNIGHT'][:] = parCat['compRetrievedTauNight'][0, :]
1479 
1480  inSuperStar = np.zeros(parCat['superstarSize'][0, :], dtype='f8')
1481  inSuperStar[:, :, :, :] = parCat['superstar'][0, :].reshape(inSuperStar.shape)
1482 
1483  return (inParInfo, inParams, inSuperStar)
1484 
1485  def _makeFgcmOutputDatasets(self, fgcmFitCycle):
1486  """
1487  Persist FGCM datasets through the butler.
1488 
1489  Parameters
1490  ----------
1491  fgcmFitCycle: `lsst.fgcm.FgcmFitCycle`
1492  Fgcm Fit cycle object
1493  """
1494  fgcmDatasetDict = {}
1495 
1496  # Save the parameters
1497  parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays()
1498 
1499  parSchema = afwTable.Schema()
1500 
1501  comma = ','
1502  lutFilterNameString = comma.join([n.decode('utf-8')
1503  for n in parInfo['LUTFILTERNAMES'][0]])
1504  fitBandString = comma.join([n.decode('utf-8')
1505  for n in parInfo['FITBANDS'][0]])
1506 
1507  parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat,
1508  lutFilterNameString, fitBandString)
1509  parCat = self._makeParCatalog(parSchema, parInfo, pars,
1510  fgcmFitCycle.fgcmPars.parSuperStarFlat,
1511  lutFilterNameString, fitBandString)
1512 
1513  fgcmDatasetDict['fgcmFitParameters'] = parCat
1514 
1515  # Save the indices of the flagged stars
1516  # (stars that have been (a) reserved from the fit for testing and
1517  # (b) bad stars that have failed quality checks.)
1518  flagStarSchema = self._makeFlagStarSchema()
1519  flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices()
1520  flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct)
1521 
1522  fgcmDatasetDict['fgcmFlaggedStars'] = flagStarCat
1523 
1524  # Save the zeropoint information and atmospheres only if desired
1525  if self.outputZeropoints:
1526  superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
1527  zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
1528 
1529  zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
1530  zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
1531 
1532  fgcmDatasetDict['fgcmZeropoints'] = zptCat
1533 
1534  # Save atmosphere values
1535  # These are generated by the same code that generates zeropoints
1536  atmSchema = makeAtmSchema()
1537  atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
1538 
1539  fgcmDatasetDict['fgcmAtmosphereParameters'] = atmCat
1540 
1541  # Save the standard stars (if configured)
1542  if self.outputStandards:
1543  stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
1544  stdSchema = makeStdSchema(len(goodBands))
1545  stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
1546 
1547  fgcmDatasetDict['fgcmStandardStars'] = stdCat
1548 
1549  return fgcmDatasetDict
1550 
1551  def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
1552  lutFilterNameString, fitBandString):
1553  """
1554  Make the parameter persistence schema
1555 
1556  Parameters
1557  ----------
1558  parInfo: `numpy.ndarray`
1559  Parameter information returned by fgcm
1560  pars: `numpy.ndarray`
1561  Parameter values returned by fgcm
1562  parSuperStarFlat: `numpy.array`
1563  Superstar flat values returned by fgcm
1564  lutFilterNameString: `str`
1565  Combined string of all the lutFilterNames
1566  fitBandString: `str`
1567  Combined string of all the fitBands
1568 
1569  Returns
1570  -------
1571  parSchema: `afwTable.schema`
1572  """
1573 
1574  parSchema = afwTable.Schema()
1575 
1576  # parameter info section
1577  parSchema.addField('nCcd', type=np.int32, doc='Number of CCDs')
1578  parSchema.addField('lutFilterNames', type=str, doc='LUT Filter names in parameter file',
1579  size=len(lutFilterNameString))
1580  parSchema.addField('fitBands', type=str, doc='Bands that were fit',
1581  size=len(fitBandString))
1582  parSchema.addField('lnTauUnit', type=np.float64, doc='Step units for ln(AOD)')
1583  parSchema.addField('lnTauSlopeUnit', type=np.float64,
1584  doc='Step units for ln(AOD) slope')
1585  parSchema.addField('alphaUnit', type=np.float64, doc='Step units for alpha')
1586  parSchema.addField('lnPwvUnit', type=np.float64, doc='Step units for ln(pwv)')
1587  parSchema.addField('lnPwvSlopeUnit', type=np.float64,
1588  doc='Step units for ln(pwv) slope')
1589  parSchema.addField('lnPwvQuadraticUnit', type=np.float64,
1590  doc='Step units for ln(pwv) quadratic term')
1591  parSchema.addField('lnPwvGlobalUnit', type=np.float64,
1592  doc='Step units for global ln(pwv) parameters')
1593  parSchema.addField('o3Unit', type=np.float64, doc='Step units for O3')
1594  parSchema.addField('qeSysUnit', type=np.float64, doc='Step units for mirror gray')
1595  parSchema.addField('filterOffsetUnit', type=np.float64, doc='Step units for filter offset')
1596  parSchema.addField('hasExternalPwv', type=np.int32, doc='Parameters fit using external pwv')
1597  parSchema.addField('hasExternalTau', type=np.int32, doc='Parameters fit using external tau')
1598 
1599  # parameter section
1600  parSchema.addField('parAlpha', type='ArrayD', doc='Alpha parameter vector',
1601  size=pars['PARALPHA'].size)
1602  parSchema.addField('parO3', type='ArrayD', doc='O3 parameter vector',
1603  size=pars['PARO3'].size)
1604  parSchema.addField('parLnTauIntercept', type='ArrayD',
1605  doc='ln(Tau) intercept parameter vector',
1606  size=pars['PARLNTAUINTERCEPT'].size)
1607  parSchema.addField('parLnTauSlope', type='ArrayD',
1608  doc='ln(Tau) slope parameter vector',
1609  size=pars['PARLNTAUSLOPE'].size)
1610  parSchema.addField('parLnPwvIntercept', type='ArrayD', doc='ln(pwv) intercept parameter vector',
1611  size=pars['PARLNPWVINTERCEPT'].size)
1612  parSchema.addField('parLnPwvSlope', type='ArrayD', doc='ln(pwv) slope parameter vector',
1613  size=pars['PARLNPWVSLOPE'].size)
1614  parSchema.addField('parLnPwvQuadratic', type='ArrayD', doc='ln(pwv) quadratic parameter vector',
1615  size=pars['PARLNPWVQUADRATIC'].size)
1616  parSchema.addField('parQeSysIntercept', type='ArrayD', doc='Mirror gray intercept parameter vector',
1617  size=pars['PARQESYSINTERCEPT'].size)
1618  parSchema.addField('compQeSysSlope', type='ArrayD', doc='Mirror gray slope parameter vector',
1619  size=pars[0]['COMPQESYSSLOPE'].size)
1620  parSchema.addField('parFilterOffset', type='ArrayD', doc='Filter offset parameter vector',
1621  size=pars['PARFILTEROFFSET'].size)
1622  parSchema.addField('parFilterOffsetFitFlag', type='ArrayI', doc='Filter offset parameter fit flag',
1623  size=pars['PARFILTEROFFSETFITFLAG'].size)
1624  parSchema.addField('parRetrievedLnPwvScale', type=np.float64,
1625  doc='Global scale for retrieved ln(pwv)')
1626  parSchema.addField('parRetrievedLnPwvOffset', type=np.float64,
1627  doc='Global offset for retrieved ln(pwv)')
1628  parSchema.addField('parRetrievedLnPwvNightlyOffset', type='ArrayD',
1629  doc='Nightly offset for retrieved ln(pwv)',
1630  size=pars['PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size)
1631  parSchema.addField('compAbsThroughput', type='ArrayD',
1632  doc='Absolute throughput (relative to transmission curves)',
1633  size=pars['COMPABSTHROUGHPUT'].size)
1634  parSchema.addField('compRefOffset', type='ArrayD',
1635  doc='Offset between reference stars and calibrated stars',
1636  size=pars['COMPREFOFFSET'].size)
1637  parSchema.addField('compRefSigma', type='ArrayD',
1638  doc='Width of reference star/calibrated star distribution',
1639  size=pars['COMPREFSIGMA'].size)
1640  parSchema.addField('compMirrorChromaticity', type='ArrayD',
1641  doc='Computed mirror chromaticity terms',
1642  size=pars['COMPMIRRORCHROMATICITY'].size)
1643  parSchema.addField('mirrorChromaticityPivot', type='ArrayD',
1644  doc='Mirror chromaticity pivot mjd',
1645  size=pars['MIRRORCHROMATICITYPIVOT'].size)
1646  parSchema.addField('compMedianSedSlope', type='ArrayD',
1647  doc='Computed median SED slope (per band)',
1648  size=pars['COMPMEDIANSEDSLOPE'].size)
1649  parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot',
1650  size=pars['COMPAPERCORRPIVOT'].size)
1651  parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope',
1652  size=pars['COMPAPERCORRSLOPE'].size)
1653  parSchema.addField('compAperCorrSlopeErr', type='ArrayD', doc='Aperture correction slope error',
1654  size=pars['COMPAPERCORRSLOPEERR'].size)
1655  parSchema.addField('compAperCorrRange', type='ArrayD', doc='Aperture correction range',
1656  size=pars['COMPAPERCORRRANGE'].size)
1657  parSchema.addField('compModelErrExptimePivot', type='ArrayD', doc='Model error exptime pivot',
1658  size=pars['COMPMODELERREXPTIMEPIVOT'].size)
1659  parSchema.addField('compModelErrFwhmPivot', type='ArrayD', doc='Model error fwhm pivot',
1660  size=pars['COMPMODELERRFWHMPIVOT'].size)
1661  parSchema.addField('compModelErrSkyPivot', type='ArrayD', doc='Model error sky pivot',
1662  size=pars['COMPMODELERRSKYPIVOT'].size)
1663  parSchema.addField('compModelErrPars', type='ArrayD', doc='Model error parameters',
1664  size=pars['COMPMODELERRPARS'].size)
1665  parSchema.addField('compExpGray', type='ArrayD', doc='Computed exposure gray',
1666  size=pars['COMPEXPGRAY'].size)
1667  parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance',
1668  size=pars['COMPVARGRAY'].size)
1669  parSchema.addField('compExpDeltaMagBkg', type='ArrayD',
1670  doc='Computed exposure offset due to background',
1671  size=pars['COMPEXPDELTAMAGBKG'].size)
1672  parSchema.addField('compNGoodStarPerExp', type='ArrayI',
1673  doc='Computed number of good stars per exposure',
1674  size=pars['COMPNGOODSTARPEREXP'].size)
1675  parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)',
1676  size=pars['COMPSIGFGCM'].size)
1677  parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)',
1678  size=pars['COMPSIGMACAL'].size)
1679  parSchema.addField('compRetrievedLnPwv', type='ArrayD', doc='Retrieved ln(pwv) (smoothed)',
1680  size=pars['COMPRETRIEVEDLNPWV'].size)
1681  parSchema.addField('compRetrievedLnPwvRaw', type='ArrayD', doc='Retrieved ln(pwv) (raw)',
1682  size=pars['COMPRETRIEVEDLNPWVRAW'].size)
1683  parSchema.addField('compRetrievedLnPwvFlag', type='ArrayI', doc='Retrieved ln(pwv) Flag',
1684  size=pars['COMPRETRIEVEDLNPWVFLAG'].size)
1685  parSchema.addField('compRetrievedTauNight', type='ArrayD', doc='Retrieved tau (per night)',
1686  size=pars['COMPRETRIEVEDTAUNIGHT'].size)
1687  # superstarflat section
1688  parSchema.addField('superstarSize', type='ArrayI', doc='Superstar matrix size',
1689  size=4)
1690  parSchema.addField('superstar', type='ArrayD', doc='Superstar matrix (flattened)',
1691  size=parSuperStarFlat.size)
1692 
1693  return parSchema
1694 
1695  def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
1696  lutFilterNameString, fitBandString):
1697  """
1698  Make the FGCM parameter catalog for persistence
1699 
1700  Parameters
1701  ----------
1702  parSchema: `lsst.afw.table.Schema`
1703  Parameter catalog schema
1704  pars: `numpy.ndarray`
1705  FGCM parameters to put into parCat
1706  parSuperStarFlat: `numpy.array`
1707  FGCM superstar flat array to put into parCat
1708  lutFilterNameString: `str`
1709  Combined string of all the lutFilterNames
1710  fitBandString: `str`
1711  Combined string of all the fitBands
1712 
1713  Returns
1714  -------
1715  parCat: `afwTable.BasicCatalog`
1716  Atmosphere and instrumental model parameter catalog for persistence
1717  """
1718 
1719  parCat = afwTable.BaseCatalog(parSchema)
1720  parCat.reserve(1)
1721 
1722  # The parameter catalog just has one row, with many columns for all the
1723  # atmosphere and instrument fit parameters
1724  rec = parCat.addNew()
1725 
1726  # info section
1727  rec['nCcd'] = parInfo['NCCD']
1728  rec['lutFilterNames'] = lutFilterNameString
1729  rec['fitBands'] = fitBandString
1730  # note these are not currently supported here.
1731  rec['hasExternalPwv'] = 0
1732  rec['hasExternalTau'] = 0
1733 
1734  # parameter section
1735 
1736  scalarNames = ['parRetrievedLnPwvScale', 'parRetrievedLnPwvOffset']
1737 
1738  arrNames = ['parAlpha', 'parO3', 'parLnTauIntercept', 'parLnTauSlope',
1739  'parLnPwvIntercept', 'parLnPwvSlope', 'parLnPwvQuadratic',
1740  'parQeSysIntercept', 'compQeSysSlope',
1741  'parRetrievedLnPwvNightlyOffset', 'compAperCorrPivot',
1742  'parFilterOffset', 'parFilterOffsetFitFlag',
1743  'compAbsThroughput', 'compRefOffset', 'compRefSigma',
1744  'compMirrorChromaticity', 'mirrorChromaticityPivot',
1745  'compAperCorrSlope', 'compAperCorrSlopeErr', 'compAperCorrRange',
1746  'compModelErrExptimePivot', 'compModelErrFwhmPivot',
1747  'compModelErrSkyPivot', 'compModelErrPars',
1748  'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm',
1749  'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope',
1750  'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
1751  'compRetrievedTauNight']
1752 
1753  for scalarName in scalarNames:
1754  rec[scalarName] = pars[scalarName.upper()]
1755 
1756  for arrName in arrNames:
1757  rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:]
1758 
1759  # superstar section
1760  rec['superstarSize'][:] = parSuperStarFlat.shape
1761  rec['superstar'][:] = parSuperStarFlat.flatten()
1762 
1763  return parCat
1764 
1765  def _makeFlagStarSchema(self):
1766  """
1767  Make the flagged-stars schema
1768 
1769  Returns
1770  -------
1771  flagStarSchema: `lsst.afw.table.Schema`
1772  """
1773 
1774  flagStarSchema = afwTable.Schema()
1775 
1776  flagStarSchema.addField('objId', type=np.int32, doc='FGCM object id')
1777  flagStarSchema.addField('objFlag', type=np.int32, doc='FGCM object flag')
1778 
1779  return flagStarSchema
1780 
1781  def _makeFlagStarCat(self, flagStarSchema, flagStarStruct):
1782  """
1783  Make the flagged star catalog for persistence
1784 
1785  Parameters
1786  ----------
1787  flagStarSchema: `lsst.afw.table.Schema`
1788  Flagged star schema
1789  flagStarStruct: `numpy.ndarray`
1790  Flagged star structure from fgcm
1791 
1792  Returns
1793  -------
1794  flagStarCat: `lsst.afw.table.BaseCatalog`
1795  Flagged star catalog for persistence
1796  """
1797 
1798  flagStarCat = afwTable.BaseCatalog(flagStarSchema)
1799  flagStarCat.resize(flagStarStruct.size)
1800 
1801  flagStarCat['objId'][:] = flagStarStruct['OBJID']
1802  flagStarCat['objFlag'][:] = flagStarStruct['OBJFLAG']
1803 
1804  return flagStarCat
table::Key< int > type
Definition: Detector.cc:163
Defines the fields and offsets for a table.
Definition: Schema.h:51
def extractReferenceMags(refStars, bands, filterMap)
Definition: utilities.py:887
def makeStdSchema(nBands)
Definition: utilities.py:742
def makeAtmCat(atmSchema, atmStruct)
Definition: utilities.py:709
def makeConfigDict(config, log, camera, maxIter, resetFitParameters, outputZeropoints, lutFilterNames, tract=None)
Definition: utilities.py:52
def translateFgcmLut(lutCat, physicalFilterMap)
Definition: utilities.py:220
def makeZptCat(zptSchema, zpStruct)
Definition: utilities.py:630
def makeStdCat(stdSchema, stdStruct, goodBands)
Definition: utilities.py:774
def makeZptSchema(superStarChebyshevSize, zptChebyshevSize)
Definition: utilities.py:544
def translateVisitCatalog(visitCat)
Definition: utilities.py:350
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603