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-23762: Check bad amps in LSSTCam are being tracked in the defects file #34

Merged
merged 3 commits into from
May 1, 2020
Merged

DM-23762: Check bad amps in LSSTCam are being tracked in the defects file #34

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f816009
Handle bad amps when plotting and persisting
plazas Apr 28, 2020
a2f23f7
Address comments by reviewer
plazas Apr 30, 2020
7c429ec
Replace list-comps
plazas May 1, 2020
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
165 changes: 95 additions & 70 deletions python/lsst/cp/pipe/ptc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -968,6 +968,9 @@ def errFunc(p, x, y):
msg = (f"\nSERIOUS: Not enough data points ({len(meanVecFinal)}) compared to the number of"
f"parameters of the PTC model({len(parsIniPtc)}). Setting {ampName} to BAD.")
self.log.warn(msg)
# The first and second parameters of initial fit are discarded (bias and gain)
# for the final NL coefficients
lenNonLinPars = self.config.polynomialFitDegreeNonLinearity - 1
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If anything, shouldn't this be +1 not -1? If I fit an x^2 i.e. 2nd order poly, I'd expect 3 coeffs, not 1, right?

dataset.badAmps.append(ampName)
dataset.gain[ampName] = np.nan
dataset.gainErr[ampName] = np.nan
Expand All @@ -978,6 +981,11 @@ def errFunc(p, x, y):
dataset.nonLinearityResiduals[ampName] = np.nan
dataset.fractionalNonLinearityResiduals[ampName] = np.nan
dataset.coefficientLinearizeSquared[ampName] = np.nan
dataset.ptcFitPars[ampName] = np.nan
dataset.ptcFitParsError[ampName] = np.nan
dataset.ptcFitReducedChiSquared[ampName] = np.nan
dataset.coefficientsLinearizePolynomial[ampName] = [np.nan]*lenNonLinPars
tableArray[i, :] = [np.nan]*self.config.maxAduForLookupTableLinearizer
continue

# Fit the PTC
Expand Down Expand Up @@ -1088,54 +1096,61 @@ def _plotPtc(self, dataset, ptcFitType, pdfPages):
pars, parsErr = dataset.ptcFitPars[amp], dataset.ptcFitParsError[amp]

if ptcFitType == 'ASTIERAPPROXIMATION':
ptcA00, ptcA00error = pars[0], parsErr[0]
ptcGain, ptcGainError = pars[1], parsErr[1]
ptcNoise = np.sqrt(np.fabs(pars[2]))
ptcNoiseError = 0.5*(parsErr[2]/np.fabs(pars[2]))*np.sqrt(np.fabs(pars[2]))
stringLegend = (f"a00: {ptcA00:.2e}+/-{ptcA00error:.2e} 1/e"
f"\n Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN"
f"\n Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n")
if len(meanVecFinal):
ptcA00, ptcA00error = pars[0], parsErr[0]
ptcGain, ptcGainError = pars[1], parsErr[1]
ptcNoise = np.sqrt(np.fabs(pars[2]))
ptcNoiseError = 0.5*(parsErr[2]/np.fabs(pars[2]))*np.sqrt(np.fabs(pars[2]))
stringLegend = (f"a00: {ptcA00:.2e}+/-{ptcA00error:.2e} 1/e"
f"\n Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN"
f"\n Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n")

if ptcFitType == 'POLYNOMIAL':
ptcGain, ptcGainError = 1./pars[1], np.fabs(1./pars[1])*(parsErr[1]/pars[1])
ptcNoise = np.sqrt(np.fabs(pars[0]))*ptcGain
ptcNoiseError = (0.5*(parsErr[0]/np.fabs(pars[0]))*(np.sqrt(np.fabs(pars[0]))))*ptcGain
stringLegend = (f"Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n"
f"Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN \n")

minMeanVecFinal = np.min(meanVecFinal)
maxMeanVecFinal = np.max(meanVecFinal)
meanVecFit = np.linspace(minMeanVecFinal, maxMeanVecFinal, 100*len(meanVecFinal))
minMeanVecOriginal = np.min(meanVecOriginal)
maxMeanVecOriginal = np.max(meanVecOriginal)
deltaXlim = maxMeanVecOriginal - minMeanVecOriginal

a.plot(meanVecFit, ptcFunc(pars, meanVecFit), color='red')
a.plot(meanVecFinal, pars[0] + pars[1]*meanVecFinal, color='green', linestyle='--')
a.scatter(meanVecFinal, varVecFinal, c='blue', marker='o', s=markerSize)
a.scatter(meanVecOutliers, varVecOutliers, c='magenta', marker='s', s=markerSize)
a.set_xlabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
a.set_xticks(meanVecOriginal)
a.set_ylabel(r'Variance (DN$^2$)', fontsize=labelFontSize)
a.tick_params(labelsize=11)
a.text(0.03, 0.8, stringLegend, transform=a.transAxes, fontsize=legendFontSize)
a.set_xscale('linear', fontsize=labelFontSize)
a.set_yscale('linear', fontsize=labelFontSize)
a.set_title(amp, fontsize=titleFontSize)
a.set_xlim([minMeanVecOriginal - 0.2*deltaXlim, maxMeanVecOriginal + 0.2*deltaXlim])

# Same, but in log-scale
a2.plot(meanVecFit, ptcFunc(pars, meanVecFit), color='red')
a2.scatter(meanVecFinal, varVecFinal, c='blue', marker='o', s=markerSize)
a2.scatter(meanVecOutliers, varVecOutliers, c='magenta', marker='s', s=markerSize)
a2.set_xlabel(r'Mean Signal ($\mu$, DN)', fontsize=labelFontSize)
a2.set_ylabel(r'Variance (DN$^2$)', fontsize=labelFontSize)
a2.tick_params(labelsize=11)
a2.text(0.03, 0.8, stringLegend, transform=a2.transAxes, fontsize=legendFontSize)
a2.set_xscale('log')
a2.set_yscale('log')
a2.set_title(amp, fontsize=titleFontSize)
a2.set_xlim([minMeanVecOriginal, maxMeanVecOriginal])
if len(meanVecFinal):
ptcGain, ptcGainError = 1./pars[1], np.fabs(1./pars[1])*(parsErr[1]/pars[1])
ptcNoise = np.sqrt(np.fabs(pars[0]))*ptcGain
ptcNoiseError = (0.5*(parsErr[0]/np.fabs(pars[0]))*(np.sqrt(np.fabs(pars[0]))))*ptcGain
stringLegend = (f"Noise: {ptcNoise:.4}+/-{ptcNoiseError:.2e} e \n"
f"Gain: {ptcGain:.4}+/-{ptcGainError:.2e} e/DN \n")

a.set_xlabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
a.set_ylabel(r'Variance (DN$^2$)', fontsize=labelFontSize)
a.tick_params(labelsize=11)
a.set_xscale('linear', fontsize=labelFontSize)
a.set_yscale('linear', fontsize=labelFontSize)

a2.set_xlabel(r'Mean Signal ($\mu$, DN)', fontsize=labelFontSize)
a2.set_ylabel(r'Variance (DN$^2$)', fontsize=labelFontSize)
a2.tick_params(labelsize=11)
a2.set_xscale('log')
a2.set_yscale('log')

if len(meanVecFinal): # Empty if the whole amp is bad, for example.
minMeanVecFinal = np.min(meanVecFinal)
maxMeanVecFinal = np.max(meanVecFinal)
meanVecFit = np.linspace(minMeanVecFinal, maxMeanVecFinal, 100*len(meanVecFinal))
minMeanVecOriginal = np.min(meanVecOriginal)
maxMeanVecOriginal = np.max(meanVecOriginal)
deltaXlim = maxMeanVecOriginal - minMeanVecOriginal

a.plot(meanVecFit, ptcFunc(pars, meanVecFit), color='red')
a.plot(meanVecFinal, pars[0] + pars[1]*meanVecFinal, color='green', linestyle='--')
a.scatter(meanVecFinal, varVecFinal, c='blue', marker='o', s=markerSize)
a.scatter(meanVecOutliers, varVecOutliers, c='magenta', marker='s', s=markerSize)
a.text(0.03, 0.8, stringLegend, transform=a.transAxes, fontsize=legendFontSize)
a.set_title(amp, fontsize=titleFontSize)
a.set_xlim([minMeanVecOriginal - 0.2*deltaXlim, maxMeanVecOriginal + 0.2*deltaXlim])

# Same, but in log-scale
a2.plot(meanVecFit, ptcFunc(pars, meanVecFit), color='red')
a2.scatter(meanVecFinal, varVecFinal, c='blue', marker='o', s=markerSize)
a2.scatter(meanVecOutliers, varVecOutliers, c='magenta', marker='s', s=markerSize)
a2.text(0.03, 0.8, stringLegend, transform=a2.transAxes, fontsize=legendFontSize)
a2.set_title(amp, fontsize=titleFontSize)
a2.set_xlim([minMeanVecOriginal, maxMeanVecOriginal])
else:
a.set_title(f"{amp} (BAD)", fontsize=titleFontSize)
a2.set_title(f"{amp} (BAD)", fontsize=titleFontSize)

f.suptitle(f"PTC \n Fit: " + stringTitle, fontsize=20)
pdfPages.savefig(f)
Expand All @@ -1148,67 +1163,77 @@ def _plotPtc(self, dataset, ptcFitType, pdfPages):
meanVecFinal = np.array(dataset.rawMeans[amp])[dataset.visitMask[amp]]
timeVecFinal = np.array(dataset.rawExpTimes[amp])[dataset.visitMask[amp]]

pars, parsErr = dataset.nonLinearity[amp], dataset.nonLinearityError[amp]
k0, k0Error = pars[0], parsErr[0]
k1, k1Error = pars[1], parsErr[1]
k2, k2Error = pars[2], parsErr[2]
stringLegend = (f"k0: {k0:.4}+/-{k0Error:.2e} DN\n k1: {k1:.4}+/-{k1Error:.2e} DN/t"
f"\n k2: {k2:.2e}+/-{k2Error:.2e} DN/t^2 \n")
a.scatter(timeVecFinal, meanVecFinal)
a.plot(timeVecFinal, self.funcPolynomial(pars, timeVecFinal), color='red')
a.set_xlabel('Time (sec)', fontsize=labelFontSize)
a.set_xticks(timeVecFinal)
a.set_ylabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
a.tick_params(labelsize=labelFontSize)
a.text(0.03, 0.75, stringLegend, transform=a.transAxes, fontsize=legendFontSize)
a.set_xscale('linear', fontsize=labelFontSize)
a.set_yscale('linear', fontsize=labelFontSize)
a.set_title(amp, fontsize=titleFontSize)

if len(meanVecFinal):
pars, parsErr = dataset.nonLinearity[amp], dataset.nonLinearityError[amp]
k0, k0Error = pars[0], parsErr[0]
k1, k1Error = pars[1], parsErr[1]
k2, k2Error = pars[2], parsErr[2]
stringLegend = (f"k0: {k0:.4}+/-{k0Error:.2e} DN\n k1: {k1:.4}+/-{k1Error:.2e} DN/t"
f"\n k2: {k2:.2e}+/-{k2Error:.2e} DN/t^2 \n")
a.scatter(timeVecFinal, meanVecFinal)
a.plot(timeVecFinal, self.funcPolynomial(pars, timeVecFinal), color='red')
a.text(0.03, 0.75, stringLegend, transform=a.transAxes, fontsize=legendFontSize)
a.set_title(f"{amp}", fontsize=titleFontSize)
else:
a.set_title(f"{amp} (BAD)", fontsize=titleFontSize)

f.suptitle("Linearity \n Fit: Polynomial (degree: %g)"
% (self.config.polynomialFitDegreeNonLinearity),
fontsize=supTitleFontSize)
pdfPages.savefig()
pdfPages.savefig(f)

# Plot linearity residual
f, ax = plt.subplots(nrows=4, ncols=4, sharex='col', sharey='row', figsize=(13, 10))
for i, (amp, a) in enumerate(zip(dataset.ampNames, ax.flatten())):
meanVecFinal = np.array(dataset.rawMeans[amp])[dataset.visitMask[amp]]
linRes = np.array(dataset.nonLinearityResiduals[amp])

a.scatter(meanVecFinal, linRes)
a.axhline(y=0, color='k')
a.axvline(x=timeVecFinal[self.config.linResidualTimeIndex], color='g', linestyle='--')
a.set_xlabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
a.set_xticks(meanVecFinal)
a.set_ylabel('LR (%)', fontsize=labelFontSize)
a.tick_params(labelsize=labelFontSize)
a.set_xscale('linear', fontsize=labelFontSize)
a.set_yscale('linear', fontsize=labelFontSize)
a.set_title(amp, fontsize=titleFontSize)
if len(meanVecFinal):
a.axvline(x=timeVecFinal[self.config.linResidualTimeIndex], color='g', linestyle='--')
a.scatter(meanVecFinal, linRes)
a.set_title(f"{amp}", fontsize=titleFontSize)
a.axhline(y=0, color='k')
a.tick_params(labelsize=labelFontSize)
else:
a.set_title(f"{amp} (BAD)", fontsize=titleFontSize)

f.suptitle(r"Linearity Residual: $100\times(1 - \mu_{\rm{ref}}/t_{\rm{ref}})/(\mu / t))$" + "\n" +
r"$t_{\rm{ref}}$: " + f"{timeVecFinal[2]} s", fontsize=supTitleFontSize)
pdfPages.savefig()
pdfPages.savefig(f)

# Plot fractional non-linearity residual (w.r.t linear part of polynomial fit)
f, ax = plt.subplots(nrows=4, ncols=4, sharex='col', sharey='row', figsize=(13, 10))
for i, (amp, a) in enumerate(zip(dataset.ampNames, ax.flatten())):
meanVecFinal = np.array(dataset.rawMeans[amp])[dataset.visitMask[amp]]
fracLinRes = np.array(dataset.fractionalNonLinearityResiduals[amp])
a.scatter(meanVecFinal, fracLinRes, c='g')

a.axhline(y=0, color='k')
a.axvline(x=0, color='k', linestyle='-')
a.set_xlabel(r'Mean signal ($\mu$, DN)', fontsize=labelFontSize)
a.set_xticks(meanVecFinal)
a.set_ylabel('Fractional nonlinearity (%)', fontsize=labelFontSize)
a.tick_params(labelsize=labelFontSize)
a.set_xscale('linear', fontsize=labelFontSize)
a.set_yscale('linear', fontsize=labelFontSize)
a.set_title(amp, fontsize=titleFontSize)

if len(meanVecFinal):
a.scatter(meanVecFinal, fracLinRes, c='g')
a.set_title(amp, fontsize=titleFontSize)
else:
a.set_title(f"{amp} (BAD)", fontsize=titleFontSize)

f.suptitle(r"Fractional NL residual" + "\n" +
r"$100\times \frac{(k_0 + k_1\times Time - \mu)}{k_0 + k_1\times Time}$",
r"$100\times \frac{(k_0 + k_1*Time-\mu)}{k_0+k_1*Time}$",
fontsize=supTitleFontSize)
pdfPages.savefig()
pdfPages.savefig(f)

return