Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-35687: Update weights in least squares fits in PTC task #146

Merged
merged 5 commits into from
Aug 9, 2022
Merged

DM-35687: Update weights in least squares fits in PTC task #146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d1fb96f
Delete redundant comment
plazas Aug 5, 2022
f315c3c
Add more comments
plazas Aug 5, 2022
1af47b7
Update weights of fits
plazas Aug 5, 2022
367848d
Fix flake8 errors
plazas Aug 5, 2022
222f87e
Define covariance weights for fits
plazas Aug 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 3 additions & 5 deletions python/lsst/cp/pipe/ptc/cpExtractPtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -345,8 +345,6 @@ def run(self, inputExp, inputDims, taskMetadata):
self.config.maximumRangeCovariancesAstier)
for ampNumber, amp in enumerate(detector):
ampName = amp.getName()
# covAstier: [(i, j, var (cov[0,0]), cov, npix) for
# (i,j) in {maxLag, maxLag}^2]
if self.config.detectorMeasurementRegion == 'AMP':
region = amp.getBBox()
elif self.config.detectorMeasurementRegion == 'FULL':
Expand All @@ -366,7 +364,6 @@ def run(self, inputExp, inputDims, taskMetadata):
# [(i, j, var (cov[0,0]), cov, npix) for (i,j) in
# {maxLag, maxLag}^2].
muDiff, varDiff, covAstier = self.measureMeanVarCov(im1Area, im2Area, imStatsCtrl, mu1, mu2)

# Estimate the gain from the flat pair
if self.config.doGain:
gain = self.getGainFromFlatPair(im1Area, im2Area, imStatsCtrl, mu1, mu2,
Expand Down Expand Up @@ -404,6 +401,7 @@ def run(self, inputExp, inputDims, taskMetadata):
if covAstier is not None:
# Turn the tuples with the measured information
# into covariance arrays.
# covrow: (i, j, var (cov[0,0]), cov, npix)
tupleRows = [(muDiff, varDiff) + covRow + (ampNumber, expTime,
ampName) for covRow in covAstier]
tempStructArray = np.array(tupleRows, dtype=tags)
Expand Down Expand Up @@ -480,7 +478,7 @@ def makeCovArray(self, inputTuple, maxRangeFromTuple):
cov : `numpy.array`
Covariance arrays, indexed by mean signal mu.
vCov : `numpy.array`
Variance arrays, indexed by mean signal mu.
Variance of the [co]variance arrays, indexed by mean signal mu.
muVals : `numpy.array`
List of mean signal values.
"""
Expand Down Expand Up @@ -611,8 +609,8 @@ def measureMeanVarCov(self, im1Area, im2Area, imStatsCtrl, mu1, mu2):
tempSize = np.array(np.log(s)/np.log(2.)).astype(int)
fftSize = np.array(2**(tempSize+1)).astype(int)
fftShape = (fftSize[0], fftSize[1])

c = CovFastFourierTransform(diffIm.image.array, w, fftShape, maxRangeCov)
# np.sum(w) is the same as npix[0][0] returned in covDiffAstier
covDiffAstier = c.reportCovFastFourierTransform(maxRangeCov)

# Compare Cov[0,0] and afwMath.VARIANCECLIP covDiffAstier[0]
Expand Down
13 changes: 10 additions & 3 deletions python/lsst/cp/pipe/ptc/cpSolvePtcTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -912,6 +912,8 @@ def errFunc(p, x, y):
dataset.expIdMask[ampName] &= mask # bitwise_and if there is already a mask
else:
dataset.expIdMask[ampName] = mask
# In case there was a previous mask stored
mask = dataset.expIdMask[ampName]
parsIniPtc = pars
meanVecFinal = meanVecOriginal[mask]
varVecFinal = varVecOriginal[mask]
Expand All @@ -927,15 +929,20 @@ def errFunc(p, x, y):
# Fill entries with NaNs
self.fillBadAmp(dataset, ptcFitType, ampName)
continue
# Fit the PTC
# Fit the PTC.
# The variance of the variance is Var(v)=2*v^2/Npix. This is
# already calculated in `makeCovArray` of CpPtcExtract.
# dataset.covariancesSqrtWeights[ampName][:,0,0]
# has 1/sqrt(Var(v)).
weightsY = dataset.covariancesSqrtWeights[ampName][:, 0, 0][mask]
if self.config.doFitBootstrap:
parsFit, parsFitErr, reducedChiSqPtc = fitBootstrap(parsIniPtc, meanVecFinal,
varVecFinal, ptcFunc,
weightsY=1./np.sqrt(varVecFinal))
weightsY=weightsY)
else:
parsFit, parsFitErr, reducedChiSqPtc = fitLeastSq(parsIniPtc, meanVecFinal,
varVecFinal, ptcFunc,
weightsY=1./np.sqrt(varVecFinal))
weightsY=weightsY)
dataset.ptcFitPars[ampName] = parsFit
dataset.ptcFitParsError[ampName] = parsFitErr
dataset.ptcFitChiSq[ampName] = reducedChiSqPtc
Expand Down
30 changes: 23 additions & 7 deletions tests/test_ptc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ def setUp(self):
self.k2NonLinearity = -5e-6
# quadratic signal-chain non-linearity
muVec = self.flux*self.timeVec + self.k2NonLinearity*self.timeVec**2
self.gain = 1.5 # e-/ADU
self.gain = 0.75 # e-/ADU
self.c1 = 1./self.gain
self.noiseSq = 2*self.gain # 7.5 (e-)^2
self.a00 = -1.2e-6
Expand All @@ -107,10 +107,11 @@ def setUp(self):

self.ampNames = [amp.getName() for amp in self.flatExp1.getDetector().getAmplifiers()]
self.dataset = PhotonTransferCurveDataset(self.ampNames, " ") # pack raw data for fitting

self.covariancesSqrtWeights = {}
for ampName in self.ampNames: # just the expTimes and means here - vars vary per function
self.dataset.rawExpTimes[ampName] = self.timeVec
self.dataset.rawMeans[ampName] = muVec
self.covariancesSqrtWeights[ampName] = []

# ISR metadata
self.metadataContents = TaskMetadata()
Expand Down Expand Up @@ -139,7 +140,7 @@ def test_covAstier(self):
solveConfig.ptcFitType = 'FULLCOVARIANCE'
solveTask = cpPipe.ptc.PhotonTransferCurveSolveTask(config=solveConfig)

inputGain = 0.75
inputGain = self.gain

muStandard, varStandard = {}, {}
expDict = {}
Expand Down Expand Up @@ -175,7 +176,7 @@ def test_covAstier(self):
self.assertAlmostEqual(v1/v2, 1.0, places=1)

def ptcFitAndCheckPtc(self, order=None, fitType=None, doTableArray=False, doFitBootstrap=False):
localDataset = copy.copy(self.dataset)
localDataset = copy.deepcopy(self.dataset)
localDataset.ptcFitType = fitType
configSolve = copy.copy(self.defaultConfigSolve)
configLin = cpPipe.linearity.LinearitySolveTask.ConfigClass()
Expand Down Expand Up @@ -209,8 +210,24 @@ def ptcFitAndCheckPtc(self, order=None, fitType=None, doTableArray=False, doFitB
else:
RuntimeError("Enter a fit function type: 'POLYNOMIAL' or 'EXPAPPROXIMATION'")

# Initialize mask and covariance weights that will be used in fits.
# Covariance weights values empirically determined from one of
# the cases in test_covAstier.
matrixSize = localDataset.covMatrixSide
maskLength = len(localDataset.rawMeans[ampName])
for ampName in self.ampNames:
localDataset.expIdMask[ampName] = np.repeat(True, len(localDataset.rawMeans[ampName]))
localDataset.expIdMask[ampName] = np.repeat(True, maskLength)
localDataset.covariancesSqrtWeights[ampName] = np.repeat(np.ones((matrixSize, matrixSize)),
maskLength).reshape((maskLength,
matrixSize,
matrixSize))
localDataset.covariancesSqrtWeights[ampName][:, 0, 0] = [0.07980188, 0.01339653, 0.0073118,
0.00502802, 0.00383132, 0.00309475,
0.00259572, 0.00223528, 0.00196273,
0.00174943, 0.00157794, 0.00143707,
0.00131929, 0.00121935, 0.0011334,
0.00105893, 0.00099357, 0.0009358,
0.00088439, 0.00083833]
configLin.maxLookupTableAdu = 200000 # Max ADU in input mock flats
configLin.maxLinearAdu = 100000
configLin.minLinearAdu = 50000
Expand Down Expand Up @@ -440,12 +457,11 @@ def runGetGainFromFlatPair(self, correctionType='NONE'):

resultsExtract = extractTask.run(inputExp=expDict, inputDims=expIds,
taskMetadata=[self.metadataContents])

for exposurePair in resultsExtract.outputCovariances:
for ampName in self.ampNames:
if exposurePair.gain[ampName] is np.nan:
continue
self.assertAlmostEqual(exposurePair.gain[ampName], inputGain, delta=0.16)
self.assertAlmostEqual(exposurePair.gain[ampName], inputGain, delta=0.04)

def test_getGainFromFlatPair(self):
for gainCorrectionType in ['NONE', 'SIMPLE', 'FULL', ]:
Expand Down