brighterFatterKernel.py

Go to the documentation of this file.

# This file is part of ip_isr.
#
# 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/>.
#
"""Brighter Fatter Kernel calibration definition."""
 
 
__all__ = ['BrighterFatterKernel']
 
 
import numpy as np
from astropy.table import Table
import lsst.afw.math as afwMath
from . import IsrCalib
 
 
class BrighterFatterKernel(IsrCalib):
    """Calibration of brighter-fatter kernels for an instrument.
 
    ampKernels are the kernels for each amplifier in a detector, as
    generated by having level == 'AMP'
 
    detectorKernel is the kernel generated for a detector as a
    whole, as generated by having level == 'DETECTOR'
 
    makeDetectorKernelFromAmpwiseKernels is a method to generate the
    kernel for a detector, constructed by averaging together the
    ampwise kernels in the detector.  The existing application code is
    only defined for kernels with level == 'DETECTOR', so this method
    is used if the supplied kernel was built with level == 'AMP'.
 
    Parameters
    ----------
    level : `str`
        Level the kernels will be generated for.
    log : `logging.Logger`, optional
        Log to write messages to.
    **kwargs :
        Parameters to pass to parent constructor.
 
    Notes
    -----
    TODO: DM-35260
        Document what is stored in the BFK calibration.
 
    Version 1.1 adds the `expIdMask` property, and substitutes
    `means` and `variances` for `rawMeans` and `rawVariances`
    from the PTC dataset.
    """
    _OBSTYPE = 'bfk'
    _SCHEMA = 'Brighter-fatter kernel'
    _VERSION = 1.1
 
    def __init__(self, camera=None, level=None, **kwargs):
        self.levellevel = level
 
        # Things inherited from the PTC
        self.expIdMaskexpIdMask = dict()
        self.rawMeansrawMeans = dict()
        self.rawVariancesrawVariances = dict()
        self.rawXcorrsrawXcorrs = dict()
        self.badAmpsbadAmps = list()
        self.shapeshape = (17, 17)
        self.gaingain = dict()
        self.noisenoise = dict()
 
        # Things calculated from the PTC
        self.meanXcorrsmeanXcorrs = dict()
        self.validvalid = dict()
 
        # Things that are used downstream
        self.ampKernelsampKernels = dict()
        self.detKernelsdetKernels = dict()
 
        super().__init__(**kwargs)
 
        if camera:
            self.initFromCamerainitFromCamera(camera, detectorId=kwargs.get('detectorId', None))
 
        self.requiredAttributes.update(['level', 'expIdMask', 'rawMeans', 'rawVariances', 'rawXcorrs',
                                        'badAmps', 'gain', 'noise', 'meanXcorrs', 'valid',
                                        'ampKernels', 'detKernels'])
 
    def updateMetadata(self, setDate=False, **kwargs):
        """Update calibration metadata.
 
        This calls the base class's method after ensuring the required
        calibration keywords will be saved.
 
        Parameters
        ----------
        setDate : `bool`, optional
            Update the CALIBDATE fields in the metadata to the current
            time. Defaults to False.
        kwargs :
            Other keyword parameters to set in the metadata.
        """
        kwargs['LEVEL'] = self.levellevel
        kwargs['KERNEL_DX'] = self.shapeshape[0]
        kwargs['KERNEL_DY'] = self.shapeshape[1]
 
        super().updateMetadata(setDate=setDate, **kwargs)
 
    def initFromCamera(self, camera, detectorId=None):
        """Initialize kernel structure from camera.
 
        Parameters
        ----------
        camera : `lsst.afw.cameraGeom.Camera`
            Camera to use to define geometry.
        detectorId : `int`, optional
            Index of the detector to generate.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel`
            The initialized calibration.
 
        Raises
        ------
        RuntimeError :
            Raised if no detectorId is supplied for a calibration with
            level='AMP'.
        """
        self._instrument_instrument = camera.getName()
 
        if detectorId is not None:
            detector = camera[detectorId]
            self._detectorId_detectorId = detectorId
            self._detectorName_detectorName = detector.getName()
            self._detectorSerial_detectorSerial = detector.getSerial()
 
        if self.levellevel == 'AMP':
            if detectorId is None:
                raise RuntimeError("A detectorId must be supplied if level='AMP'.")
 
            self.badAmpsbadAmps = []
 
            for amp in detector:
                ampName = amp.getName()
                self.expIdMaskexpIdMask[ampName] = []
                self.rawMeansrawMeans[ampName] = []
                self.rawVariancesrawVariances[ampName] = []
                self.rawXcorrsrawXcorrs[ampName] = []
                self.gaingain[ampName] = amp.getGain()
                self.noisenoise[ampName] = amp.getReadNoise()
                self.meanXcorrsmeanXcorrs[ampName] = []
                self.ampKernelsampKernels[ampName] = []
                self.validvalid[ampName] = []
        elif self.levellevel == 'DETECTOR':
            if detectorId is None:
                for det in camera:
                    detName = det.getName()
                    self.detKernelsdetKernels[detName] = []
            else:
                self.detKernelsdetKernels[self._detectorName_detectorName] = []
 
        return self
 
    def getLengths(self):
        """Return the set of lengths needed for reshaping components.
 
        Returns
        -------
        kernelLength : `int`
            Product of the elements of self.shapeshape.
        smallLength : `int`
            Size of an untiled covariance.
        nObs : `int`
            Number of observation pairs used in the kernel.
        """
        kernelLength = self.shapeshape[0] * self.shapeshape[1]
        smallLength = int((self.shapeshape[0] - 1)*(self.shapeshape[1] - 1)/4)
        if self.levellevel == 'AMP':
            nObservations = set([len(self.rawMeansrawMeans[amp]) for amp in self.rawMeansrawMeans])
            if len(nObservations) != 1:
                raise RuntimeError("Inconsistent number of observations found.")
            nObs = nObservations.pop()
        else:
            nObs = 0
 
        return (kernelLength, smallLength, nObs)
 
    @classmethod
    def fromDict(cls, dictionary):
        """Construct a calibration from a dictionary of properties.
 
        Parameters
        ----------
        dictionary : `dict`
            Dictionary of properties.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel
            Constructed calibration.
 
        Raises
        ------
        RuntimeError :
            Raised if the supplied dictionary is for a different
            calibration.
            Raised if the version of the supplied dictionary is 1.0.
        """
        calib = cls()
 
        if calib._OBSTYPE != (found := dictionary['metadata']['OBSTYPE']):
            raise RuntimeError(f"Incorrect brighter-fatter kernel supplied.  Expected {calib._OBSTYPE}, "
                               f"found {found}")
        calib.setMetadata(dictionary['metadata'])
        calib.calibInfoFromDict(dictionary)
 
        calib.level = dictionary['metadata'].get('LEVEL', 'AMP')
        calib.shape = (dictionary['metadata'].get('KERNEL_DX', 0),
                       dictionary['metadata'].get('KERNEL_DY', 0))
 
        calibVersion = dictionary['metadata']['bfk_VERSION']
        if calibVersion == 1.0:
            calib.log.warning("Old Version of brighter-fatter kernel found. Current version: "
                              f"{calib._VERSION}. The new attribute 'expIdMask' will be "
                              "populated with 'True' values, and the new attributes 'rawMeans'"
                              "and 'rawVariances' will be populated with the masked 'means'."
                              "and 'variances' values."
                              )
            # use 'means', because 'expIdMask' does not exist.
            calib.expIdMask = {amp: np.repeat(True, len(dictionary['means'][amp])) for amp in
                               dictionary['means']}
            calib.rawMeans = {amp: np.array(dictionary['means'][amp]) for amp in dictionary['means']}
            calib.rawVariances = {amp: np.array(dictionary['variances'][amp]) for amp in
                                  dictionary['variances']}
        elif calibVersion == 1.1:
            calib.expIdMask = {amp: np.array(dictionary['expIdMask'][amp]) for amp in dictionary['expIdMask']}
            calib.rawMeans = {amp: np.array(dictionary['rawMeans'][amp]) for amp in dictionary['rawMeans']}
            calib.rawVariances = {amp: np.array(dictionary['rawVariances'][amp]) for amp in
                                  dictionary['rawVariances']}
        else:
            raise RuntimeError(f"Unknown version for brighter-fatter kernel: {calibVersion}")
 
        # Lengths for reshape:
        _, smallLength, nObs = calib.getLengths()
        smallShapeSide = int(np.sqrt(smallLength))
 
        calib.rawXcorrs = {amp: np.array(dictionary['rawXcorrs'][amp]).reshape((nObs,
                                                                                smallShapeSide,
                                                                                smallShapeSide))
                           for amp in dictionary['rawXcorrs']}
 
        calib.gain = dictionary['gain']
        calib.noise = dictionary['noise']
 
        calib.meanXcorrs = {amp: np.array(dictionary['meanXcorrs'][amp]).reshape(calib.shape)
                            for amp in dictionary['rawXcorrs']}
        calib.ampKernels = {amp: np.array(dictionary['ampKernels'][amp]).reshape(calib.shape)
                            for amp in dictionary['ampKernels']}
        calib.valid = {amp: bool(value) for amp, value in dictionary['valid'].items()}
        calib.badAmps = [amp for amp, valid in dictionary['valid'].items() if valid is False]
 
        calib.detKernels = {det: np.array(dictionary['detKernels'][det]).reshape(calib.shape)
                            for det in dictionary['detKernels']}
 
        calib.updateMetadata()
        return calib
 
    def toDict(self):
        """Return a dictionary containing the calibration properties.
 
        The dictionary should be able to be round-tripped through
        `fromDict`.
 
        Returns
        -------
        dictionary : `dict`
            Dictionary of properties.
        """
        self.updateMetadataupdateMetadata()
 
        outDict = {}
        metadata = self.getMetadata()
        outDict['metadata'] = metadata
 
        # Lengths for ravel:
        kernelLength, smallLength, nObs = self.getLengthsgetLengths()
 
        outDict['expIdMask'] = {amp: np.array(self.expIdMaskexpIdMask[amp]).tolist() for amp in self.expIdMaskexpIdMask}
        outDict['rawMeans'] = {amp: np.array(self.rawMeansrawMeans[amp]).tolist() for amp in self.rawMeansrawMeans}
        outDict['rawVariances'] = {amp: np.array(self.rawVariancesrawVariances[amp]).tolist() for amp in
                                   self.rawVariancesrawVariances}
        outDict['rawXcorrs'] = {amp: np.array(self.rawXcorrsrawXcorrs[amp]).reshape(nObs*smallLength).tolist()
                                for amp in self.rawXcorrsrawXcorrs}
        outDict['badAmps'] = self.badAmpsbadAmps
        outDict['gain'] = self.gaingain
        outDict['noise'] = self.noisenoise
 
        outDict['meanXcorrs'] = {amp: self.meanXcorrsmeanXcorrs[amp].reshape(kernelLength).tolist()
                                 for amp in self.meanXcorrsmeanXcorrs}
        outDict['ampKernels'] = {amp: self.ampKernelsampKernels[amp].reshape(kernelLength).tolist()
                                 for amp in self.ampKernelsampKernels}
        outDict['valid'] = self.validvalid
 
        outDict['detKernels'] = {det: self.detKernelsdetKernels[det].reshape(kernelLength).tolist()
                                 for det in self.detKernelsdetKernels}
        return outDict
 
    @classmethod
    def fromTable(cls, tableList):
        """Construct calibration from a list of tables.
 
        This method uses the `fromDict` method to create the
        calibration, after constructing an appropriate dictionary from
        the input tables.
 
        Parameters
        ----------
        tableList : `list` [`astropy.table.Table`]
            List of tables to use to construct the brighter-fatter
            calibration.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel`
            The calibration defined in the tables.
        """
        ampTable = tableList[0]
 
        metadata = ampTable.meta
        inDict = dict()
        inDict['metadata'] = metadata
 
        amps = ampTable['AMPLIFIER']
 
        # Determine version for expected values.  The ``fromDict``
        # method can unpack either, but the appropriate fields need to
        # be supplied.
        calibVersion = metadata['bfk_VERSION']
 
        if calibVersion == 1.0:
            # We expect to find ``means`` and ``variances`` for this
            # case, and will construct an ``expIdMask`` from these
            # parameters in the ``fromDict`` method.
            rawMeanList = ampTable['MEANS']
            rawVarianceList = ampTable['VARIANCES']
 
            inDict['means'] = {amp: mean for amp, mean in zip(amps, rawMeanList)}
            inDict['variances'] = {amp: var for amp, var in zip(amps, rawVarianceList)}
        elif calibVersion == 1.1:
            # This will have ``rawMeans`` and ``rawVariances``, which
            # are filtered via the ``expIdMask`` fields.
            expIdMaskList = ampTable['EXP_ID_MASK']
            rawMeanList = ampTable['RAW_MEANS']
            rawVarianceList = ampTable['RAW_VARIANCES']
 
            inDict['expIdMask'] = {amp: mask for amp, mask in zip(amps, expIdMaskList)}
            inDict['rawMeans'] = {amp: mean for amp, mean in zip(amps, rawMeanList)}
            inDict['rawVariances'] = {amp: var for amp, var in zip(amps, rawVarianceList)}
        else:
            raise RuntimeError(f"Unknown version for brighter-fatter kernel: {calibVersion}")
 
        rawXcorrs = ampTable['RAW_XCORRS']
        gainList = ampTable['GAIN']
        noiseList = ampTable['NOISE']
 
        meanXcorrs = ampTable['MEAN_XCORRS']
        ampKernels = ampTable['KERNEL']
        validList = ampTable['VALID']
 
        inDict['rawXcorrs'] = {amp: kernel for amp, kernel in zip(amps, rawXcorrs)}
        inDict['gain'] = {amp: gain for amp, gain in zip(amps, gainList)}
        inDict['noise'] = {amp: noise for amp, noise in zip(amps, noiseList)}
        inDict['meanXcorrs'] = {amp: kernel for amp, kernel in zip(amps, meanXcorrs)}
        inDict['ampKernels'] = {amp: kernel for amp, kernel in zip(amps, ampKernels)}
        inDict['valid'] = {amp: bool(valid) for amp, valid in zip(amps, validList)}
 
        inDict['badAmps'] = [amp for amp, valid in inDict['valid'].items() if valid is False]
 
        if len(tableList) > 1:
            detTable = tableList[1]
            inDict['detKernels'] = {det: kernel for det, kernel
                                    in zip(detTable['DETECTOR'], detTable['KERNEL'])}
        else:
            inDict['detKernels'] = {}
 
        return cls.fromDictfromDict(inDict)
 
    def toTable(self):
        """Construct a list of tables containing the information in this
        calibration.
 
        The list of tables should create an identical calibration
        after being passed to this class's fromTable method.
 
        Returns
        -------
        tableList : `list` [`lsst.afw.table.Table`]
            List of tables containing the crosstalk calibration
            information.
 
        """
        tableList = []
        self.updateMetadataupdateMetadata()
 
        # Lengths
        kernelLength, smallLength, nObs = self.getLengthsgetLengths()
 
        ampList = []
        expIdMaskList = []
        rawMeanList = []
        rawVarianceList = []
        rawXcorrs = []
        gainList = []
        noiseList = []
 
        meanXcorrsList = []
        kernelList = []
        validList = []
 
        if self.levellevel == 'AMP':
            for amp in self.rawMeansrawMeans.keys():
                ampList.append(amp)
                expIdMaskList.append(self.expIdMaskexpIdMask[amp])
                rawMeanList.append(self.rawMeansrawMeans[amp])
                rawVarianceList.append(self.rawVariancesrawVariances[amp])
                rawXcorrs.append(np.array(self.rawXcorrsrawXcorrs[amp]).reshape(nObs*smallLength).tolist())
                gainList.append(self.gaingain[amp])
                noiseList.append(self.noisenoise[amp])
 
                meanXcorrsList.append(self.meanXcorrsmeanXcorrs[amp].reshape(kernelLength).tolist())
                kernelList.append(self.ampKernelsampKernels[amp].reshape(kernelLength).tolist())
                validList.append(int(self.validvalid[amp] and not (amp in self.badAmpsbadAmps)))
 
        ampTable = Table({'AMPLIFIER': ampList,
                          'EXP_ID_MASK': expIdMaskList,
                          'RAW_MEANS': rawMeanList,
                          'RAW_VARIANCES': rawVarianceList,
                          'RAW_XCORRS': rawXcorrs,
                          'GAIN': gainList,
                          'NOISE': noiseList,
                          'MEAN_XCORRS': meanXcorrsList,
                          'KERNEL': kernelList,
                          'VALID': validList,
                          })
 
        ampTable.meta = self.getMetadata().toDict()
        tableList.append(ampTable)
 
        if len(self.detKernelsdetKernels):
            detList = []
            kernelList = []
            for det in self.detKernelsdetKernels.keys():
                detList.append(det)
                kernelList.append(self.detKernelsdetKernels[det].reshape(kernelLength).tolist())
 
            detTable = Table({'DETECTOR': detList,
                              'KERNEL': kernelList})
            detTable.meta = self.getMetadata().toDict()
            tableList.append(detTable)
 
        return tableList
 
    # Implementation methods
    def makeDetectorKernelFromAmpwiseKernels(self, detectorName, ampsToExclude=[]):
        """Average the amplifier level kernels to create a detector level
        kernel.
        """
        inKernels = np.array([self.ampKernelsampKernels[amp] for amp in
                              self.ampKernelsampKernels if amp not in ampsToExclude])
        averagingList = np.transpose(inKernels)
        avgKernel = np.zeros_like(inKernels[0])
        sctrl = afwMath.StatisticsControl()
        sctrl.setNumSigmaClip(5.0)
        for i in range(np.shape(avgKernel)[0]):
            for j in range(np.shape(avgKernel)[1]):
                avgKernel[i, j] = afwMath.makeStatistics(averagingList[i, j],
                                                         afwMath.MEANCLIP, sctrl).getValue()
 
        self.detKernelsdetKernels[detectorName] = avgKernel
 
    def replaceDetectorKernelWithAmpKernel(self, ampName, detectorName):
        self.detKernel[detectorName] = self.ampKernel[ampName]