Skip to content

Commit

Permalink
Merge pull request #276 from lsst/tickets/DM-46509
Browse files Browse the repository at this point in the history
DM-46509: Add FULLCOVRIANCE_NO_B fit option to PTC solver
  • Loading branch information
Alex-Broughton authored Nov 19, 2024
2 parents 1be69a5 + b7a62c3 commit 8882459
Show file tree
Hide file tree
Showing 4 changed files with 220 additions and 220 deletions.
4 changes: 4 additions & 0 deletions pipelines/LSSTComCam/cpCti.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
description: LSSTCam cti calibration construction
instrument: lsst.obs.lsst.LsstComCam
imports:
- location: $CP_PIPE_DIR/pipelines/_ingredients/cpCtiLSST.yaml
110 changes: 55 additions & 55 deletions python/lsst/cp/pipe/ptc/cpPtcSolve.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,8 @@ class PhotonTransferCurveSolveConfig(pipeBase.PipelineTaskConfig,
allowed={
"POLYNOMIAL": "n-degree polynomial (use 'polynomialFitDegree' to set 'n').",
"EXPAPPROXIMATION": "Approximation in Astier+19 (Eq. 16).",
"FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)"
"FULLCOVARIANCE_NO_B": "Full covariances model in Astier+19 (Eq. 15)",
"FULLCOVARIANCE": "Full covariances model in Astier+19 (Eq. 20)",
}
)
minMeanSignal = pexConfig.DictField(
Expand Down Expand Up @@ -399,7 +400,7 @@ def run(self, inputCovariances, camera=None, detId=0):
index = np.argsort(detectorMeans)
datasetPtc.sort(index)

if self.config.ptcFitType == "FULLCOVARIANCE":
if self.config.ptcFitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]:
# Fit the measured covariances vs mean signal to
# the Astier+19 full model (Eq. 20). Before that
# do a preliminary fit to the variance (C_00) vs mean
Expand All @@ -419,7 +420,7 @@ def run(self, inputCovariances, camera=None, detId=0):
# previous fit.
for ampName in datasetPtc.ampNames:
datasetPtc.expIdMask[ampName] = tempDatasetPtc.expIdMask[ampName]
datasetPtc.fitType = "FULLCOVARIANCE"
datasetPtc.fitType = self.config.ptcFitType
datasetPtc = self.fitMeasurementsToModel(datasetPtc)
# The other options are: self.config.ptcFitType in
# ("EXPAPPROXIMATION", "POLYNOMIAL")
Expand Down Expand Up @@ -485,7 +486,7 @@ def fitMeasurementsToModel(self, dataset):
`PhotonTransferCurveDatase`.
"""
fitType = dataset.ptcFitType
if fitType in ["FULLCOVARIANCE", ]:
if fitType in ["FULLCOVARIANCE", "FULLCOVARIANCE_NO_B"]:
# This model uses the full covariance matrix in the fit.
# The PTC is technically defined as variance vs signal,
# with variance = Cov_00
Expand Down Expand Up @@ -560,10 +561,7 @@ def fitDataFullCovariance(self, dataset):
dataset.covariancesSqrtWeights[ampName] = listNanMatrix
dataset.aMatrix[ampName] = nanMatrixFit
dataset.bMatrix[ampName] = nanMatrixFit
dataset.covariancesModelNoB[ampName] = listNanMatrixFit
dataset.aMatrixNoB[ampName] = nanMatrixFit
dataset.noiseMatrix[ampName] = nanMatrixFit
dataset.noiseMatrixNoB[ampName] = nanMatrixFit

dataset.expIdMask[ampName] = np.repeat(False, lenInputTimes)
dataset.gain[ampName] = np.nan
Expand Down Expand Up @@ -608,31 +606,45 @@ def fitDataFullCovariance(self, dataset):
covSqrtWeightsAtAmpForFitMasked
)

# Fit full model (Eq. 20 of Astier+19) and same model with
# b=0 (c=0 in this code).
# Fit full model (Eq. 20 of Astier+19)
pInit = np.concatenate((a0.ravel(), c0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)
functionsDict = {'fullModel': self.funcFullCovarianceModel,
'fullModelNoB': self.funcFullCovarianceModelNoB}
fitResults = {'fullModel': {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []},
'fullModelNoB': {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}}
for key in functionsDict:
params, paramsErr, _ = fitLeastSq(pInit, muAtAmpMasked,
covAtAmpForFitMasked.ravel(), functionsDict[key],
weightsY=covSqrtWeightsAtAmpForFitMasked.ravel())
a = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
gain = params[-1]

fitResults[key]['a'] = a
fitResults[key]['c'] = c
fitResults[key]['noiseMatrix'] = noiseMatrix
fitResults[key]['gain'] = gain
fitResults[key]['paramsErr'] = paramsErr

# Initialize empty results dictionary
fitResults = {'a': [], 'c': [], 'noiseMatrix': [], 'gain': [], 'paramsErr': []}

# Pick the correct full covariance model function
ptcModel = self.funcFullCovarianceModel
if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
ptcModel = self.funcFullCovarianceModelNoB
pInit = np.concatenate((a0.ravel(), noiseMatrix0.ravel(), np.array(gain0)), axis=None)

params, paramsErr, _ = fitLeastSq(
pInit,
muAtAmpMasked,
covAtAmpForFitMasked.ravel(),
ptcModel,
weightsY=covSqrtWeightsAtAmpForFitMasked.ravel(),
)

if dataset.ptcFitType == "FULLCOVARIANCE_NO_B":
zeros = np.zeros_like(params[:lenParams])
params = np.insert(params, lenParams, zeros)
paramsErr = np.insert(paramsErr, lenParams, zeros)

a = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
c = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
gain = params[-1]

fitResults['a'] = a
fitResults['c'] = c
fitResults['noiseMatrix'] = noiseMatrix
fitResults['gain'] = gain
fitResults['paramsErr'] = paramsErr

# Put the information in the PTC dataset

# Not used when ptcFitType is 'FULLCOVARIANCE'
# Not used when ptcFitType is 'FULLCOVARIANCE*'
dataset.ptcFitPars[ampName] = np.array([np.nan])
dataset.ptcFitParsError[ampName] = np.array([np.nan])
dataset.ptcFitChiSq[ampName] = np.nan
Expand All @@ -646,33 +658,24 @@ def fitDataFullCovariance(self, dataset):
# masked amps.
dataset.covariancesModel[ampName] = self.evalCovModel(
muAtAmp,
fitResults['fullModel']['a'],
fitResults['fullModel']['c'],
fitResults['fullModel']['noiseMatrix'],
fitResults['fullModel']['gain']
fitResults['a'],
fitResults['c'],
fitResults['noiseMatrix'],
fitResults['gain'],
setBtoZero=(dataset.ptcFitType == "FULLCOVARIANCE_NO_B"),
)
dataset.covariancesSqrtWeights[ampName] = covSqrtWeightsAtAmp
dataset.aMatrix[ampName] = fitResults['fullModel']['a']
dataset.bMatrix[ampName] = fitResults['fullModel']['c']/fitResults['fullModel']['a']
dataset.covariancesModelNoB[ampName] = self.evalCovModel(
muAtAmp,
fitResults['fullModelNoB']['a'],
fitResults['fullModelNoB']['c'],
fitResults['fullModelNoB']['noiseMatrix'],
fitResults['fullModelNoB']['gain'],
setBtoZero=True
)
dataset.aMatrixNoB[ampName] = fitResults['fullModelNoB']['a']
dataset.gain[ampName] = fitResults['fullModel']['gain']
dataset.gainErr[ampName] = fitResults['fullModel']['paramsErr'][-1]
readoutNoiseSquared = fitResults['fullModel']['noiseMatrix'][0][0]
dataset.aMatrix[ampName] = fitResults['a']
dataset.bMatrix[ampName] = fitResults['c']/fitResults['a']
dataset.gain[ampName] = fitResults['gain']
dataset.gainErr[ampName] = fitResults['paramsErr'][-1]
readoutNoiseSquared = fitResults['noiseMatrix'][0][0]
readoutNoise = np.sqrt(np.fabs(readoutNoiseSquared))
dataset.noise[ampName] = readoutNoise
readoutNoiseSquaredSigma = fitResults['fullModel']['paramsErr'][2*lenParams]
readoutNoiseSquaredSigma = fitResults['paramsErr'][2*lenParams]
noiseErr = 0.5*(readoutNoiseSquaredSigma/np.fabs(readoutNoiseSquared))*readoutNoise
dataset.noiseErr[ampName] = noiseErr
dataset.noiseMatrix[ampName] = fitResults['fullModel']['noiseMatrix']
dataset.noiseMatrixNoB[ampName] = fitResults['fullModelNoB']['noiseMatrix']
dataset.noiseMatrix[ampName] = fitResults['noiseMatrix']

dataset.finalVars[ampName] = covAtAmp[:, 0, 0].copy()
dataset.finalVars[ampName][~maskAtAmp] = np.nan
Expand Down Expand Up @@ -779,7 +782,7 @@ def funcFullCovarianceModelNoB(self, params, x):
Parameters
----------
params : `list`
Parameters of the model: aMatrix, CMatrix, noiseMatrix,
Parameters of the model: aMatrix, noiseMatrix,
gain (e/adu).
x : `numpy.array`, (N,)
Signal mu (adu)
Expand All @@ -792,8 +795,8 @@ def funcFullCovarianceModelNoB(self, params, x):
matrixSideFit = self.config.maximumRangeCovariancesAstierFullCovFit
lenParams = matrixSideFit*matrixSideFit
aMatrix = params[:lenParams].reshape((matrixSideFit, matrixSideFit))
cMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
noiseMatrix = params[2*lenParams:3*lenParams].reshape((matrixSideFit, matrixSideFit))
cMatrix = np.zeros_like(aMatrix)
noiseMatrix = params[lenParams:2*lenParams].reshape((matrixSideFit, matrixSideFit))
gain = params[-1]

return self.evalCovModel(x, aMatrix, cMatrix, noiseMatrix, gain, setBtoZero=True).flatten()
Expand Down Expand Up @@ -1090,10 +1093,7 @@ def fitPtc(self, dataset):
dataset.covariancesModel[amp] = listNanMatrixFit
dataset.aMatrix[amp] = nanMatrixFit
dataset.bMatrix[amp] = nanMatrixFit
dataset.covariancesModelNoB[amp] = listNanMatrixFit
dataset.aMatrixNoB[amp] = nanMatrixFit
dataset.noiseMatrix[amp] = nanMatrixFit
dataset.noiseMatrixNoB[amp] = nanMatrixFit

def errFunc(p, x, y):
return ptcFunc(p, x) - y
Expand Down
2 changes: 1 addition & 1 deletion tests/test_defects.py
Original file line number Diff line number Diff line change
Expand Up @@ -549,7 +549,7 @@ def test_maskBlocks_y_out_of_order_dm38103(self):
"""A test for maskBlocksIfIntermitentBadPixelsInColumn, y out of order.
This test is a variant of
notest_maskBlocks_every_other_pixel_bad_greater_than_threshold with
test_maskBlocks_every_other_pixel_bad_greater_than_threshold with
an extra out-of-y-order bad pixel to trigger DM-38103.
"""
expectedDefects = [Box2I(corner=Point2I(50, 110), dimensions=Extent2I(1, 31))]
Expand Down
Loading

0 comments on commit 8882459

Please sign in to comment.