Coverage for python / lsst / fgcmcal / fgcmCalibrateTractBase.py: 14%
161 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:32 +0000
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 fgcmStarObservationCat
294 del fgcmRefCat
296 fgcmFitCycle.setLUT(fgcmLut)
297 fgcmFitCycle.setStars(fgcmStars, fgcmPars)
298 fgcmFitCycle.setPars(fgcmPars)
299 fgcmFitCycle.finishSetup()
301 converged = False
302 cycleNumber = 0
304 previousReservedRawRepeatability = np.zeros(fgcmPars.nBands) + 1000.0
305 previousParInfo = None
306 previousParams = None
307 previousSuperStar = None
309 while (not converged and cycleNumber < self.config.maxFitCycles):
311 fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
313 if cycleNumber > 0:
314 # Use parameters from previous cycle
315 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
316 fgcmExpInfo,
317 previousParInfo,
318 previousParams,
319 previousSuperStar)
321 expGrayPhotometricCutDict = fgcmFitCycle.fgcmConfig.expGrayPhotometricCutDict
322 expGrayHighCutDict = fgcmFitCycle.fgcmConfig.expGrayHighCutDict
323 for i, key in enumerate(expGrayPhotometricCutDict.keys()):
324 expGrayPhotometricCutDict[key] = float(fgcmFitCycle.updatedPhotometricCut[i])
325 expGrayHighCutDict[key] = float(fgcmFitCycle.updatedHighCut[i])
327 fgcmFitCycle.updateConfigNextCycle(
328 cycleNumber,
329 resetParameters=True,
330 outputStandards=False,
331 outputZeropoints=False,
332 freezeStdAtmosphere=False,
333 expGrayPhotometricCutDict=expGrayPhotometricCutDict,
334 expGrayHighCutDict=expGrayHighCutDict,
335 )
337 fgcmFitCycle.fgcmStars.reloadStarMagnitudes()
338 fgcmFitCycle.fgcmStars.computeAllNobs(fgcmPars)
340 fgcmFitCycle.setPars(fgcmPars)
341 fgcmFitCycle.finishReset()
343 fgcmFitCycle.run()
345 # Grab the parameters for the next cycle
346 previousParInfo, previousParams = fgcmFitCycle.fgcmPars.parsToArrays()
347 previousSuperStar = fgcmFitCycle.fgcmPars.parSuperStarFlat.copy()
349 self.log.info("Raw repeatability after cycle number %d is:" % (cycleNumber))
350 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
351 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
352 continue
353 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
354 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
356 # Check for convergence
357 if np.all((previousReservedRawRepeatability
358 - fgcmFitCycle.fgcmPars.compReservedRawRepeatability)
359 < self.config.convergenceTolerance):
360 self.log.info("Raw repeatability has converged after cycle number %d." % (cycleNumber))
361 converged = True
362 else:
363 previousReservedRawRepeatability[:] = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
364 self.log.info("Setting exposure gray photometricity cuts to:")
365 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
366 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
367 continue
368 cut = fgcmFitCycle.updatedPhotometricCut[i] * 1000.0
369 self.log.info(" Band %s, photometricity cut: %.2f mmag" % (band, cut))
371 cycleNumber += 1
373 # Log warning if not converged
374 if not converged:
375 self.log.warning("Maximum number of fit cycles exceeded (%d) without convergence.", cycleNumber)
377 # Do final clean-up iteration
378 expGrayPhotometricCutDict = fgcmFitCycle.fgcmConfig.expGrayPhotometricCutDict
379 expGrayHighCutDict = fgcmFitCycle.fgcmConfig.expGrayHighCutDict
380 for i, key in enumerate(expGrayPhotometricCutDict.keys()):
381 expGrayPhotometricCutDict[key] = float(fgcmFitCycle.updatedPhotometricCut[i])
382 expGrayHighCutDict[key] = float(fgcmFitCycle.updatedHighCut[i])
384 fgcmFitCycle.updateConfigNextCycle(
385 cycleNumber,
386 maxIter=0,
387 resetParameters=False,
388 outputStandards=True,
389 outputZeropoints=True,
390 freezeStdAtmosphere=False,
391 expGrayPhotometricCutDict=expGrayPhotometricCutDict,
392 expGrayHighCutDict=expGrayHighCutDict,
393 )
395 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
396 fgcmExpInfo,
397 previousParInfo,
398 previousParams,
399 previousSuperStar)
401 fgcmFitCycle.fgcmStars.reloadStarMagnitudes()
402 fgcmFitCycle.fgcmStars.computeAllNobs(fgcmPars)
404 fgcmFitCycle.setPars(fgcmPars)
405 fgcmFitCycle.finishReset()
407 self.log.info("Running final clean-up fit cycle...")
408 fgcmFitCycle.run()
410 self.log.info("Raw repeatability after clean-up cycle is:")
411 for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
412 if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
413 continue
414 rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
415 self.log.info(" Band %s, repeatability: %.2f mmag" % (band, rep))
417 # Do the outputs. Need to keep track of tract.
419 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
420 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
422 zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
423 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
425 atmSchema = makeAtmSchema()
426 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
428 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
429 stdSchema = makeStdSchema(len(goodBands))
430 stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
432 outStruct = self.fgcmOutputProducts.generateTractOutputProducts(handleDict,
433 tract,
434 visitCat,
435 zptCat, atmCat, stdCat,
436 self.config.fgcmBuildStars)
438 outStruct.repeatability = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
440 fgcmFitCycle.freeSharedMemory()
442 return outStruct