DM-25114: Update internal fgcm math to always use 64-bit floats

fgcm/_version.py

@@ -1,5 +1,5 @@
 from __future__ import division, absolute_import, print_function
 
-__version__ = '3.2.0'
+__version__ = '3.3.0'
 
 __version_info__ = __version__.split('.')

fgcm/fgcmChisq.py

@@ -590,7 +590,10 @@ def _magWorker(self, goodStarsAndObs):
         qeSysGO = self.fgcmPars.expQESys[obsExpIndexGO]
         filterOffsetGO = self.fgcmPars.expFilterOffset[obsExpIndexGO]
 
-        obsMagGO = obsMagADU[goodObs] + 2.5*np.log10(I0GO) + qeSysGO + filterOffsetGO
+        # Explicitly update obsMagADU to float64 (internally is 32-bit)
+        # I0GO, qeSysGO, filterOffsetGO are 64 bit
+        obsMagGO = obsMagADU[goodObs].astype(np.float64) + \
+            2.5*np.log10(I0GO) + qeSysGO + filterOffsetGO
 
         if (self.fgcmGray is not None):
             # We want to apply the "CCD Gray Crunch"
@@ -619,15 +622,15 @@ def _magWorker(self, goodStarsAndObs):
                 obsMagGO[ok] += ccdGray[obsExpIndexGO[ok], obsCCDIndexGO[ok]]
 
         # Compute the sub-selected error-squared, using model error when available
-        obsMagErr2GO = obsMagADUModelErr[goodObs]**2.
+        obsMagErr2GO = obsMagADUModelErr[goodObs].astype(np.float64)**2.
 
         if (self.computeSEDSlopes):
             # first, compute mean mags (code same as below.  FIXME: consolidate, but how?)
 
             # make temp vars.  With memory overhead
 
-            wtSum = np.zeros_like(objMagStdMean,dtype='f8')
-            objMagStdMeanTemp = np.zeros_like(objMagStdMean)
+            wtSum = np.zeros_like(objMagStdMean, dtype='f8')
+            objMagStdMeanTemp = np.zeros_like(objMagStdMean, dtype='f8')
 
             add_at_2d(wtSum,
                    (obsObjIDIndexGO,obsBandIndexGO),
@@ -691,8 +694,8 @@ def _magWorker(self, goodStarsAndObs):
         #  indexing in the np.add.at() more difficult
 
         wtSum = np.zeros_like(objMagStdMean,dtype='f8')
-        objMagStdMeanTemp = np.zeros_like(objMagStdMean)
-        objMagStdMeanNoChromTemp = np.zeros_like(objMagStdMeanNoChrom)
+        objMagStdMeanTemp = np.zeros_like(objMagStdMean, dtype='f8')
+        objMagStdMeanNoChromTemp = np.zeros_like(objMagStdMeanNoChrom, dtype='f8')
 
         add_at_2d(wtSum,
                (obsObjIDIndexGO,obsBandIndexGO),
@@ -806,7 +809,7 @@ def _chisqWorker(self, goodStarsAndObs):
                                        lutIndicesGO) / I0GO
 
         # Compute the sub-selected error-squared, using model error when available
-        obsMagErr2GO = obsMagADUModelErr[goodObs]**2.
+        obsMagErr2GO = obsMagADUModelErr[goodObs].astype(np.float64)**2.
 
         obsMagStdLock.acquire()
 

fgcm/fgcmConfig.py

@@ -733,7 +733,7 @@ def saveConfigForNextCycle(self,fileName,parFile,flagStarFile):
         # These need to be converted to lists of floats
         for i, b in enumerate(self.bands):
             configDict['expGrayPhotometricCutDict'][b] = float(self.expGrayPhotometricCut[i])
-            configDict['expGrayHighCut'][b] = float(self.expGrayHighCut[i])
+            configDict['expGrayHighCutDict'][b] = float(self.expGrayHighCut[i])
 
         with open(fileName,'w') as f:
             yaml.dump(configDict, stream=f)

fgcm/fgcmFitCycle.py

@@ -106,7 +106,7 @@ def runWithFits(self):
         self._setupWithFits()
         self.run()
 
-    def setStars(self, fgcmStars):
+    def setStars(self, fgcmStars, fgcmPars):
         """
         Record the star information.  This is a separate method to allow
          for memory to be cleared.
@@ -115,10 +115,13 @@ def setStars(self, fgcmStars):
         ----------
         fgcmStars: FgcmStars
            Object with star information
+        fgcmPars: FgcmParameters
+           Parameter information
         """
         # this has to be done outside for memory issues
 
         self.fgcmStars = fgcmStars
+        self.fgcmStars.prepStars(fgcmPars)
 
     def setLUT(self, fgcmLUT):
         """
@@ -170,6 +173,7 @@ def _setupWithFits(self):
         # Read in the stars
         self.fgcmStars = FgcmStars(self.fgcmConfig)
         self.fgcmStars.loadStarsFromFits(self.fgcmPars, computeNobs=True)
+        self.fgcmStars.prepStars(self.fgcmPars)
 
         self.finishSetup()
 
@@ -560,7 +564,7 @@ def run(self):
                 self.fgcmStars.saveStdStars(outStarFile, self.fgcmPars)
 
             # Auto-update photometric cuts
-            for i, b in enumerate(self.bands):
+            for i, b in enumerate(self.fgcmPars.bands):
                 self.fgcmConfig.expGrayPhotometricCutDict[b] = self.updatedPhotometricCut[i]
                 self.fgcmConfig.expGrayHighCutDict[b] = self.updatedHighCut[i]
 

fgcm/fgcmLUT.py

@@ -791,28 +791,28 @@ def __init__(self, indexVals, lutFlat, lutDerivFlat, stdVals, sedLUT=None, filte
         sizeTuple = (len(self.filterNames), self.pwv.size, self.o3.size,
                      self.tau.size, self.alpha.size, self.zenith.size, self.nCCDStep)
 
-        self.lutI0Handle = snmm.createArray(sizeTuple,dtype='f4')
+        self.lutI0Handle = snmm.createArray(sizeTuple,dtype='f8')
         snmm.getArray(self.lutI0Handle)[:, :, :, :, :, :, :] = \
             lutFlat['I0'].reshape(origSizeTuple)[usedLutFilterMark, :, :, :, :, :, :]
 
-        self.lutI1Handle = snmm.createArray(sizeTuple,dtype='f4')
+        self.lutI1Handle = snmm.createArray(sizeTuple,dtype='f8')
         snmm.getArray(self.lutI1Handle)[:, :, :, :, :, :, :] = \
             lutFlat['I1'].reshape(origSizeTuple)[usedLutFilterMark, :, :, :, :, :, :]
 
         # and read in the derivatives
 
         # create shared memory
-        self.lutDLnPwvHandle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDO3Handle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDLnTauHandle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDAlphaHandle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDSecZenithHandle = snmm.createArray(sizeTuple,dtype='f4')
-
-        self.lutDLnPwvI1Handle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDO3I1Handle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDLnTauI1Handle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDAlphaI1Handle = snmm.createArray(sizeTuple,dtype='f4')
-        self.lutDSecZenithI1Handle = snmm.createArray(sizeTuple,dtype='f4')
+        self.lutDLnPwvHandle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDO3Handle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDLnTauHandle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDAlphaHandle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDSecZenithHandle = snmm.createArray(sizeTuple,dtype='f8')
+
+        self.lutDLnPwvI1Handle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDO3I1Handle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDLnTauI1Handle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDAlphaI1Handle = snmm.createArray(sizeTuple,dtype='f8')
+        self.lutDSecZenithI1Handle = snmm.createArray(sizeTuple,dtype='f8')
 
         snmm.getArray(self.lutDLnPwvHandle)[:, :, :, :, :, :, :] = \
             lutDerivFlat['D_LNPWV'].reshape(origSizeTuple)[usedLutFilterMark, :, :, :, :, :, :]

fgcm/fgcmParameters.py

@@ -278,25 +278,25 @@ def _initializeNewParameters(self, expInfo, fgcmLUT):
         self._loadEpochAndWashInfo()
 
         # and make the new parameter arrays
-        self.parAlpha = np.zeros(self.campaignNights.size,dtype=np.float32) + fgcmLUT.alphaStd
-        self.parO3 = np.zeros(self.campaignNights.size,dtype=np.float32) + fgcmLUT.o3Std
-        self.parLnTauIntercept = np.zeros(self.campaignNights.size,dtype=np.float32) + fgcmLUT.lnTauStd
-        self.parLnTauSlope = np.zeros(self.campaignNights.size,dtype=np.float32)
+        self.parAlpha = np.zeros(self.campaignNights.size,dtype=np.float64) + fgcmLUT.alphaStd
+        self.parO3 = np.zeros(self.campaignNights.size,dtype=np.float64) + fgcmLUT.o3Std
+        self.parLnTauIntercept = np.zeros(self.campaignNights.size,dtype=np.float64) + fgcmLUT.lnTauStd
+        self.parLnTauSlope = np.zeros(self.campaignNights.size,dtype=np.float64)
         # these we will always have, won't always fit
-        self.parLnPwvIntercept = np.zeros(self.campaignNights.size, dtype=np.float32) + fgcmLUT.lnPwvStd
-        self.parLnPwvSlope = np.zeros(self.campaignNights.size, dtype=np.float32)
-        self.parLnPwvQuadratic = np.zeros(self.campaignNights.size, dtype=np.float32)
+        self.parLnPwvIntercept = np.zeros(self.campaignNights.size, dtype=np.float64) + fgcmLUT.lnPwvStd
+        self.parLnPwvSlope = np.zeros(self.campaignNights.size, dtype=np.float64)
+        self.parLnPwvQuadratic = np.zeros(self.campaignNights.size, dtype=np.float64)
 
         # parameters with per-epoch values
-        self.parSuperStarFlat = np.zeros((self.nEpochs,self.nLUTFilter,self.nCCD,self.superStarNPar),dtype=np.float64)
+        self.parSuperStarFlat = np.zeros((self.nEpochs,self.nLUTFilter,self.nCCD,self.superStarNPar),dtype=np.float32)
         # The first term should be 1.0 with new flux units
         self.parSuperStarFlat[:, :, :, 0] = 1.0
 
         # parameters with per-wash values
 
         # We always have these parameters, even if we don't fit them
-        self.parQESysIntercept = np.zeros((self.nBands, self.nWashIntervals), dtype=np.float32)
-        self.compQESysSlope = np.zeros((self.nBands, self.nWashIntervals), dtype=np.float32)
+        self.parQESysIntercept = np.zeros((self.nBands, self.nWashIntervals), dtype=np.float64)
+        self.compQESysSlope = np.zeros((self.nBands, self.nWashIntervals), dtype=np.float64)
         self.compQESysSlopeApplied = np.zeros_like(self.compQESysSlope)
 
         # parameters for "absolute" offsets (and relative between filters)
@@ -890,7 +890,7 @@ def parsToArrays(self):
                ('LUTFILTERNAMES', 'a%d' % (maxFilterLen), (len(self.lutFilterNames), )),
                ('BANDS', 'a%d' % (maxBandLen), (len(self.bands), )),
                ('FITBANDS', 'a%d' % (maxBandLen), (len(self.fitBands), )),
-               ('NOTFITBANDS', 'a%d' % (maxBandLen), (len(self.notFitBands), )),
+               # ('NOTFITBANDS', 'a%d' % (maxBandLen), (len(self.notFitBands), )),
                ('LNTAUUNIT', 'f8'),
                ('LNTAUSLOPEUNIT', 'f8'),
                ('ALPHAUNIT', 'f8'),
@@ -915,7 +915,7 @@ def parsToArrays(self):
         parInfo['LUTFILTERNAMES'] = self.lutFilterNames
         parInfo['BANDS'] = self.bands
         parInfo['FITBANDS'] = self.fitBands
-        parInfo['NOTFITBANDS'] = self.notFitBands
+        # parInfo['NOTFITBANDS'] = self.notFitBands
 
         parInfo['HASEXTERNALPWV'] = self.hasExternalPwv
         if (self.hasExternalPwv):
@@ -1052,11 +1052,11 @@ def loadExternalPwv(self, externalPwvDeltaT):
         pwvIndex = np.clip(np.searchsorted(pwvTable['MJD'],self.expMJD),0,pwvTable.size-1)
         # this will be True or False...
         self.externalPwvFlag[:] = (np.abs(pwvTable['MJD'][pwvIndex] - self.expMJD) < externalPwvDeltaT)
-        self.externalLnPwv = np.zeros(self.nExp,dtype=np.float32)
+        self.externalLnPwv = np.zeros(self.nExp,dtype=np.float64)
         self.externalLnPwv[self.externalPwvFlag] = np.log(np.clip(pwvTable['PWV'][pwvIndex[self.externalPwvFlag]], self.pwvRange[0], self.pwvRange[1]))
 
         # and new PWV scaling pars!
-        self.parExternalLnPwvOffset = np.zeros(self.nCampaignNights,dtype=np.float32)
+        self.parExternalLnPwvOffset = np.zeros(self.nCampaignNights,dtype=np.float64)
         self.parExternalLnPwvScale = 1.0
 
         match, = np.where(self.externalPwvFlag)
@@ -1382,8 +1382,8 @@ def getParBounds(self, fitterUnits=False):
         else:
             units = np.ones(self.nFitPars)
 
-        parLow = np.zeros(self.nFitPars,dtype=np.float32)
-        parHigh = np.zeros(self.nFitPars,dtype=np.float32)
+        parLow = np.zeros(self.nFitPars,dtype=np.float64)
+        parHigh = np.zeros(self.nFitPars,dtype=np.float64)
 
         if not self.useRetrievedPwv:
             u = units[self.parLnPwvInterceptLoc: