Merge pull request #39 from lsst/tickets/DM-42285

DM-42285: Update to support numpy 1.25+ and various optimizations.

fgcm/_version.py

@@ -1,3 +1,3 @@
-__version__ = '3.10.3'
+__version__ = '3.11.0'
 
 __version_info__ = __version__.split('.')

fgcm/fgcmBrightObs.py

@@ -217,7 +217,7 @@ def _worker(self,goodStarsAndObs):
         # find the brightest (minmag) object at each index
         np.fmin.at(objMagStdMeanTemp,
                    (obsObjIDIndexGO, obsBandIndexGO),
-                   obsMagStdGO)
+                   obsMagStdGO.astype(objMagStdMeanTemp.dtype))
 
         # now which observations are bright *enough* to consider?
         brightEnoughGO, = np.where((obsMagStdGO -
@@ -236,11 +236,11 @@ def _worker(self,goodStarsAndObs):
         np.add.at(wtSum,
                   (obsObjIDIndexGO[brightEnoughGO],
                    obsBandIndexGO[brightEnoughGO]),
-                  1./obsMagErr2GOBE)
+                  1./obsMagErr2GOBE.astype(wtSum.dtype))
         np.add.at(objMagStdMeanTemp,
                   (obsObjIDIndexGO[brightEnoughGO],
                    obsBandIndexGO[brightEnoughGO]),
-                  obsMagStdGO[brightEnoughGO]/obsMagErr2GOBE)
+                  (obsMagStdGO[brightEnoughGO]/obsMagErr2GOBE).astype(objMagStdMeanTemp.dtype))
         np.add.at(objNGoodObsTemp,
                   (obsObjIDIndexGO[brightEnoughGO],
                    obsBandIndexGO[brightEnoughGO]),

fgcm/fgcmChromaticity.py

@@ -5,6 +5,7 @@
 import time
 import matplotlib.pyplot as plt
 import scipy.optimize as optimize
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -226,10 +227,10 @@ def __call__(self, fitPars, returnExpGray=False):
 
             np.add.at(expGrayColorSplit[:, k],
                       self.obsExpIndexGO[sel],
-                      EGrayGOS / self.EGrayErr2GO[sel])
+                      (EGrayGOS / self.EGrayErr2GO[sel]).astype(expGrayColorSplit.dtype))
             np.add.at(expGrayColorSplitWt[:, k],
                       self.obsExpIndexGO[sel],
-                      1. / self.EGrayErr2GO[sel])
+                      (1. / self.EGrayErr2GO[sel]).astype(expGrayColorSplitWt.dtype))
 
             ok, = np.where(expGrayColorSplitWt[:, k] > 0.0)
             expGrayColorSplit[ok, k] /= expGrayColorSplitWt[ok, k]
@@ -240,14 +241,14 @@ def __call__(self, fitPars, returnExpGray=False):
 
             np.add.at(dExpGraydC0ColorSplit[:, :, k],
                       (self.obsExpIndexGO[sel], self.fgcmPars.expCoatingIndex[self.obsExpIndexGO[sel]]),
-                      dDeltadC0GOS / self.EGrayErr2GO[sel])
+                      (dDeltadC0GOS / self.EGrayErr2GO[sel]).astype(dExpGraydC0ColorSplit.dtype))
 
             for i in range(c0s.size):
                 dExpGraydC0ColorSplit[ok, i, k] /= expGrayColorSplitWt[ok, k]
 
             np.add.at(dExpGraydC1ColorSplit[:, k],
                       self.obsExpIndexGO[sel],
-                      dDeltadC0GOS * self.deltaTGO[sel] / self.EGrayErr2GO[sel])
+                      (dDeltadC0GOS * self.deltaTGO[sel] / self.EGrayErr2GO[sel]).astype(dExpGraydC1ColorSplit.dtype))
 
             dExpGraydC1ColorSplit[ok, k] /= expGrayColorSplitWt[ok, k]
 
@@ -266,7 +267,7 @@ def __call__(self, fitPars, returnExpGray=False):
                   (2.0 * (expGrayColorSplitWt[ok, 0] + expGrayColorSplitWt[ok, 1]) *
                    (expGrayColorSplit[ok, 0] - expGrayColorSplit[ok, 1]) *
                    (dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 0] -
-                    dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 1])))
+                    dExpGraydC0ColorSplit[ok, self.fgcmPars.expCoatingIndex[ok], 1])).astype(deriv.dtype))
 
         deriv /= (ok.size - len(fitPars))
 
@@ -354,7 +355,7 @@ def computeCCDChromaticity(self):
         ccdFilterHash = (obsLUTFilterIndexGO.astype(np.int64)*(self.fgcmPars.nCCD + 1) +
                          obsCCDIndexGO.astype(np.int64))
 
-        h, rev = esutil.stat.histogram(ccdFilterHash, rev=True)
+        h, rev = histogram_rev_sorted(ccdFilterHash)
 
         # Make a simple cut here.
         use, = np.where(h >= 10)
@@ -382,7 +383,7 @@ def computeCCDChromaticity(self):
                 self.fgcmLog.warning("Found out-of-bounds value for chromaticity for filter %s, detector %d." % (self.fgcmPars.lutFilterNames[self.fInd], self.cInd + self.ccdStartIndex))
                 continue
 
-            self.fgcmPars.compCCDChromaticity[self.cInd, self.fInd] = res.x
+            self.fgcmPars.compCCDChromaticity[self.cInd, self.fInd] = res.x[0]
             cWasFit[self.cInd, self.fInd] = True
 
         # And make plots if necessary.

fgcm/fgcmConnectivity.py

@@ -4,6 +4,7 @@
 import matplotlib.pyplot as plt
 
 from .fgcmUtilities import objFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
 class FgcmConnectivity(object):
@@ -76,8 +77,11 @@ def plotConnectivity(self):
                                 ((objFlag[obsObjIDIndex] & mask) == 0))
 
             # Split into nights with a histogram
-            h, rev = esutil.stat.histogram(self.fgcmPars.expNightIndex[obsExpIndex[goodObs]],
-                                           min=0, max=self.fgcmPars.nCampaignNights, rev=True)
+            h, rev = histogram_rev_sorted(
+                self.fgcmPars.expNightIndex[obsExpIndex[goodObs]],
+                min=0,
+                max=self.fgcmPars.nCampaignNights,
+            )
 
             done = False
             while not done:

fgcm/fgcmDeltaAper.py

@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 
 from .fgcmUtilities import dataBinner
+from .fgcmUtilities import histogram_rev_sorted
 
 import multiprocessing
 
@@ -89,7 +90,7 @@ def computeDeltaAperExposures(self):
         self.fgcmPars.compMedDeltaAper[:] = self.illegalValue
         self.fgcmPars.compEpsilon[:] = self.illegalValue
 
-        h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True, min=0)
+        h, rev = histogram_rev_sorted(obsExpIndex[goodObs], min=0)
         expIndices, = np.where(h >= self.minStarPerExp)
 
         for expIndex in expIndices:
@@ -248,7 +249,11 @@ def computeEpsilonMap(self):
         # Do the mapping
         self.fgcmLog.info("Computing delta-aper epsilon spatial map at nside %d" % self.deltaAperFitSpatialNside)
         ipring = hpg.angle_to_pixel(self.deltaAperFitSpatialNside, objRA, objDec, nest=False)
-        h, rev = esutil.stat.histogram(ipring, min=0, max=hpg.nside_to_npixel(self.deltaAperFitSpatialNside) - 1, rev=True)
+        h, rev = histogram_rev_sorted(
+            ipring,
+            min=0,
+            max=hpg.nside_to_npixel(self.deltaAperFitSpatialNside) - 1,
+        )
 
         offsetMap = np.zeros(h.size, dtype=[('nstar_fit', 'i4', (self.fgcmStars.nBands, )),
                                             ('epsilon', 'f4', (self.fgcmStars.nBands, ))])
@@ -365,7 +370,7 @@ def computeEpsilonPerCcd(self):
 
         filterCcdHash = ccdIndexGO*(self.fgcmPars.nLUTFilter + 1) + lutFilterIndexGO
 
-        h, rev = esutil.stat.histogram(filterCcdHash, rev=True)
+        h, rev = histogram_rev_sorted(filterCcdHash)
 
         # Arbitrary minimum number here
         gdHash, = np.where(h > 10)
@@ -392,7 +397,7 @@ def computeEpsilonPerCcd(self):
 
             xyBinHash = xBin[i1a]*(self.deltaAperFitPerCcdNy + 1) + yBin[i1a]
 
-            h2, rev2 = esutil.stat.histogram(xyBinHash, rev=True)
+            h2, rev2 = histogram_rev_sorted(xyBinHash)
 
             gdHash2, = np.where(h2 > 10)
             for j in gdHash2:
@@ -517,10 +522,10 @@ def _starWorker(self, goodStarsAndObs):
 
         np.add.at(objDeltaAperMeanTemp,
                   (obsObjIDIndex[goodObs], obsBandIndex[goodObs]),
-                  (obsDeltaAper[goodObs] - self.fgcmPars.compMedDeltaAper[obsExpIndex[goodObs]])/obsMagErr2GO)
+                  ((obsDeltaAper[goodObs] - self.fgcmPars.compMedDeltaAper[obsExpIndex[goodObs]])/obsMagErr2GO).astype(objDeltaAperMeanTemp.dtype))
         np.add.at(wtSum,
                   (obsObjIDIndex[goodObs], obsBandIndex[goodObs]),
-                  1./obsMagErr2GO)
+                  (1./obsMagErr2GO).astype(wtSum.dtype))
 
         gd = np.where(wtSum > 0.0)
 

fgcm/fgcmExposureSelector.py

@@ -4,6 +4,7 @@
 import esutil
 
 from .fgcmUtilities import expFlagDict, logFlaggedExposuresPerBand, checkFlaggedExposuresPerBand
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -143,8 +144,11 @@ def selectCalibratableNights(self):
         self.fgcmLog.info('Flagging exposures on %d bad nights with too few photometric exposures.' % (badNights.size))
 
         # and we need to use *all* the exposures to flag bad nights
-        h,rev = esutil.stat.histogram(self.fgcmPars.expNightIndex,min=0,
-                                      max=self.fgcmPars.nCampaignNights-1,rev=True)
+        h, rev = histogram_rev_sorted(
+            self.fgcmPars.expNightIndex,
+            min=0,
+            max=self.fgcmPars.nCampaignNights - 1,
+        )
 
         for badNight in badNights:
             i1a=rev[rev[badNight]:rev[badNight+1]]

fgcm/fgcmGray.py

@@ -12,6 +12,7 @@
 from .fgcmUtilities import Cheb2dField
 from .fgcmUtilities import computeDeltaRA
 from .fgcmUtilities import expFlagDict
+from .fgcmUtilities import histogram_rev_sorted
 
 from .sharedNumpyMemManager import SharedNumpyMemManager as snmm
 
@@ -210,10 +211,10 @@ def computeExpGrayForInitialSelection(self):
 
         np.add.at(expGrayForInitialSelection,
                   obsExpIndex[goodObs],
-                  EGray[goodObs])
+                  (EGray[goodObs]).astype(expGrayForInitialSelection.dtype))
         np.add.at(expGrayRMSForInitialSelection,
                   obsExpIndex[goodObs],
-                  EGray[goodObs]**2.)
+                  (EGray[goodObs]).astype(expGrayRMSForInitialSelection.dtype)**2.)
         np.add.at(expNGoodStarForInitialSelection,
                   obsExpIndex[goodObs],
                   1)
@@ -375,7 +376,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
 
         np.add.at(ccdGrayWt,
                   (obsExpIndex[goodObs],obsCCDIndex[goodObs]),
-                  1./EGrayErr2GO)
+                  (1./EGrayErr2GO).astype(ccdGrayWt.dtype))
         np.add.at(ccdNGoodStars,
                   (obsExpIndex[goodObs],obsCCDIndex[goodObs]),
                   1)
@@ -385,7 +386,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
                               obsBandIndex[goodObs]])
         np.add.at(ccdDeltaStd,
                   (obsExpIndex[goodObs], obsCCDIndex[goodObs]),
-                  obsDeltaStdGO/EGrayErr2GO)
+                  (obsDeltaStdGO/EGrayErr2GO).astype(ccdDeltaStd.dtype))
 
         gd = np.where((ccdNGoodStars >= 3) & (ccdGrayWt > 0.0))
         ccdDeltaStd[gd] /= ccdGrayWt[gd]
@@ -403,7 +404,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             FGrayGO = 10.**(EGrayGO/(-2.5))
             FGrayErrGO = (np.log(10.)/2.5)*np.sqrt(EGrayErr2GO)*FGrayGO
 
-            h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True)
+            h, rev = histogram_rev_sorted(obsExpIndex[goodObs])
 
             use, = np.where(h >= 3)
             for i in use:
@@ -462,10 +463,10 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
 
                     np.add.at(ccdGrayTemp,
                               obsCCDIndex[goodObs[i1a]],
-                              EGrayGO[i1a]/EGrayErr2GO[i1a])
+                              (EGrayGO[i1a]/EGrayErr2GO[i1a]).astype(ccdGrayTemp.dtype))
                     np.add.at(ccdGrayRMSTemp,
                               obsCCDIndex[goodObs[i1a]],
-                              EGrayGO[i1a]**2./EGrayErr2GO[i1a])
+                              (EGrayGO[i1a]**2./EGrayErr2GO[i1a]).astype(ccdGrayRMSTemp.dtype))
 
                     gdTemp, = np.where((ccdNGoodStars[eInd, :] >= 3) &
                                        (ccdGrayWt[eInd, :] > 0.0) &
@@ -494,7 +495,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
                     ccdGrayEvalNStars = np.zeros(self.fgcmPars.nCCD, dtype=np.int32)
                     np.add.at(ccdGrayEval,
                               obsCCDIndex[goodObs[i1a]],
-                              ccdGrayEvalStars)
+                              ccdGrayEvalStars.astype(ccdGrayEval.dtype))
                     np.add.at(ccdGrayEvalNStars,
                               obsCCDIndex[goodObs[i1a]],
                               1)
@@ -535,10 +536,10 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             # we can do all of this at once.
             np.add.at(ccdGray,
                       (obsExpIndex[goodObs],obsCCDIndex[goodObs]),
-                      EGrayGO/EGrayErr2GO)
+                      (EGrayGO/EGrayErr2GO).astype(ccdGray.dtype))
             np.add.at(ccdGrayRMS,
                       (obsExpIndex[goodObs],obsCCDIndex[goodObs]),
-                      EGrayGO**2./EGrayErr2GO)
+                      (EGrayGO**2./EGrayErr2GO).astype(ccdGrayRMS.dtype))
 
             # need at least 3 or else computation can blow up
             gd = np.where((ccdNGoodStars >= 3) & (ccdGrayWt > 0.0) & (ccdGrayRMS > 0.0))
@@ -574,7 +575,7 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
             expCcdHash = (obsExpIndex[goodObs]*(self.fgcmPars.nCCD + 1) +
                           obsCCDIndex[goodObs])
 
-            h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+            h, rev = histogram_rev_sorted(expCcdHash)
 
             # Anything with 2 or fewer stars will be marked bad
             use, = np.where(h >= 3)
@@ -678,25 +679,25 @@ def computeCCDAndExpGray(self, onlyObsErr=False):
 
         np.add.at(expGrayWt,
                   goodCCD[0],
-                  1./ccdGrayErr[goodCCD]**2.)
+                  (1./ccdGrayErr[goodCCD]**2.).astype(ccdGrayWt.dtype))
         np.add.at(expGray,
                   goodCCD[0],
-                  ccdGray[goodCCD]/ccdGrayErr[goodCCD]**2.)
+                  (ccdGray[goodCCD]/ccdGrayErr[goodCCD]**2.).astype(expGray.dtype))
         np.add.at(expGrayRMS,
                   goodCCD[0],
-                  ccdGray[goodCCD]**2./ccdGrayErr[goodCCD]**2.)
+                  (ccdGray[goodCCD]**2./ccdGrayErr[goodCCD]**2.).astype(expGrayRMS.dtype))
         np.add.at(expNGoodCCDs,
                   goodCCD[0],
                   1)
         np.add.at(expNGoodTilings,
                   goodCCD[0],
-                  ccdNGoodTilings[goodCCD])
+                  (ccdNGoodTilings[goodCCD]).astype(expNGoodTilings.dtype))
         np.add.at(expNGoodStars,
                   goodCCD[0],
-                  ccdNGoodStars[goodCCD])
+                  (ccdNGoodStars[goodCCD]).astype(expNGoodStars.dtype))
         np.add.at(expDeltaStd,
                   goodCCD[0],
-                  ccdDeltaStd[goodCCD]/ccdGrayErr[goodCCD]**2.)
+                  (ccdDeltaStd[goodCCD]/ccdGrayErr[goodCCD]**2.).astype(expDeltaStd.dtype))
 
         if self.fgcmPars.nCCD >= 3:
             # Regular mode, when we have a multi-detector camera.
@@ -811,13 +812,13 @@ def computeExpGrayColorSplit(self):
 
             np.add.at(expGrayColorSplit[:, c],
                       obsExpIndex[goodObs[use]],
-                      EGrayGO[use] / EGrayErr2GO[use])
+                      (EGrayGO[use] / EGrayErr2GO[use]).astype(expGrayColorSplit.dtype))
             np.add.at(expGrayWtColorSplit[:, c],
                       obsExpIndex[goodObs[use]],
-                      1. / EGrayErr2GO[use])
+                      (1. / EGrayErr2GO[use]).astype(expGrayWtColorSplit.dtype))
             np.add.at(expGrayRMSColorSplit[:, c],
                       obsExpIndex[goodObs[use]],
-                      EGrayGO[use]**2. / EGrayErr2GO[use])
+                      (EGrayGO[use]**2. / EGrayErr2GO[use]).astype(expGrayRMSColorSplit.dtype))
             np.add.at(expGrayNGoodStarsColorSplit[:, c],
                       obsExpIndex[goodObs[use]],
                       1)
@@ -1206,7 +1207,7 @@ def computeCCDAndExpDeltaMagBkg(self):
         goodObs, = np.where(obsFlag == 0)
 
         # Do the exposures first
-        h, rev = esutil.stat.histogram(obsExpIndex[goodObs], rev=True)
+        h, rev = histogram_rev_sorted(obsExpIndex[goodObs])
 
         expDeltaMagBkg[:] = self.illegalValue
 
@@ -1228,7 +1229,7 @@ def computeCCDAndExpDeltaMagBkg(self):
         # Do the exp/ccd second
         expCcdHash = obsExpIndex[goodObs]*(self.fgcmPars.nCCD + 1) + obsCCDIndex[goodObs]
 
-        h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+        h, rev = histogram_rev_sorted(expCcdHash)
 
         # We need at least 3 for a median, and of those from the percentile...
         use, = np.where(h > int(3./self.deltaMagBkgOffsetPercentile))
@@ -1289,7 +1290,7 @@ def computeExposureReferenceOffsets(self):
             obsYGO = snmm.getArray(self.fgcmStars.obsYHandle)[goodObs]
             expCcdHash = (obsExpIndexGO[ok] * (self.fgcmPars.nCCD + 1) +
                           obsCCDIndexGO[ok])
-            h, rev = esutil.stat.histogram(expCcdHash, rev=True)
+            h, rev = histogram_rev_sorted(expCcdHash)
             use, = np.where(h > 0)
             for i in use:
                 i1a = rev[rev[i]: rev[i + 1]]
@@ -1318,7 +1319,7 @@ def computeExposureReferenceOffsets(self):
 
         # And then this can be split per exposure.
 
-        h, rev = esutil.stat.histogram(obsExpIndexGO[goodRefObsGO], rev=True)
+        h, rev = histogram_rev_sorted(obsExpIndexGO[goodRefObsGO])
 
         use, = np.where(h >= self.minStarPerExp)