jointcal.py

Go to the documentation of this file.

# This file is part of jointcal.
#
# 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/>.

import collections
import numpy as np
import astropy.units as u

import lsst.geom
import lsst.utils
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
from lsst.afw.image import fluxErrFromABMagErr
import lsst.pex.exceptions as pexExceptions
import lsst.afw.table
import lsst.log
import lsst.meas.algorithms
from lsst.pipe.tasks.colorterms import ColortermLibrary
from lsst.verify import Job, Measurement

from lsst.meas.algorithms import LoadIndexedReferenceObjectsTask, ReferenceSourceSelectorTask
from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry

from .dataIds import PerTractCcdDataIdContainer

import lsst.jointcal
from lsst.jointcal import MinimizeResult

__all__ = ["JointcalConfig", "JointcalRunner", "JointcalTask"]

Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))


# TODO: move this to MeasurementSet in lsst.verify per DM-12655.
def add_measurement(job, name, value):
    meas = Measurement(job.metrics[name], value)
    job.measurements.insert(meas)


class JointcalRunner(pipeBase.ButlerInitializedTaskRunner):
    """Subclass of TaskRunner for jointcalTask

    jointcalTask.runDataRef() takes a number of arguments, one of which is a list of dataRefs
    extracted from the command line (whereas most CmdLineTasks' runDataRef methods take
    single dataRef, are are called repeatedly). This class transforms the processed
    arguments generated by the ArgumentParser into the arguments expected by
    Jointcal.runDataRef().

    See pipeBase.TaskRunner for more information.
    """

    @staticmethod
    def getTargetList(parsedCmd, **kwargs):
        """
        Return a list of tuples per tract, each containing (dataRefs, kwargs).

        Jointcal operates on lists of dataRefs simultaneously.
        """
        kwargs['profile_jointcal'] = parsedCmd.profile_jointcal
        kwargs['butler'] = parsedCmd.butler

        # organize data IDs by tract
        refListDict = {}
        for ref in parsedCmd.id.refList:
            refListDict.setdefault(ref.dataId["tract"], []).append(ref)
        # we call runDataRef() once with each tract
        result = [(refListDict[tract], kwargs) for tract in sorted(refListDict.keys())]
        return result

    def __call__(self, args):
        """
        Parameters
        ----------
        args
            Arguments for Task.runDataRef()

        Returns
        -------
        pipe.base.Struct
            if self.doReturnResults is False:

            - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.

            if self.doReturnResults is True:

            - ``result``: the result of calling jointcal.runDataRef()
            - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.
        """
        exitStatus = 0  # exit status for shell

        # NOTE: cannot call self.makeTask because that assumes args[0] is a single dataRef.
        dataRefList, kwargs = args
        butler = kwargs.pop('butler')
        task = self.TaskClass(config=self.config, log=self.log, butler=butler)
        result = None
        try:
            result = task.runDataRef(dataRefList, **kwargs)
            exitStatus = result.exitStatus
            job_path = butler.get('verify_job_filename')
            result.job.write(job_path[0])
        except Exception as e:  # catch everything, sort it out later.
            if self.doRaise:
                raise e
            else:
                exitStatus = 1
                eName = type(e).__name__
                tract = dataRefList[0].dataId['tract']
                task.log.fatal("Failed processing tract %s, %s: %s", tract, eName, e)

        # Put the butler back into kwargs for the other Tasks.
        kwargs['butler'] = butler
        if self.doReturnResults:
            return pipeBase.Struct(result=result, exitStatus=exitStatus)
        else:
            return pipeBase.Struct(exitStatus=exitStatus)


class JointcalConfig(pexConfig.Config):
    """Configuration for JointcalTask"""

    doAstrometry = pexConfig.Field(
        doc="Fit astrometry and write the fitted result.",
        dtype=bool,
        default=True
    )
    doPhotometry = pexConfig.Field(
        doc="Fit photometry and write the fitted result.",
        dtype=bool,
        default=True
    )
    coaddName = pexConfig.Field(
        doc="Type of coadd, typically deep or goodSeeing",
        dtype=str,
        default="deep"
    )
    positionErrorPedestal = pexConfig.Field(
        doc="Systematic term to apply to the measured position error (pixels)",
        dtype=float,
        default=0.02,
    )
    photometryErrorPedestal = pexConfig.Field(
        doc="Systematic term to apply to the measured error on flux or magnitude as a "
        "fraction of source flux or magnitude delta (e.g. 0.05 is 5% of flux or +50 millimag).",
        dtype=float,
        default=0.0,
    )
    # TODO: DM-6885 matchCut should be an geom.Angle
    matchCut = pexConfig.Field(
        doc="Matching radius between fitted and reference stars (arcseconds)",
        dtype=float,
        default=3.0,
    )
    minMeasurements = pexConfig.Field(
        doc="Minimum number of associated measured stars for a fitted star to be included in the fit",
        dtype=int,
        default=2,
    )
    minMeasuredStarsPerCcd = pexConfig.Field(
        doc="Minimum number of measuredStars per ccdImage before printing warnings",
        dtype=int,
        default=100,
    )
    minRefStarsPerCcd = pexConfig.Field(
        doc="Minimum number of measuredStars per ccdImage before printing warnings",
        dtype=int,
        default=30,
    )
    allowLineSearch = pexConfig.Field(
        doc="Allow a line search during minimization, if it is reasonable for the model"
        " (models with a significant non-linear component, e.g. constrainedPhotometry).",
        dtype=bool,
        default=False
    )
    astrometrySimpleOrder = pexConfig.Field(
        doc="Polynomial order for fitting the simple astrometry model.",
        dtype=int,
        default=3,
    )
    astrometryChipOrder = pexConfig.Field(
        doc="Order of the per-chip transform for the constrained astrometry model.",
        dtype=int,
        default=1,
    )
    astrometryVisitOrder = pexConfig.Field(
        doc="Order of the per-visit transform for the constrained astrometry model.",
        dtype=int,
        default=5,
    )
    useInputWcs = pexConfig.Field(
        doc="Use the input calexp WCSs to initialize a SimpleAstrometryModel.",
        dtype=bool,
        default=True,
    )
    astrometryModel = pexConfig.ChoiceField(
        doc="Type of model to fit to astrometry",
        dtype=str,
        default="constrained",
        allowed={"simple": "One polynomial per ccd",
                 "constrained": "One polynomial per ccd, and one polynomial per visit"}
    )
    photometryModel = pexConfig.ChoiceField(
        doc="Type of model to fit to photometry",
        dtype=str,
        default="constrainedMagnitude",
        allowed={"simpleFlux": "One constant zeropoint per ccd and visit, fitting in flux space.",
                 "constrainedFlux": "Constrained zeropoint per ccd, and one polynomial per visit,"
                 " fitting in flux space.",
                 "simpleMagnitude": "One constant zeropoint per ccd and visit,"
                 " fitting in magnitude space.",
                 "constrainedMagnitude": "Constrained zeropoint per ccd, and one polynomial per visit,"
                 " fitting in magnitude space.",
                 }
    )
    applyColorTerms = pexConfig.Field(
        doc="Apply photometric color terms to reference stars?"
            "Requires that colorterms be set to a ColortermLibrary",
        dtype=bool,
        default=False
    )
    colorterms = pexConfig.ConfigField(
        doc="Library of photometric reference catalog name to color term dict.",
        dtype=ColortermLibrary,
    )
    photometryVisitOrder = pexConfig.Field(
        doc="Order of the per-visit polynomial transform for the constrained photometry model.",
        dtype=int,
        default=7,
    )
    photometryDoRankUpdate = pexConfig.Field(
        doc="Do the rank update step during minimization. "
        "Skipping this can help deal with models that are too non-linear.",
        dtype=bool,
        default=True,
    )
    astrometryDoRankUpdate = pexConfig.Field(
        doc="Do the rank update step during minimization (should not change the astrometry fit). "
        "Skipping this can help deal with models that are too non-linear.",
        dtype=bool,
        default=True,
    )
    outlierRejectSigma = pexConfig.Field(
        doc="How many sigma to reject outliers at during minimization.",
        dtype=float,
        default=5.0,
    )
    maxPhotometrySteps = pexConfig.Field(
        doc="Maximum number of minimize iterations to take when fitting photometry.",
        dtype=int,
        default=20,
    )
    maxAstrometrySteps = pexConfig.Field(
        doc="Maximum number of minimize iterations to take when fitting photometry.",
        dtype=int,
        default=20,
    )
    astrometryRefObjLoader = pexConfig.ConfigurableField(
        target=LoadIndexedReferenceObjectsTask,
        doc="Reference object loader for astrometric fit",
    )
    photometryRefObjLoader = pexConfig.ConfigurableField(
        target=LoadIndexedReferenceObjectsTask,
        doc="Reference object loader for photometric fit",
    )
    sourceSelector = sourceSelectorRegistry.makeField(
        doc="How to select sources for cross-matching",
        default="astrometry"
    )
    astrometryReferenceSelector = pexConfig.ConfigurableField(
        target=ReferenceSourceSelectorTask,
        doc="How to down-select the loaded astrometry reference catalog.",
    )
    photometryReferenceSelector = pexConfig.ConfigurableField(
        target=ReferenceSourceSelectorTask,
        doc="How to down-select the loaded photometry reference catalog.",
    )
    astrometryReferenceErr = pexConfig.Field(
        doc="Uncertainty on reference catalog coordinates [mas] to use in place of the `coord_*Err` fields."
            " If None, then raise an exception if the reference catalog is missing coordinate errors."
            " If specified, overrides any existing `coord_*Err` values.",
        dtype=float,
        default=None,
        optional=True
    )
    writeInitMatrix = pexConfig.Field(
        dtype=bool,
        doc="Write the pre/post-initialization Hessian and gradient to text files, for debugging."
            "The output files will be of the form 'astrometry_preinit-mat.txt', in the current directory."
            "Note that these files are the dense versions of the matrix, and so may be very large.",
        default=False
    )
    writeChi2FilesInitialFinal = pexConfig.Field(
        dtype=bool,
        doc="Write .csv files containing the contributions to chi2 for the initialization and final fit.",
        default=False
    )
    writeChi2FilesOuterLoop = pexConfig.Field(
        dtype=bool,
        doc="Write .csv files containing the contributions to chi2 for the outer fit loop.",
        default=False
    )
    sourceFluxType = pexConfig.Field(
        dtype=str,
        doc="Source flux field to use in source selection and to get fluxes from the catalog.",
        default='Calib'
    )

    def validate(self):
        super().validate()
        if self.doPhotometry and self.applyColorTerms and len(self.colorterms.data) == 0:
            msg = "applyColorTerms=True requires the `colorterms` field be set to a ColortermLibrary."
            raise pexConfig.FieldValidationError(JointcalConfig.colorterms, self, msg)
        if self.doAstrometry and not self.doPhotometry and self.applyColorTerms:
            msg = ("Only doing astrometry, but Colorterms are not applied for astrometry;"
                   "applyColorTerms=True will be ignored.")
            lsst.log.warn(msg)

    def setDefaults(self):
        # Use science source selector which can filter on extendedness, SNR, and whether blended
        self.sourceSelector.name = 'science'
        # Use only stars because aperture fluxes of galaxies are biased and depend on seeing
        self.sourceSelector['science'].doUnresolved = True
        # with dependable signal to noise ratio.
        self.sourceSelector['science'].doSignalToNoise = True
        # Min SNR must be > 0 because jointcal cannot handle negative fluxes,
        # and S/N > 10 to use sources that are not too faint, and thus better measured.
        self.sourceSelector['science'].signalToNoise.minimum = 10.
        # Base SNR on CalibFlux because that is the flux jointcal that fits and must be positive
        fluxField = f"slot_{self.sourceFluxType}Flux_instFlux"
        self.sourceSelector['science'].signalToNoise.fluxField = fluxField
        self.sourceSelector['science'].signalToNoise.errField = fluxField + "Err"
        # Do not trust blended sources' aperture fluxes which also depend on seeing
        self.sourceSelector['science'].doIsolated = True
        # Do not trust either flux or centroid measurements with flags,
        # chosen from the usual QA flags for stars)
        self.sourceSelector['science'].doFlags = True
        badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_saturated',
                    'base_PixelFlags_flag_interpolatedCenter', 'base_SdssCentroid_flag',
                    'base_PsfFlux_flag', 'base_PixelFlags_flag_suspectCenter']
        self.sourceSelector['science'].flags.bad = badFlags


class JointcalTask(pipeBase.CmdLineTask):
    """Jointly astrometrically and photometrically calibrate a group of images."""

    ConfigClass = JointcalConfig
    RunnerClass = JointcalRunner
    _DefaultName = "jointcal"

    def __init__(self, butler=None, profile_jointcal=False, **kwargs):
        """
        Instantiate a JointcalTask.

        Parameters
        ----------
        butler : `lsst.daf.persistence.Butler`
            The butler is passed to the refObjLoader constructor in case it is
            needed. Ignored if the refObjLoader argument provides a loader directly.
            Used to initialize the astrometry and photometry refObjLoaders.
        profile_jointcal : `bool`
            Set to True to profile different stages of this jointcal run.
        """
        pipeBase.CmdLineTask.__init__(self, **kwargs)
        self.profile_jointcal = profile_jointcal
        self.makeSubtask("sourceSelector")
        if self.config.doAstrometry:
            self.makeSubtask('astrometryRefObjLoader', butler=butler)
            self.makeSubtask("astrometryReferenceSelector")
        else:
            self.astrometryRefObjLoader = None
        if self.config.doPhotometry:
            self.makeSubtask('photometryRefObjLoader', butler=butler)
            self.makeSubtask("photometryReferenceSelector")
        else:
            self.photometryRefObjLoader = None

        # To hold various computed metrics for use by tests
        self.job = Job.load_metrics_package(subset='jointcal')

    # We don't currently need to persist the metadata.
    # If we do in the future, we will have to add appropriate dataset templates
    # to each obs package (the metadata template should look like `jointcal_wcs`).
    def _getMetadataName(self):
        return None

    @classmethod
    def _makeArgumentParser(cls):
        """Create an argument parser"""
        parser = pipeBase.ArgumentParser(name=cls._DefaultName)
        parser.add_argument("--profile_jointcal", default=False, action="store_true",
                            help="Profile steps of jointcal separately.")
        parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=6789 ccd=0..9",
                               ContainerClass=PerTractCcdDataIdContainer)
        return parser

    def _build_ccdImage(self, dataRef, associations, jointcalControl):
        """
        Extract the necessary things from this dataRef to add a new ccdImage.

        Parameters
        ----------
        dataRef : `lsst.daf.persistence.ButlerDataRef`
            DataRef to extract info from.
        associations : `lsst.jointcal.Associations`
            Object to add the info to, to construct a new CcdImage
        jointcalControl : `jointcal.JointcalControl`
            Control object for associations management

        Returns
        ------
        namedtuple
            ``wcs``
                The TAN WCS of this image, read from the calexp
                (`lsst.afw.geom.SkyWcs`).
            ``key``
                A key to identify this dataRef by its visit and ccd ids
                (`namedtuple`).
            ``filter``
                This calexp's filter (`str`).
        """
        if "visit" in dataRef.dataId.keys():
            visit = dataRef.dataId["visit"]
        else:
            visit = dataRef.getButler().queryMetadata("calexp", ("visit"), dataRef.dataId)[0]

        src = dataRef.get("src", flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS, immediate=True)

        visitInfo = dataRef.get('calexp_visitInfo')
        detector = dataRef.get('calexp_detector')
        ccdId = detector.getId()
        photoCalib = dataRef.get('calexp_photoCalib')
        tanWcs = dataRef.get('calexp_wcs')
        bbox = dataRef.get('calexp_bbox')
        filt = dataRef.get('calexp_filter')
        filterName = filt.getName()

        goodSrc = self.sourceSelector.run(src)

        if len(goodSrc.sourceCat) == 0:
            self.log.warn("No sources selected in visit %s ccd %s", visit, ccdId)
        else:
            self.log.info("%d sources selected in visit %d ccd %d", len(goodSrc.sourceCat), visit, ccdId)
        associations.createCcdImage(goodSrc.sourceCat,
                                    tanWcs,
                                    visitInfo,
                                    bbox,
                                    filterName,
                                    photoCalib,
                                    detector,
                                    visit,
                                    ccdId,
                                    jointcalControl)

        Result = collections.namedtuple('Result_from_build_CcdImage', ('wcs', 'key', 'filter'))
        Key = collections.namedtuple('Key', ('visit', 'ccd'))
        return Result(tanWcs, Key(visit, ccdId), filterName)

    @pipeBase.timeMethod
    def runDataRef(self, dataRefs, profile_jointcal=False):
        """
        Jointly calibrate the astrometry and photometry across a set of images.

        Parameters
        ----------
        dataRefs : `list` of `lsst.daf.persistence.ButlerDataRef`
            List of data references to the exposures to be fit.
        profile_jointcal : `bool`
            Profile the individual steps of jointcal.

        Returns
        -------
        result : `lsst.pipe.base.Struct`
            Struct of metadata from the fit, containing:

            ``dataRefs``
                The provided data references that were fit (with updated WCSs)
            ``oldWcsList``
                The original WCS from each dataRef
            ``metrics``
                Dictionary of internally-computed metrics for testing/validation.
        """
        if len(dataRefs) == 0:
            raise ValueError('Need a non-empty list of data references!')

        exitStatus = 0  # exit status for shell

        sourceFluxField = "slot_%sFlux" % (self.config.sourceFluxType,)
        jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
        associations = lsst.jointcal.Associations()

        visit_ccd_to_dataRef = {}
        oldWcsList = []
        filters = []
        load_cat_prof_file = 'jointcal_build_ccdImage.prof' if profile_jointcal else ''
        with pipeBase.cmdLineTask.profile(load_cat_prof_file):
            # We need the bounding-box of the focal plane for photometry visit models.
            # NOTE: we only need to read it once, because its the same for all exposures of a camera.
            camera = dataRefs[0].get('camera', immediate=True)
            self.focalPlaneBBox = camera.getFpBBox()
            for ref in dataRefs:
                result = self._build_ccdImage(ref, associations, jointcalControl)
                oldWcsList.append(result.wcs)
                visit_ccd_to_dataRef[result.key] = ref
                filters.append(result.filter)
        filters = collections.Counter(filters)

        associations.computeCommonTangentPoint()

        boundingCircle = associations.computeBoundingCircle()
        center = lsst.geom.SpherePoint(boundingCircle.getCenter())
        radius = lsst.geom.Angle(boundingCircle.getOpeningAngle().asRadians(), lsst.geom.radians)

        # Determine a default filter associated with the catalog. See DM-9093
        defaultFilter = filters.most_common(1)[0][0]
        self.log.debug("Using %s band for reference flux", defaultFilter)

        # TODO: need a better way to get the tract.
        tract = dataRefs[0].dataId['tract']

        if self.config.doAstrometry:
            astrometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
                                                      name="astrometry",
                                                      refObjLoader=self.astrometryRefObjLoader,
                                                      referenceSelector=self.astrometryReferenceSelector,
                                                      fit_function=self._fit_astrometry,
                                                      profile_jointcal=profile_jointcal,
                                                      tract=tract)
            self._write_astrometry_results(associations, astrometry.model, visit_ccd_to_dataRef)
        else:
            astrometry = Astrometry(None, None, None)

        if self.config.doPhotometry:
            photometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
                                                      name="photometry",
                                                      refObjLoader=self.photometryRefObjLoader,
                                                      referenceSelector=self.photometryReferenceSelector,
                                                      fit_function=self._fit_photometry,
                                                      profile_jointcal=profile_jointcal,
                                                      tract=tract,
                                                      filters=filters,
                                                      reject_bad_fluxes=True)
            self._write_photometry_results(associations, photometry.model, visit_ccd_to_dataRef)
        else:
            photometry = Photometry(None, None)

        return pipeBase.Struct(dataRefs=dataRefs,
                               oldWcsList=oldWcsList,
                               job=self.job,
                               astrometryRefObjLoader=self.astrometryRefObjLoader,
                               photometryRefObjLoader=self.photometryRefObjLoader,
                               defaultFilter=defaultFilter,
                               exitStatus=exitStatus)

    def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
                                filters=[],
                                tract="", profile_jointcal=False, match_cut=3.0,
                                reject_bad_fluxes=False, *,
                                name="", refObjLoader=None, referenceSelector=None,
                                fit_function=None):
        """Load reference catalog, perform the fit, and return the result.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        defaultFilter : `str`
            filter to load from reference catalog.
        center : `lsst.geom.SpherePoint`
            ICRS center of field to load from reference catalog.
        radius : `lsst.geom.Angle`
            On-sky radius to load from reference catalog.
        name : `str`
            Name of thing being fit: "astrometry" or "photometry".
        refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
            Reference object loader to use to load a reference catalog.
        referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
            Selector to use to pick objects from the loaded reference catalog.
        fit_function : callable
            Function to call to perform fit (takes Associations object).
        filters : `list` [`str`], optional
            List of filters to load from the reference catalog.
        tract : `str`, optional
            Name of tract currently being fit.
        profile_jointcal : `bool`, optional
            Separately profile the fitting step.
        match_cut : `float`, optional
            Radius in arcseconds to find cross-catalog matches to during
            associations.associateCatalogs.
        reject_bad_fluxes : `bool`, optional
            Reject refCat sources with NaN/inf flux or NaN/0 fluxErr.

        Returns
        -------
        result : `Photometry` or `Astrometry`
            Result of `fit_function()`
        """
        self.log.info("====== Now processing %s...", name)
        # TODO: this should not print "trying to invert a singular transformation:"
        # if it does that, something's not right about the WCS...
        associations.associateCatalogs(match_cut)
        add_measurement(self.job, 'jointcal.associated_%s_fittedStars' % name,
                        associations.fittedStarListSize())

        applyColorterms = False if name.lower() == "astrometry" else self.config.applyColorTerms
        refCat, fluxField = self._load_reference_catalog(refObjLoader, referenceSelector,
                                                         center, radius, defaultFilter,
                                                         applyColorterms=applyColorterms)

        if self.config.astrometryReferenceErr is None:
            refCoordErr = float('nan')
        else:
            refCoordErr = self.config.astrometryReferenceErr

        associations.collectRefStars(refCat,
                                     self.config.matchCut*lsst.geom.arcseconds,
                                     fluxField,
                                     refCoordinateErr=refCoordErr,
                                     rejectBadFluxes=reject_bad_fluxes)
        add_measurement(self.job, 'jointcal.collected_%s_refStars' % name,
                        associations.refStarListSize())

        associations.prepareFittedStars(self.config.minMeasurements)

        self._check_star_lists(associations, name)
        add_measurement(self.job, 'jointcal.selected_%s_refStars' % name,
                        associations.nFittedStarsWithAssociatedRefStar())
        add_measurement(self.job, 'jointcal.selected_%s_fittedStars' % name,
                        associations.fittedStarListSize())
        add_measurement(self.job, 'jointcal.selected_%s_ccdImages' % name,
                        associations.nCcdImagesValidForFit())

        load_cat_prof_file = 'jointcal_fit_%s.prof'%name if profile_jointcal else ''
        dataName = "{}_{}".format(tract, defaultFilter)
        with pipeBase.cmdLineTask.profile(load_cat_prof_file):
            result = fit_function(associations, dataName)
        # TODO DM-12446: turn this into a "butler save" somehow.
        # Save reference and measurement chi2 contributions for this data
        if self.config.writeChi2FilesInitialFinal:
            baseName = f"{name}_final_chi2-{dataName}"
            result.fit.saveChi2Contributions(baseName+"{type}")

        return result

    def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterName,
                                applyColorterms=False):
        """Load the necessary reference catalog sources, convert fluxes to
        correct units, and apply color term corrections if requested.

        Parameters
        ----------
        refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
            The reference catalog loader to use to get the data.
        referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
            Source selector to apply to loaded reference catalog.
        center : `lsst.geom.SpherePoint`
            The center around which to load sources.
        radius : `lsst.geom.Angle`
            The radius around ``center`` to load sources in.
        filterName : `str`
            The name of the camera filter to load fluxes for.
        applyColorterms : `bool`
            Apply colorterm corrections to the refcat for ``filterName``?

        Returns
        -------
        refCat : `lsst.afw.table.SimpleCatalog`
            The loaded reference catalog.
        fluxField : `str`
            The name of the reference catalog flux field appropriate for ``filterName``.
        """
        skyCircle = refObjLoader.loadSkyCircle(center,
                                               radius,
                                               filterName)

        selected = referenceSelector.run(skyCircle.refCat)
        # Need memory contiguity to get reference filters as a vector.
        if not selected.sourceCat.isContiguous():
            refCat = selected.sourceCat.copy(deep=True)
        else:
            refCat = selected.sourceCat

        if self.config.astrometryReferenceErr is None and 'coord_raErr' not in refCat.schema:
            msg = ("Reference catalog does not contain coordinate errors, "
                   "and config.astrometryReferenceErr not supplied.")
            raise pexConfig.FieldValidationError(JointcalConfig.astrometryReferenceErr,
                                                 self.config,
                                                 msg)

        if self.config.astrometryReferenceErr is not None and 'coord_raErr' in refCat.schema:
            self.log.warn("Overriding reference catalog coordinate errors with %f/coordinate [mas]",
                          self.config.astrometryReferenceErr)

        if applyColorterms:
            try:
                refCatName = refObjLoader.ref_dataset_name
            except AttributeError:
                # NOTE: we need this try:except: block in place until we've completely removed a.net support.
                raise RuntimeError("Cannot perform colorterm corrections with a.net refcats.")
            self.log.info("Applying color terms for filterName=%r reference catalog=%s",
                          filterName, refCatName)
            colorterm = self.config.colorterms.getColorterm(
                filterName=filterName, photoCatName=refCatName, doRaise=True)

            refMag, refMagErr = colorterm.getCorrectedMagnitudes(refCat, filterName)
            refCat[skyCircle.fluxField] = u.Magnitude(refMag, u.ABmag).to_value(u.nJy)
            # TODO: I didn't want to use this, but I'll deal with it in DM-16903
            refCat[skyCircle.fluxField+'Err'] = fluxErrFromABMagErr(refMagErr, refMag) * 1e9

        return refCat, skyCircle.fluxField

    def _check_star_lists(self, associations, name):
        # TODO: these should be len(blah), but we need this properly wrapped first.
        if associations.nCcdImagesValidForFit() == 0:
            raise RuntimeError('No images in the ccdImageList!')
        if associations.fittedStarListSize() == 0:
            raise RuntimeError('No stars in the {} fittedStarList!'.format(name))
        if associations.refStarListSize() == 0:
            raise RuntimeError('No stars in the {} reference star list!'.format(name))

    def _logChi2AndValidate(self, associations, fit, model, chi2Label="Model",
                            writeChi2Name=None):
        """Compute chi2, log it, validate the model, and return chi2.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        fit : `lsst.jointcal.FitterBase`
            The fitter to use for minimization.
        model : `lsst.jointcal.Model`
            The model being fit.
        chi2Label : str, optional
            Label to describe the chi2 (e.g. "Initialized", "Final").
        writeChi2Name : `str`, optional
            Filename prefix to write the chi2 contributions to.
            Do not supply an extension: an appropriate one will be added.

        Returns
        -------
        chi2: `lsst.jointcal.Chi2Accumulator`
            The chi2 object for the current fitter and model.

        Raises
        ------
        FloatingPointError
            Raised if chi2 is infinite or NaN.
        ValueError
            Raised if the model is not valid.
        """
        if writeChi2Name is not None:
            fit.saveChi2Contributions(writeChi2Name+"{type}")
            self.log.info("Wrote chi2 contributions files: %s", writeChi2Name)

        chi2 = fit.computeChi2()
        self.log.info("%s %s", chi2Label, chi2)
        self._check_stars(associations)
        if not np.isfinite(chi2.chi2):
            raise FloatingPointError(f'{chi2Label} chi2 is invalid: {chi2}')
        if not model.validate(associations.getCcdImageList(), chi2.ndof):
            raise ValueError("Model is not valid: check log messages for warnings.")
        return chi2

    def _fit_photometry(self, associations, dataName=None):
        """
        Fit the photometric data.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        dataName : `str`
            Name of the data being processed (e.g. "1234_HSC-Y"), for
            identifying debugging files.

        Returns
        -------
        fit_result : `namedtuple`
            fit : `lsst.jointcal.PhotometryFit`
                The photometric fitter used to perform the fit.
            model : `lsst.jointcal.PhotometryModel`
                The photometric model that was fit.
        """
        self.log.info("=== Starting photometric fitting...")

        # TODO: should use pex.config.RegistryField here (see DM-9195)
        if self.config.photometryModel == "constrainedFlux":
            model = lsst.jointcal.ConstrainedFluxModel(associations.getCcdImageList(),
                                                       self.focalPlaneBBox,
                                                       visitOrder=self.config.photometryVisitOrder,
                                                       errorPedestal=self.config.photometryErrorPedestal)
            # potentially nonlinear problem, so we may need a line search to converge.
            doLineSearch = self.config.allowLineSearch
        elif self.config.photometryModel == "constrainedMagnitude":
            model = lsst.jointcal.ConstrainedMagnitudeModel(associations.getCcdImageList(),
                                                            self.focalPlaneBBox,
                                                            visitOrder=self.config.photometryVisitOrder,
                                                            errorPedestal=self.config.photometryErrorPedestal)
            # potentially nonlinear problem, so we may need a line search to converge.
            doLineSearch = self.config.allowLineSearch
        elif self.config.photometryModel == "simpleFlux":
            model = lsst.jointcal.SimpleFluxModel(associations.getCcdImageList(),
                                                  errorPedestal=self.config.photometryErrorPedestal)
            doLineSearch = False  # purely linear in model parameters, so no line search needed
        elif self.config.photometryModel == "simpleMagnitude":
            model = lsst.jointcal.SimpleMagnitudeModel(associations.getCcdImageList(),
                                                       errorPedestal=self.config.photometryErrorPedestal)
            doLineSearch = False  # purely linear in model parameters, so no line search needed

        fit = lsst.jointcal.PhotometryFit(associations, model)
        # TODO DM-12446: turn this into a "butler save" somehow.
        # Save reference and measurement chi2 contributions for this data
        if self.config.writeChi2FilesInitialFinal:
            baseName = f"photometry_initial_chi2-{dataName}"
        else:
            baseName = None
        self._logChi2AndValidate(associations, fit, model, "Initialized", writeChi2Name=baseName)

        def getChi2Name(whatToFit):
            if self.config.writeChi2FilesOuterLoop:
                return f"photometry_init-%s_chi2-{dataName}" % whatToFit
            else:
                return None

        # The constrained model needs the visit transform fit first; the chip
        # transform is initialized from the singleFrame PhotoCalib, so it's close.
        dumpMatrixFile = "photometry_preinit" if self.config.writeInitMatrix else ""
        if self.config.photometryModel.startswith("constrained"):
            # no line search: should be purely (or nearly) linear,
            # and we want a large step size to initialize with.
            fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
            self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("ModelVisit"))
            dumpMatrixFile = ""  # so we don't redo the output on the next step

        fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
        self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Model"))

        fit.minimize("Fluxes")  # no line search: always purely linear.
        self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Fluxes"))

        fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
        self._logChi2AndValidate(associations, fit, model, "Fit prepared",
                                 writeChi2Name=getChi2Name("ModelFluxes"))

        model.freezeErrorTransform()
        self.log.debug("Photometry error scales are frozen.")

        chi2 = self._iterate_fit(associations,
                                 fit,
                                 self.config.maxPhotometrySteps,
                                 "photometry",
                                 "Model Fluxes",
                                 doRankUpdate=self.config.photometryDoRankUpdate,
                                 doLineSearch=doLineSearch,
                                 dataName=dataName)

        add_measurement(self.job, 'jointcal.photometry_final_chi2', chi2.chi2)
        add_measurement(self.job, 'jointcal.photometry_final_ndof', chi2.ndof)
        return Photometry(fit, model)

    def _fit_astrometry(self, associations, dataName=None):
        """
        Fit the astrometric data.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        dataName : `str`
            Name of the data being processed (e.g. "1234_HSC-Y"), for
            identifying debugging files.

        Returns
        -------
        fit_result : `namedtuple`
            fit : `lsst.jointcal.AstrometryFit`
                The astrometric fitter used to perform the fit.
            model : `lsst.jointcal.AstrometryModel`
                The astrometric model that was fit.
            sky_to_tan_projection : `lsst.jointcal.ProjectionHandler`
                The model for the sky to tangent plane projection that was used in the fit.
        """

        self.log.info("=== Starting astrometric fitting...")

        associations.deprojectFittedStars()

        # NOTE: need to return sky_to_tan_projection so that it doesn't get garbage collected.
        # TODO: could we package sky_to_tan_projection and model together so we don't have to manage
        # them so carefully?
        sky_to_tan_projection = lsst.jointcal.OneTPPerVisitHandler(associations.getCcdImageList())

        if self.config.astrometryModel == "constrained":
            model = lsst.jointcal.ConstrainedAstrometryModel(associations.getCcdImageList(),
                                                             sky_to_tan_projection,
                                                             chipOrder=self.config.astrometryChipOrder,
                                                             visitOrder=self.config.astrometryVisitOrder)
        elif self.config.astrometryModel == "simple":
            model = lsst.jointcal.SimpleAstrometryModel(associations.getCcdImageList(),
                                                        sky_to_tan_projection,
                                                        self.config.useInputWcs,
                                                        nNotFit=0,
                                                        order=self.config.astrometrySimpleOrder)

        fit = lsst.jointcal.AstrometryFit(associations, model, self.config.positionErrorPedestal)
        # TODO DM-12446: turn this into a "butler save" somehow.
        # Save reference and measurement chi2 contributions for this data
        if self.config.writeChi2FilesInitialFinal:
            baseName = f"astrometry_initial_chi2-{dataName}"
        else:
            baseName = None
        self._logChi2AndValidate(associations, fit, model, "Initial", writeChi2Name=baseName)

        def getChi2Name(whatToFit):
            if self.config.writeChi2FilesOuterLoop:
                return f"astrometry_init-%s_chi2-{dataName}" % whatToFit
            else:
                return None

        dumpMatrixFile = "astrometry_preinit" if self.config.writeInitMatrix else ""
        # The constrained model needs the visit transform fit first; the chip
        # transform is initialized from the detector's cameraGeom, so it's close.
        if self.config.astrometryModel == "constrained":
            fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
            self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("DistortionsVisit"))
            dumpMatrixFile = ""  # so we don't redo the output on the next step

        fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
        self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Distortions"))

        fit.minimize("Positions")
        self._logChi2AndValidate(associations, fit, model, writeChi2Name=getChi2Name("Positions"))

        fit.minimize("Distortions Positions")
        self._logChi2AndValidate(associations, fit, model, "Fit prepared",
                                 writeChi2Name=getChi2Name("DistortionsPositions"))

        chi2 = self._iterate_fit(associations,
                                 fit,
                                 self.config.maxAstrometrySteps,
                                 "astrometry",
                                 "Distortions Positions",
                                 doRankUpdate=self.config.astrometryDoRankUpdate,
                                 dataName=dataName)

        add_measurement(self.job, 'jointcal.astrometry_final_chi2', chi2.chi2)
        add_measurement(self.job, 'jointcal.astrometry_final_ndof', chi2.ndof)

        return Astrometry(fit, model, sky_to_tan_projection)

    def _check_stars(self, associations):
        """Count measured and reference stars per ccd and warn/log them."""
        for ccdImage in associations.getCcdImageList():
            nMeasuredStars, nRefStars = ccdImage.countStars()
            self.log.debug("ccdImage %s has %s measured and %s reference stars",
                           ccdImage.getName(), nMeasuredStars, nRefStars)
            if nMeasuredStars < self.config.minMeasuredStarsPerCcd:
                self.log.warn("ccdImage %s has only %s measuredStars (desired %s)",
                              ccdImage.getName(), nMeasuredStars, self.config.minMeasuredStarsPerCcd)
            if nRefStars < self.config.minRefStarsPerCcd:
                self.log.warn("ccdImage %s has only %s RefStars (desired %s)",
                              ccdImage.getName(), nRefStars, self.config.minRefStarsPerCcd)

    def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit,
                     dataName="",
                     doRankUpdate=True,
                     doLineSearch=False):
        """Run fitter.minimize up to max_steps times, returning the final chi2.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        fitter : `lsst.jointcal.FitterBase`
            The fitter to use for minimization.
        max_steps : `int`
            Maximum number of steps to run outlier rejection before declaring
            convergence failure.
        name : {'photometry' or 'astrometry'}
            What type of data are we fitting (for logs and debugging files).
        whatToFit : `str`
            Passed to ``fitter.minimize()`` to define the parameters to fit.
        dataName : `str`, optional
            Descriptive name for this dataset (e.g. tract and filter),
            for debugging.
        doRankUpdate : `bool`, optional
            Do an Eigen rank update during minimization, or recompute the full
            matrix and gradient?
        doLineSearch : `bool`, optional
            Do a line search for the optimum step during minimization?

        Returns
        -------
        chi2: `lsst.jointcal.Chi2Statistic`
            The final chi2 after the fit converges, or is forced to end.

        Raises
        ------
        FloatingPointError
            Raised if the fitter fails with a non-finite value.
        RuntimeError
            Raised if the fitter fails for some other reason;
            log messages will provide further details.
        """
        dumpMatrixFile = "%s_postinit" % name if self.config.writeInitMatrix else ""
        for i in range(max_steps):
            if self.config.writeChi2FilesOuterLoop:
                writeChi2Name = f"{name}_iterate_{i}_chi2-{dataName}"
            else:
                writeChi2Name = None
            result = fitter.minimize(whatToFit,
                                     self.config.outlierRejectSigma,
                                     doRankUpdate=doRankUpdate,
                                     doLineSearch=doLineSearch,
                                     dumpMatrixFile=dumpMatrixFile)
            dumpMatrixFile = ""  # clear it so we don't write the matrix again.
            chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(),
                                            writeChi2Name=writeChi2Name)

            if result == MinimizeResult.Converged:
                if doRankUpdate:
                    self.log.debug("fit has converged - no more outliers - redo minimization "
                                   "one more time in case we have lost accuracy in rank update.")
                    # Redo minimization one more time in case we have lost accuracy in rank update
                    result = fitter.minimize(whatToFit, self.config.outlierRejectSigma)
                    chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(), "Fit completed")

                # log a message for a large final chi2, TODO: DM-15247 for something better
                if chi2.chi2/chi2.ndof >= 4.0:
                    self.log.error("Potentially bad fit: High chi-squared/ndof.")

                break
            elif result == MinimizeResult.Chi2Increased:
                self.log.warn("still some outliers but chi2 increases - retry")
            elif result == MinimizeResult.NonFinite:
                filename = "{}_failure-nonfinite_chi2-{}.csv".format(name, dataName)
                # TODO DM-12446: turn this into a "butler save" somehow.
                fitter.saveChi2Contributions(filename)
                msg = "Nonfinite value in chi2 minimization, cannot complete fit. Dumped star tables to: {}"
                raise FloatingPointError(msg.format(filename))
            elif result == MinimizeResult.Failed:
                raise RuntimeError("Chi2 minimization failure, cannot complete fit.")
            else:
                raise RuntimeError("Unxepected return code from minimize().")
        else:
            self.log.error("%s failed to converge after %d steps"%(name, max_steps))

        return chi2

    def _write_astrometry_results(self, associations, model, visit_ccd_to_dataRef):
        """
        Write the fitted astrometric results to a new 'jointcal_wcs' dataRef.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        model : `lsst.jointcal.AstrometryModel`
            The astrometric model that was fit.
        visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
            Dict of ccdImage identifiers to dataRefs that were fit.
        """

        ccdImageList = associations.getCcdImageList()
        for ccdImage in ccdImageList:
            # TODO: there must be a better way to identify this ccdImage than a visit,ccd pair?
            ccd = ccdImage.ccdId
            visit = ccdImage.visit
            dataRef = visit_ccd_to_dataRef[(visit, ccd)]
            self.log.info("Updating WCS for visit: %d, ccd: %d", visit, ccd)
            skyWcs = model.makeSkyWcs(ccdImage)
            try:
                dataRef.put(skyWcs, 'jointcal_wcs')
            except pexExceptions.Exception as e:
                self.log.fatal('Failed to write updated Wcs: %s', str(e))
                raise e

    def _write_photometry_results(self, associations, model, visit_ccd_to_dataRef):
        """
        Write the fitted photometric results to a new 'jointcal_photoCalib' dataRef.

        Parameters
        ----------
        associations : `lsst.jointcal.Associations`
            The star/reference star associations to fit.
        model : `lsst.jointcal.PhotometryModel`
            The photoometric model that was fit.
        visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
            Dict of ccdImage identifiers to dataRefs that were fit.
        """

        ccdImageList = associations.getCcdImageList()
        for ccdImage in ccdImageList:
            # TODO: there must be a better way to identify this ccdImage than a visit,ccd pair?
            ccd = ccdImage.ccdId
            visit = ccdImage.visit
            dataRef = visit_ccd_to_dataRef[(visit, ccd)]
            self.log.info("Updating PhotoCalib for visit: %d, ccd: %d", visit, ccd)
            photoCalib = model.toPhotoCalib(ccdImage)
            try:
                dataRef.put(photoCalib, 'jointcal_photoCalib')
            except pexExceptions.Exception as e:
                self.log.fatal('Failed to write updated PhotoCalib: %s', str(e))
                raise e