fgcmCalibrateTractBase.py

Go to the documentation of this file.

# This file is part of fgcmcal.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""Base class for running fgcmcal on a single tract using src tables
or sourceTable_visit tables.
"""
 
import sys
import traceback
 
import numpy as np
 
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
 
from .fgcmBuildStars import FgcmBuildStarsTask, FgcmBuildStarsConfig
from .fgcmFitCycle import FgcmFitCycleConfig
from .fgcmOutputProducts import FgcmOutputProductsTask
from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog
from .utilities import computeCcdOffsets, computeApertureRadiusFromDataRef, extractReferenceMags
from .utilities import makeZptSchema, makeZptCat
from .utilities import makeAtmSchema, makeAtmCat
from .utilities import makeStdSchema, makeStdCat
 
import fgcm
 
__all__ = ['FgcmCalibrateTractConfigBase', 'FgcmCalibrateTractBaseTask', 'FgcmCalibrateTractRunner']
 
 
class FgcmCalibrateTractConfigBase(pexConfig.Config):
    """Config for FgcmCalibrateTract"""
 
    fgcmBuildStars = pexConfig.ConfigurableField(
        target=FgcmBuildStarsTask,
        doc="Task to load and match stars for fgcm",
    )
    fgcmFitCycle = pexConfig.ConfigField(
        dtype=FgcmFitCycleConfig,
        doc="Config to run a single fgcm fit cycle",
    )
    fgcmOutputProducts = pexConfig.ConfigurableField(
        target=FgcmOutputProductsTask,
        doc="Task to output fgcm products",
    )
    convergenceTolerance = pexConfig.Field(
        doc="Tolerance on repeatability convergence (per band)",
        dtype=float,
        default=0.005,
    )
    maxFitCycles = pexConfig.Field(
        doc="Maximum number of fit cycles",
        dtype=int,
        default=5,
    )
    doDebuggingPlots = pexConfig.Field(
        doc="Make plots for debugging purposes?",
        dtype=bool,
        default=False,
    )
 
    def setDefaults(self):
        pexConfig.Config.setDefaults(self)
 
        self.fgcmFitCycle.quietMode = True
        self.fgcmOutputProducts.doReferenceCalibration = False
        self.fgcmOutputProducts.doRefcatOutput = False
        self.fgcmOutputProducts.cycleNumber = 0
        self.fgcmOutputProducts.photoCal.applyColorTerms = False
 
    def validate(self):
        super().validate()
 
        for band in self.fgcmFitCycle.bands:
            if not self.fgcmFitCycle.useRepeatabilityForExpGrayCutsDict[band]:
                msg = 'Must set useRepeatabilityForExpGrayCutsDict[band]=True for all bands'
                raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict,
                                                     self, msg)
 
 
class FgcmCalibrateTractRunner(pipeBase.ButlerInitializedTaskRunner):
    """Subclass of TaskRunner for FgcmCalibrateTractTask
 
    fgcmCalibrateTractTask.run() takes a number of arguments, one of which is
    the butler (for persistence and mapper data), and a list of dataRefs
    extracted from the command line.  This task runs on a constrained set
    of dataRefs, typically a single tract.
    This class transforms the process arguments generated by the ArgumentParser
    into the arguments expected by FgcmCalibrateTractTask.run().
    This runner does not use any parallelization.
    """
 
    @staticmethod
    def getTargetList(parsedCmd):
        """
        Return a list with one element: a tuple with the butler and
        list of dataRefs.
        """
        return [(parsedCmd.butler, parsedCmd.id.refList)]
 
    def __call__(self, args):
        """
        Parameters
        ----------
        args: `tuple` with (butler, dataRefList)
 
        Returns
        -------
        exitStatus: `list` with `lsst.pipe.base.Struct`
           exitStatus (0: success; 1: failure)
           May also contain results if `self.doReturnResults` is `True`.
        """
        butler, dataRefList = args
 
        task = self.TaskClass(config=self.config, log=self.log)
 
        exitStatus = 0
        if self.doRaise:
            results = task.runDataRef(butler, dataRefList)
        else:
            try:
                results = task.runDataRef(butler, dataRefList)
            except Exception as e:
                exitStatus = 1
                task.log.fatal("Failed: %s" % e)
                if not isinstance(e, pipeBase.TaskError):
                    traceback.print_exc(file=sys.stderr)
 
        task.writeMetadata(butler)
 
        if self.doReturnResults:
            return [pipeBase.Struct(exitStatus=exitStatus,
                                    results=results)]
        else:
            return [pipeBase.Struct(exitStatus=exitStatus)]
 
    def run(self, parsedCmd):
        """
        Run the task, with no multiprocessing
 
        Parameters
        ----------
        parsedCmd: `lsst.pipe.base.ArgumentParser` parsed command line
        """
 
        resultList = []
 
        if self.precall(parsedCmd):
            targetList = self.getTargetList(parsedCmd)
            resultList = self(targetList[0])
 
        return resultList
 
 
class FgcmCalibrateTractBaseTask(pipeBase.CmdLineTask):
    """
    Base class to calibrate a single tract using fgcmcal
    """
    def __init__(self, butler=None, **kwargs):
        """
        Instantiate an `FgcmCalibrateTractTask`.
 
        Parameters
        ----------
        butler : `lsst.daf.persistence.Butler`
        """
 
        pipeBase.CmdLineTask.__init__(self, **kwargs)
 
    # no saving of metadata for now
    def _getMetadataName(self):
        return None
 
    @pipeBase.timeMethod
    def runDataRef(self, butler, dataRefs):
        """
        Run full FGCM calibration on a single tract, including building star list,
        fitting multiple cycles, and making outputs.
 
        Parameters
        ----------
        butler:  `lsst.daf.persistence.Butler`
        dataRefs: `list` of `lsst.daf.persistence.ButlerDataRef`
           Data references for the input visits.
           These may be either per-ccd "src" or per-visit"sourceTable_visit"
           references.
 
        Raises
        ------
        RuntimeError: Raised if `config.fgcmBuildStars.doReferenceMatches` is
           not True, or if fgcmLookUpTable is not available, or if
           doSubtractLocalBackground is True and aperture radius cannot be
           determined.
        """
        datasetType = dataRefs[0].butlerSubset.datasetType
        self.log.info("Running with %d %s dataRefs" % (len(dataRefs), datasetType))
 
        if not butler.datasetExists('fgcmLookUpTable'):
            raise RuntimeError("Must run FgcmCalibrateTract with an fgcmLookUpTable")
 
        if not self.config.fgcmBuildStars.doReferenceMatches:
            raise RuntimeError("Must run FgcmCalibrateTract with fgcmBuildStars.doReferenceMatches")
        if isinstance(self.config.fgcmBuildStars, FgcmBuildStarsConfig):
            if self.config.fgcmBuildStars.checkAllCcds:
                raise RuntimeError("Cannot run FgcmCalibrateTract with "
                                   "fgcmBuildStars.checkAllCcds set to True")
 
        self.makeSubtask("fgcmBuildStars", butler=butler)
        self.makeSubtask("fgcmOutputProducts", butler=butler)
 
        # Compute the aperture radius if necessary.  This is useful to do now before
        # any heavy lifting has happened (fail early).
        calibFluxApertureRadius = None
        if self.config.fgcmBuildStars.doSubtractLocalBackground:
            try:
                field = self.config.fgcmBuildStars.instFluxField
                calibFluxApertureRadius = computeApertureRadiusFromDataRef(dataRefs[0],
                                                                           field)
            except RuntimeError:
                raise RuntimeError("Could not determine aperture radius from %s. "
                                   "Cannot use doSubtractLocalBackground." %
                                   (field))
 
        # Run the build stars tasks
        tract = int(dataRefs[0].dataId['tract'])
        self.log.info("Running on tract %d" % (tract))
 
        # Note that we will need visitCat at the end of the procedure for the outputs
        groupedDataRefs = self.fgcmBuildStars.findAndGroupDataRefs(butler, dataRefs)
        camera = butler.get('camera')
        visitCat = self.fgcmBuildStars.fgcmMakeVisitCatalog(camera, groupedDataRefs)
        rad = calibFluxApertureRadius
        fgcmStarObservationCat = self.fgcmBuildStars.fgcmMakeAllStarObservations(groupedDataRefs,
                                                                                 visitCat,
                                                                                 calibFluxApertureRadius=rad)
 
        fgcmStarIdCat, fgcmStarIndicesCat, fgcmRefCat = \
            self.fgcmBuildStars.fgcmMatchStars(butler,
                                               visitCat,
                                               fgcmStarObservationCat)
 
        # Load the LUT
        lutCat = butler.get('fgcmLookUpTable')
        fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat,
                                                         dict(self.config.fgcmFitCycle.filterMap))
        del lutCat
 
        # Translate the visit catalog into fgcm format
        fgcmExpInfo = translateVisitCatalog(visitCat)
 
        camera = butler.get('camera')
        configDict = makeConfigDict(self.config.fgcmFitCycle, self.log, camera,
                                    self.config.fgcmFitCycle.maxIterBeforeFinalCycle,
                                    True, False, tract=tract)
        # Turn off plotting in tract mode
        configDict['doPlots'] = False
 
        # Use the first orientation.
        # TODO: DM-21215 will generalize to arbitrary camera orientations
        ccdOffsets = computeCcdOffsets(camera, fgcmExpInfo['TELROT'][0])
        del camera
 
        # Set up the fit cycle task
 
        noFitsDict = {'lutIndex': lutIndexVals,
                      'lutStd': lutStd,
                      'expInfo': fgcmExpInfo,
                      'ccdOffsets': ccdOffsets}
 
        fgcmFitCycle = fgcm.FgcmFitCycle(configDict, useFits=False,
                                         noFitsDict=noFitsDict, noOutput=True)
 
        # We determine the conversion from the native units (typically radians) to
        # degrees for the first star.  This allows us to treat coord_ra/coord_dec as
        # numpy arrays rather than Angles, which would we approximately 600x slower.
        conv = fgcmStarObservationCat[0]['ra'].asDegrees() / float(fgcmStarObservationCat[0]['ra'])
 
        # To load the stars, we need an initial parameter object
        fgcmPars = fgcm.FgcmParameters.newParsWithArrays(fgcmFitCycle.fgcmConfig,
                                                         fgcmLut,
                                                         fgcmExpInfo)
 
        # Match star observations to visits
        # Only those star observations that match visits from fgcmExpInfo['VISIT'] will
        # actually be transferred into fgcm using the indexing below.
 
        obsIndex = fgcmStarIndicesCat['obsIndex']
        visitIndex = np.searchsorted(fgcmExpInfo['VISIT'],
                                     fgcmStarObservationCat['visit'][obsIndex])
 
        refMag, refMagErr = extractReferenceMags(fgcmRefCat,
                                                 self.config.fgcmFitCycle.bands,
                                                 self.config.fgcmFitCycle.filterMap)
        refId = fgcmRefCat['fgcm_id'][:]
 
        fgcmStars = fgcm.FgcmStars(fgcmFitCycle.fgcmConfig)
        fgcmStars.loadStars(fgcmPars,
                            fgcmStarObservationCat['visit'][obsIndex],
                            fgcmStarObservationCat['ccd'][obsIndex],
                            fgcmStarObservationCat['ra'][obsIndex] * conv,
                            fgcmStarObservationCat['dec'][obsIndex] * conv,
                            fgcmStarObservationCat['instMag'][obsIndex],
                            fgcmStarObservationCat['instMagErr'][obsIndex],
                            fgcmExpInfo['FILTERNAME'][visitIndex],
                            fgcmStarIdCat['fgcm_id'][:],
                            fgcmStarIdCat['ra'][:],
                            fgcmStarIdCat['dec'][:],
                            fgcmStarIdCat['obsArrIndex'][:],
                            fgcmStarIdCat['nObs'][:],
                            obsX=fgcmStarObservationCat['x'][obsIndex],
                            obsY=fgcmStarObservationCat['y'][obsIndex],
                            obsDeltaMagBkg=fgcmStarObservationCat['deltaMagBkg'][obsIndex],
                            psfCandidate=fgcmStarObservationCat['psf_candidate'][obsIndex],
                            refID=refId,
                            refMag=refMag,
                            refMagErr=refMagErr,
                            flagID=None,
                            flagFlag=None,
                            computeNobs=True)
 
        # Clear out some memory
        del fgcmStarIdCat
        del fgcmStarIndicesCat
        del fgcmRefCat
 
        fgcmFitCycle.setLUT(fgcmLut)
        fgcmFitCycle.setStars(fgcmStars, fgcmPars)
 
        converged = False
        cycleNumber = 0
 
        previousReservedRawRepeatability = np.zeros(fgcmPars.nBands) + 1000.0
        previousParInfo = None
        previousParams = None
        previousSuperStar = None
 
        while (not converged and cycleNumber < self.config.maxFitCycles):
 
            fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
 
            if cycleNumber > 0:
                # Use parameters from previous cycle
                fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
                                                                  fgcmExpInfo,
                                                                  previousParInfo,
                                                                  previousParams,
                                                                  previousSuperStar)
                # We need to reset the star magnitudes and errors for the next
                # cycle
                fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
                                                            fgcmStarObservationCat['instMagErr'][obsIndex])
                fgcmFitCycle.initialCycle = False
 
            fgcmFitCycle.setPars(fgcmPars)
            fgcmFitCycle.finishSetup()
 
            fgcmFitCycle.run()
 
            # Grab the parameters for the next cycle
            previousParInfo, previousParams = fgcmFitCycle.fgcmPars.parsToArrays()
            previousSuperStar = fgcmFitCycle.fgcmPars.parSuperStarFlat.copy()
 
            self.log.info("Raw repeatability after cycle number %d is:" % (cycleNumber))
            for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
                if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
                    continue
                rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
                self.log.info("  Band %s, repeatability: %.2f mmag" % (band, rep))
 
            # Check for convergence
            if np.alltrue((previousReservedRawRepeatability -
                           fgcmFitCycle.fgcmPars.compReservedRawRepeatability) <
                          self.config.convergenceTolerance):
                self.log.info("Raw repeatability has converged after cycle number %d." % (cycleNumber))
                converged = True
            else:
                fgcmFitCycle.fgcmConfig.expGrayPhotometricCut[:] = fgcmFitCycle.updatedPhotometricCut
                fgcmFitCycle.fgcmConfig.expGrayHighCut[:] = fgcmFitCycle.updatedHighCut
                fgcmFitCycle.fgcmConfig.precomputeSuperStarInitialCycle = False
                fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
                previousReservedRawRepeatability[:] = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
                self.log.info("Setting exposure gray photometricity cuts to:")
                for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
                    if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
                        continue
                    cut = fgcmFitCycle.updatedPhotometricCut[i] * 1000.0
                    self.log.info("  Band %s, photometricity cut: %.2f mmag" % (band, cut))
 
            cycleNumber += 1
 
        # Log warning if not converged
        if not converged:
            self.log.warn("Maximum number of fit cycles exceeded (%d) without convergence." % (cycleNumber))
 
        # Do final clean-up iteration
        fgcmFitCycle.fgcmConfig.freezeStdAtmosphere = False
        fgcmFitCycle.fgcmConfig.resetParameters = False
        fgcmFitCycle.fgcmConfig.maxIter = 0
        fgcmFitCycle.fgcmConfig.outputZeropoints = True
        fgcmFitCycle.fgcmConfig.outputStandards = True
        fgcmFitCycle.fgcmConfig.doPlots = self.config.doDebuggingPlots
        fgcmFitCycle.fgcmConfig.updateCycleNumber(cycleNumber)
        fgcmFitCycle.initialCycle = False
 
        fgcmPars = fgcm.FgcmParameters.loadParsWithArrays(fgcmFitCycle.fgcmConfig,
                                                          fgcmExpInfo,
                                                          previousParInfo,
                                                          previousParams,
                                                          previousSuperStar)
        fgcmFitCycle.fgcmStars.reloadStarMagnitudes(fgcmStarObservationCat['instMag'][obsIndex],
                                                    fgcmStarObservationCat['instMagErr'][obsIndex])
        fgcmFitCycle.setPars(fgcmPars)
        fgcmFitCycle.finishSetup()
 
        self.log.info("Running final clean-up fit cycle...")
        fgcmFitCycle.run()
 
        self.log.info("Raw repeatability after clean-up cycle is:")
        for i, band in enumerate(fgcmFitCycle.fgcmPars.bands):
            if not fgcmFitCycle.fgcmPars.hasExposuresInBand[i]:
                continue
            rep = fgcmFitCycle.fgcmPars.compReservedRawRepeatability[i] * 1000.0
            self.log.info("  Band %s, repeatability: %.2f mmag" % (band, rep))
 
        # Do the outputs.  Need to keep track of tract.
 
        superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1]
        zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1]
 
        zptSchema = makeZptSchema(superStarChebSize, zptChebSize)
        zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct)
 
        atmSchema = makeAtmSchema()
        atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct)
 
        stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars)
        stdSchema = makeStdSchema(len(goodBands))
        stdCat = makeStdCat(stdSchema, stdStruct, goodBands)
 
        outStruct = self.fgcmOutputProducts.generateTractOutputProducts(butler, tract,
                                                                        visitCat,
                                                                        zptCat, atmCat, stdCat,
                                                                        self.config.fgcmBuildStars,
                                                                        self.config.fgcmFitCycle)
        outStruct.repeatability = fgcmFitCycle.fgcmPars.compReservedRawRepeatability
 
        return outStruct