DM-31880: Update to version 3.9.0 with new reference star plots and cuts.

fgcm/_version.py

@@ -1,3 +1,3 @@
-__version__ = '3.8.5'
+__version__ = '3.9.1'
 
 __version_info__ = __version__.split('.')

fgcm/fgcmChisq.py

@@ -67,6 +67,7 @@ def __init__(self,fgcmConfig,fgcmPars,fgcmStars,fgcmLUT):
         self.ccdGraySubCCD = fgcmConfig.ccdGraySubCCD
         self.useRefStarsWithInstrument = fgcmConfig.useRefStarsWithInstrument
         self.instrumentParsPerBand = fgcmConfig.instrumentParsPerBand
+        self.useExposureReferenceOffset = fgcmConfig.useExposureReferenceOffset
         self.saveParsForDebugging = fgcmConfig.saveParsForDebugging
         self.quietMode = fgcmConfig.quietMode
 
@@ -620,6 +621,11 @@ def _magWorker(self, goodStarsAndObs):
                 # Regular non-sub-ccd
                 obsMagGO[ok] += ccdGray[obsExpIndexGO[ok], obsCCDIndexGO[ok]]
 
+            if self.useExposureReferenceOffset:
+                # Apply the reference offsets as well.
+                ok, = np.where(self.fgcmPars.compExpRefOffset[obsExpIndexGO] > self.illegalValue)
+                obsMagGO[ok] += self.fgcmPars.compExpRefOffset[obsExpIndexGO[ok]]
+
         # Compute the sub-selected error-squared, using model error when available
         obsMagErr2GO = obsMagADUModelErr[goodObs].astype(np.float64)**2.
 
@@ -861,8 +867,9 @@ def _chisqWorker(self, goodStarsAndObs):
                 # Get good observations of reference stars
                 # This must be two steps because we first need the indices to
                 # avoid out-of-bounds
+                mask = objFlagDict['REFSTAR_OUTLIER'] | objFlagDict['REFSTAR_BAD_COLOR']
                 goodRefObsGO, = np.where((objRefIDIndex[obsObjIDIndexGO] >= 0) &
-                                         ((objFlag[obsObjIDIndexGO] & objFlagDict['REFSTAR_OUTLIER']) == 0))
+                                         ((objFlag[obsObjIDIndexGO] & mask) == 0))
 
                 # And check that these are all quality observations
                 tempUse, = np.where((objMagStdMean[obsObjIDIndexGO[goodRefObsGO],

fgcm/fgcmConfig.py

@@ -145,6 +145,7 @@ class FgcmConfig(object):
     sedBoundaryTermDict = ConfigField(dict, required=True)
     sedTermDict = ConfigField(dict, required=True)
     starColorCuts = ConfigField(list, required=True)
+    refStarColorCuts = ConfigField(list, required=False, default=[])
     quantityCuts = ConfigField(list, default=[])
     cycleNumber = ConfigField(int, default=0)
     outfileBase = ConfigField(str, required=True)
@@ -187,6 +188,7 @@ class FgcmConfig(object):
     refStarOutlierNSig = ConfigField(float, default=4.0)
     applyRefStarColorCuts = ConfigField(bool, default=True)
     useRefStarsWithInstrument = ConfigField(bool, default=True)
+    useExposureReferenceOffset = ConfigField(bool, default=False)
 
     mapNSide = ConfigField(int, default=256)
     nStarPerRun = ConfigField(int, default=200000)
@@ -570,15 +572,29 @@ def __init__(self, configDict, lutIndex, lutStd, expInfo, checkFiles=False, noOu
         #  note that self.starColorCuts is a copy so that we don't overwrite.
 
         for cCut in self.starColorCuts:
-            if (not isinstance(cCut[0],int)) :
-                if (cCut[0] not in self.bands):
+            if (not isinstance(cCut[0], int)) :
+                cCut[0] = cCut[0].strip()
+                if cCut[0] not in self.bands:
                     raise ValueError("starColorCut band %s not in list of bands!" % (cCut[0]))
                 cCut[0] = list(self.bands).index(cCut[0])
-            if (not isinstance(cCut[1],int)) :
-                if (cCut[1] not in self.bands):
+            if (not isinstance(cCut[1], int)) :
+                cCut[1] = cCut[1].strip()
+                if cCut[1] not in self.bands:
                     raise ValueError("starColorCut band %s not in list of bands!" % (cCut[1]))
                 cCut[1] = list(self.bands).index(cCut[1])
 
+        for cCut in self.refStarColorCuts:
+            if not isinstance(cCut[0], int):
+                cCut[0] = cCut[0].strip()
+                if cCut[0] not in self.bands:
+                    raise ValueError("refStarColorCut band %s not in list of bands!" % (cCut[0]))
+                cCut[0] = list(self.bands).index(cCut[0])
+            if not isinstance(cCut[1], int):
+                cCut[1] = cCut[1].strip()
+                if cCut[1] not in self.bands:
+                    raise ValueError("refStarColorCut band %s not in list of bands!" % (cCut[1]))
+                cCut[1] = list(self.bands).index(cCut[1])
+
         # Check for input aperture corrections.
         if self.aperCorrFitNBins == 0 and np.any(self.aperCorrInputSlopes == self.illegalValue):
             self.fgcmLog.warning("Aperture corrections will not be fit; strongly recommend setting aperCorrInputSlopeDict")

fgcm/fgcmDeltaAper.py

@@ -116,7 +116,8 @@ def computeDeltaAperExposures(self):
 
             fit, _ = self._fitEpsilonWithDataBinner(mag[ok], deltaAper[ok])
 
-            self.fgcmPars.compEpsilon[expIndex] = self._normalizeEpsilon(fit)
+            if fit is not None:
+                self.fgcmPars.compEpsilon[expIndex] = self._normalizeEpsilon(fit)
 
     def computeDeltaAperStars(self, debug=False):
         """
@@ -197,15 +198,16 @@ def computeDeltaAperStars(self, debug=False):
 
             fit, bin_struct = self._fitEpsilonWithDataBinner(mag_std[r], delta[r])
 
-            globalEpsilon[i] = self._normalizeEpsilon(fit)
-            globalOffset[i] = fit[1]
+            if fit is not None:
+                globalEpsilon[i] = self._normalizeEpsilon(fit)
+                globalOffset[i] = fit[1]
 
             # Store the value of njy_per_arcsec2
             self.fgcmPars.compGlobalEpsilon[i] = globalEpsilon[i]
             self.fgcmLog.info('Global background offset in %s band: %.5f nJy/arcsec2' %
                               (band, globalEpsilon[i]))
 
-            if self.plotPath is not None:
+            if self.plotPath is not None and fit is not None:
                 # Do plots
 
                 st = np.argsort(mag_std[r])
@@ -270,8 +272,9 @@ def computeEpsilonMap(self):
 
                 fit, _ = self._fitEpsilonWithDataBinner(mag_std, delta)
 
-                offsetMap['nstar_fit'][i, j] = len(mag_std)
-                offsetMap['epsilon'][i, j] = self._normalizeEpsilon(fit)
+                if fit is not None:
+                    offsetMap['nstar_fit'][i, j] = len(mag_std)
+                    offsetMap['epsilon'][i, j] = self._normalizeEpsilon(fit)
 
         # Store the offsetmap in njy_per_arcsec2
         self.fgcmPars.compEpsilonMap[:, :] = offsetMap['epsilon']
@@ -601,6 +604,9 @@ def _fitEpsilonWithDataBinner(self, mag, delta_aper, binsize=0.2):
         bin_struct = dataBinner(mag, delta_aper, binsize, [mag_min, mag_max])
         u, = np.where(bin_struct['Y_ERR'] > 0.0)
 
+        if u.size < 5:
+            return None, None
+
         bin_flux = 10.**((bin_struct['X'] - self.njyZp)/(-2.5))
         fit = np.polyfit((2.5/np.log(10.0))/bin_flux[u],
                          bin_struct['Y'][u],

fgcm/fgcmFitCycle.py

@@ -581,6 +581,10 @@ def run(self):
             sigRef = FgcmSigmaRef(self.fgcmConfig, self.fgcmPars, self.fgcmStars)
             sigRef.computeSigmaRef()
 
+            self.fgcmStars.plotRefStarColorTermResiduals(self.fgcmPars)
+
+            self.fgcmGray.computeExposureReferenceOffsets()
+
         # Make Zeropoints
         # We always want to compute these because of the plots
         # In the future we might want to streamline if something is bogging down.
@@ -674,7 +678,7 @@ def _doFit(self, doPlots=True, ignoreRef=False, maxIter=None):
         computeAbsThroughput = self.fgcmStars.hasRefstars
 
         try:
-            fun = optimize.optimize.MemoizeJac(self.fgcmChisq)
+            fun = optimize._optimize.MemoizeJac(self.fgcmChisq)
             jac = fun.derivative
 
             res = optimize.minimize(fun,

fgcm/fgcmGray.py

@@ -11,6 +11,7 @@
 from .fgcmUtilities import histoGauss
 from .fgcmUtilities import Cheb2dField
 from .fgcmUtilities import computeDeltaRA
+from .fgcmUtilities import expFlagDict
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -283,7 +284,8 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             raise ValueError("Must run FgcmChisq to compute magStd before computeCCDAndExpGray")
 
         startTime = time.time()
-        self.fgcmLog.debug('Computing CCDGray and ExpGray.')
+        if not self.quietMode:
+            self.fgcmLog.info('Computing CCDGray and ExpGray.')
 
         # Note: this computes the gray values for all exposures, good and bad
 
@@ -1227,6 +1229,129 @@ def computeCCDAndExpDeltaMagBkg(self):
             self.fgcmLog.info('Computed ccdDeltaBkg and expDeltaBkg in %.2f seconds.' %
                               (time.time() - startTime))
 
+    def computeExposureReferenceOffsets(self):
+        """Compute exposure reference offsets.
+
+        This method computes per-exposure offsets between the calibrated stars and the reference
+        stars for an end-of-calibration "fixup".
+        """
+        if not self.fgcmStars.hasRefstars:
+            # Nothing to do here
+            return
+
+        obsObjIDIndex = snmm.getArray(self.fgcmStars.obsObjIDIndexHandle)
+        obsMagStd = snmm.getArray(self.fgcmStars.obsMagStdHandle)
+        obsBandIndex = snmm.getArray(self.fgcmStars.obsBandIndexHandle)
+        obsExpIndex = snmm.getArray(self.fgcmStars.obsExpIndexHandle)
+        obsCCDIndex = snmm.getArray(self.fgcmStars.obsCCDHandle) - self.ccdStartIndex
+
+        objRefIDIndex = snmm.getArray(self.fgcmStars.objRefIDIndexHandle)
+        refMag = snmm.getArray(self.fgcmStars.refMagHandle)
+
+        goodStars = self.fgcmStars.getGoodStarIndices(checkMinObs=True, removeRefstarOutliers=True, removeRefstarBadcols=True)
+        _, goodObs = self.fgcmStars.getGoodObsIndices(goodStars, expFlag=self.fgcmPars.expFlag, checkBadMag=True)
+
+        # Add in the gray values
+        obsExpIndexGO = obsExpIndex[goodObs]
+        obsCCDIndexGO = obsCCDIndex[goodObs]
+        obsBandIndexGO = obsBandIndex[goodObs]
+        obsMagStdGO = obsMagStd[goodObs]
+
+        ccdGray = snmm.getArray(self.ccdGrayHandle)
+        if np.any(self.ccdGraySubCCD):
+            ccdGraySubCCDPars = snmm.getArray(self.ccdGraySubCCDParsHandle)
+
+        ok, = np.where(ccdGray[obsExpIndexGO, obsCCDIndexGO] > self.illegalValue)
+
+        if np.any(self.ccdGraySubCCD):
+            obsXGO = snmm.getArray(self.fgcmStars.obsXHandle)[goodObs]
+            obsYGO = snmm.getArray(self.fgcmStars.obsYHandle)[goodObs]
+            expCcdHash = (obsExpIndexGO[ok] * (self.fgcmPars.nCCD + 1) +
+                          obsCCDIndexGO[ok])
+            h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+            use, = np.where(h > 0)
+            for i in use:
+                i1a = rev[rev[i]: rev[i + 1]]
+                eInd = obsExpIndexGO[ok[i1a[0]]]
+                cInd = obsCCDIndexGO[ok[i1a[0]]]
+                field = Cheb2dField(self.deltaMapperDefault['x_size'][cInd],
+                                    self.deltaMapperDefault['y_size'][cInd],
+                                    ccdGraySubCCDPars[eInd, cInd, :])
+                fluxScale = field.evaluate(obsXGO[ok[i1a]], obsYGO[ok[i1a]])
+                obsMagStdGO[ok[i1a]] += -2.5 * np.log10(np.clip(fluxScale, 0.1, None))
+        else:
+            # Regular non-sub-ccd
+            obsMagStdGO[ok] += ccdGray[obsExpIndexGO[ok], obsCCDIndexGO[ok]]
+
+        goodRefObsGO, = np.where(objRefIDIndex[obsObjIDIndex[goodObs]] >= 0)
+
+        obsUse, = np.where((obsMagStd[goodObs[goodRefObsGO]] < 90.0) &
+                           (refMag[objRefIDIndex[obsObjIDIndex[goodObs[goodRefObsGO]]],
+                                   obsBandIndex[goodObs[goodRefObsGO]]] < 90.0))
+
+        goodRefObsGO = goodRefObsGO[obsUse]
+
+        EGrayGRO = (refMag[objRefIDIndex[obsObjIDIndex[goodObs[goodRefObsGO]]],
+                           obsBandIndex[goodObs[goodRefObsGO]]] -
+                    obsMagStdGO[goodRefObsGO])
+
+        # And then this can be split per exposure.
+
+        h, rev = esutil.stat.histogram(obsExpIndexGO[goodRefObsGO], rev=True)
+
+        use, = np.where(h >= self.minStarPerExp)
+
+        self.fgcmPars.compExpRefOffset[:] = self.illegalValue
+        for i in use:
+            i1a = rev[rev[i]: rev[i + 1]]
+
+            eInd = obsExpIndexGO[goodRefObsGO[i1a[0]]]
+            bInd = obsBandIndexGO[goodRefObsGO[i1a[0]]]
+
+            self.fgcmPars.compExpRefOffset[eInd] = np.median(EGrayGRO[i1a])
+
+        # Do plots if necessary.
+        if self.plotPath is None:
+            return
+
+        rejectMask = (expFlagDict['TOO_FEW_STARS'] |
+                      expFlagDict['NO_STARS'] |
+                      expFlagDict['BAND_NOT_IN_LUT'])
+
+        expUse, = np.where((self.fgcmPars.expFlag & rejectMask) == 0)
+
+        for i in range(self.fgcmPars.nBands):
+            inBand, = np.where((self.fgcmPars.expBandIndex[expUse] == i) &
+                               (self.fgcmPars.compExpRefOffset[expUse] > self.illegalValue))
+
+            if inBand.size == 0:
+                continue
+
+            fig = plt.figure(1, figsize=(8, 6))
+            fig.clf()
+            ax = fig.add_subplot(111)
+
+            coeff = histoGauss(ax, self.fgcmPars.compExpRefOffset[expUse[inBand]]*1000.0)
+            coeff[1] /= 1000.0
+            coeff[2] /= 1000.0
+
+            ax.tick_params(axis='both', which='major', labelsize=14)
+            ax.locator_params(axis='x', nbins=5)
+
+            text=r'$(%s)$' % (self.fgcmPars.bands[i]) + '\n' + \
+                r'$\mathrm{Cycle\ %d}$' % (self.cycleNumber) + '\n' + \
+                r'$\mu = %.2f$' % (coeff[1]*1000.0) + '\n' + \
+                r'$\sigma = %.2f$' % (coeff[2]*1000.0)
+
+            ax.annotate(text, (0.95, 0.93), xycoords='axes fraction', ha='right', va='top', fontsize=16)
+            ax.set_xlabel(r'$\mathrm{EXP}^{\mathrm{ref}}\,(\mathrm{mmag})$', fontsize=16)
+            ax.set_ylabel(r'# of Exposures',fontsize=14)
+
+            fig.savefig('%s/%s_expref_%s.png' % (self.plotPath,
+                                                 self.outfileBaseWithCycle,
+                                                 self.fgcmPars.bands[i]))
+            plt.close(fig)
+
     def __getstate__(self):
         # Don't try to pickle the logger.
 

fgcm/fgcmParameters.py

@@ -387,6 +387,9 @@ def _initializeNewParameters(self, expInfo, fgcmLUT):
         self.compVarGray = np.zeros(self.nExp,dtype='f8')
         self.compNGoodStarPerExp = np.zeros(self.nExp,dtype='i4')
 
+        # Per-Exposure Reference offset
+        self.compExpRefOffset = np.zeros(self.nExp, dtype='f8')
+
         # We also have the median delta aperture and epsilon per exposure
         self.compMedDeltaAper = np.zeros(self.nExp, dtype='f8')
         self.compEpsilon = np.zeros(self.nExp, dtype='f8')
@@ -504,6 +507,11 @@ def _loadOldParameters(self, expInfo, inParInfo, inParams, inSuperStar):
         self.compVarGray = np.atleast_1d(inParams['COMPVARGRAY'][0])
         self.compNGoodStarPerExp = np.atleast_1d(inParams['COMPNGOODSTARPEREXP'][0])
 
+        try:
+            self.compExpRefOffset = np.atleast_1d(inParams['COMPEXPREFOFFSET'][0])
+        except:
+            self.compExpRefOffset = np.zeros(self.nExp, dtype='f8')
+
         npix = hp.nside2npix(self.deltaAperFitSpatialNside)
         try:
             self.compMedDeltaAper = np.atleast_1d(inParams['COMPMEDDELTAAPER'][0])
@@ -1010,6 +1018,7 @@ def parsToArrays(self):
                ('COMPMODELERRPARS', 'f8', (self.compModelErrPars.size, )),
                ('COMPEXPGRAY', 'f8', (self.compExpGray.size, )),
                ('COMPVARGRAY', 'f8', (self.compVarGray.size, )),
+               ('COMPEXPREFOFFSET', 'f8', (self.compExpRefOffset.size, )),
                ('COMPMEDDELTAAPER', 'f8', (self.compMedDeltaAper.size, )),
                ('COMPEPSILON', 'f8', (self.compEpsilon.size, )),
                ('COMPGLOBALEPSILON', 'f4', (self.compGlobalEpsilon.size, )),
@@ -1081,6 +1090,8 @@ def parsToArrays(self):
         pars['COMPVARGRAY'][:] = self.compVarGray
         pars['COMPNGOODSTARPEREXP'][:] = self.compNGoodStarPerExp
 
+        pars['COMPEXPREFOFFSET'][:] = self.compExpRefOffset
+
         pars['COMPMEDDELTAAPER'][:] = self.compMedDeltaAper
         pars['COMPEPSILON'][:] = self.compEpsilon
         pars['COMPGLOBALEPSILON'][:] = self.compGlobalEpsilon

fgcm/fgcmSigmaRef.py

@@ -66,8 +66,10 @@ def computeSigmaRef(self):
 
         # Select only stars that have reference magnitudes
         # and that are not flagged as outliers
+        # We allow all the reference colors in to get color range statistics.
+        mask = objFlagDict['REFSTAR_OUTLIER']
         use, = np.where((objRefIDIndex[goodStars] >= 0) &
-                        ((objFlag[goodStars] & objFlagDict['REFSTAR_OUTLIER']) == 0))
+                        ((objFlag[goodStars] & mask) == 0))
         goodRefStars = goodStars[use]
 
         # We need to have a branch of "small-number" and "large number" of reference stars