lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask Class Reference

Inheritance diagram for lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask:

Public Member Functions
def	runQuantum (self, butlerQC, inputRefs, outputRefs)
def	run (self, inputPtc, dummy, camera, inputDims)
def	averageCorrelations (self, xCorrList, name)
def	quadraticCorrelations (self, xCorrList, fluxList, name)
def	successiveOverRelax (self, source, maxIter=None, eLevel=None)

Static Public Attributes
	ConfigClass = BrighterFatterKernelSolveConfig

Detailed Description

Measure appropriate Brighter-Fatter Kernel from the PTC dataset.

Definition at line 142 of file makeBrighterFatterKernel.py.

Member Function Documentation

averageCorrelations()

def lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.averageCorrelations	(		self,
			xCorrList,
			name
	)

Average input correlations.

Parameters
----------
xCorrList : `list` [`numpy.array`]
    List of cross-correlations.
name : `str`
    Name for log messages.

Returns
-------
meanXcorr : `numpy.array`
    The averaged cross-correlation.

Definition at line 327 of file makeBrighterFatterKernel.py.

    def averageCorrelations(self, xCorrList, name):
        """Average input correlations.
 
        Parameters
        ----------
        xCorrList : `list` [`numpy.array`]
            List of cross-correlations.
        name : `str`
            Name for log messages.
 
        Returns
        -------
        meanXcorr : `numpy.array`
            The averaged cross-correlation.
        """
        meanXcorr = np.zeros_like(xCorrList[0])
        xCorrList = np.transpose(xCorrList)
        sctrl = afwMath.StatisticsControl()
        sctrl.setNumSigmaClip(self.config.nSigmaClip)
        for i in range(np.shape(meanXcorr)[0]):
            for j in range(np.shape(meanXcorr)[1]):
                meanXcorr[i, j] = afwMath.makeStatistics(xCorrList[i, j],
                                                         afwMath.MEANCLIP, sctrl).getValue()
 
        # To match previous definitions, pad by one element.
        meanXcorr = np.pad(meanXcorr, ((1, 1)))
 
        return meanXcorr
 

quadraticCorrelations()

def lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.quadraticCorrelations	(		self,
			xCorrList,
			fluxList,
			name
	)

Measure a quadratic correlation model.

Parameters
----------
xCorrList : `list` [`numpy.array`]
    List of cross-correlations.
fluxList : `numpy.array`
    Associated list of fluxes.
name : `str`
    Name for log messages.

Returns
-------
meanXcorr : `numpy.array`
    The averaged cross-correlation.

Definition at line 356 of file makeBrighterFatterKernel.py.

    def quadraticCorrelations(self, xCorrList, fluxList, name):
        """Measure a quadratic correlation model.
 
        Parameters
        ----------
        xCorrList : `list` [`numpy.array`]
            List of cross-correlations.
        fluxList : `numpy.array`
            Associated list of fluxes.
        name : `str`
            Name for log messages.
 
        Returns
        -------
        meanXcorr : `numpy.array`
            The averaged cross-correlation.
        """
        meanXcorr = np.zeros_like(xCorrList[0])
        fluxList = np.square(fluxList)
        xCorrList = np.array(xCorrList)
 
        for i in range(np.shape(meanXcorr)[0]):
            for j in range(np.shape(meanXcorr)[1]):
                # Fit corrlation_i(x, y) = a0 + a1 * (flux_i)^2 The
                # i,j indices are inverted to apply the transposition,
                # as is done in the averaging case.
                linearFit, linearFitErr, chiSq, weights = irlsFit([0.0, 1e-4], fluxList,
                                                                  xCorrList[:, j, i], funcPolynomial)
                meanXcorr[i, j] = linearFit[1]  # Discard the intercept.
                self.log.debug("Quad fit meanXcorr[%d,%d] = %g", i, j, linearFit[1])
 
        # To match previous definitions, pad by one element.
        meanXcorr = np.pad(meanXcorr, ((1, 1)))
 
        return meanXcorr
 

run()

def lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.run	(		self,
			inputPtc,
			dummy,
			camera,
			inputDims
	)

Combine covariance information from PTC into brighter-fatter kernels.

Parameters
----------
inputPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
    PTC data containing per-amplifier covariance measurements.
dummy : `lsst.afw.image.Exposure`
    The exposure used to select the appropriate PTC dataset.
    In almost all circumstances, one of the input exposures
    used to generate the PTC dataset is the best option.
camera : `lsst.afw.cameraGeom.Camera`
    Camera to use for camera geometry information.
inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
    DataIds to use to populate the output calibration.

Returns
-------
results : `lsst.pipe.base.Struct`
    The resulst struct containing:

    ``outputBfk`` : `lsst.ip.isr.BrighterFatterKernel`
        Resulting Brighter-Fatter Kernel.

Definition at line 169 of file makeBrighterFatterKernel.py.

    def run(self, inputPtc, dummy, camera, inputDims):
        """Combine covariance information from PTC into brighter-fatter kernels.
 
        Parameters
        ----------
        inputPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
            PTC data containing per-amplifier covariance measurements.
        dummy : `lsst.afw.image.Exposure`
            The exposure used to select the appropriate PTC dataset.
            In almost all circumstances, one of the input exposures
            used to generate the PTC dataset is the best option.
        camera : `lsst.afw.cameraGeom.Camera`
            Camera to use for camera geometry information.
        inputDims : `lsst.daf.butler.DataCoordinate` or `dict`
            DataIds to use to populate the output calibration.
 
        Returns
        -------
        results : `lsst.pipe.base.Struct`
            The resulst struct containing:
 
            ``outputBfk`` : `lsst.ip.isr.BrighterFatterKernel`
                Resulting Brighter-Fatter Kernel.
        """
        if len(dummy) == 0:
            self.log.warn("No dummy exposure found.")
 
        detector = camera[inputDims['detector']]
        detName = detector.getName()
 
        if self.config.level == 'DETECTOR':
            detectorCorrList = list()
 
        bfk = BrighterFatterKernel(camera=camera, detectorId=detector.getId(), level=self.config.level)
        bfk.means = inputPtc.finalMeans  # ADU
        bfk.variances = inputPtc.finalVars  # ADU^2
        # Use the PTC covariances as the cross-correlations.  These
        # are scaled before the kernel is generated, which performs
        # the conversion.
        bfk.rawXcorrs = inputPtc.covariances  # ADU^2
        bfk.badAmps = inputPtc.badAmps
        bfk.shape = (inputPtc.covMatrixSide*2 + 1, inputPtc.covMatrixSide*2 + 1)
        bfk.gain = inputPtc.gain
        bfk.noise = inputPtc.noise
        bfk.meanXcorrs = dict()
        bfk.valid = dict()
 
        for amp in detector:
            ampName = amp.getName()
            mask = np.array(inputPtc.expIdMask[ampName], dtype=bool)
 
            gain = bfk.gain[ampName]
            fluxes = np.array(bfk.means[ampName])[mask]
            variances = np.array(bfk.variances[ampName])[mask]
            xCorrList = [np.array(xcorr) for xcorr in bfk.rawXcorrs[ampName]]
            xCorrList = np.array(xCorrList)[mask]
 
            if gain <= 0:
                # We've received very bad data.
                self.log.warn("Impossible gain recieved from PTC for %s: %f.  Skipping amplifier.",
                              ampName, gain)
                bfk.meanXcorrs[ampName] = np.zeros(bfk.shape)
                bfk.ampKernels[ampName] = np.zeros(bfk.shape)
                bfk.valid[ampName] = False
                continue
 
            fluxes = np.array([flux*gain for flux in fluxes])  # Now in e^-
            variances = np.array([variance*gain*gain for variance in variances])  # Now in e^2-
 
            # This should duplicate Coulton et al. 2017 Equation 22-29 (arxiv:1711.06273)
            scaledCorrList = list()
            for xcorrNum, (xcorr, flux, var) in enumerate(zip(xCorrList, fluxes, variances), 1):
                q = np.array(xcorr) * gain * gain  # xcorr now in e^-
                q *= 2.0  # Remove factor of 1/2 applied in PTC.
                self.log.info("Amp: %s %d/%d Flux: %f  Var: %f  Q(0,0): %g  Q(1,0): %g  Q(0,1): %g",
                              ampName, xcorrNum, len(xCorrList), flux, var, q[0][0], q[1][0], q[0][1])
 
                # Normalize by the flux, which removes the (0,0)
                # component attributable to Poisson noise.  This
                # contains the two "t I delta(x - x')" terms in Coulton et al. 2017 equation 29
                q[0][0] -= 2.0*(flux)
 
                if q[0][0] > 0.0:
                    self.log.warn("Amp: %s %d skipped due to value of (variance-mean)=%f",
                                  ampName, xcorrNum, q[0][0])
                    continue
 
                # This removes the "t (I_a^2 + I_b^2)" factor in Coulton et al. 2017 equation 29.
                q /= -2.0*(flux**2)
                scaled = self._tileArray(q)
 
                xcorrCheck = np.abs(np.sum(scaled))/np.sum(np.abs(scaled))
                if (xcorrCheck > self.config.xcorrCheckRejectLevel) or not (np.isfinite(xcorrCheck)):
                    self.log.warn("Amp: %s %d skipped due to value of triangle-inequality sum %f",
                                  ampName, xcorrNum, xcorrCheck)
                    continue
 
                scaledCorrList.append(scaled)
                self.log.info("Amp: %s %d/%d  Final: %g  XcorrCheck: %f",
                              ampName, xcorrNum, len(xCorrList), q[0][0], xcorrCheck)
 
            if len(scaledCorrList) == 0:
                self.log.warn("Amp: %s All inputs rejected for amp!", ampName)
                bfk.meanXcorrs[ampName] = np.zeros(bfk.shape)
                bfk.ampKernels[ampName] = np.zeros(bfk.shape)
                bfk.valid[ampName] = False
                continue
 
            if self.config.useAmatrix:
                # Use the aMatrix, ignoring the meanXcorr generated above.
                preKernel = np.pad(self._tileArray(np.array(inputPtc.aMatrix[ampName])), ((1, 1)))
            elif self.config.correlationQuadraticFit:
                # Use a quadratic fit to the correlations as a function of flux.
                preKernel = self.quadraticCorrelations(scaledCorrList, fluxes, f"Amp: {ampName}")
            else:
                # Use a simple average of the measured correlations.
                preKernel = self.averageCorrelations(scaledCorrList, f"Amp: {ampName}")
 
            center = int((bfk.shape[0] - 1) / 2)
 
            if self.config.forceZeroSum:
                totalSum = np.sum(preKernel)
 
                if self.config.correlationModelRadius < (preKernel.shape[0] - 1) / 2:
                    # Assume a correlation model of Corr(r) = -preFactor * r^(2 * slope)
                    preFactor = np.sqrt(preKernel[center, center + 1] * preKernel[center + 1, center])
                    slopeFactor = 2.0 * np.abs(self.config.correlationModelSlope)
                    totalSum += 2.0*np.pi*(preFactor / (slopeFactor*(center + 0.5))**slopeFactor)
 
                preKernel[center, center] -= totalSum
                self.log.info("%s Zero-Sum Scale: %g", ampName, totalSum)
 
            finalSum = np.sum(preKernel)
            bfk.meanXcorrs[ampName] = preKernel
 
            postKernel = self.successiveOverRelax(preKernel)
            bfk.ampKernels[ampName] = postKernel
            if self.config.level == 'DETECTOR':
                detectorCorrList.extend(scaledCorrList)
            bfk.valid[ampName] = True
            self.log.info("Amp: %s Sum: %g  Center Info Pre: %g  Post: %g",
                          ampName, finalSum, preKernel[center, center], postKernel[center, center])
 
        # Assemble a detector kernel?
        if self.config.level == 'DETECTOR':
            preKernel = self.averageCorrelations(detectorCorrList, f"Det: {detName}")
            finalSum = np.sum(preKernel)
            center = int((bfk.shape[0] - 1) / 2)
 
            postKernel = self.successiveOverRelax(preKernel)
            bfk.detKernels[detName] = postKernel
            self.log.info("Det: %s Sum: %g  Center Info Pre: %g  Post: %g",
                          detName, finalSum, preKernel[center, center], postKernel[center, center])
 
        return pipeBase.Struct(
            outputBFK=bfk,
        )
 

runQuantum()

def lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.runQuantum	(		self,
			butlerQC,
			inputRefs,
			outputRefs
	)

Ensure that the input and output dimensions are passed along.

Parameters
----------
butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
    Butler to operate on.
inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
    Input data refs to load.
ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
    Output data refs to persist.

Definition at line 149 of file makeBrighterFatterKernel.py.

    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        """Ensure that the input and output dimensions are passed along.
 
        Parameters
        ----------
        butlerQC : `lsst.daf.butler.butlerQuantumContext.ButlerQuantumContext`
            Butler to operate on.
        inputRefs : `lsst.pipe.base.connections.InputQuantizedConnection`
            Input data refs to load.
        ouptutRefs : `lsst.pipe.base.connections.OutputQuantizedConnection`
            Output data refs to persist.
        """
        inputs = butlerQC.get(inputRefs)
 
        # Use the dimensions to set calib/provenance information.
        inputs['inputDims'] = inputRefs.inputPtc.dataId.byName()
 
        outputs = self.run(**inputs)
        butlerQC.put(outputs, outputRefs)
 

successiveOverRelax()

def lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.successiveOverRelax	(		self,
			source,
			maxIter = None,
			eLevel = None
	)

An implementation of the successive over relaxation (SOR) method.

A numerical method for solving a system of linear equations
with faster convergence than the Gauss-Seidel method.

Parameters:
-----------
source : `numpy.ndarray`
    The input array.
maxIter : `int`, optional
    Maximum number of iterations to attempt before aborting.
eLevel : `float`, optional
    The target error level at which we deem convergence to have
    occurred.

Returns:
--------
output : `numpy.ndarray`
    The solution.

Definition at line 432 of file makeBrighterFatterKernel.py.

    def successiveOverRelax(self, source, maxIter=None, eLevel=None):
        """An implementation of the successive over relaxation (SOR) method.
 
        A numerical method for solving a system of linear equations
        with faster convergence than the Gauss-Seidel method.
 
        Parameters:
        -----------
        source : `numpy.ndarray`
            The input array.
        maxIter : `int`, optional
            Maximum number of iterations to attempt before aborting.
        eLevel : `float`, optional
            The target error level at which we deem convergence to have
            occurred.
 
        Returns:
        --------
        output : `numpy.ndarray`
            The solution.
        """
        if not maxIter:
            maxIter = self.config.maxIterSuccessiveOverRelaxation
        if not eLevel:
            eLevel = self.config.eLevelSuccessiveOverRelaxation
 
        assert source.shape[0] == source.shape[1], "Input array must be square"
        # initialize, and set boundary conditions
        func = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
        resid = np.zeros([source.shape[0] + 2, source.shape[1] + 2])
        rhoSpe = np.cos(np.pi/source.shape[0])  # Here a square grid is assumed
 
        # Calculate the initial error
        for i in range(1, func.shape[0] - 1):
            for j in range(1, func.shape[1] - 1):
                resid[i, j] = (func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
                               + func[i + 1, j] - 4*func[i, j] - source[i - 1, j - 1])
        inError = np.sum(np.abs(resid))
 
        # Iterate until convergence
        # We perform two sweeps per cycle,
        # updating 'odd' and 'even' points separately
        nIter = 0
        omega = 1.0
        dx = 1.0
        while nIter < maxIter*2:
            outError = 0
            if nIter%2 == 0:
                for i in range(1, func.shape[0] - 1, 2):
                    for j in range(1, func.shape[1] - 1, 2):
                        resid[i, j] = float(func[i, j-1] + func[i, j + 1] + func[i - 1, j]
                                            + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
                        func[i, j] += omega*resid[i, j]*.25
                for i in range(2, func.shape[0] - 1, 2):
                    for j in range(2, func.shape[1] - 1, 2):
                        resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
                                            + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
                        func[i, j] += omega*resid[i, j]*.25
            else:
                for i in range(1, func.shape[0] - 1, 2):
                    for j in range(2, func.shape[1] - 1, 2):
                        resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
                                            + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
                        func[i, j] += omega*resid[i, j]*.25
                for i in range(2, func.shape[0] - 1, 2):
                    for j in range(1, func.shape[1] - 1, 2):
                        resid[i, j] = float(func[i, j - 1] + func[i, j + 1] + func[i - 1, j]
                                            + func[i + 1, j] - 4.0*func[i, j] - dx*dx*source[i - 1, j - 1])
                        func[i, j] += omega*resid[i, j]*.25
            outError = np.sum(np.abs(resid))
            if outError < inError*eLevel:
                break
            if nIter == 0:
                omega = 1.0/(1 - rhoSpe*rhoSpe/2.0)
            else:
                omega = 1.0/(1 - rhoSpe*rhoSpe*omega/4.0)
            nIter += 1
 
        if nIter >= maxIter*2:
            self.log.warn("Failure: SuccessiveOverRelaxation did not converge in %s iterations."
                          "\noutError: %s, inError: %s," % (nIter//2, outError, inError*eLevel))
        else:
            self.log.info("Success: SuccessiveOverRelaxation converged in %s iterations."
                          "\noutError: %s, inError: %s", nIter//2, outError, inError*eLevel)
        return func[1: -1, 1: -1]

Member Data Documentation

ConfigClass

lsst.cp.pipe.makeBrighterFatterKernel.BrighterFatterKernelSolveTask.ConfigClass = BrighterFatterKernelSolveConfig

static

Definition at line 146 of file makeBrighterFatterKernel.py.

---

The documentation for this class was generated from the following file:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.6.0/Linux64/cp_pipe/22.0.1-18-gb17765a+2264247a6b/python/lsst/cp/pipe/makeBrighterFatterKernel.py