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.

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

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

Merged
merged 5 commits into from
Aug 9, 2022
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
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
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
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