Coverage for python/lsst/fgcmcal/fgcmCalibrateTractBase.py: 14%
157 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-09 04:36 -0700
1# This file is part of fgcmcal.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""Base class for running fgcmcal on a single tract using src tables
22or sourceTable_visit tables.
23"""
25import abc
27import numpy as np
29import lsst.pex.config as pexConfig
30import lsst.pipe.base as pipeBase
32from .fgcmBuildStarsTable import FgcmBuildStarsTableTask
33from .fgcmFitCycle import FgcmFitCycleConfig
34from .fgcmOutputProducts import FgcmOutputProductsTask
35from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
36from .utilities import computeApertureRadiusFromName, extractReferenceMags
37from .utilities import makeZptSchema, makeZptCat
38from .utilities import makeAtmSchema, makeAtmCat
39from .utilities import makeStdSchema, makeStdCat
40from .focalPlaneProjector import FocalPlaneProjector
42import fgcm
44__all__ = ['FgcmCalibrateTractConfigBase', 'FgcmCalibrateTractBaseTask']
47class FgcmCalibrateTractConfigBase(pexConfig.Config):
48 """Config for FgcmCalibrateTract"""
50 fgcmBuildStars = pexConfig.ConfigurableField(
51 target=FgcmBuildStarsTableTask,
52 doc="Task to load and match stars for fgcm",
53 )
54 fgcmFitCycle = pexConfig.ConfigField(
55 dtype=FgcmFitCycleConfig,
56 doc="Config to run a single fgcm fit cycle",
57 )
58 fgcmOutputProducts = pexConfig.ConfigurableField(
59 target=FgcmOutputProductsTask,
60 doc="Task to output fgcm products",
61 )
62 convergenceTolerance = pexConfig.Field(
63 doc="Tolerance on repeatability convergence (per band)",
64 dtype=float,
65 default=0.005,
66 )
67 maxFitCycles = pexConfig.Field(
68 doc="Maximum number of fit cycles",
69 dtype=int,
70 default=5,
71 )
72 doDebuggingPlots = pexConfig.Field(
73 doc="Make plots for debugging purposes?",
74 dtype=bool,
75 default=False,
76 )
78 def setDefaults(self):
79 pexConfig.Config.setDefaults(self)
81 self.fgcmFitCycle.quietMode = True
82 self.fgcmFitCycle.doPlots = False
83 self.fgcmOutputProducts.doReferenceCalibration = False
84 self.fgcmOutputProducts.photoCal.applyColorTerms = False
86 def validate(self):
87 super().validate()
89 for band in self.fgcmFitCycle.bands:
90 if not self.fgcmFitCycle.useRepeatabilityForExpGrayCutsDict[band]:
91 msg = 'Must set useRepeatabilityForExpGrayCutsDict[band]=True for all bands'
92 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
93 self, msg)
96class FgcmCalibrateTractBaseTask(pipeBase.PipelineTask, abc.ABC):
97 """Base class to calibrate a single tract using fgcmcal
98 """
99 def __init__(self, initInputs=None, **kwargs):
100 super().__init__(**kwargs)
101 self.makeSubtask("fgcmBuildStars", initInputs=initInputs)
102 self.makeSubtask("fgcmOutputProducts")
104 def run(self, handleDict, tract,
105 buildStarsRefObjLoader=None, returnCatalogs=True):
106 """Run the calibrations for a single tract with fgcm.
108 Parameters
109 ----------
110 handleDict : `dict`
111 All handles are `lsst.daf.butler.DeferredDatasetHandle`
112 handle dictionary with the following keys. Note that all
113 keys need not be set based on config parameters.
115 ``"camera"``
116 Camera object (`lsst.afw.cameraGeom.Camera`)
117 ``"source_catalogs"``
118 `list` of handles for input source catalogs.
119 ``"sourceSchema"``
120 Schema for the source catalogs.
121 ``"fgcmLookUpTable"``
122 handle for the FGCM look-up table.
123 ``"calexps"``
124 `list` of handles for the input calexps
125 ``"fgcmPhotoCalibs"``
126 `dict` of output photoCalib handles. Key is
127 (tract, visit, detector).
128 Present if doZeropointOutput is True.
129 ``"fgcmTransmissionAtmospheres"``
130 `dict` of output atmosphere transmission handles.
131 Key is (tract, visit).
132 Present if doAtmosphereOutput is True.
133 tract : `int`
134 Tract number
135 buildStarsRefObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`, optional
136 Reference object loader object for fgcmBuildStars.
137 returnCatalogs : `bool`, optional
138 Return photoCalibs as per-visit exposure catalogs.
140 Returns
141 -------
142 outstruct : `lsst.pipe.base.Struct`
143 Output structure with keys:
145 offsets : `np.ndarray`
146 Final reference offsets, per band.
147 repeatability : `np.ndarray`
148 Raw fgcm repeatability for bright stars, per band.
149 atmospheres : `generator` [(`int`, `lsst.afw.image.TransmissionCurve`)]
150 Generator that returns (visit, transmissionCurve) tuples.
151 photoCalibs : `generator` [(`int`, `int`, `str`, `lsst.afw.image.PhotoCalib`)]
152 Generator that returns (visit, ccd, filtername, photoCalib) tuples.
153 (returned if returnCatalogs is False).
154 photoCalibCatalogs : `generator` [(`int`, `lsst.afw.table.ExposureCatalog`)]
155 Generator that returns (visit, exposureCatalog) tuples.
156 (returned if returnCatalogs is True).
157 """
158 self.log.info("Running on tract %d", (tract))
160 # Compute the aperture radius if necessary. This is useful to do now before
161 # any heavy lifting has happened (fail early).
162 calibFluxApertureRadius = None
163 if self.config.fgcmBuildStars.doSubtractLocalBackground:
164 try:
165 field = self.config.fgcmBuildStars.instFluxField
166 calibFluxApertureRadius = computeApertureRadiusFromName(field)
167 except RuntimeError:
168 raise RuntimeError("Could not determine aperture radius from %s. "
169 "Cannot use doSubtractLocalBackground." %
170 (field))
172 # Run the build stars tasks
174 # Note that we will need visitCat at the end of the procedure for the outputs
175 groupedHandles = self.fgcmBuildStars._groupHandles(handleDict['sourceTableHandleDict'],
176 handleDict['visitSummaryHandleDict'])
178 lutCat = handleDict["fgcmLookUpTable"].get()
179 camera = handleDict["camera"]
180 if len(camera) == lutCat[0]["nCcd"]:
181 useScienceDetectors = False
182 else:
183 # If the LUT has a different number of detectors than
184 # the camera, then we only want to use science detectors
185 # in the focal plane projector.
186 useScienceDetectors = True
187 del lutCat
189 visitCat = self.fgcmBuildStars.fgcmMakeVisitCatalog(
190 camera,
191 groupedHandles,
192 useScienceDetectors=useScienceDetectors,
193 )
194 rad = calibFluxApertureRadius
195 fgcmStarObservationCat = self.fgcmBuildStars.fgcmMakeAllStarObservations(groupedHandles,
196 visitCat,
197 handleDict['sourceSchema'],
198 handleDict['camera'],
199 calibFluxApertureRadius=rad)
201 if self.fgcmBuildStars.config.doReferenceMatches:
202 lutHandle = handleDict['fgcmLookUpTable']
203 self.fgcmBuildStars.makeSubtask("fgcmLoadReferenceCatalog",
204 refObjLoader=buildStarsRefObjLoader,
205 refCatName=self.fgcmBuildStars.config.connections.refCat)
206 else:
207 lutHandle = None
209 fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat = \
210 self.fgcmBuildStars.fgcmMatchStars(visitCat,
211 fgcmStarObservationCat,
212 lutHandle=lutHandle)
214 # Load the LUT
215 lutCat = handleDict['fgcmLookUpTable'].get()
216 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
217 dict(self.config.fgcmFitCycle.physicalFilterMap))
218 del lutCat
220 # Translate the visit catalog into fgcm format
221 fgcmExpInfo = translateVisitCatalog(visitCat)
223 configDict = makeConfigDict(self.config.fgcmFitCycle, self.log, handleDict['camera'],
224 self.config.fgcmFitCycle.maxIterBeforeFinalCycle,
225 True, False, lutIndexVals[0]['FILTERNAMES'],
226 tract=tract)
228 focalPlaneProjector = FocalPlaneProjector(handleDict['camera'],
229 self.config.fgcmFitCycle.defaultCameraOrientation)
231 # Set up the fit cycle task
233 noFitsDict = {'lutIndex': lutIndexVals,
234 'lutStd': lutStd,
235 'expInfo': fgcmExpInfo,
236 'focalPlaneProjector': focalPlaneProjector}
238 fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
239 noFitsDict=noFitsDict, noOutput=True)
241 # We determine the conversion from the native units (typically radians) to
242 # degrees for the first star. This allows us to treat coord_ra/coord_dec as
243 # numpy arrays rather than Angles, which would we approximately 600x slower.
244 conv = fgcmStarObservationCat[0]['ra'].asDegrees() / float(fgcmStarObservationCat[0]['ra'])
246 # To load the stars, we need an initial parameter object
247 fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
248 fgcmLut,
249 fgcmExpInfo)
251 # Match star observations to visits
252 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
253 # actually be transferred into fgcm using the indexing below.
255 obsIndex = fgcmStarIndicesCat['obsIndex']
256 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'],
257 fgcmStarObservationCat['visit'][obsIndex])
259 refMag, refMagErr = extractReferenceMags(fgcmRefCat,
260 self.config.fgcmFitCycle.bands,
261 self.config.fgcmFitCycle.physicalFilterMap)
262 refId = fgcmRefCat['fgcm_id'][:]
264 fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
265 fgcmStars.loadStars(fgcmPars,
266 fgcmStarObservationCat['visit'][obsIndex],
267 fgcmStarObservationCat['ccd'][obsIndex],
268 fgcmStarObservationCat['ra'][obsIndex] * conv,
269 fgcmStarObservationCat['dec'][obsIndex] * conv,
270 fgcmStarObservationCat['instMag'][obsIndex],
271 fgcmStarObservationCat['instMagErr'][obsIndex],
272 fgcmExpInfo['FILTERNAME'][visitIndex],
273 fgcmStarIdCat['fgcm_id'][:],
274 fgcmStarIdCat['ra'][:],
275 fgcmStarIdCat['dec'][:],
276 fgcmStarIdCat['obsArrIndex'][:],
277 fgcmStarIdCat['nObs'][:],
278 obsX=fgcmStarObservationCat['x'][obsIndex],
279 obsY=fgcmStarObservationCat['y'][obsIndex],
280 obsDeltaMagBkg=fgcmStarObservationCat['deltaMagBkg'][obsIndex],
281 obsDeltaAper=fgcmStarObservationCat['deltaMagAper'][obsIndex],
282 psfCandidate=fgcmStarObservationCat['psf_candidate'][obsIndex],
283 refID=refId,
284 refMag=refMag,
285 refMagErr=refMagErr,
286 flagID=None,
287 flagFlag=None,
288 computeNobs=True)
290 # Clear out some memory
291 del fgcmStarIdCat
292 del fgcmStarIndicesCat
293 del fgcmRefCat
295 fgcmFitCycle.setLUT(fgcmLut)
296 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
298 converged = False
299 cycleNumber = 0
301 previousReservedRawRepeatability = np.zeros(fgcmPars.nBands) + 1000.0
302 previousParInfo = None
303 previousParams = None
304 previousSuperStar = None
306 while (not converged and cycleNumber < self.config.maxFitCycles):
308 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
310 if cycleNumber > 0:
311 # Use parameters from previous cycle
312 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
313 fgcmExpInfo,
314 previousParInfo,
315 previousParams,
316 previousSuperStar)
317 # We need to reset the star magnitudes and errors for the next
318 # cycle
319 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
320 fgcmStarObservationCat['instMagErr'][obsIndex])
321 fgcmFitCycle.initialCycle = False
323 fgcmFitCycle.setPars(fgcmPars)
324 fgcmFitCycle.finishSetup()
326 fgcmFitCycle.run()
328 # Grab the parameters for the next cycle
329 previousParInfo, previousParams = fgcmFitCycle.fgcmPars.parsToArrays()
330 previousSuperStar = fgcmFitCycle.fgcmPars.parSuperStarFlat.copy()
332 self.log.info("Raw repeatability after cycle number %d is:" % (cycleNumber))
333 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
334 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
335 continue
336 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
337 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
339 # Check for convergence
340 if np.all((previousReservedRawRepeatability
341 - fgcmFitCycle.fgcmPars.compReservedRawRepeatability)
342 < self.config.convergenceTolerance):
343 self.log.info("Raw repeatability has converged after cycle number %d." % (cycleNumber))
344 converged = True
345 else:
346 fgcmFitCycle.fgcmConfig.expGrayPhotometricCut[:] = fgcmFitCycle.updatedPhotometricCut
347 fgcmFitCycle.fgcmConfig.expGrayHighCut[:] = fgcmFitCycle.updatedHighCut
348 fgcmFitCycle.fgcmConfig.precomputeSuperStarInitialCycle = False
349 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
350 previousReservedRawRepeatability[:] = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
351 self.log.info("Setting exposure gray photometricity cuts to:")
352 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
353 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
354 continue
355 cut = fgcmFitCycle.updatedPhotometricCut[i] * 1000.0
356 self.log.info(" Band %s, photometricity cut: %.2f mmag" % (band, cut))
358 cycleNumber += 1
360 # Log warning if not converged
361 if not converged:
362 self.log.warning("Maximum number of fit cycles exceeded (%d) without convergence.", cycleNumber)
364 # Do final clean-up iteration
365 fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
366 fgcmFitCycle.fgcmConfig.resetParameters = False
367 fgcmFitCycle.fgcmConfig.maxIter = 0
368 fgcmFitCycle.fgcmConfig.outputZeropoints = True
369 fgcmFitCycle.fgcmConfig.outputStandards = True
370 fgcmFitCycle.fgcmConfig.doPlots = self.config.doDebuggingPlots
371 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
372 fgcmFitCycle.initialCycle = False
374 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
375 fgcmExpInfo,
376 previousParInfo,
377 previousParams,
378 previousSuperStar)
379 fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
380 fgcmStarObservationCat['instMagErr'][obsIndex])
381 fgcmFitCycle.setPars(fgcmPars)
382 fgcmFitCycle.finishSetup()
384 self.log.info("Running final clean-up fit cycle...")
385 fgcmFitCycle.run()
387 self.log.info("Raw repeatability after clean-up cycle is:")
388 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
389 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
390 continue
391 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
392 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
394 # Do the outputs. Need to keep track of tract.
396 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
397 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
399 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
400 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
402 atmSchema = makeAtmSchema()
403 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
405 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
406 stdSchema = makeStdSchema(len(goodBands))
407 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
409 outStruct = self.fgcmOutputProducts.generateTractOutputProducts(handleDict,
410 tract,
411 visitCat,
412 zptCat, atmCat, stdCat,
413 self.config.fgcmBuildStars)
415 outStruct.repeatability = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
417 fgcmFitCycle.freeSharedMemory()
419 return outStruct