Merge pull request #28 from erykoff/cterm_plots

Add reference star color term residual plots.

fgcm/_version.py

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

fgcm/fgcmFitCycle.py

@@ -581,6 +581,8 @@ def run(self):
             sigRef = FgcmSigmaRef(self.fgcmConfig, self.fgcmPars, self.fgcmStars)
             sigRef.computeSigmaRef()
 
+            self.fgcmStars.plotRefStarColorTermResiduals(self.fgcmPars)
+
         # 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 +676,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/fgcmStars.py

@@ -89,6 +89,7 @@ def __init__(self,fgcmConfig):
         self.bandNotFitIndex = fgcmConfig.bandNotFitIndex
         self.lutFilterNames = fgcmConfig.lutFilterNames
         self.filterToBand = fgcmConfig.filterToBand
+        self.colorSplitBands = fgcmConfig.colorSplitBands
         self.colorSplitIndices = fgcmConfig.colorSplitIndices
 
         self.superStarSubCCD = fgcmConfig.superStarSubCCD
@@ -2013,6 +2014,73 @@ def retrieveStdStarCatalog(self, fgcmPars):
 
         return outCat, goodBandNames
 
+    def plotRefStarColorTermResiduals(self, fgcmPars):
+        """
+        Plot reference star color-term residuals.
+
+        Parameters
+        ----------
+        fgcmPars : `fgcm.FgcmParameters`
+        """
+        if not self.hasRefstars:
+            self.fgcmLog.info("No reference stars for color term residual plots.")
+
+        if self.plotPath is None:
+            return
+
+        objMagStdMean = snmm.getArray(self.objMagStdMeanHandle)
+        objNGoodObs = snmm.getArray(self.objNGoodObsHandle)
+        objFlag = snmm.getArray(self.objFlagHandle)
+
+        objRefIDIndex = snmm.getArray(self.objRefIDIndexHandle)
+        refMag = snmm.getArray(self.refMagHandle)
+
+        goodStars = self.getGoodStarIndices(includeReserve=True, checkMinObs=True)
+
+        # Select only stars that have reference magnitudes
+        use, = np.where(objRefIDIndex[goodStars] >= 0)
+        goodRefStars = goodStars[use]
+
+        # Compute "g-i" based on the configured colorSplitIndices
+        gmiGRS = (objMagStdMean[goodRefStars, self.colorSplitIndices[0]] -
+                  objMagStdMean[goodRefStars, self.colorSplitIndices[1]])
+
+        okColor, = np.where((objMagStdMean[goodRefStars, self.colorSplitIndices[0]] < 90.0) &
+                            (objMagStdMean[goodRefStars, self.colorSplitIndices[1]] < 90.0))
+
+        for bandIndex, band in enumerate(self.bands):
+            if not fgcmPars.hasExposuresInBand[bandIndex]:
+                continue
+
+            fig = plt.figure(figsize=(8, 6))
+            fig.clf()
+            ax = fig.add_subplot(111)
+
+            refUse, = np.where((refMag[objRefIDIndex[goodRefStars[okColor]], bandIndex] < 90.0) &
+                               (objMagStdMean[goodRefStars[okColor], bandIndex] < 90.0))
+            refUse = okColor[refUse]
+
+            delta = (objMagStdMean[goodRefStars[refUse], bandIndex] -
+                     refMag[objRefIDIndex[goodRefStars[refUse]], bandIndex])
+
+            st = np.argsort(delta)
+            ylow = delta[st[int(0.02*refUse.size)]]
+            yhigh = delta[st[int(0.98*refUse.size)]]
+            st = np.argsort(gmiGRS[refUse])
+            xlow = gmiGRS[refUse[st[int(0.02*refUse.size)]]]
+            xhigh = gmiGRS[refUse[st[int(0.98*refUse.size)]]]
+
+            ax.hexbin(gmiGRS[refUse], delta, bins='log', extent=[xlow, xhigh, ylow, yhigh])
+            ax.set_title('%s band' % (band))
+            ax.set_xlabel('%s - %s' % (self.colorSplitBands[0], self.colorSplitBands[1]))
+            ax.set_ylabel('%s_std - %s_ref' % (band, band))
+
+            fig.tight_layout()
+            fig.savefig('%s/%s_refresidvscol_%s.png' % (self.plotPath,
+                                                        self.outfileBaseWithCycle,
+                                                        band))
+            plt.close(fig)
+
     def __getstate__(self):
         # Don't try to pickle the logger.