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.
 
    """
    _OBSTYPE = 'bfk'
    _SCHEMA = 'Brighter-fatter kernel'
    _VERSION = 1.0
 
    def __init__(self, camera=None, level=None, **kwargs):
        self.levellevel = level
 
        # Things inherited from the PTC
        self.meansmeans = dict()
        self.variancesvariances = 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', 'means', 'variances', '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.meansmeans[ampName] = []
                self.variancesvariances[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.shape.
        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.meansmeans[amp]) for amp in self.meansmeans])
            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.
        """
        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))
 
        calib.means = {amp: np.array(dictionary['means'][amp]) for amp in dictionary['means']}
        calib.variances = {amp: np.array(dictionary['variances'][amp]) for amp in dictionary['variances']}
 
        # 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['means'] = {amp: np.array(self.meansmeans[amp]).tolist() for amp in self.meansmeans}
        outDict['variances'] = {amp: np.array(self.variancesvariances[amp]).tolist() for amp in self.variancesvariances}
        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']
 
        meanList = ampTable['MEANS']
        varianceList = ampTable['VARIANCES']
 
        rawXcorrs = ampTable['RAW_XCORRS']
        gainList = ampTable['GAIN']
        noiseList = ampTable['NOISE']
 
        meanXcorrs = ampTable['MEAN_XCORRS']
        ampKernels = ampTable['KERNEL']
        validList = ampTable['VALID']
 
        inDict['means'] = {amp: mean for amp, mean in zip(amps, meanList)}
        inDict['variances'] = {amp: var for amp, var in zip(amps, varianceList)}
        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 = []
        meanList = []
        varianceList = []
        rawXcorrs = []
        gainList = []
        noiseList = []
 
        meanXcorrsList = []
        kernelList = []
        validList = []
 
        if self.levellevel == 'AMP':
            for amp in self.meansmeans.keys():
                ampList.append(amp)
                meanList.append(self.meansmeans[amp])
                varianceList.append(self.variancesvariances[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,
                          'MEANS': meanList,
                          'VARIANCES': varianceList,
                          '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]