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-36944: Update to fgcm 3.9.6 for look-up table floating point fix. #33

Merged
merged 4 commits into from
Nov 9, 2022
Merged

DM-36944: Update to fgcm 3.9.6 for look-up table floating point fix. #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d10776c
Update look-up table generation code to use float64 internally.
erykoff Nov 9, 2022
41aebc7
Update version.
erykoff Nov 9, 2022
af8ef35
Merge pull request #35 from erykoff/DM-36944_lutfloats
erykoff Nov 9, 2022
a0b1f14
Merge branch 'master' into tickets/DM-36944
erykoff Nov 9, 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
2 changes: 1 addition & 1 deletion fgcm/_version.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
__version__ = '3.9.5'
__version__ = '3.9.6'

__version_info__ = __version__.split('.')
54 changes: 27 additions & 27 deletions fgcm/fgcmAtmosphereTable.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,33 +246,33 @@ def loadTable(self):
parStruct = pyfits.getdata(self.atmosphereTableFile, ext=('PARS', 1))


self.elevation = parStruct['ELEVATION'][0]
self.pmbElevation = parStruct['PMBELEVATION'][0]
self.pmbStd = parStruct['PMBSTD'][0]
self.pwvStd = parStruct['PWVSTD'][0]
self.o3Std = parStruct['O3STD'][0]
self.tauStd = parStruct['TAUSTD'][0]
self.alphaStd = parStruct['ALPHASTD'][0]
self.secZenithStd = parStruct['AIRMASSSTD'][0]
self.lambdaRange = parStruct['LAMBDARANGE'][0]
self.lambdaStep = parStruct['LAMBDASTEP'][0]
self.lambdaNorm = parStruct['LAMBDANORM'][0]

self.atmLambda = parStruct['ATMLAMBDA'][0]
self.atmStdTrans = parStruct['ATMSTDTRANS'][0]

self.pmb = parStruct['PMB'][0]
self.pmbDelta = parStruct['PMBDELTA'][0]
self.lnPwv = parStruct['LNPWV'][0]
self.lnPwvDelta = parStruct['LNPWVDELTA'][0]
self.o3 = parStruct['O3'][0]
self.o3Delta = parStruct['O3DELTA'][0]
self.lnTau = parStruct['LNTAU'][0]
self.lnTauDelta = parStruct['LNTAUDELTA'][0]
self.alpha = parStruct['ALPHA'][0]
self.alphaDelta = parStruct['ALPHADELTA'][0]
self.secZenith = parStruct['SECZENITH'][0]
self.secZenithDelta = parStruct['SECZENITHDELTA'][0]
self.elevation = parStruct['ELEVATION'][0].astype(np.float64)
self.pmbElevation = parStruct['PMBELEVATION'][0].astype(np.float64)
self.pmbStd = parStruct['PMBSTD'][0].astype(np.float64)
self.pwvStd = parStruct['PWVSTD'][0].astype(np.float64)
self.o3Std = parStruct['O3STD'][0].astype(np.float64)
self.tauStd = parStruct['TAUSTD'][0].astype(np.float64)
self.alphaStd = parStruct['ALPHASTD'][0].astype(np.float64)
self.secZenithStd = parStruct['AIRMASSSTD'][0].astype(np.float64)
self.lambdaRange = parStruct['LAMBDARANGE'][0].astype(np.float64)
self.lambdaStep = parStruct['LAMBDASTEP'][0].astype(np.float64)
self.lambdaNorm = parStruct['LAMBDANORM'][0].astype(np.float64)

self.atmLambda = parStruct['ATMLAMBDA'][0].astype(np.float64)
self.atmStdTrans = parStruct['ATMSTDTRANS'][0].astype(np.float64)

self.pmb = parStruct['PMB'][0].astype(np.float64)
self.pmbDelta = parStruct['PMBDELTA'][0].astype(np.float64)
self.lnPwv = parStruct['LNPWV'][0].astype(np.float64)
self.lnPwvDelta = parStruct['LNPWVDELTA'][0].astype(np.float64)
self.o3 = parStruct['O3'][0].astype(np.float64)
self.o3Delta = parStruct['O3DELTA'][0].astype(np.float64)
self.lnTau = parStruct['LNTAU'][0].astype(np.float64)
self.lnTauDelta = parStruct['LNTAUDELTA'][0].astype(np.float64)
self.alpha = parStruct['ALPHA'][0].astype(np.float64)
self.alphaDelta = parStruct['ALPHADELTA'][0].astype(np.float64)
self.secZenith = parStruct['SECZENITH'][0].astype(np.float64)
self.secZenithDelta = parStruct['SECZENITHDELTA'][0].astype(np.float64)

if fits_package == 'fitsio':
self.pwvAtmTable = fitsio.read(self.atmosphereTableFile, ext='PWVATM')
Expand Down
37 changes: 18 additions & 19 deletions fgcm/fgcmLUT.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,9 +162,9 @@ def setThroughputs(self, throughputDict):
except:
raise ValueError("Wavelength LAMBDA not found for filter %s in throughputDict!" % (filterName))

tput = np.zeros(lam.size, dtype=[('LAMBDA','f4'),
('THROUGHPUT_AVG','f4'),
('THROUGHPUT_CCD','f4',self.nCCD)])
tput = np.zeros(lam.size, dtype=[('LAMBDA', 'f8'),
('THROUGHPUT_AVG', 'f8'),
('THROUGHPUT_CCD', 'f8', self.nCCD)])
tput['LAMBDA'][:] = lam
for ccdIndex in range(self.nCCD):
try:
Expand Down Expand Up @@ -282,7 +282,7 @@ def makeLUT(self):
self.atmStdTrans = self.atmosphereTable.atmStdTrans

self.pmbElevation = self.atmosphereTable.pmbElevation

erykoff marked this conversation as resolved.
Show resolved Hide resolved
self.pmb = self.atmosphereTable.pmb
self.pmbDelta = self.atmosphereTable.pmbDelta
pmbPlus = np.append(self.pmb, self.pmb[-1] + self.pmbDelta)
Expand Down Expand Up @@ -328,9 +328,9 @@ def makeLUT(self):
for i in range(len(self.filterNames)):
inLam = self.inThroughputs[i]['LAMBDA']

tput = np.zeros(self.atmLambda.size, dtype=[('LAMBDA','f4'),
('THROUGHPUT_AVG','f4'),
('THROUGHPUT_CCD','f4',self.nCCD)])
tput = np.zeros(self.atmLambda.size, dtype=[('LAMBDA', 'f8'),
('THROUGHPUT_AVG', 'f8'),
('THROUGHPUT_CCD', 'f8', self.nCCD)])
tput['LAMBDA'][:] = self.atmLambda

for ccdIndex in range(self.nCCD):
Expand All @@ -347,16 +347,16 @@ def makeLUT(self):
self.fgcmLog.info("Computing lambdaB")
self.lambdaB = np.zeros(len(self.filterNames))
for i in range(len(self.filterNames)):
num = integrate.simps(self.atmLambda * self.throughputs[i]['THROUGHPUT_AVG'] / self.atmLambda, self.atmLambda)
denom = integrate.simps(self.throughputs[i]['THROUGHPUT_AVG'] / self.atmLambda, self.atmLambda)
num = integrate.simpson(self.atmLambda * self.throughputs[i]['THROUGHPUT_AVG'] / self.atmLambda, self.atmLambda)
denom = integrate.simpson(self.throughputs[i]['THROUGHPUT_AVG'] / self.atmLambda, self.atmLambda)
self.lambdaB[i] = num / denom
self.fgcmLog.info("Filter: %s, lambdaB = %.3f" % (self.filterNames[i], self.lambdaB[i]))

self.fgcmLog.info("Computing lambdaStdFilter")
self.lambdaStdFilter = np.zeros(len(self.filterNames))
for i in range(len(self.filterNames)):
num = integrate.simps(self.atmLambda * self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
denom = integrate.simps(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
num = integrate.simpson(self.atmLambda * self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
denom = integrate.simpson(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
self.lambdaStdFilter[i] = num / denom
self.fgcmLog.info("Filter: %s, lambdaStdFilter = %.3f" % (self.filterNames[i],self.lambdaStdFilter[i]))

Expand All @@ -378,9 +378,9 @@ def makeLUT(self):
self.I2Std = np.zeros(len(self.filterNames))

for i in range(len(self.filterNames)):
self.I0Std[i] = integrate.simps(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
self.I1Std[i] = integrate.simps(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans * (self.atmLambda - self.lambdaStd[i]) / self.atmLambda, self.atmLambda)
self.I2Std[i] = integrate.simps(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans * (self.atmLambda - self.lambdaStd[i])**2. / self.atmLambda, self.atmLambda)
self.I0Std[i] = integrate.simpson(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans / self.atmLambda, self.atmLambda)
self.I1Std[i] = integrate.simpson(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans * (self.atmLambda - self.lambdaStd[i]) / self.atmLambda, self.atmLambda)
self.I2Std[i] = integrate.simpson(self.throughputs[i]['THROUGHPUT_AVG'] * self.atmStdTrans * (self.atmLambda - self.lambdaStd[i])**2. / self.atmLambda, self.atmLambda)

self.I10Std = self.I1Std / self.I0Std

Expand Down Expand Up @@ -422,8 +422,8 @@ def makeLUT(self):
else:
Sb = self.throughputs[i]['THROUGHPUT_CCD'][:,p] * o2AtmTable[o,:] * rayleighAtmTable[o,:] * pwvAtmTable[j,o,:] * o3AtmTable[k,o,:] * aerosolTauLambda

lutPlus['I0'][i,j,k,m,n,o,p] = integrate.simps(Sb / self.atmLambda, self.atmLambda)
lutPlus['I1'][i,j,k,m,n,o,p] = integrate.simps(Sb * (self.atmLambda - self.lambdaStd[i]) / self.atmLambda, self.atmLambda)
lutPlus['I0'][i,j,k,m,n,o,p] = integrate.simpson(Sb / self.atmLambda, self.atmLambda)
lutPlus['I1'][i,j,k,m,n,o,p] = integrate.simpson(Sb * (self.atmLambda - self.lambdaStd[i]) / self.atmLambda, self.atmLambda)

# and create the LUT (not plus)
self.lut = np.zeros((len(self.filterNames),
Expand Down Expand Up @@ -540,7 +540,6 @@ def makeLUT(self):
self.secZenithDelta)
)


if (self.makeSeds):
# and the SED LUT
self.fgcmLog.info("Building SED LUT")
Expand Down Expand Up @@ -596,8 +595,8 @@ def makeLUT(self):

# compute synthetic mags
for j in range(len(self.filterNames)):
num = integrate.simps(fnu * self.throughputs[j]['THROUGHPUT_AVG'][:] * self.atmStdTrans / self.atmLambda, self.atmLambda)
denom = integrate.simps(self.throughputs[j]['THROUGHPUT_AVG'][:] * self.atmStdTrans / self.atmLambda, self.atmLambda)
num = integrate.simpson(fnu * self.throughputs[j]['THROUGHPUT_AVG'][:] * self.atmStdTrans / self.atmLambda, self.atmLambda)
denom = integrate.simpson(self.throughputs[j]['THROUGHPUT_AVG'][:] * self.atmStdTrans / self.atmLambda, self.atmLambda)

self.sedLUT['SYNTHMAG'][i,j] = -2.5*np.log10(num/denom)

Expand Down